热方程的数值解

Question

在用解析法求解热方程之后，我试图用显式欧拉法数值求解它。我得到了后续离散化，其中T是温度。我的代码是：

#include <string.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define N 100

double T[N+1][N];

int main(){
    int i,j;
    double dt=1./N;
    double dz=1./N;
    double b=43351./94400;
    for (i=0;i<N+2;i++){
        T[i][0]=b;
        T[i][N-1]=b;
    }
    
    for (j=0;j<N+1;j++){
        T[0][j+1]=b;
    }
    for (i=0;i<N+1;i++){
        for (j=1;j<N;j++){
            T[i+1][j] = (dt/pow(dz, 2))*(T[i][j+1] - 2*T[i][j] + T[i][j-1]) + dt + T[i][j];
            
        }
    }

    FILE* output;
    
    output = fopen("numerica.txt", "w");
    
    for (i=0;i<N+2;i++){
        for (j=0;j<N+1;j++){
            fprintf(output, "%lf\t", T[i][j]);
        }
        fprintf(output,"\n");
    }
    fclose(output);
    return 0;
        
}

我要做的是创建一个保存函数中所有值的N+1xN矩阵。编译之后，我有一个无限的.txt文件。有人帮忙吗？

Answer 1

输出文件应该有102行，除非程序崩溃，因为您的循环是错误的。

$ gcc -O -Wall t.c 
t.c: In function ‘main’:
t.c:35:13: warning: iteration 100 invokes undefined behavior [-Waggressive-loop-optimizations]
   35 |             fprintf(output, "%lf\t", T[i][j]);
      |             ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
t.c:34:19: note: within this loop
   34 |         for (j=0;j<N+1;j++){
      |                  ~^~~~
t.c:24:46: warning: iteration 98 invokes undefined behavior [-Waggressive-loop-optimizations]
   24 |             T[i+1][j] = (dt/pow(dz, 2))*(T[i][j+1] - 2*T[i][j] + T[i][j-1]) + dt + T[i][j];
      |                                          ~~~~^~~~~
t.c:23:19: note: within this loop
   23 |         for (j=1;j<N;j++){
      |                   ^
t.c:20:18: warning: iteration 99 invokes undefined behavior [-Waggressive-loop-optimizations]
   20 |         T[0][j+1]=b;
      |         ~~~~~~~~~^~
t.c:19:15: note: within this loop
   19 |     for (j=0;j<N+1;j++){
      |              ~^~~~
t.c:15:16: warning: iteration 101 invokes undefined behavior [-Waggressive-loop-optimizations]
   15 |         T[i][0]=b;
      |         ~~~~~~~^~
t.c:14:15: note: within this loop
   14 |     for (i=0;i<N+2;i++){
      |              ~^~~~

Answer 2

正如vm666所指出的，编译器注意到循环读取过数组T的末尾，这是未定义的行为。编译器可以以未定义的行为执行任何操作。在这种情况下--在循环的某个迭代中--编译器可能假设迭代101从未发生。这是对未定义行为的一种允许的解释，因为它(令人惊讶)通常是正确的。然后，编译器可以推断，如果迭代101从未发生，某种神奇的力量(比如调用exit)必须在迭代101之前停止循环，所以它永远不能超过迭代101，所以它不需要检查它是否超过了迭代101，而且如果程序从未真正检查并且循环是无限的，那么程序可能会更快一些。

狡猾的编译器。

如果这是实际原因，您可以通过不从数组外部打印数字来修复它。您似乎意外地在i和j循环中添加了另一个+1，所以只需再次取消它。