今回はある四角形の形をしたパルス派の進行を数値計算してみる。両端は固定端なので、パルス派は正負逆転するはずである。

#include <stdio.h>

#define N 100           //系の長さ
#define tmax 30         //時間的最大繰り返し回数
#define c 1.0           //波の速度
#define dt 1.0          //時間の刻み幅
#define dx 1.0          //系の刻み幅

int main()
{
        int t,i;
        double y[3][N];

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

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

                if(i>10 && i<20) {
                        y[0][i]=10.0;
                }
                else {
                        y[0][i]=0.0;
                }

                fprintf(output,"%d %d %f\n",0,i,y[0][i]);

        }
        fprintf(output,"\n");

/*初期状態から次の状態を計算する*/
        for(i=1;i<N-1;i++) {
                y[1][i]=y[0][i]+c*c/2.0*dt*dt/(dx*dx)*(y[0][i+1]+y[0][i-1]-2.0*y[0][i]);
        }
        y[1][0]=y[1][N-1]=0.0;          //境界条件
        for(i=0;i<N;i++) {
                fprintf(output,"%d %d %f\n",1,i,y[1][i]);
        }
        fprintf(output,"\n");

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

                for(i=1;i<N-1;i++) {
                        y[2][i]=2.0*y[1][i]-y[0][i]+c*c*dt*dt/(dx*dx)*(y[1][i+1]+y[1][i-1]-2.0*y[1][i]);
                }
                y[2][0]=y[2][N-1]=0.0;          //境界条件

                for(i=0;i<N;i++) {
                        fprintf(output,"%d %d %f\n",t,i,y[2][i]);
                }
                fprintf(output,"\n");

                for(i=0;i<N;i++) {
                        y[0][i]=y[1][i];
                        y[1][i]=y[2][i];
                }

        }

        fclose(output);

        return 0;
}

下の画像は、0 から 100 まである横軸は系の位置を表し、もう一方の横軸は時間の経過を表す。縦軸はパルス波の高さである。初期状態のパルスは、左右に半分ずつの高さになり別れ、一方では固定端で反射して、正負逆転していることがわかる。

下の画像は、上の画像のパルス波の高さを色で塗り分けた 2次元の図である。縦軸は波の位置を、横軸は時間経過を表す。