采用时间分裂法求解耦合微分方程组的C代码问题咨询

安全阿岳

2026-9-24

问题背景

我编写了用于求解耦合微分方程组的C代码，目标方程如下：

du/dt - Dd²u/dx² - v = 0
dv/dt - Dd²v/dx² + u = 0
求解采用时间分裂法，将方程拆分为二阶x导数项A₁、含u和v的交叉项A₂两部分，计算时使用两个0.5dt的半步长处理A₁，1个全步长dt处理A₂。目前代码运行后结果可以收敛，但数值解不符合预期，希望获得问题排查的方向指引。另外我使用gcc编译该代码时，会弹出与double main(void)相关的警告，也想了解该警告的产生原因。
代码内容如下：

/******************************************************************************

                            Online C Compiler.
                Code, Compile, Run and Debug C program online.
Write your code in this editor and press "Run" button to compile and execute it.

*******************************************************************************/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>

#define PI  3.141592

void read_input(double *D, double *L, int *nx, double *t_F);

double main(void) {
  /******************************/
  /* Declarations of parameters */
  /******************************/
  /* Number of grid points */                  
  int nx;
  /* Length of domain */
  double L;
  /* Equation coefficients */
  double D;
  /* Length of time to run simulation. */
  double t_F;

  /* Read in from file; */
  read_input(&D, &L, &nx, &t_F);

  /* Grid spacing */
  double dx = L/nx;
  double invdx2 = 1.0/(dx*dx);      
  /* Time step */
  double dt = 0.25/invdx2; // changed to 0.25/dx^2 to satisfy the stability condition

  /************************************************/
  /* Solution Storage at Current / Next time step */
  /************************************************/
  double *uc, *un, *vc, *vn;
  /* Time splitting solutions */
  double *uts1, *uts2, *vts1, *vts2;
  /* Derivative used in finite difference */
  double deriv;

  /* Allocate memory according to size of nx */
  uc = malloc(nx  * sizeof(double));
  un = malloc(nx * sizeof(double));
  vc = malloc(nx * sizeof(double));
  vn = malloc(nx * sizeof(double));
  uts1 = malloc(nx * sizeof(double));
  uts2 = malloc(nx * sizeof(double));
  vts1 = malloc(nx * sizeof(double));
  vts2 = malloc(nx * sizeof(double));

  /* Check the allocation pointers */
  if (uc==NULL||un==NULL||vc==NULL||vn==NULL||uts1==NULL||
  uts2==NULL||vts1==NULL||vts2==NULL) {
    printf("Memory allocation failed\n");
    return 1;
  }
  
  int k;
  double x;
  /* Current time */
  double ctime;

  /* Initialise arrays */
  for(k = 0; k < nx; k++) {
    x = k*dx;
    uc[k]  = 1.0 + sin(2.0*PI*x/L);
    vc[k]  = 0.0;
    /* Set other arrays to 0 */
    uts1[k] = 0; uts2[k] = 0;
    vts1[k] = 0; vts2[k] = 0;
  }

  /* Loop over timesteps */ 
  while (ctime < t_F){
  
    /* Rotation factors for time-splitting scheme. */
    double cfac = cos(dt); //changed from 2*dt to dt
    double sfac = sin(dt);

    /* First substep for diffusion equation, A_1 */ 
    for (k = 0; k < nx; k++) {
      x = k*dx;
      /* Diffusion at half time step. */
      deriv = (uc[k-1] + uc[k+1] - 2*uc[k])*invdx2 ;
      uts1[k] = uc[k] + (D * deriv + vc[k])* 0.5*dt; //
      deriv = (vc[k-1] + vc[k+1] - 2*vc[k])*invdx2;
      vts1[k] = vc[k] + (D * deriv - uc[k]) * 0.5*dt;
    }

    /* Second substep for decay/growth terms, A_2 */
    for (k = 0; k < nx; k++) {
      x = k*dx;
      /* Apply rotation matrix to u and v, */
      uts2[k] = cfac*uts1[k] + sfac*vts1[k];
      vts2[k] = -sfac*uts1[k] + cfac*vts1[k];
    }

    /* Third substep for diffusion terms, A_1 */
    for (k = 0; k < nx; k++) {
      x = k*dx;
      deriv = (uts2[k-1] + uts2[k+1] - 2*uts2[k])*invdx2;
      un[k] = uts2[k] + (D * deriv + vts2[k]) * 0.5*dt;
      deriv = (vts2[k-1] + vts2[k+1] - 2*vts2[k])*invdx2;
      vn[k] = vts2[k] + (D * deriv - uts2[k]) * 0.5*dt;
    }
       
    /* Copy next values at timestep to u, v arrays. */
    memcpy(uc,un, sizeof(double) * nx);
    memcpy(vc,vn, sizeof(double) * nx);

    /* Increment time. */   
    ctime += dt;
    for (k = 0; k < nx; k++ ) {
      x = k*dx;
      printf("%g %g %g %g\n",ctime,x,uc[k],vc[k]);
    }
    
  }
  /* Free allocated memory */
  free(uc); free(un);
  free(vc); free(vn);
  free(uts1); free(uts2);
  free(vts1); free(vts2);

  return 0;
}

// The lines below don't contain any bugs! Don't modify them
void read_input(double *D, double *L, int *nx, double *t_F) {
   FILE *infile;
   if(!(infile=fopen("input.txt","r"))) {
       printf("Error opening file\n");
       exit(1);
   }
   if(4!=fscanf(infile,"%lf %lf %d %lf",D,L,nx,t_F)) {
       printf("Error reading parameters from file\n");
       exit(1);
   }
   fclose(infile);
}

编译警告原因

C语言标准明确规定main函数的返回值必须为int类型，你声明为double main(void)不符合标准要求，因此gcc会触发返回值类型不匹配的警告。只需将返回值类型改为int main(void)即可消除该警告。

数值解问题排查方向

未初始化变量问题：变量ctime仅声明未赋值初始值，直接在while (ctime < t_F)中使用会读取到内存中的随机垃圾值，导致循环执行次数完全不符合预期，需在进入循环前添加ctime = 0.0;进行初始化。
数组越界访问问题：有限差分计算时，当k=0访问k-1=-1、k=nx-1访问k+1=nx均属于数组越界，读取的是非法内存的随机值。你需要先明确边界条件（如周期性边界、固定值边界、Neumann边界等），对边界点的差分计算做特殊处理，比如周期性边界可以将k-1替换为(k-1 + nx) % nx，k+1替换为(k+1) % nx。
时间分裂算子拆分错误：按照你的设计，A₁仅包含扩散项、A₂仅包含u/v交叉项，但现有代码在A₁的两个半步计算中都混入了交叉项（+vc[k]、-uc[k]），相当于交叉项被多计算了一次，完全偏离了原定的算子拆分逻辑。正确的A₁半步计算应当只保留扩散项：uts1[k] = uc[k] + D * deriv * 0.5*dt;，vts1[k] = vc[k] + D * deriv * 0.5*dt;，交叉项仅在A₂的全步中处理。
时间步稳定性问题：现有时间步dt = 0.25/invdx2 = 0.25*dx²未包含扩散系数D，显式扩散格式的稳定条件为dt <= dx²/(2D)，如果D的取值不等于1，当前dt可能不满足稳定条件，导致数值结果失真。
精度问题：自定义的PI精度仅为3.141592，会导致初始条件的正弦函数计算存在误差，建议改为更高精度的定义#define PI 3.141592653589793，或者直接使用math.h中定义的M_PI（部分编译器需要定义_USE_MATH_DEFINES后引入math.h才可使用）。

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