ルンゲクッタ法

考え方は修正オイラー法と同じです。よって詳細な説明を省略しますが、ルンゲクッタ法ではテイラー展開において4次の項まで利用します。やや求めるのが大変ですが、以下のような式が導けます。
ルンゲクッタ法ではこの式を用いて数値計算を行います。また、プログラムのソースコードを書く際に、考えやすいように上記のような形に書いてあります。修正オイラー法等と比較すると計算回数は多くなりますが、非常に精度が高いです。

ここでは単振動の問題をルンゲクッタ法を使って解いてみます。質量の値とばね定数の値を1とすると運動方程式は
となりますが次のように変形します。
vは速度です。常微分方程式が2つ出てきましたが難しいことはありません。それぞれにルンゲクッタ法を適用して解いていきます。

#include <stdio.h>

double f1(double t,double x,double v);
double f2(double t,double x,double v);

int main()
{
        double t,x,v,dt,tmax;
        double k1[2],k2[2],k3[2],k4[2];

        x=1.0;                  //位置の初期値
        v=0.0;                  //速度の初期値
        dt=0.01;                //刻み幅
        tmax=500.0;             //繰り返し最大回数

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

        for(t=0;t<tmax;t+=dt) {
                k1[0]=dt*f1(t,x,v);
                k1[1]=dt*f2(t,x,v);

                k2[0]=dt*f1(t+dt/2.0,x+k1[0]/2.0,v+k1[1]/2.0);
                k2[1]=dt*f2(t+dt/2.0,x+k1[0]/2.0,v+k1[1]/2.0);

                k3[0]=dt*f1(t+dt/2.0,x+k2[0]/2.0,v+k2[1]/2.0);
                k3[1]=dt*f2(t+dt/2.0,x+k2[0]/2.0,v+k2[1]/2.0);

                k4[0]=dt*f1(t+dt,x+k3[0],v+k3[1]);
                k4[1]=dt*f2(t+dt,x+k3[0],v+k3[1]);

                x=x+(k1[0]+2.0*k2[0]+2.0*k3[0]+k4[0])/6.0;
                v=v+(k1[1]+2.0*k2[1]+2.0*k3[1]+k4[1])/6.0;

                fprintf(output,"%f %f %f\n",t,x,v);
        }

        fclose(output);

        return 0;
}

double f1(double t,double x,double v)
{
        return v;
}

double f2(double t,double x,double v)
{
        return (-x);
}

下の画像は数値計算の結果をGNUPLOTで表示したものです。横軸が位置xを表し、縦軸が速度vを表しています。単振動において、外部からの力や摩擦や抵抗がなければ、位置と速度の最大値と最小値は変わらないから、一定の円を描きます。ルンゲクッタ法で解いた結果はそのようになっています。オイラー法と比べると、その精度の高さはひと目でわかりますが、修正オイラー法と比べると、GNUPLOT出力結果からは判断が難しいですが、数値データの出力結果を見るとルンゲクッタ法の精度の高さがよくわかります。