显式欧拉法求解热方程时生成无限大TXT文件的问题求助

核心问题：数组越界访问 #

你定义的数组是double T[N+1][N];，意味着：

• 第一维（时间步）的有效索引范围是 0 到 N（共 N+1 个元素）
• 第二维（空间步）的有效索引范围是 0 到 N-1（共 N 个元素）

但你的多个循环都超出了这个范围，直接导致内存破坏，程序行为异常（包括生成异常大的输出文件）。同时输出换行的代码写法错误，也会加剧格式异常。

逐点修复说明 #

1. 边界初始化循环越界 #

原代码循环范围超出数组索引上限，修正后：

// 初始化所有时间步的空间边界（z=0和z=1）
for (i=0;i<=N;i++){  // i范围0到N，覆盖T的全部时间步
    T[i][0] = b;
    T[i][N-1] = b;
}

// 初始化初始时间步的所有空间点
for (j=0;j<N;j++){  // j范围0到N-1，覆盖T的全部空间步
    T[0][j] = b;
}

2. 数值计算循环越界 #

原代码的时间步循环会访问超出数组范围的索引，修正后：

// 显式欧拉法计算后续时间步
for (i=0;i<N;i++){  // i从0到N-1，计算i+1到N（共N次迭代）
    for (j=1;j<N-1;j++){  // j从1到N-2，避开已固定的边界点
        T[i+1][j] = (dt/pow(dz, 2))*(T[i][j+1] - 2*T[i][j] + T[i][j-1]) + dt + T[i][j];
    }
}

3. 文件输出循环越界与换行错误 #

原代码的输出循环同样越界，且换行符写法错误，修正后：

// 输出所有有效数据：N+1行，N列
for (i=0;i<=N;i++){
    for (j=0;j<N;j++){
        fprintf(output, "%lf\t", T[i][j]);
    }
    fprintf(output,"\n");  // 使用标准换行转义符
}

修复后的完整代码 #

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

// 按需求定义N+1行（时间步0到N）、N列（空间步0到N-1）的矩阵
double T[N+1][N];

int main(){
    int i,j;
    double dt=1./N;
    double dz=1./N;
    double b=43351./94400;
    // 预计算扩散系数，避免重复调用pow函数
    double alpha = dt / (dz * dz);
    
    // 初始化空间边界
    for (i=0;i<=N;i++){
        T[i][0] = b;
        T[i][N-1] = b;
    }
    
    // 初始化初始时间步
    for (j=0;j<N;j++){
        T[0][j] = b;
    }
    
    // 显式欧拉法迭代计算
    for (i=0;i<N;i++){
        for (j=1;j<N-1;j++){
            T[i+1][j] = alpha*(T[i][j+1] - 2*T[i][j] + T[i][j-1]) + dt + T[i][j];
        }
    }

    FILE* output = fopen("numerica.txt", "w");
    if (!output) {  // 增加文件打开失败的错误判断
        perror("Failed to open file");
        return 1;
    }
    
    // 输出结果
    for (i=0;i<=N;i++){
        for (j=0;j<N;j++){
            fprintf(output, "%lf\t", T[i][j]);
        }
        fprintf(output,"\n");
    }
    fclose(output);
    return 0;
}

额外优化与注意事项 #

• 稳定性检查：显式欧拉法求解热方程需满足CFL条件alpha <= 0.5，你当前参数下alpha=1，会导致数值解发散。建议调整时间步，比如设dt = 0.4/(N*N)，确保稳定性。
• 边界条件确认：当前代码将所有初始值和边界值设为b，需根据你的实际物理场景调整。

内容的提问来源于stack exchange，提问作者ErikLAndre