計算物理学２第 11 回：微分方程式の解法
ver. 2018/7/6
本講義の最後の課題として微分方程式を取り扱います。今回はその中でも常微分方程式で初期値が与えられている場合の数値解の求め方を解説します。

1 Euler 法
まず解きたい方程式は一階常微分方程式でdy= f (x, y),dxy(x0 ) = y0(1)dyh + O(h2 ) = y(x) + f (x, y)h + O(h2 )dx(2)の形で与えられているとしましょう。

y(x) の Taylor 展開はy(x + h) = y(x) +で与えられます。Euler 法では x + h での y の値を x での関数の値を用いて y(x) + f (x, y)h で近似する方法です。初期値の与えられている x = x0 から計算を始め、まずは x0 での y の値を y0 とセットします。続いて、x1 = x0 + h での y の値 y1 はy1 = y(x0 + h) = y0 + f (x0 , y0 )h(3)と x0 と y0 を使って求まります。続いて x2 = x0 + 2h での値は y2 = y(x2 ) = y(x0 + 2h) を x1 = x0 + h で展開して、

y2 = y1 +dy
h = y1 + f (x1 , y1 )hdx x=x1(4)となり、x1 と y1 の値を用いて計算することができます。このように一般に xn , yn を用いて xn+1 = xn + hでの yn+1 を求めることができます。

Euler 法は最も簡単に実装することができますが、Taylor 展開の 1 次までしか考慮しないため、精度がよくありません。また、同じ計算を何度も繰り返すため、初期値の点から離れるほど誤差が大きくなっていく問題があります。

2 Runge-Kutta 法
2.1 2 次の Runge-Kutta 法
初期値のある常微分方程式の問題は xn , yn = y(xn ) の値を用いて xn+1 = xn + h での y の値 yn+1 =y(xn+1 ) を求める問題となります。この厳密解は∫ xn+1∫ xn+1dyyn+1 = yn +f (x, y(x))dx
dx = yn +dx
xnxn(5)の積分によって与えられます。Euler 法では式 (3) や (4) を用いて x の値を進めていきますのでこの積分を∫ xn+1f (x, y(x))dx = f (xn , yn )h + O(h2 )xn1(6)

と近似たことに対応します。積分の評価を改良することによって、Euler 法よりも精度のよい計算を行うことができます。これがルンゲクッタ (Runge-Kutta) 法です。Euler 法では f (x, y) の (x, y) の値は積分の左側の点 xn とそこでの yn = f (xn ) を積分区間 [xn , xn+1 ] で通して使っていました。xn の代わりに、被積分関数の値を積分区間の中点で評価された点を用いることで精度の向上が期待できます。積分区間の中点の x はxn + h/2 で与えられます。y の中点は yn+1 が不明なので正確には与えられませんが、Euler 法を用いて推測します。Euler 法を用いると yn+1 はyn+1 = yn + hf (xn , yn )(7)となりますので、k1 ≡ hf (xn , yn ) とおいて、y の中点は yn + k1 /2 と書きます。これらの変数を被積分関数に代入して評価したものが 2 次の Runge-Kutta 法です。式をまとめると

k1 =hf (xn , yn ),hk1, yn + ),22yn+1 =yn + k2 + O(h3 )k2 =hf (xn +(8)
(9)
(10)

2 次の Runge-Kutta 法では h2 のオーダーまで考慮されているため、Euler 法と比べて精度がよくなっています。プログラムへの実装も簡単で、式 (8) から式 (10) を順番に実行することで計算できます。

2.2 4 次の Runge-Kutta 法
広く使われているのは 4 次の Runge-Kutta 法です。

k1 =hf (xn , yn ),(11)hk1, yn + ),
(12)
22hk2
k3 =hf (xn + , yn + ),(13)
22k4 =hf (xn + h, yn + k3 ),(14)1(15)
yn+1 =yn + (k1 + 2k2 + 2k3 + k4 ) + O(h5 )6
はじめの２つは 2 次の Runge-Kutta と同じです。まず Euler 法を行って y の変化分を計算し (k1 )、それを用
k2 =hf (xn +いて中点での値を用いて y の変化分を計算します (k2 )。続いて y の評価を先ほどの中点での変化分による評価での値に更新して y の変化分を計算し (k3 )、これを用いて右端の点での値を計算します (k4 )。これらの 4点を重み付きで足し合わせたものが 4 次の Runge-Kutta の公式となります。
最後の足し合わせは Simpson の公式に近いものになっています。特に f (x, y) が y によらない場合、つまり f (x, y) = g(x) の形であれば、∫ xn+1xn[()]hhg(x)dx =g(xn ) + 4g xn ++ g(xn+1 ) + O(h5 )62
1= (k1 + 2k2 + 2k3 + k4 ) + O(h5 )6(16)となり、Simpson の公式と一致します。

Runge-Kutta 法の導出は複雑ですが実装は 2 次の場合と同様に簡単に行えるためよく使われます。y の値を更新するステップとして式 (11) から式 (15) を順番にプログラムの中に組み込むことによって実装できます。

2

3 二階常微分方程式
続いて二階常微分方程式を考えます。

d2 y=fdx2(dyx, y,dx)(17)

これも Euler 法や Runge-Kutta 法を用いて解くことができます。二階常微分方程式は連立一階常微分方程式の形に変換します。dy/dx = z とおくと、式 (17) は

dy=z,dxdz=f (x, y, z),dx(18)(19)となります。これれは y および z の一階微分のみが含まれている連立一階微分方程式です。
それではこれを 4 次の Runge-Kutta 法で解くアルゴリズムを考えてみます。ある x = xn での y = yn , z =zn の値が既知であるとします。これらをもとに x = xn+1 = xn + h での yn+1 と zn+1 の値を計算します。4次の Runge-Kutta の式により

k1 =hzn ,()l1
k2 =h zn +,2)
(l2,k3 =h zn +2
k4 =h (zn + l3 ) ,1yn+1 =yn + (k1 + 2k2 + 2k3 + k4 ),6

l1 =hf (xn , yn , zn ),k1l1h, yn + , zn + ),222hk2l2
l3 =hf (xn + , yn + , zn + ),222
l4 =hf (xn + h, yn + k3 , zn + l3 ),1zn+1 =zn + (l1 + 2l2 + 2l3 + l4 ),6l2 =hf (xn +(20)
(21)
(22)
(23)
(24)と書けます。上から一行ずつ順番に実行することで連立微分方程式の場合でも Runge-Kutta 法の適用が可能です。また、二階のみならず高階の常微分方程式は同様にして新しい変数を導入することにより連立一階微分方程式の形に変換して解くことができます。

4 演習問題
微分方程式

dy= cos x,dxy(0) = 0について(39) ステップサイズ h = 0.1 として Euler 法で解を求めよ。x > 0 のみでよい。
(40) h = 0.1 として 4 次の Runge-Kutta 法で解を求めよ。x > 0 のみでよい。
速度に比例する空気抵抗がある場合の質点 m の鉛直方向の運動方程式mz̈ = −mg − κż

3を考える。これは連立一階微分方程式ż =v,v̇ = − g −κvmとみなすことができる。

(41) g = 9.8, κ/m = 0.5, z(t = 0) = 0.0, v(t = 0) = 100 として z(t) および v(t) を図示せよ。

4