微分方程式の数値計算（オイラー法）-数学アルゴリズム演習ノート-

数値計算による物理シミュレーションで出てくる微分方程式は、積分など解析的な方法では解くのが困難、もしくは解けない事があります。しかし、微分方程式とは言ってみれば変数の「変化」のしかたを表した方程式ですので、単純にその方程式にしたがって変数を変化させてその結果を積み重ねていけば（数値積分）、大体の様子はわかるはずですね。さらにそのように数値計算で解くのは、「計算機」たるコンピュータの最も得意とする分野とも言えるでしょう。

例えば、dy/dx=2x、x=1でy=1という微分方程式があったとします。これは積分すれば解析的にも簡単に解けます（y=x²)が、初期値を入れた後、y をdy/dx=2x に従って変化させていっても大体の値は求まるわけです。まず、x=1でのyは1で傾きはその時点でのx*2の値,すなわち2です。これを使ってx=1.1でのyを求めてみましょう。

すでに傾きがわかっているので、xの増加分に傾きをかけてやればyの増加分が求まります。後はそれをy に加えれば良いわけです。yの変化量を計算すると、2*0.1=0.2ですから、それにyの元の値1を加えると1.2になります。実際の値（微分方程式を積分して解析的に解いた値）は、(1.1)²=1.21 ですので、わりと近い値と言えるでしょう。

オイラー法とは、このように微分方程式で表された変数の変化率を代入していって変数の値を求め、微分方程式を解く数値計算のアルゴリズムです。このように数値を入れながら数値計算で解いていく方法だと、数値を代入できる微分方程式の形に表された問題なら、どんな問題でも（とりあえず）解く事ができます。ただし、これで求まる解はあくまでも近似値であり、条件によっては誤差がかなり大きくなるので注意してください。
特に物理シミュレーションなどの科学計算では、ある程度の精度で誤差の評価を行えなければ計算結果もほとんど無意味なものになりかねません。さらに、数値計算の解は数値であって一般解や特殊解のような「式」は得られないので、解の「意味」をとらえにくい面もあるかもしれません。

では、例としてオイラー法で微分方程式 dy/dx=axy（a は任意の定数）を解くプログラムを作ってみましょう。初期値とaの値を入れたら後はオイラー法でひたすら数値積分していきます。その手順は、

となります。例えば、x が０から６までの値を幅ｈで求めるなら

for (x=0;x<6.0;x+=h) {

    dy=a*x*y; /* y の変化率を求める */
    y+=dy*h;  /* 変化率にx の幅をかけた変化量をy に加える */

}

という感じになるでしょう。これでx, yが求まるのでグラフを描いてみる事にします。幅を変えると誤差がどうなるかも調べてみたいので、幅も選べるようにしましょうか。
それから、解析的に解いた値（実際の値）も欲しいですね。簡単に解ける式ですから、解いておきましょう。

・微分方程式 dy/dx=axy の解

となります。これで、計算の準備が整いましたので、後は計算するだけです。ただ、絶対値が邪魔なので、今回はa の値とx, yの初期値をすべて正の値に限定しましょう。そうすればyは常に正ですから絶対値を外せます。

a, y0（x=0でのy）, h（幅）を指定したらG o ボタンをクリックしてください。h を変えながら描いてみると幅によってどれだけ誤差が出てくるかわかると思います。緑色の線が解析的に解いたグラフ（y=ｅ^(1/2ax²+C),C=logy0）で赤が数値計算の結果です。幅が0.02 だとほぼ一致しますが、0.1 ではかなりずれてきますね。