みその計算物理学
ホーム はじめに リンク集
ルンゲクッタ法の数値計算例(C言語)
  • プログラムソース (1次元の単振動を解く)
    #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出力結果

    前のオイラー法と修正オイラー法に比べて、繰り返し最大回数を多くしたにもかかわらず、単振動の位相空間の特徴である一定の円が得られた。かなりの回数を繰り返しても、一定の円からずれることはない。非常にルンゲクッタ法は実用的である。