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

これで最終回です。全体的にはメモ程度ですが、何かの足がかりになった人がいれば幸いです。

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

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

時間・空間共に中心差分スキームを単純に適用するとリチャードソン(Richardson)スキーム(CTCSスキームとも)が得られます。ここで、空間に出てくる現在のタイムステップの項を、その前後の平均で置き換えると、次のデュフォート・フランケル(DuFort-Frankel)スキームが得られます。(これは蛙飛び(リープフロッグ・Leap-frog)とも呼ばれます。)

無条件不安定で精度も二次の陽的スキームですがはじめに一段階FTCSスキーム等を入れてやらなければなりません。

以下、最初に一段階FTCSスキームを入れたプログラム例です(C言語)。

#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,*u2;

  ul_int i,j,k;

  FILE *fp;

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

  scanf("%llu",&dt);

  printf("type FTCS 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));

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

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

  fclose(fp);

  for(i=0;i<layer;i++) *(u2+i)=*(u+i);

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

  }

  dt*=tstep;

  printf("Dufort time is %llu sec\n",dt);

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

  scanf("%llu",&tstep);

  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)=(nu/(dx*dx)*(*(c_vec+j+1)-*(u2+j)+*(c_vec+j-1))+*(u2+j)/(2*dt))/(1.0/(2*dt)+nu/(dx*dx));

    *u=*(u+1);

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

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

  }

    

  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;

}

以上、お付き合いありがとうございました。