微分方程式の数値計算（ルンゲ・クッタ法）-数学アルゴリズム演習ノート-

今回は、前回のオイラー法よりも高精度に微分方程式の数値解を求められるルンゲ・クッタ法を試してみます。ルンゲ・クッタ法は、数値計算のときに「中間点」を使って精度を上げるもので、計算も簡単なため広く使われているようです。

微分方程式dy/dxがx, y の関数になっている（dy/dx=f(x, y)）なら、その微分方程式をルンゲ・クッタ法で解くアルゴリズムは以下の通りです。ここでは、計算する刻みの幅をh、y の増加量をdy とし、x, y は現時点のx ,yを使います。

h/2 という中間点も使って計算したk1 ,k2 ,k3 ,k4という中間結果に重みをかけてdy を求めるのが特徴で、dy を求めたらx にはh を、y にはdy を加えて次のx ,yを計算していきます。なぜそれで精度が向上するかは....良くわからん（^^;。
まあ、中間点を計算する事で、誤差を修正している、らしいのですが、k2で一度中間点を計算した後で、さらにその値を使ってk3という中間点の値を再計算するのがポイントなのでしょうか？
これを使って前回の微分方程式 dy/dx=axy を計算するには、f(x, y)=axyとしてやれば良いわけですね。では、さっそくプログラム例です。

上で見たようにルンゲ・クッタ法で微分方程式を解いてx ,yを求めていく部分は、

for (x=0;x<6.0;x+=h) { /* ルンゲ・クッタ法で解いて表示 */

    k1=f(x,y)*h;
    k2=f(x+h/2,y+k1/2)*h;
    k3=f(x+h/2,y+k2/2)*h;
    k4=f(x+h,y+k3)*h;

    dy=(k1+2*k2+2*k3+k4)/6;
    y+=dy;

}

double f(double x,double y) { /* f(x,y) */

    return a*x*y;

}

と定義しておきます。前回からの修正はこれだけ。で、精度は....解析的に解いた曲線と見事に一致していますね。完全に重なってしまい、一本の線にしか見えません。0.4 でも完全に一致していますから、オイラー法に比べれば大幅に精度が向上しているようです（もちろん誤差はあるはずですが、このグラフでは識別できない）。

a, y0（x=0 でのy）, h（幅）を指定したらG o ボタンをクリックしてください。緑色の線が解析的に解いたグラフ（y=ｅ^(1/2ax²+C),C=logy0）で赤がルンゲ・クッタ法による数値計算の結果です。