10. 常微分方程式 (1) 概要, 1段法, 多段法, 使われなくなった解法
目次
常微分方程式 (1): 概要, 1段法, 多段法, 使われなくなった解法
常微分方程式 (2): 誤差の推定, 補外法, ステップ幅の自動調節, 密出力, 2階常微分方程式, 遅延微分方程式
常微分方程式 (3): 解法の安定性, スティフな方程式, 微分代数方程式
常微分方程式 (4): 実用ルーチンのベンチマーク (実用ルーチンの選択, スティフではない問題, スティフな問題, 微分代数方程式)
常微分方程式の初期値問題について
初期値問題
1階の常微分方程式は次の形に書くことができる.
y' = f(t, y(t))
ここで, y(t) は未知の関数である. 常微分方程式の初期値問題(initial value problem)は, 初期点 t0 における関数値 y(t0) が与えられたとき, 最終点 tf における値 y(tf)を求める問題である. 初期点における関数値を与える次式は初期条件(initial condition)とよばれる.
y(t0) = y0
1階の連立方程式
1階の常微分方程式の初期値問題は, スカラー y をベクトル y に変えてそのまま連立方程式に拡張することができる.
y' = f(t, y) y(t0) = y0
例えば, 2変数(ベクトルの要素数が2)の場合について書き下すと次のようになる.
y1' = f1(t, y1, y2) y2' = f2(t, y1, y2)
この場合, 初期条件も2つ必要になる.
y1(t0) = y1,0 y2(t0) = y2,0
高階の常微分方程式
2階の常微分方程式は1階の連立常微分方程式の形にして扱うことができる. 例えば, 次の方程式を考える.
y'' = f(t, y, y')
これは, y1 = y, y2 = y’ とおけば, 次のように1階の連立常微分方程式になる.
y1' = y2 y2' = f(t, y1, y2)
この場合, 2つの初期条件が必要になる.
y1(t0) = y(t0) = y0 y2(t0) = y'(t0) = y'0
同様にしてさらに高階の常微分方程式を解くこともできる.
常微分方程式の初期値問題の数値解法
数値解法の分類
変数 t を離散化して, t0, t1, … , tn, … と表す. これらの点の間隔をステップ幅あるいはステップ幅(step size)といい, hn で表す.
hn = tn+1 - tn
ステップ幅は必ずしも一定でなくてよいが, n によらず一定の場合には単に h と表すことがある. 未知関数 y(t) の tn における値 y(tn) の近似値を yn と表すことにする.
差分法(微分を差分で近似する方法)により数値解法を導いてみる. 代表的な差分の求め方には, 次のようなものがある.
- 前進差分: (yn+1 – yn) / h
- 後退差分: (yn – yn-1) / h
- 中心差分: (yn+1 – yn-1) / 2h
まず, 前進差分を使って微分を近似する.
y'n = f(tn, yn) ≒ (yn+1 - yn) / h
これより次の近似式が得られる.
yn+1 = yn + hf(tn, yn)
この式を使い, 初期値 y0 から始めて, y1, y2, … と次々に各点における近似解を求めることができることがわかる. これをオイラー法 (Euler method) という.
後退差分を使った場合には, 同様にして次式により近似解を求めることができる.
yn+1 = yn + hf(tn+1, yn+1)
これを後退オイラー法 (backward Euler method) という.
後退オイラー法では右辺の f() の中に yn+1 が含まれており, yn+1 を求めるために方程式を解かなければならない(一般的には非線形方程式になる). このような解法を陰解法 (陰的公式, implicit formula) という. これに対して, オイラー法の場合には方程式を解くことなく直接 yn+1 を計算することができた. そのような解法を陽解法 (陽的公式, explicit formula) という.
次に, 中心差分を使った場合には, 次式により近似解を求めることができる.
yn+1 = yn-1 + 2hf(tn, yn)
これは中点則とよばれる.
中点則では, yn+1 の計算には yn の他に過去の値 yn-1 も必要で, 2つの値を使うので2段法とよばれ, 多段法 (多段階法, multistep method) とよばれる解法のひとつになる. これに対して, オイラー法と後退オイラー法は yn だけしか使わないので, 1段法 (1段階法, one-step method) とよばれる解法に分類される.
まとめると, 常微分方程式の初期値問題の解法を大きく分類すると1段法と多段法があり, それぞれに陽解法と陰解法がある. ここでは, 1段法の陽解法, および, 多段法について説明する. 1段法の陰解法についてはスティフな方程式の項で説明する.
数値解法の次数
関数 y(t) (真の解) を tn においてテイラー展開すると次のようになる.
y(tn+1) = y(tn) + hy'(tn) + (1/2)h2y''(tn) + (1/6)h3y'''(tn) + ... + (1/p)hpy(p)(tn) + ... = y(tn) + hf(tn, y(tn)) + (1/2)h2f'(tn, y(tn))+ (1/6)h3f''(tn, y(tn)) + ... + (1/p)hpf(p-1)(tn, y(tn)) + ...
数値解法による近似解 yn+1 の右辺が上のテイラー展開 y(tn+1) のp次の項まで一致するとき, すなわち次式が成り立つときにその数値解法はp次の方法であるという.
yn+1 - y(tn+1) = O(hp+1)
オイラー法はテイラー展開の1次(h)の項までで打ち切ったものに等しいから明らかに1次の方法である. また, 後退オイラー法は1次, 中点則は2次の方法であることを確かめることができる.
数値実験 (1)
1段法
テイラー法
関数 y(t) の tn におけるテイラー展開を再度次に示す.
y(tn+1) = y(tn) + hf(tn, y(tn)) + (1/2)h2f'(tn, y(tn))+ (1/6)h3f''(tn, y(tn)) + ... + (1/p)hpf(p-1)(tn, y(tn)) + ...
オイラー法はテイラー展開を1次の項までで打ち切って近似したものに相当した. そこで, より高次の項まで使った近似式を使えばオイラー法より精度のよい, すなわち, より次数の高い解法を得ることができると考えられる. これをテイラー法という.
上式でf'(t, y(t))とf”(t, y(t))は具体的には次のように計算する.
f'(t, y(t)) = ft(t, y(t)) + fy(t, y(t))・f(t, y(t)) f''(t, y(t)) = ftt(t, y(t)) + 2fty(t, y(t))・f(t, y(t)) + fyy(t, y(t))・f2(t, y(t)) + ft(t, y(t))・fy(t, y(t)) + fy2(t, y(t))・f(t, y(t))
2次の項を計算するのでさえ結構複雑なことがわかる. そのため, f()の高階の微分が簡単に計算できる関数でない限りこの方法を実際に適用するのは難しい.
数値実験 (2)
2次のルンゲ・クッタ法
ルンゲ・クッタ法 (Runge-Kutta method) は, f() の微分を計算することなく高次のテイラー展開に相当する精度を得る1段法である. いくつかのtとyについてf(t, y)を計算しそれらの重み付き平均を作ってテイラー展開のp次の項まで一致させることによりp次の公式を作る.
2次のルンゲ・クッタ法の一般形は次式のとおりである. a, b1, b2, c はパラメータである.
k1 = f(tn, yn) k2 = f(tn + ch, yn + hak1) yn+1 = yn + h(b1k1 + b2k2)
k2を (tn, yn)のまわりでテイラー展開すると次のようになる.
k2 = fn + chft,n + ak1fy,n + ...
これを上の式に代入して, 次式が得られる.
yn+1 = yn + (b1 + b2)hfn + b2h2(cft,n + afy,nfn + ... )
この式とテイラー展開の式を比較して, hとh2の係数が等しいとすると次式が得られる.
b1 + b2 = 1 b2c = 1/2 b2a = 1/2
未知数が4個で式が3本なので自由度が1ある. よく使われるものとしては, c = 1 としたものをホイン法(Heun method)とよぶ. また, c = 1/2 としたものを修正オイラー法とよぶことがある (2次のルンゲ・クッタ法とはあまりよばれないが, これらの総称ととらえてよい).
ホイン法の計算式は次のとおり.
k1 = f(tn, yn) k2 = f(tn + h, yn + hk1) yn+1 = yn + h(k1 + k2)/2
修正オイラー法の計算式は次のとおり.
k1 = f(tn, yn) k2 = f(tn + h/2, yn + hk1/2) yn+1 = yn + hk2
数値実験 (3)
ルンゲ・クッタ法の一般形
上と同様に3次, 4次・・・と高次のルンゲ・クッタ法の公式を作ることができ, 多くの公式が提案されている.
ルンゲ・クッタ法は一般的に次のように表される.
yn+1 = yn + hΣbiki (Σ は i = 1, 2, ... , s) ki = f(tn + cih, yn + hΣaijkj) (Σ は j = 1, 2, ... , s) (i = 1, 2, ... , s)
s は段数(stage)とよばれ, 1ステップ進むのに必要な f(t, y) の計算回数を表す(1段法や多段法の”段”とは異なる). aij, bi および ci はパラメータで, bi および ci は次式を満たす.
ci = Σaij (Σ は j = 1, 2, ... , s) (i = 1, 2, ... , s) Σbi = 1 (Σ は i = 1, 2, ... , s)
ルンゲ・クッタ法のパラメータには次のような表記法がある (ブッチャー配列という).
c1 | a11 a12 ... a1s c2 | a21 a22 ... a2s c3 | a31 a32 ... a3s . | . . cs | as1 as2 ... ass ---------------------------- | b1 b2 ... bs
aijの値によって3つの場合に分けられる.
(1) aij = 0 (i ≤ j) の場合 (右上三角部分と対角部分が0)
陽解法. (単に「ルンゲ・クッタ法」あるいは「陽的(explicit)ルンゲ・クッタ法」という)
c1 = 0 である.
ブッチャー配列ではaijの右上三角部分と対角部分を書かない(空白にする).
(2) aij = 0 (i < j) の場合 (右上三角部分が0)
半陰解法. (半陰的(semi-implicit)ルンゲ・クッタ法という)
(3) aij ≠ 0 (i < j) の場合 (右上三角部分に0でないものがある)
陰解法. (陰的(implicit)ルンゲ・クッタ法という)
ここでは(1)についてのみ説明する. (2), (3)については別途スティフな方程式の項で触れる.
上で示したホイン法と修正オイラー法をブッチャー配列で表すと次のようになる.
ホイン法
0 | 1 | 1 ---------------- | 1/2 1/2
修正オイラー法
0 | 1/2 | 1/2 ------------- | 0 1
(陽的)ルンゲ・クッタ法ではできるだけ高次の公式になるようにパラメータを決めるが, s段ならばs次になるわけではなく, s段の公式が達成しうる最大の次数pは次のようになることがわかっている.
段数 s | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 |
最大次数 p | 1 | 2 | 3 | 4 | 4 | 5 | 6 | 6 | 7 | 7 | 8 |
4次のルンゲ・クッタ法
4段4次の公式を作ろうとすると, パラメータが12個で自由度が4なので様々な公式を作ることが可能である. その中で最良のものを選ぶための多くの研究が行われた. 中でもクッタが選んだ以下の2組の係数が有名である.
古典的ルンゲ・クッタ法
単に「ルンゲ・クッタ法」というとこの公式を指す.
0 | 1/2 | 1/2 1/2 | 0 1/2 1 | 0 0 1 ----------------------- | 1/6 1/3 1/3 1/6
計算式は次のとおり.
k1 = f(tn, yn) k2 = f(tn + h/2, yn + hk1/2) k3 = f(tn + h/2, yn + hk2/2) k4 = f(tn + h, yn + hk3) yn+1 = yn + h(k1 + 2k2 + 2k3 + k4)/6
クッタの3/8公式
0 | 1/3 | 1/3 2/3 | -1/3 1 1 | 1 -1 1 ----------------------- | 1/8 3/8 3/8 1/8
計算式は次のとおり.
k1 = f(tn, yn) k2 = f(tn + h/2, yn + hk1/3) k3 = f(tn + h/2, yn + h(-k1/3 + k2)) k4 = f(, tn + h, yn + h(k1 - k2 + k3)) yn+1 = yn + h(k1 + 3k2 + 3k3 + k4)/8
ルンゲ・クッタ法一族において最もよく使われているのは4次の古典的ルンゲ・クッタ法である. いくつかある4段4次の公式間での精度の差はほとんどないが, この公式は係数aijに0が多く, すべての係数が正であるなどの理由により他のものより使いやすい. また, 5次の公式は6段以上でないと作れないため, これ以上高次のものを求めると複雑になるためこれが普及したのかもしれない.
数値実験 (4)
多段法
1段法は yn+1 の計算に現在の値 yn のみを使う解法であったが, 過去の値(それまでに計算済の値) yn-1, yn-2, … も使った方がよい結果が得られそうに思われる. このアイデアを実現したのが多段法である.
1階の常微分方程式の積分形を考える.
y(tn+1) = y(tn) + ∫ f(t, y(t)) dt [tn, tn+1]
ここで, 被積分関数f()を補間多項式で近似して右辺の積分を求める方法を一般にアダムス法(Adams method)という.
アダムス・バシュフォース法
y(tn+1) と y(tn)は近似値 yn+1 と yn を使い, 右辺の積分は (ti, fi) (i = n – k +1, …, n) を通る補間多項式 p(t) で f(t, y(t)) を置き換えると, 次のように k段の陽解法の公式が得られ, アダムス・バシュフォース法 (Adams-Bashforth method) とよばれる. ただし, fi = f(ti, yi) を表し, yi (i = n – k + 1, …, n) が計算済なので fi はすぐに求めることができる. この公式は, k次の公式になる.
yn+1 = yn + ∫ p(t) dt [tn, tn+1] = yn + hΣγj▽jfn (Σは j = 0, ... , k-1)
ただし,
▽0fn = fn, ▽j+1fn = ▽jfn - ▽jfn-1 (後退差分) γj = (-1)j ∫ C(-s, j) ds [0, 1]
である. ただし, C(n, k)は二項係数である. なお, γj は直接計算するよりも次の漸化式を用いて求めるとよい.
γ0 = 1, γj + (1/2)γj-1 + (1/3)γj-2 + ... + (1/(j+1))γ0 = 1
j = 0 ~ 7 のときの値は次のとおりである.
j | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
γj | 1 | 1/2 | 5/12 | 3/8 | 251/720 | 95/288 | 19087/60480 | 5257/17280 |
この方法を図に示すと次のようになる. 区間 [tn, tn+1] の積分を計算するために, 区間 [tn-k+1, tn] の補間データを使用した補間関数 p(t) を使う. 図で, 黒丸は補間点, 色つきの部分は積分範囲を表す.
k = 1 ~ 5 について, fn, fn-1, fn-2, … について整理すると次式のようになる.
k = 1 : yn+1 = yn + hfn k = 2 : yn+1 = yn + h/2(3fn - fn-1) k = 3 : yn+1 = yn + h/12(23fn - 16fn-1 + 5fn-2) k = 4 : yn+1 = yn + h/24(55fn - 59fn-1 + 37fn-2 - 9fn-3) k = 5 : yn+1 = yn + h/720(1901fn - 2774fn-1 + 2616fn-2 - 1274fn-3 + 251fn-4)
k = 1 の場合は1段法とみることができオイラー法に一致する.
なお, アダムス・バシュフォース法はk次である.
アダムス・モルトン法
アダムス・バシュフォース法では補間関数を使って外挿により積分を求めたが, 精度を上げるために tn+1を加えて区間 [tn-k+1, tn+1] の補間データを使う, すなわち, (ti, fi) (i = n – k + 1, …, n, n + 1) を通る補間多項式 p*(t) を使い, 外挿することを避けたのがアダムス・モルトン法 (Adams-Moulton method) である. こちらの方が精度がよさそうであるが, その代わり右辺に yn+1 が現れ陰解法になる. 式は次のとおりである. この公式は, k+1次の公式になる.
yn+1 = yn + ∫ p*(x) dx [tn, tn+1] = yn + hΣγ*j▽jfn+1 (Σは j = 0, ... , k)
ただし,
▽0fn = fn, ▽j+1fn = ▽jfn - ▽jfn-1 (後退差分) γ*j = (-1)j ∫ C(-s+1, j) ds [0, 1]
である. ただし, C(n, k)は二項係数である. なお, γ*j は次の漸化式を用いて計算するとよい.
γ*0 = 1, γ*j + (1/2)γ*j-1 + (1/3)γ*j-2 + ... + (1/(j+1))γ*0 = 0
j = 0 ~ 7 のときの値は次のとおりである.
j | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
γ*j | 1 | -1/2 | -1/12 | -1/24 | -19/720 | -3/160 | -863/60480 | -275/24192 |
この方法を図に示すと次のようになる. 下図で, 黒丸は補間点, 色つきの部分は積分範囲を表す. 積分範囲は補間区間に含まれており, 外挿されていない.
k = 0 ~ 5 について, fn+1, fn, fn-1, … について整理すると次式のようになる.
k = 0 : yn+1 = yn + hfn+1 k = 1 : yn+1 = yn + h/2(fn+1 + fn) k = 2 : yn+1 = yn + h/12(5fn+1 + 8fn - fn-1) k = 3 : yn+1 = yn + h/24(9fn+1 + 19fn - 5fn-1 + fn-2) k = 4 : yn+1 = yn + h/720(251fn+1 + 646fn - 264fn-1 + 106fn-2 - 19fn-3) k = 5 : yn+1 = yn + h/1440 (475fn+1 + 1427fn - 798fn-1 + 482fn-2 - 173fn-3 + 27fn-4)
k = 0 の場合は1段法とみることができ後退オイラー法に一致する. k = 1 の場合は台形則とよばれる公式に一致する.
なお, アダムス・モルトン法はk+1次である.
数値実験 (5)
予測子・修正子法
多段法においては, 陽解法は計算が容易であるが満足のいく精度が得られない. 一方, 陰解法は精度はよいが一般的には非線形方程式を毎回解かなければならないため計算が難しい.
そこで広く使われている方法が予測子・修正子法(predictor-corrector method)である. まず予測子(predictor)とよばれる陽解法によりyn+1の第1近似を求める. 次にその値を使って修正子(corrector)とよばれる陰解法を1回~数回適用しこれを改良する方法である.
予測子と修正子の組み合わせにはバリエーションがあるが, 予測子にアダムス・バシュフォース法, 修正子にアダムス・モルトン法を使った予測子・修正子法をアダムス・バシュフォース・モルトン法(あるいは単にアダムス法)とよび, 広く使われている.
予測子・修正子法の計算手順は, (1) 予測子を適用して Yn+1を求める(P), (2) それを使って関数値 fn+1を計算する(E), (3) 修正子を適用して Yn+1を更新する(C), (4) 更新した Yn+1を使って fn+1を更新する(E) となる. これを PECE と表すことにすると, 最後の E を省略して PEC とする方法, 修正子の適用を 2回繰り返して PECECE とする方法など, バリエーションがある.
数値実験 (6)
使われなくなった解法
よく知られているが現在では使われなくなった解法についてここで触れておく.
ルンゲ・クッタ・ジル法
ルンゲ・クッタ・ジル法(Runge-Kutta-Gill method)は係数が次のように表される4次のルンゲ・クッタ法である.
0 | 1/2 | 1/2 1/2 | (√2-1)/2 (2-√2)/2 1 | 0 (-√2)/2 1+(√2)/2 --------------------------------------- | 1/6 (2-√2)/6 (2+√2)/6 1/6
Gillの提案は係数の違いだけを見てもよくわからないが, その特長は次の2点である.
(1) 普通の4次のルンゲ・クッタ法では計算のために4個の変数(連立方程式の場合には配列)が必要なのに対し, この方法では係数を工夫して3個の変数(配列)だけで計算できるようになっていてその手順も示されている. [メモリの節約]
(2) 丸め誤差を減らすための手順が含まれている(これは他のルンゲ・クッタ法にも適用可能). [精度の向上]
メモリが少なくおそらく単精度計算が主流であった初期の計算機におけるプログラムに適した方法だったと思われる. かつて広く使われていたが, メモリが豊富で倍精度計算があたりまえの現代では必要性が薄れ使われなくなった.
ここで, (2)について少し説明する.
常微分方程式の基本的な数値解法では各ステップにおいて関数値の更新差分(ここでは kと表す)を計算し, 次のようにして関数値を更新する.
yn+1 = yn + k
この計算を計算機の浮動小数データ(ワード長は有限)で行うと下図のようになる.
ここで, k をワード長いっぱい正しく求めたとしても, この更新に使用されるのは色のついた一部だけで残りの部分は捨てられてしまう(丸め誤差という). そこで, この捨てられる部分を別の場所に累積しておき, 色のついた部分に頭が出るほどに貯まるたびにkに加えて関数値の更新に使用することにより丸め誤差を減らす方法が考案された. この方法がルンゲ・クッタ・ジル法に使われている.
数値実験 (7)
ミルン法
アダムス法と似ているが, 積分範囲が異なる次の積分形を考える.
y(tn+1) = y(tn-1) + ∫ f(t, y(t)) dt [tn-1, tn+1]
アダムス法と同様に (ti, fi) (i = n – k + 1, …, n) を通る補間多項式p(t)で積分を近似すると, 次の陽解法の公式が得られる.
yn+1 = yn-1 + hΣκj▽jfn (Σは j = 0, ... , k-1)
ただし,
κj = (-1)j ∫ C(-s, j) ds [-1, 1]
また, (ti, fi) (i = n – k +1, …, n+1) を通る補間多項式 p*(x) で積分を近似すると, 次の陰解法の公式が得られる.
yn+1 = yn-1 + hΣκ*j▽jfn+1 (Σは j = 0, ... , k)
ただし,
κ*j = (-1)j ∫ C(-s+1, j) ds [-1, 1]
予測子・修正子法の一種にミルン法(Milne method)とよばれる方法があり, 上の陰解法の式において k = 2 とした次の式がその修正子として使われている.
yn+1 = yn-1 + h/3 (fn+1 + 4fn + fn-1)
ミルン法の予測子には, 積分範囲をさらに [tn-3, tn+1] に広げた積分形より導かれる次式が使われている.
yn+1 = yn-3 + h/3 (8fn - 4fn-1 + 8fn-2)
ミルン法は4次の公式で, コンピュータが普及する以前の手計算の時代によく使われていた. しかし, コンピュータで使われるようになると安定性の面で問題があることがわかり, 現在では使われなくなった.
なお, 安定性の問題を解決するために, ミルン法の予測子と組み合わせて使う安定な修正子として次のハミング(Hamming) の修正子が提案されている.
yn+1 = 9/8 yn - 1/8 yn-2 + h/8 (3fn+1 + 6fn - 3fn-1)
数値実験 (8)
数値実験 (9)
参考文献
[1] 名取亮「数値解析とその応用」(1990) コロナ社[2] 三井斌友「常微分方程式の数値解法」(2003) 岩波書店
[3] E. Hairer, 他「常微分方程式の数値解法 I」(2007) シュプリンガー・ジャパン