1次元波動方程式を解く

プログラム

解法には差分法を用います。差分法とは、微分方程式を漸化式に近似し、離散的に解く解法です。
精度は多少落ちますが、比較的計算しやすく、直感的にもわかりやすい方法で、偏微分方程式は、だいたい差分法で解けます。

解くべき式は、

\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}

です。普通の波動方程式ですが、初期条件によっては紐を両端で持って中央を上下させるような運動です。

ここでu(x,t)は波の変位を表す関数であり、cは波速度です。この方程式を差分法によって解くために、空間方向と時間方向にグリッドを作り、各点の変位を数値で求めます。具体的には、以下のように差分をとります。

u(i,j+1) = u(i,j) + c^2dt/dx^2(u(i+1,j) - 2u(i,j) + u(i-1,j))

ここでu(i,j)はグリッド上の(i,j)の位置の変位を表しています。この式は、(i,j)周りの3点の変位を使って(i,j+1)の変位を求めるものであり、波動方程式を離散化したものです。この差分方程式は、時間ステップごとに解くことで波動方程式を数値的に解くことができます。

このコードでは、初期条件として与えられた関数f(x)をsin(2πx)としています。空間方向のグリッド数nx、空間方向の刻み幅dx、時間方向の刻み幅dt、波速度c、計算の終了時刻tmaxなどのパラメータも定義されています。初期条件の設定後、時間ステップのループに入り、先ほどの差分方程式を用いて、時間ステップごとに変位u(x,t)を計算しています。最後に計算結果を出力しています。

 u[i] = func(x, 0.0);

の部分は、初期条件を表しています。t=０における、各点の初期位置をsin2πxで表現しています。
（ちなみにxは離散の実数なので、sin2πxは常に0ではないことに注意。）

#include <stdio.h>
#include <math.h>

#define PI 3.14159265358979323846

double func(double x, double t) {
    // 初期条件として与えられた関数 f(x)
    return sin(2*PI*x);
}

int main() {
    int nx = 100;          // 空間方向のグリッド数
    double dx = 0.1;       // 空間方向の刻み幅
    double dt = 0.001;     // 時間方向の刻み幅
    double c = 1.0;        // 波速度
    double tmax = 1.0;     // 計算の終了時刻
    double t = 0.0;        // 現在の時間
    double u[nx];          // u(x, t)の配列
    double u_new[nx];      // 時間ステップt+dtでのu(x)の配列

    // 初期条件の設定
    for (int i = 0; i < nx; i++) {
        double x = i*dx;
        u[i] = func(x, 0.0);
    }


    // タイムステップのループ
    while (t < tmax) {
        // 端点以外の空間方向の差分を計算
        for (int i = 1; i < nx-1; i++) {
            u_new[i] = u[i] + c*c*dt/dx/dx*(u[i+1] - 2*u[i] + u[i-1]);
        }
        // 端点の差分を計算
        u_new[0] = u[0] + c*c*dt/dx/dx*(u[1] - 2*u[0] + u[nx-1]);
        u_new[nx-1] = u[nx-1] + c*c*dt/dx/dx*(u[0] - 2*u[nx-1] + u[nx-2]);
        // u(x, t+dt)を更新
        for (int i = 0; i < nx; i++) {
            u[i] = u_new[i];
        }
        // 時間を更新
        t += dt;
    }

    // 結果の出力
    for (int i = 0; i < nx; i++) {
        double x = i*dx;
        printf("%f %f\n", x, u[i]);
    }

    return 0;
}