ルンゲクッタ法

となる微分方程式の解を数値的に解く方法を考えてみる。

いま、x=x0におけるｙの値y0がわかっているとき、y(x0+h)のx0におけるテイラー展開は次のようになる。
　　　（1）

ここで、この式をhのn次の項まで打 ち切ると、

　　(2)

を用いてx1=x0+hのときのy1を求めることができる。
このときの誤差は、O(hn+1)となり、hのn+1乗に比例した誤差となる。すると、hが小さいほど、誤差も小さくなるわけである。

　さて、(2)式において、n=1までで計算を打ち切ったも のが先週のオイラー法になる。
当然,こ の方法ではhの値が大きいと誤差も大きくなるのでできるだけ小さいhを求めた方が精度よく計算でき る。
が、そのときには計算回数が大きくなってしまう。
一方、大きいn次まで採用すると、その誤差は小さくなるが、一回の計算に必要な処理が大きくなってしま う。
また、n次 の導関数の算出が出来ない場合もある。
そこで、x0とx1の中間のxにおける 値を計算しそれらの重み付き平均を用いることで、
テイラー展開した場合とおなじ次数までをカバーする方法がルンゲクッタ法です。

比較的簡単な場合のルンゲクッタ 法－2次 の場合



でy(x0+h)を算出するとする。
ここで、φはx0とx1の中間のxにおけるf(x)の値k1とk2の重み付き平均である。
k2に2変数関数のテイラー展開をして 2次以上を無視すると、次のようになる。



これを代入すると、

　　(3)

となる。
また、2次 までのテイラー展開の式は(1)式より

（4）

となる。
(3)式と(4)式が等しくなる条件は、w1+w2=1, w2a=w2b=1/2 である。
このとき、数値計算的なf(x,y)の2変数導関数fx(x,y)、fy(x,y)を直接算出する代わりに
f(x,y)の重み付き平均で求めることができる。
ここで、w2=1の場合、

となるが、これが2次のルンゲクッタ公式（修正オイラー法）である。


まとめ：

  1. 数値計算で微分方程式を解くときは、テイラー展開を利用する。
  2. n次の導関数を求めて計算する変わりに、xnとxn+1の中間値ｘにおける導関数を複数回計算 し、
    その重み付き平均を用いる。
  3. 誤差は、hの次数乗に比例する。
4次のルンゲクッタ法

一階の常微分方程式

dy/dx=f(x,y)

において、xがx0からhだけ増加した点x1=x0+hにおける値y1を以下の計算で求めます。
k1=hf(x0,y0)
k2=hf(x0+h/2,y0+k1/2)
k3=hf(x0+h/2,y0+k2/2)
k4=hf(x0+h,y0+k3)
k=(k1+2k2+2k3+k4)/6

y1=y0+k
これらの式の意味は図１に示す通りです。


図を捕捉すると、
①点(x0,y0)で傾きf(x0,y0)でhだけ進むとyはk1だけ変化することになる。
②この①の線上の点(x0+h/2,y0+k1/2)で傾きを計算し、(x0,y0)からhだけ進むとyはk2だけ変化することになる。
③同じように(x0,y0)から傾き(x0+h/2,y0+k2/2)でhだけ進むとk3の変化となる。
④(x0,y0)から傾き(x0+h,y0+k3)でhだけ進むとk4の変化となる。
⑤以上の計算結果の加重平均として点(x1,y1)を得る。

これがちょうど、x=x0におけるy(x0+h)の4次までのテーラー展開に一致します。


例題１
dy/dx=sin(x)+cos(y) を初期条件(x0,y0)=(0,0) のもとでRunge-Kuta法によりx=piまで
30等分して計算する。

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

double f(double x,double y)
{
    return sin(x)+cos(y);
}

void main()
{
    int i=0;
    double h,x,y,k1,k2,k3,k4,pi=3.141592;

    h = pi/30.0;
    x=0.0;
    y = 0.0;
    printf(" i x y\n");
    printf("%2d %lf %lf\n",i,x,y);

    for(x = 0 ; x <= pi ; x+=h){
        i+=1;
        k1 = h*f(x,y);
        k2 = h*f(x+h*0.5,y+k1*0.5); /*2で割るより、0.5掛ける方が計算が速い*/
        k3 = h*f(x+h*0.5,y+k2*0.5);
        k4 = h*f(x+h,y+k3);
        y = y + (k1+2.0*k2+2.0*k3+k4)/6.0;
        printf("%2d %lf %lf\n",i,x+h,y);
    }

}


