2次元波動方程式の数値計算例(C言語)

プログラムソース

ここでは波動方程式(PDF形式)に書いてある2次元波動方程式を数値的に求めるアルゴリズムの通りに数値計算してみる。初期条件は四角形の形をした波とした。

#include <stdio.h>

#define N 100           //系の長さ
#define c 1.0           //波の速度
#define dt 0.1          //時間の刻み幅
#define dd 2.0          //系の刻み幅

int main()
{
        int t,i,j,tmax;
        double z[3][N][N];

        FILE *output;
        output=fopen("output.data","w");

        printf("時間的最大繰り返し回数を整数で３以上の値を入力してください。\n");
        scanf("%d",&tmax);

/*波の初期条件の決定*/
        for(i=0;i<N;i++) {
                for(j=0;j<N;j++) {

                        if(i>40 && i<60 && j>40 && j<60) {
                                z[0][i][j]=20.0;
                        }
                        else {
                                z[0][i][j]=0.0;
                        }
                }
        }

/*初期状態から次の状態を計算する*/
        for(i=1;i<N-1;i++) {
                for(j=1;j<N-1;j++) {
                        z[1][i][j]=z[0][i][j]+c*c/2.0*dt*dt/(dd*dd)*(z[0][i+1][j]+z[0][i-1][j]+z[0][i][j+1]+z[0][i][j-1]-4.0*z[0][i][j]);
                }
        }
        for(i=0;i<N;i++) {
                z[1][i][0]=z[1][i][N-1]=z[1][0][i]=z[1][N-1][i]=0.0;    //境界条件
        }

/*以降の波の状態を計算する*/
        for(t=2;t<tmax;t++) {

                for(i=1;i<N-1;i++) {
                        for(j=1;j<N-1;j++) {
                                z[2][i][j]=2.0*z[1][i][j]-z[0][i][j]+c*c*dt*dt/(dd*dd)*(z[1][i+1][j]+z[1][i-1][j]+z[1][i][j+1]+z[1][i][j-1]-4.0*z[1][i][j]);
                        }
                }

                for(i=0;i<N;i++) {
                        z[2][i][0]=z[2][i][N-1]=z[2][0][i]=z[2][N-1][i]=0.0;    //境界条件
                }

/*次の計算のために配列の数値を入れかえる。ここで過去の情報は失われる。*/
                for(i=0;i<N;i++) {
                        for(j=0;j<N;j++) {
                                z[0][i][j]=z[1][i][j];
                                z[1][i][j]=z[2][i][j];
                        }
                }
        }

/*最後の結果をファイルに書き込む*/
        for(i=0;i<N;i++) {
                for(j=0;j<N;j++) {
                        fprintf(output,"%d %d %f\n",i,j,z[2][i][j]);
                }
                fprintf(output,"\n");
        }
        fprintf(output,"\n");

        fclose(output);

        return 0;
}

GNUPLOT出力結果

プログラムを実行すると時間的最大繰り返し回数の入力を求められるので、結果を見たい時間経過に応じてその回数を入力する。以下の画像は、それぞれの時間的繰り返し回数別の結果である。

回数	結果
3
100
500
1000
2000
3000
4000
5000
10000

さて、これらの結果は正しいであろうか。