パソコンで数値計算

オイラー法

常微分方程式をテイラー展開を用いて数値的に解けるような形にする

常微分方程式の基本形
   (1)
において座標点を考え、におけるを求めたいとします。hは微小な値とします。そこでテイラー展開
   (2)
を用いて式(1)と(2)からhの1次までの項を考えるとすると
   (3)
という形の式が導けます。この式を利用して数値計算を行っていきますが、テイラー展開の1次までしかとってないので誤差が大きく精度が低いです。

プログラムソース

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

#include <stdio.h>

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

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

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

/*初期値*/
        x=1.0;
        v=0.0;
        dt=0.01;
        tmax=100;

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

                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を表しています。単振動において、外部からの力や摩擦や抵抗がなければ、位置と速度の最大値と最小値は変わらないから、一定の円を描くはずですが、オイラー法で解いた結果はそうなっていません。この原因は誤差によるものです。



2009/04/19 更新