#include <stdio.h>

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

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

        dt=0.1;                 //刻み幅設定
        x=1.0;                  //位置初期値
        v=0.0;                  //速度初期値
        tmax=100;               //繰り返し最大回数

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

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

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

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

        fclose(output);

        return 0;
}

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

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

  

下の画像を見ると、わずかだが、まだ誤差の影響が見られる。ここでは繰り返し回数が少ないからわかりにくいが、繰り返し回数を多くすれば、一定の円からずれていく。しかし、オイラー法よりは誤差の影響がずいぶんと小さくなったといえる。この方法よりももっと優れたルンゲクッタ法というものがある。常微分方程式を解く時は、ルンゲクッタ法を使うのがよい。