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

所要により、熱拡散方程式の数値解法を扱うことになり、それを調べてみて、日本語でまとまったものを見つけるのに苦労したので、自分のところで3回にわたってまとめてみました。必要だったのが1次元塩分拡散だったので、数値計算するのに必要な拡散係数などは塩分濃度のものに合わせていますが、それ以外は特段の変更を伴わなくとも一般に適用できるかと思います。

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

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

これを解くスキームの一つにFTCS(Forward Time Centered Space)スキームがあります。これは、次のようなスキームです。

時間に対して前進差分を、空間に対して中央差分を用いたスキームで、陽的です。

増幅因子を計算するとわかりますが、このスキームはタイムステップが空間ステップの2乗に比例する関数よりも小さくならなければ安定ではありません。これが、スキームの問題点です。タイムステップを速くするための方法として、次回はCrank-Nicholsonスキームを、その次は、この一連の記事を書くための元となったDuFort-Frankelスキームを紹介します。

以下に、C言語(C99)によるFTCSスキームの実装例(塩分濃度拡散)を示しておきます。

#include<stdio.h>

#include<stdlib.h>

typedef unsigned long long int ul_int;

typedef long double ldouble;

int main(void){

  ldouble dx,c;

  const ldouble nu=1.47E-9;

  ldouble length;

  ul_int layer,tstep,dt;

  char filename[128];

  ldouble *u,*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));

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

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

  fclose(fp);

  c=(nu*dt)/(double)(dx*dx);

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

    for(j=0;j<layer;j++) *(c_vec+j)=*(u+j);

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

    *u=*(u+1);

    *(u+layer-1)=*(u+layer-2);

  }

    

  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);

  return 0;

}