10. 常微分方程式 (1)
常微分方程式は科学技術計算においてしばしば現れるが、解析的に解けることは少なく数値解法に頼ることが多い。
常微分方程式の初期値問題
初期値問題
1階の常微分方程式は次の形に書くことができる。
y' = f(y, t)
ここで、y(t) は未知の関数である。常微分方程式の初期値問題(initial value problem)は、初期点 t0 における関数値 y(t0) が与えられたとき、最終点 tf における値 y(tf)を求める問題である。初期点における関数値を与える次式は初期条件(initial condition)とよばれる。
y(t0) = y0
1階の連立方程式
1階の常微分方程式の初期値問題は、スカラー y をベクトル y に代えてそのまま連立方程式に拡張することができる。
y' = f(y, t)
y(t0) = y0
例えば、2変数の場合について書き下すと次のようになる。
y1' = f1(y1, y2, t)
y2' = f2(y1, y2, t)
この場合、初期条件も2つ必要になる。
y1(t0) = y1,0
y2(t0) = y2,0
高階の常微分方程式
高階の常微分方程式、すなわち、2階以上の微分を含んだ常微分方程式は、1階の連立常微分方程式の形にして扱うことができる。例えば、次の方程式を考える。
y''' = f(y, y', y'', t)
これは、y1 = y, y2 = y’, y3 = y” とおけば、次のように1階の連立常微分方程式になる。
y1' = y2
y2' = y3
y3' = f(y1, y2, y3, t)
この場合、3つの初期条件が必要になる。
y1(t0) = y(t0) = y0
y2(t0) = y'(t0) = y'0
y3(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
前進差分を使って微分を近似すると次式が得られる。これを使い、初期値 y0 から始めて、y1, y2, … と次々に各点における近似解を求めることができる。
yn+1 = yn + hf(yn, tn)
これをオイラー法 (Euler method) という。
次に、後退差分を使った場合には、同様にして次式により近似解を求めることができる。
yn+1 = yn + hf(yn+1, tn+1)
これを後退オイラー法 (backward Euler method) という。この式を使うと、右辺の f の中に yn+1 が含まれており、yn+1 を求めるために一般的には非線形方程式を解かなければならない。このような解法を陰解法 (陰的公式、implicit formula) という。これに対し、上に示したオイラー法の場合には方程式を解くことなく直接 yn+1 を計算することができた。そのような解法を陽解法 (陽的公式、explicit formula) という。
次に、中心差分を使った場合には、次式により近似解を求めることができる。
yn+1 = yn-1 + 2hf(yn, tn)
これを中点法というが、yn+1 の計算に yn の他に過去の値 yn-1 の2つを使うので2段法とよばれ、多段法 (多段階法、multistep method) とよばれる解法の一種である。これに対して、オイラー法と後退オイラー法は yn だけしか使わないので、1段法 (1段階法、one-step method) とよばれる。
ここまで基本的な解法を例としてあげたが、常微分方程式の初期値問題の解法を大きく分類すると1段法と多段法があり、別の分類の仕方によると陽解法と陰解法がある。ここでは、1段法の陽解法および多段法について説明する。
1段法
もっとも基本的な1段法はオイラー法であるが、多くのより高精度な解法が開発されている。
テイラー法 (Taylor method)
関数 y を tn においてテイラー展開すると次のようになる。
y(tn+1) = y(tn) + hy'(tn) + (1/2)h2y''(tn) + (1/6)h3y'''(tn) + ...
式を見比べると、オイラー法はテイラー展開を1次の項までで打ち切って近似したものであることがわかる。従って、より高次の項まで使って近似すれば精度よく解を求めることがきると期待される。しかし、2次以上の項を計算するためには ｆの微分が必要となり汎用性に乏しいため、この方法はあまり使われない。
ただし、f の高階の微分が簡単に計算できる関数に限れば有効であるため、特殊関数の計算などに使われることがある。
ルンゲ・クッタ法 (RK法) (Runge-Kutta method)
RK法は、f の微分を計算することなく高次のテイラー展開に相当する精度を得られる1段法である。パラメータの決め方により多くの種類があり、一般形は次のように表される。
yn+1 = yn + hΣbiki (Σ は i = 1, 2, ... , s)
ki = f(yn + hΣaijkj, tn + cih) (Σ は j = 1, 2, ... , s) (i = 1, 2, ... , s)
s は段数(stage)とよばれ、1ステップ進むのに必要な f(y, t) の計算回数を表す(1段法や多段法の”段”とは異なる)。aij, bi および ci はパラメータで、bi および ci は次式を満たす。aij=0 (i<j) の場合、陽解法になる。また、c1 = 0 である。
ci = Σaij (Σ は j = 1, 2, ... , s) (i = 1, 2, ... , s)
Σbi = 1 (Σ は i = 1, 2, ... , s)
RK法のパラメータには次のような表記法がある (Butcher表という)。
0 |
c2 | a21
c3 | a31 a32
. | . .
cs | as1 as2 ... as,s-1
----------------------------
| b1 b2 ... bs-1 bs
RK法が p次であるとは、テイラー展開の p次の項まで正しく求めることができることを表す。s段公式が達成しうる最大の次数 p は次のようになる。
| s | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
| p | 1 | 2 | 3 | 4 | 4 | 5 | 6 | 6 | 7 | 7 |
低次の公式
s = 1 の場合、a11 = c1 = 0, b1 = 1 となり1段1次のRK法が得られるが、これはオイラー法に一致する。
s = 2 の場合、2段2次のRK法の公式を2種類作ることができる。
修正オイラー法 (modified Euler method)
0 |
1/2 | 1/2
-------------
| 0 1
計算式は次のとおり。
yn+1 = yn + hk2
k1 = f(yn, tn)
k2 = f(yn + hk1/2, tn + h/2)
ホイン法 (Heun method)
0 |
1 | 1
----------------
| 1/2 1/2
計算式は次のとおり。
yn+1 = yn + h(k1 + k2)/2
k1 = f(yn, tn)
k2 = f(yn + hk1, tn + h)
4次の公式 (ルンゲ・クッタ法)
低次の公式では精度が不十分で、解の精度を確保するためには刻み幅を小さくする必要があり計算時間がかかる。そこで、実用的によく使われているのは4段4次の公式である。
4段4次の公式を作ろうとすると、パラメータが12個で自由度が4なので様々な公式を作ることが可能である。その中で最良のものを選ぶための多くの研究が行われてきたが、いわゆるルンゲ・クッタ法 (古典的ルンゲ・クッタ法) の係数は次のとおりである。
0 |
1/2 | 1/2
1/2 | 0 1/2
1 | 0 0 1
-----------------------
| 1/6 1/3 1/3 1/6
計算式は次のとおり。
yn+1 = yn + h(k1 + 2k2 + 2k3 + k4)/6
k1 = f(yn, tn)
k2 = f(yn + hk1/2, tn + h/2)
k3 = f(yn + hk2/2, tn + h/2)
k4 = f(yn + hk3, tn + h)
単にルンゲ・クッタ法というとこの公式を指し、広く使われている。
多段法
1段法は yn+1 の計算に現在の値 yn のみを使う解法であったが、過去の値(それまでに計算済の値) yn-1, yn-2, … も使った方がよい結果が得られそうに思われる。このアイデアを実現したのが多段法である。
アダムス・バシュフォース法 (Adams-Bashforth method)
1階の常微分方程式の積分形を考える。
y(tn+1) = y(tn) + ∫ f(y(x), x) dx [tn, tn+1]
ここで、y(tn+1) と y(tn)は近似値 yn+1 と yn を使い、右辺の積分は (ti, fi) (i = n – k +1, …, n) を通る補間多項式 p(x) で f(y(x), x) を置き換えると、次のように k段の陽解法の公式が得られる。ただし、fi = f(yi, ti) を表し、yi (i = n – k +1, …, n) が計算済なので fi はすぐに求めることができる。この公式は、k次の公式になる。
yn+1 = yn + ∫ p(x) dx [tn, tn+1]
= yn + hΣγj▽jfn (Σは j = 0, ... , k-1)
ただし、
▽0fn = fn, ▽j+1fn = ▽jfn - ▽jfn-1 (後退差分)
γj = (-1)j ∫-sCj ds [0, 1]
である。なお、γ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 について具体的な式を示すと次のようになる。
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段法で、オイラー法に一致する。
アダムス・モルトン法 (Adams-Moulton method)
アダムス・バシュフォース法では補間関数を使って外挿により積分を求めたが、精度を上げるために tn+1を加えて区間 [tn-k+1, tn+1] の補間データを使う、すなわち、(ti, fi) (i = n – k +1, …, n, n+1) を通る補間多項式 p*(t) を使い、外挿することを避けたのがアダムス・モルトン法である。こちらの方が精度がよさそうであるが、その代わり右辺に 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 ∫-s+1Cj ds [0, 1]
である。なお、γ*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 について具体的な式を示すと次のようになる。
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 の場合も1段法であるが、こちらは台形則とよばれる公式に一致する。
予測子・修正子法 (predictor-corrector method)
多段法においては、陽解法は計算が容易であるが満足のいく精度が得られない。また、陰解法は精度はよいが、一般的には非線形方程式を毎回解かなければならないため計算が難しい。
そこで、広く使われている方法が予測子・修正子法 (PC法ともいう) である。まず陽解法の予測子(predictor)によりyn+1の第1近似を求め、次にその値を使って陰解法の修正子(corrector)を1回～数回適用しこれを改良する方法である。
多段法では、はじめに出発値 y1, …, yk-1 の値を求めておく必要があるが、そのために1段解法のRK法などを補助的に使わなければならない。この精度が悪いと以降の計算値に影響するため、使用する多段法の公式の精度に十分見合うだけの精度を持った方法を使用する必要がある。
予測子と修正子の組み合わせにはバリエーションがあるが、予測子にアダムス・バシュフォース法、修正子にアダムス・モルトン法を使った予測子・修正子法を、アダムス・バシュフォース・モルトン法 (あるいは単にアダムス法) とよび広く使われている。
使われなくなった解法
よく知られているが現在では使われなくなった解法についてここで触れておく。
ルンゲ・クッタ・ジル法 (Runge-Kutta-Gill method)
次のような4次のRK法である。
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
計算式は次のとおり。
yn+1 = yn + h(k1 + 2(1 - 1/√2)k2 + 2(1 + 1/√2)k3 + k4)/6
k1 = f(yn, tn)
k2 = f(yn + hk1/2, tn + h/2)
k3 = f(yn + h((√2 - 1)/2 k1 + (1 - 1/√2) k2), tn + h/2)
k4 = f(yn + h(-1/√2 k2 + (1 + 1/√2)k3), tn + h)
普通の4次のRK法では計算のために5個(yn, k1, k2, k3, k4)の変数(連立方程式の場合には配列)が必要なのに対し、この方法では工夫された手順により3個の変数だけで計算できるようになっている。また、丸め誤差を減らす工夫も盛り込まれていることが特長である。メモリが少なかった初期の計算機におけるプログラム、特に機械語プログラムに適しており、広く使われていた。
しかし、現在のコンピュータではメモリの心配がなくなり、また、倍精度で計算するのが普通になったため、欠点があるわけではないが特に優位性はなくなり、使われることが少なくなった。
ミルン法 (Milne method)
アダムス法と似ているが、積分範囲が異なる次の積分形を考える。
y(tn+1) = y(tn-1) + ∫ f(y(x), x) dx [tn-1, tn+1]
アダムス法と同様に (ti, fi) (i = n – k +1, …, n) を通る補間多項式 p(x) で積分を求めると、次の陽解法の公式が得られる。
yn+1 = yn-1 + hΣκj▽jfn (Σは j = 0, ... , k-1)
ただし、
κj = (-1)j ∫-sCj 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 ∫-s+1Cj ds [-1, 1]
予測子・修正子法の一種にミルン法とよばれる方法があり、上の陰解法の式において 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)
数値実験
ここまで説明した代表的な解法について、数値実験を行って比較してみる。
実験 1
最も簡単な微分方程式の例である。
dy/dt = y, t = 0 で y = 1
解析解は、y = et となる。t = 1 での真値 y = e と比較して誤差の絶対値を種々の公式についてプロットすると下図のようになった。横軸は刻み幅、縦軸は誤差 |yN – e| である。
オイラー法が1次、ホイン法と修正オイラー法が2次、ルンゲ・クッタ(RK)法と4次のアダムス法が4次収束であることが確かめられる。
目的にもよるが、実用的にはやはり4次のRK法程度の性能は欲しいところで、オイラー法やホイン法では不十分に感じる。
実験 2
同じ微分方程式の例で、3次～5次のRK法の収束を見てみる。
比較したのは、3次のRK法、3/8公式(4次)、4次のRK法、RKG (ルンゲ・クッタ・ジル法、4次)、5次のRK法 (Kutta-Nystrom法) である。それぞれ、次数なりの収束を示した。3種類の4次の公式はほとんどぴったり同じ結果を示した。ただし、きざみ幅が10-3程度よりも小さくなるとRKGの誤差が他よりも小さくなり最終的にはぴったり0になった。これは、RKGの丸め誤差軽減の工夫が効いているためと思われる。
実験 3
次に、同じ微分方程式の例でアダムス法を詳しく見てみる。
4次のアダムス法では、予測子に4次 (アダムス・バシュフォースの k=4) を使うが、修正子は文献により4次 (アダムス・モルトンの k=3) を使っている場合と5次 (アダムス・モルトンの k=4) を使っている場合がある。計算手順は、(1) 予測子を適用して Yn+1を求める(P)、(2) それを使って関数値 fn+1を計算する(E)、(3) 修正子を適用して Yn+1を更新する(C)、(4) 更新した Yn+1を使って fn+1を更新する(E) となる。これを PECE と表すことにする。ここで、最後の E を省略して PEC とする方法、修正子の適用を 2回繰り返して PECECE とする方法など、種々のバリエーションがある。また、出発値は4次の場合には初期値を含めて4個必要であるが、それを求める方法もいくつか考えられる。
ここでは、例えば、予測子が4次、修正子が5次、手順が PECE のとき、4-5 PECE と表すことにして、いくつかの組み合わせを比較してみた。なお、出発値は、予測子が4次の場合は4次のRK法、5次の場合は5次のRK法で求めた。
まず PEC と PECE であるが、この例は簡単な関数であるためあまり差が出なかったが、難しい関数で試してみると大きく差が出ることがあった。関数評価回数が 2倍になるが PECE の方が安全かもしれない。
次に修正子であるが、5次を使うと計算は面倒になるが関数評価回数が増えるわけではない。収束速度が図のように純粋な4次と5次の大体中間になるので5次の修正子を使う価値があると思われる。
出発値については、今回はRK法だけしか実験していないが、次数を合わせておいた方がよいことが確かめられた。すなわち、4次のRK法で求めた出発値を5次のアダムス法で使うと後々まで影響がおよび収束が遅くなった。また、5次のRK法で求めた出発値を4次のアダムス法で使っても、きざみ幅を小さくしていったときに途中から失速して最終的にきざみ幅の小さいところでは収束が速くなるわけではなかった。
参考文献
[1] E. Hairer、他「常微分方程式の数値解法 I」(2007) シュプリンガー・ジャパン[2] 三井斌友「岩波講座 応用数学13 [方法3] 微分方程式の数値解法 I」(1993) 岩波書店
[3] 名取亮「数値解析とその応用」(1990) コロナ社
[4] 森口繁一「数値計算工学」(1989) 岩波書店