連立常微分方程式も同様の計算で求めることができる。
具体例を用いて、２階常微分方程式を解くことを考える。

d2y/dx2=x*dy/dx+y

を初期条件(x0,y0)=(0,1)、dy/dx=1のもとで、Runge-Kutta法によりx=1.0まで
10等分して計算する。

この場合、上の式を次のように連立常微分方程式に変換して求める。
dy/dx=y' -----(a)
dy'/dx=xy'+y -----(b)

Runge-Kutta法はそれぞれの常微分方程式に対して計算すればよい。
例題2では例題１での計算を拡張し、配列によって２つの常微分方程式を計算している。
それぞれの計算をk1からk4までの加重平均で求めているのは例題１と同様である。
しかし、(b)の中に(a)が含まれている部分をyworkという変数によって逐次、(b)の計算に
組み入れている。

例題2

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

double df(int i, double x, double s[])
{
    /*    i:0 以下    dy/dx=y' の計算
        i:1 以上    dy'/dx=xy'+y　の計算 */
    if( i <= 0){
        return s[1];
    }
    else{
        return x*s[1]+s[0];
    }
}

int main(int argc,char *argv[])
{
    int i=0,j=0;
    double h,x,x0,y0,dy0,k1[2],k2[2],k3[2],k4[2],y[2],ywork[2];

    h = 0.1;
    x=0.0;
    y[0] = 1.0;        /*y[0]はy座標の計算値*/
    y[1] = 1.0;        /*y[0]はy方向の速度の計算値*/
    x0 = 0.0;
    y0 = 1.0;
    dy0 = 1.0;

    printf(" j x y0 y1\n");
    printf("%2d %lf %lf %lf\n", j, x, y[0], y[1]);

    for(j =1 ; j <= 10 ; j+=1){
        /*    i:0 以下    dy/dx=y' の計算
            i:1 以上    dy'/dx=xy'+y　の計算 */
       
        for( i=0; i<2; i++)
            k1[i]=h*df(i,x,y);
        for( i=0; i<2; i++)
            ywork[i]=y[i]+0.5*k1[i];
        for( i=0; i<2; i++)
            k2[i]=h*df(i,x+0.5*h,ywork);
        for( i=0; i<2; i++)
            ywork[i]=y[i]+0.5*k2[i];
        for( i=0; i<2; i++)
            k3[i]=h*df(i,x+0.5*h,ywork);
        for( i=0; i<2; i++)
            ywork[i]=y[i]+k3[i];
        for( i=0; i<2; i++)
            k4[i]=h*df(i,x+h,ywork);
        for( i=0; i<2; i++)
            y[i]=y[i]+(k1[i]+2.0*k2[i]+2.0*k3[i]+k4[i])/6.0;
        x=x+h;
        printf("%2d %lf %lf %lf\n", j, x, y[0], y[1]);

        x0=x;
        y0=y[0];
        dy0=y[1];
    }

    printf("Hit Any Key To End!");
    getch();

}

課題
1.　dy/dx=2y/(1+x)　をRunge-Kutta法で解きなさい。その際、計算区間を[0,1]とし、
計算間隔をh=0.1、初期値はx=0,のときのyの値を手計算で求めなさい。その際、結果は
ファイルに保存するプログラムとしなさい。

2.　減衰がある場合のバネの振動を考える。このモデルの運動方程式は、

Md2x/dt2 + Cdx/dt + Kx = 0

で与えられる。
初期条件として、M=１、C=2、K＝ 10、t=0のときx=４、Ｖ（つまり、dx/dt）＝－７として、運動方程式を
数値計算で解きなさい。
なお、この条件での厳密解は、 x = exp(-t) (4cos3t - sin3t) である。
結果を厳密解による結果とあわせてグラフで表示できるようにすること。

計算結果は、テキスト形式のファイルに保存すること。
また、数値計算と厳密解との差（誤差）も同時に保存しておくこと。