c - 有限差分一维热传导的解决方案不随着网格点的增加而减少误差

Question

我正在求解一维热传导方程。

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

void grid(int nx, double xst, double xen, double *x, double *dx)
{
  int i;

  *dx = (xen-xst)/(double)(nx-1);

  for(i=0; i<nx; i++)
    x[i] = (double)i * (*dx); // ensure x[0] == 0.0 and x[nx-1] == 1.0
}

void get_exact_solution(int nx, double time, double *x, double *Texact)
{
// calculates exact solution
  int i;

  for(i=0; i<nx; i++)
    Texact[i] = erf((x[i]-0.5)/2.0/sqrt(time));
}


void enforce_bcs(int nx, double *x, double *T)
{
  T[0] = -1.0;
  T[nx-1] = 1.0;
}

void set_initial_condition(int nx, double *x, double *T)
{
  int i;

  for(i=0; i<nx; i++)
    T[i] = tanh((x[i]-0.5)/0.1);

  enforce_bcs(nx,x,T); //ensure BCs are satisfied at t = 0
}

void get_rhs(int nx, double dx, double *x, double *T, double *rhs)
{
  int i;
  double dxsq = dx*dx;

  // compute rhs. For this problem, d2T/dx2
  for(i=1; i<nx-1; i++)
    rhs[i] = (T[i+1]+T[i-1]-2.0*T[i])/dxsq;
}

void timestep_Euler(int nx, double dt, double dx, double *x, double *T, double *rhs)
{

  int i;

  // compute rhs
  get_rhs(nx,dx,x,T,rhs);

  // (Forward) Euler scheme
  for(i=1; i<nx-1; i++)
    T[i] = T[i] + dt*rhs[i];  // T^(it+1)[i] = T^(it)[i] + dt * rhs[i];

  // set Dirichlet BCs
  enforce_bcs(nx,x,T);

}

void output_soln(int nx, int it, double tcurr, double *x, double *T)
{
  int i;
  FILE* fp;
  char fname[100];

  sprintf(fname, "T_x_%04d.dat", it);
  //printf("\n%s\n", fname);

  fp = fopen(fname, "w");
  for(i=0; i<nx; i++)
    fprintf(fp, "%lf %lf\n", x[i], T[i]);
  fclose(fp);
}

int main()
{

  int nx;
  double *x, *T, *rhs, tst, ten, xst, xen, dx, dt, tcurr, *Texact;
  int i, it, num_time_steps, it_print;
  FILE* fp;  


  nx =10; //grid points

  xst = 0, //starting and ending in space
  xen = 1;
  tst=0; // start and end in time
  ten=1


  printf("Inputs are: %d %lf %lf %lf %lf\n", nx, xst, xen, tst, ten);

  x = (double *)malloc(nx*sizeof(double));
  T = (double *)malloc(nx*sizeof(double));
  Texact = (double *)malloc(nx*sizeof(double));
  rhs = (double *)malloc(nx*sizeof(double));


  grid(nx,xst,xen,x,&dx);         // initialize the grid

  set_initial_condition(nx,x,T);  // initial condition

  // prepare for time loop
  dt = 0.001; 

  num_time_steps = (int)((ten-tst)/dt) + 1; // why add 1 to this?
  it_print = num_time_steps/10;             // write out approximately 10 intermediate results

  // start time stepping loop
  for(it=0; it<num_time_steps; it++)
  {
    tcurr = tst + (double)it * dt;

    timestep_Euler(nx,dt,dx,x,T,rhs);    // update T

    // output soln every it_print time steps
    //if(it%it_print==0)
      //output_soln(nx,it,tcurr,x,T);
  }

  get_exact_solution( nx, ten, x, Texact);

  int loop;
  double sq_diff,total,err;
  total = 0;

  for(loop = 0; loop < nx; loop++)
      printf("%f ", T[loop]);

  printf("\n");

  for(loop = 0; loop < nx; loop++)
      printf("%f ", Texact[loop]);

  for(loop = 0; loop < nx; loop++){
      sq_diff = (T[loop] - Texact[loop])*(T[loop] - Texact[loop]);
      total = total+sq_diff;
    }

  err = sqrt(total);
  printf("\n L2 Norm error = %f",err);


  printf("\n");
  free(rhs);
  free(T);
  free(x);

  return 0;
}

理想情况下，增加网格点的数量（nx在程序中）应该会减少我ten使用精确解Texact和数值解在最终时间计算的误差T。

误差应该随着小幅增加而减小，nx直到解稳定并遵循稳定性准则（冯诺依曼稳定性准则）。但在我的情况下，这并没有发生，并且错误随着nx.

我正在计算 和 之间的 L2_normT误差Texact。 Texact由函数计算get_exact_solution。

该图显示了不同情况下的问题。

请帮助纠正这个程序。

标签： cnumerical-methodseuclidean-distancevon-neumann