熱拡散方程式の数値解法について(2)

前回の続きです。今回はクランク・ニコルソン(Crank-Nicholson)スキームです。

熱拡散方程式の初期値問題

を数値的に解くことを考えます。

前回書いたFTCSスキームと、BTCSスキームを用います。BTCSスキームは、FTCSスキームと時間については同じ近似でありながら、基準を次の時刻として差分化する方法です。つまり、n+1タイムステップを基準にして、時間に対して後退差分を、空間に対して中心差分を取ります。

FTCSスキームとBTCSスキームの和より、次のクランク・ニコルソン(Crank-Nicholson)スキームが導かれます。

精度は時間・空間共に二次まで保証されるので、実用に耐えられると思います。増幅因子を計算すると、無条件安定スキームであることが言えますが、陰解法であり、並列化・ベクトルコンピューティングには難しいかもしれません。同じオーダーの精度であり、陽解法であるDuFort-Frankelスキームは次回紹介します。

以下は、やはり塩分濃度拡散の場合の、クランクニコルソンスキームの実装例です。C言語を用いています。

#include<stdio.h>

#include<stdlib.h>

typedef unsigned long long int ul_int;

typedef long double ldouble;

int main(void){

  ldouble c1,c2;

  ldouble dx;

  const ldouble nu=1.47E-9;

  ldouble length;

  ul_int layer,tstep,dt;

  char filename[128];

  ldouble *u;

  ldouble *matrix[3],*c_vec;

  ul_int i,j,k;

  FILE *fp;

  printf("type delta time(s)>>");

  scanf("%llu",&dt);

  printf("type timestep number (integer)>>");

  scanf("%llu",&tstep);

  printf("type z-length(m)(float)>>");

  scanf("%Lf",&length);

  printf("type z-layer number(integer)>>");

  scanf("%llu",&layer);

  dx=(ldouble)length/layer;

  printf("type data filename(a float number is contained per line)\n>>");

  scanf("%s",filename);

  if((fp=fopen(filename,"r"))==NULL){

    puts("Bad filename");

    return -1;

  }

  u=(ldouble *)calloc(layer,sizeof(ldouble));

  for(i=0;i<3;i++) matrix[i]=(ldouble *)calloc(layer,sizeof(ldouble));

  c_vec=(ldouble *)calloc(layer,sizeof(ldouble));

  for(i=0;i<layer;i++) fscanf(fp,"%Lf",u+i);

  fclose(fp);

  *matrix[0]=0;

  for(i=1;i<layer-1;i++) *(matrix[0]+i)=1;

  *(matrix[0]+layer-1)=-1;

  *matrix[2]=-1;

  for(i=1;i<layer-1;i++) *(matrix[2]+i)=1;

  c1=2+dx/nu*dx/dt;

  c2=4-c1;

  *matrix[1]=1;

  *(matrix[1]+layer-1)=1;

  for(i=1;i<layer-1;i++) *(matrix[1]+i)=-c1;

  for(i=1;i<layer;i++){

    *(matrix[0]+i)/=*(matrix[1]+i-1);

    *(matrix[1]+i)-=*(matrix[0]+i)*(*(matrix[2]+i-1));

  }

  for(i=0;i<tstep;i++){

    *c_vec=0;

    for(j=1;j<layer-1;j++) *(c_vec+j)=-*(u+j-1)+*(u+j)*c2-*(u+j+1);

    *(c_vec+layer-1)=0;

    for(j=1;j<layer;j++) *(c_vec+j)-=*(c_vec+j-1)* *(matrix[0]+j);

    *(u+layer-1)=*(c_vec+layer-1)/(*(matrix[1]+layer-1));

    for(j=layer-2;j>0;j--) *(u+j)=(*(c_vec+j)-*(matrix[2]+j)**(u+j+1))/(*(matrix[1]+j));

    *u=*(u+1);

  }

    

  printf("type output filename(a float number is contained per line)\n>>");

  scanf("%s",filename);

  if((fp=fopen(filename,"w"))==NULL){

    puts("Bad filename");

    return -2;

  }

  

  for(i=0;i<layer;i++) fprintf(fp,"%.8Lf\n",*(u+i));

  free(u);

  free(c_vec);

  for(i=0;i<3;i++) free(matrix[i]);

  return 0;

}