10. 常微分方程式 (3) 解法の安定性, スティフな方程式, 微分代数方程式

解法の安定性

初期値問題を解いている間に数値解が異常な値になったり発散したりすることがある. その原因は, 微分方程式自体の問題である場合と, 数値解法の問題である場合がある.

微分方程式自体の安定性

「誤差の推定」の説明図に示されるように, 微分方程式の解曲線が時間が進むにつれて離れていく場合, 初期値のわずかな差が大きく拡大する現象が生じることがある. これは不安定な微分方程式の例である.

逆に, 微分方程式の解曲線が時間が進むにつれて近づいていく場合, むしろ誤差が縮小していくので安定な微分方程式の例である.

繰り返すと, 関数がどちらの傾向を示すかは微分方程式の ∂f/∂y の符号により決まる.

数値解法の安定性

安定な微分方程式であっても, 使用する数値解法によって不安定現象を示すことがある. これは数値解法の安定性の問題である.

ステップ幅の制約

いま, 微分方程式の導関数 ∂f/∂y が一定値 λ で, λ < 0 とする(これは安定な微分方程式の初期値問題である). また, 各ステップで加えられる誤差はステップによらず一定であるとする. その場合, 数値解法の安定性は次の線形テスト問題を調べることでわかる.

  
  dy/dt = λy (t > 0, λ(定数) < 0), y(0) = y0

 
この問題の解は,

  
  y(t) = y0exp(λt)

 
となるから, t → ∞ のとき y(t) → 0 となる.

(1) オイラー法

上のテスト問題をオイラー法で解くと,

  
  yn = y0(1 + hλ)n

 
となるから, n → ∞ のとき y_n → 0 とならなければならない. そのためには,

  
  |1 + hλ| < 1

 
でなければならない. これを安定条件という.

λを複素数として z = hλ とおき, 安定条件 |R(z)| < 1 (ただし, R(z) = 1 + z) を複素平面上に図示すると次のようになる.

色付き部分 { z ∈ C; |R(z)| < 1 } を絶対安定領域という. λが実数であれば -2 < z = hλ < 0 である. 区間 [-2, 0] を絶対安定区間という. このとき, ステップ幅は 0 < h < -2/λ としなければならない.

(2) 後退オイラー法

後退オイラー法ではR(z)は次のようになる.

  
  R(z) = 1/(1 - z)

 
絶対安定領域を図示すると次のようになる.

このように, 左半平面が絶対安定領域に含まれるときにA安定であるといい, ステップ幅 h に対する制限がなくなる.

(3) ルンゲ・クッタ法

s段p次のRK法では R(z) は次のようになる.

  
  R(z) = 1 + z + z2/2! + ... + zp/p! + γp+1zp+1 + ... + γszs

 
ここで, γ_p+1, … はRK法のパラメータにより決まる.

p = s = 2, 3, 4 について絶対安定領域を図示すると次のようになる.

次数が上がると安定領域も広くなることがわかる. 絶対安定区間はそれぞれ, [-2, 0], [-2.513, 0], [-2.785, 0] である. ルンゲ・クッタ法の陽解法は A安定ではない.

(4) 多段法 (予測子・修正子法)

多段法も陽解法は A安定ではない.

多段法の陰解法では A安定な解法は2次以下である (A安定な多段法で最も高次なものは台形則である).

(5) 台形則 (陰解法, アダムス・モルトン法の k=1)

台形則ではR(z)は次のようになる.

  
  R(z) = (1 + z/2)/(1 - z/2)

 
この絶対安定領域は左半平面全部となり, 台形則はA安定である.

ステップ幅の制約 (連立方程式の場合)

連立方程式の場合, 線形テスト問題で y をベクトルに拡張して扱う.

  
  dy/dt = Jy, y(0) = y0

 
J は y の元数と同じ次元数の正方行列となる. 安定な微分方程式の初期値問題であるための条件は, J のすべての固有値 λ_i の実部が負であることである.

  
  Re(λi) < 0 (すべての i について)

 
また, z = hλ_i として, 複素平面上で z が絶対安定領域に含まれていればその h について安定である.

数値実験 (20)

  
  y' = -10y, y(0) = 1

注 – 「使われなくなった解法」で説明したミルン法の不安定性はここで説明した安定性とは別の問題で, ステップ幅を小さくしても不安定性はおさまらない.

スティフな方程式

スティフな方程式の例

スティフな方程式(硬い方程式ともいう)は化学反応や電気回路などの問題で現れる. 非常に急な変化をする解を持ち, 通常の解法では解きにくいものをいう.

次の初期値問題はスティフな方程式の例である.

  
  dy/dt = ( -2    1    ) y + ( -cos(t)             )
          ( 1998 -1999 )     ( 1999cos(t) - sin(t) )

  y(0) = ( 1 )
         ( 2 )

 
一般解は次のとおりである.

  
  y(t) = C1exp(-t) ( 1 ) + C2exp(-2000t) (   1   ) + (   0    )
                   ( 1 )                 ( -1998 )   ( cos(t) )

 
上の方程式の y の係数行列の固有値は実数で, -1 および -2000 である. 従って, 「解法の安定性」の議論から, オイラー法ならばステップ幅 h < 0.001, 4次のルンゲ・クッタ法ならば h < 0.00139 としなければ解くことができないはずである. 実際に試してみると, h をこれらの制約条件よりも大きくすると発散して解くことができないことが確かめられた.

XLPackに収録されている通常のステップ幅自動調節機能付きのルーチン Derkf, Dverk, Deabm で同じ例題を Tol = 0.01 として解いてみた. それぞれ非常に小さなステップ幅が使われながらも解が求められた. 後で示すBDF法に比べると100倍くらいの関数評価回数を要した.

スティフな方程式の定義

スティフな方程式は, 上の例のように安定性の問題により通常の解法ではステップ幅に制約があるために, 計算に非常に時間がかかったり, 条件によっては振動や発散をしたりするものである. 上の例のように, 一見急な変化をするようには見えないものであっても, 一般解に指数項を含んでいたりするとスティフな方程式となる.

次のような連立常微分方程式を考える.

  
  dy/dt = Ay + φ(t)
  y(0) = y0

 
行列Aの固有値λ_iはすべて異なりその実数部は負であるとする.

硬度比(stiffness ratio)を次のように, 実数部が最大の固有値λ_maxと実数部が最小の固有値λ_minの実数部の値の比と定義する.

  
  |λmaxの実数部|/|λminの実数部|

 
この硬度比が大きい(> 10⁴ くらい)ときにスティフな方程式ということが多いようである. 上の例はほぼスティフな方程式といえる.

しかし, 硬度比が小さくてもスティフなケースや, 硬度比が大きくても初期値によってはスティフではなく解けるケースがあり, 文献[1]では次のような定義も与えている.

非斉次項の変動量M(t)を次のように定義する.

  
  M(t) = ||φ'(t)||/||φ(t)||,  ただし φ(t) = 0 のときは M(t) = 0 とする

 
次に減衰比(damping ratio)を次のように定義する.

  
  |λmaxの実数部|/M(t)

 
このとき次のように定義する.

– 減衰比が大きな問題を, そのtの近傍で硬い系という.
– 硬度比が大きいが, 減衰比がそれよりずっと小さい問題は疑似的に硬い系という.

上の例は硬度比は2000で, 減衰比はtにより変動し1～約1999²となっている.

数値実験 (21)

スティフな方程式の解法

スティフな方程式を解くためにはステップ幅に制約がない解法, すなわち, A安定な解法を使う必要がある.

A安定な解法としては後退オイラー法や陰的台形則があるが, 実用的にはもっと高次の公式が必要であり種々の解法が研究されてきた. 以下, 代表的な解法として, 後退微分公式(BDF)および陰的ルンゲ・クッタ(IRK)法について説明する. また, ローゼンブロック法と補外法についても少し触れる.

後退微分公式(BDF)

アダムス・モルトン法では (t_i, f_i) (i = n – k +1, …, n+1) を補間データとして f を補間したが, y そのものを補間する方法がある. すなわち, (t_i, y_i) (i = n – k +1, …, n+1) を補間データとして補間多項式 q(t) を求め, それが少なくても格子点の1か所(i=n+1-r の点)で微分方程式を満たすようにする.

  
  q'(tn+1-r) = f(yn+1-r, tn+1-r)

 
ここで, r = 1 とすると陽解法になり, k = 1 と k = 2 の場合はそれぞれオイラー法と中点則になる. しかし, k ≧ 3 の場合には不安定になり微分方程式の解法として使うことはできない.

r = 0 とすると陰解法になるが k ≦ 6 で安定な解法が得られる. これは後退微分公式 (Backward differentiation formula (BDF) または ギヤ法(Gear method)) とよばれ, k 次の公式になる.

  
  Σ (1/j)▽jyn+1 = hfn+1  (Σは j = 1, ... , k)

 
ただし,

  
  ▽0fn = fn,  ▽j+1fn = ▽jfn - ▽jfn-1  (後退差分)

 
k = 1～6 について, y_n+1, y_n, y_n-1, … について整理すると次のようになる.

  
  k = 1 : yn+1 - yn = hfn+1  (後退オイラー法に一致)
  k = 2 : (3/2)yn+1 - 2yn + (1/2)yn-1 = hfn+1
  k = 3 : (11/6)yn+1 - 3yn + (3/2)yn-1 - (1/3)yn-2 = hfn+1
  k = 4 : (25/12)yn+1 - 4yn + 3yn-1 - (4/3)yn-2 + (1/4)yn-3 = hfn+1
  k = 5 : (137/60)yn+1 - 5yn + 5yn-1 - (10/3)yn-2 + (5/4)yn-3 - (1/5)yn-4 = hfn+1
  k = 6 : (147/60)yn+1 - 6yn + (15/2)yn-1 - (20/3)yn-2 + (15/4)yn-3 - (6/5)yn-4 + (1/6)yn-5 = hfn+1

 
BDFの絶対安定領域を k = 1～4 について下図に示す.

図の曲線の外側が安定領域である. k = 1 (後退オイラー法) と k = 2 は A安定であるが, それ以外では少し左半平面にはみだしている部分があり A安定ではない. k が大きくなるほど安定領域は狭くなり, k ≧ 7 では原点でも不安定になる.

スティフな方程式を解く際に, 実際には多くの場合に完全にA安定である必要はない. そこで, 条件を緩めた代わりに高次な解法としてBDFがよく使われている. 4～6次のBDFはA安定の条件を緩めたA(α)安定である.

原点を頂点として, 左側に非有界であり, 実軸に対して上にα°下にα°の扇型の領域が絶対安定領域に含まれる場合, A(α)安定であるという. A(90°)安定はA安定に等しい. 方程式の係数行列の固有値がすべてこの中に入っているならば, 実質的に A安定と同じである. 3次のBDFはA(86°)安定, 4次のBDFはA(73°)安定, 5次のBDFはA(52°)安定, 6次のBDFはA(18°)安定である.

左半平面がほとんど絶対安定領域に含まれるが, 原点付近を含まない虚軸近傍だけが絶対安定領域に含まれない場合を硬安定という. 虚軸に近いところに係数行列の固有値ない限り, これも実質的に A安定と同じである. 4～6次のBDFは硬安定である.

数値実験 (22)

陰的ルンゲ・クッタ法

s段ルンゲ・クッタ法の一般形は次のようなものであった.

  
  ki = f(tn + cih, yn + hΣaijkj) (Σ は j = 1, 2, ... , s) (i = 1, 2, ... , s)
  yn+1 = yn + hΣbiki (Σ は i = 1, 2, ... , s)

 
a_ij, b_i および c_i はパラメータで, b_i および c_i は次式を満たす.

  
  ci = Σaij (Σ は j = 1, 2, ... , s) (i = 1, 2, ... , s)
  Σbi = 1 (Σ は i = 1, 2, ... , s)

 
a_ij = 0 (i ≦ j) の場合, 陽解法になる. それ以外の場合, 陰的(implicit)ルンゲ・クッタ(IRK)法という. IRK法は A安定なものを作れるが, 1 ステップごとに連立非線形方程式を解かなければならないため計算の手間がかかるという問題がある.

ガウス法

ガウス積分則に基づく最大次数のIRK法である. s段ガウス法の次数は2s次である. c₁, c₂, …, c_s は区間[0, 1]に変数変換されたルジャンドル多項式

  
  ds/dxs (xs(x - 1)s)

 
のゼロ点で与えられる.

クンツマン-ブッチャー法とよばれる4次と6次の公式の係数を以下に示す. この方法はA安定である.

2段4次クンツマン-ブッチャー法

  
  1/2-√3/6 |    1/4     1/4-√3/6
  1/2+√3/6 | 1/4+√3/6      1/4
  ----------------------------------
            |    1/2        1/2
    

 
3段6次クンツマン-ブッチャー法

  
  1/2-√15/10 |    5/36      2/9-√15/15  5/36-√15/30
       1/2    | 5/36+√15/24      2/9     5/36-√15/24
  1/2+√15/10 | 5/36+√15/30  2/9+√15/15     5/36
  ------------------------------------------------------
              |    5/18          4/9        5/18
    

ラダウ法とロバット法

ラダウとロバットの積分則に基づくIRK法である. c₁, c₂, …, c_s は次の多項式のゼロ点で与えられる.

  
  ds-1/dxs-1 (xs(x - 1)s-1)     (Ⅰ型: ラダウ右)
  ds-1/dxs-1 (xs-1(x - 1)s)     (Ⅱ型: ラダウ左)
  ds-2/dxs-2 (xs-1(x - 1)s-1)   (Ⅲ型: ロバット)

 
Ⅰ型のラダウⅠA法(c₁ = 0 になる), Ⅱ型のラダウⅡA法(c_s = 1 になる)の次数はs段のとき2s-1次である. 両方ともA安定である. 1段ラダウⅡA法は陰的オイラー法になる.

2段3次ラダウⅡA法

  
  1/3  |  5/12  -1/12
   1   |  3/4    1/4
  ---------------------
       |  3/4    1/4

 
3段5次ラダウⅡA法

  
  (4-√6)/10 |   (88-7√6)/360    (296-169√6)/1800  (-2+3√6)/225
  (4+√6)/10 | (296+169√6)/1800    (88+7√6)/360    (-2-3√6)/225
       1     |    (16-√6)/36        (16+√6)/36           1/9
  ------------------------------------------------------------------
             |    (16-√6)/36        (16+√6)/36           1/9
    

 
Ⅲ型の方法はロバットⅢA, ⅢB, ⅢC の３つがある. s段のとき2s-2次となる. いずれもA安定である.
以下にロバットⅢA法の例を示す。

3段4次ロバットⅢA法

  
     0   |    0     0      0
    1/2  |  5/24   1/3   -1/24
     1   |   1/6   2/3    1/6
  ------------------------------
         |   1/6   2/3    1/6

 
4段6次ロバットⅢA法

  
       0     |       0             0              0             0
  (5-√5)/10 | (11+√5)/120   (25-√5)/120   (25-13√5)/120  (-1+√5)/120
  (5+√5)/10 | (11-√5)/120  (25+13√5)/120   (25+√5)/120   (-1-√5)/120
       1     |     1/12          5/12           5/12          1/12
  -------------------------------------------------------------------------
             |     1/12          5/12           5/12          1/12

数値実験 (23)

対角陰的ルンゲ・クッタ法とローゼンブロック法

IRK法ではステップごとにn*s元非線形連立方程式を解かなければならないが, 場合によっては大きな方程式になってしまうのでこの負担を軽減する工夫が検討されている.

a_ij = 0 (i < j) で a_ii ≠ 0 (少なくても 1つの i について) としたものを対角陰的(diagonally implicit)ルンゲ・クッタ(DIRK)法という. 半陰的(semi-implicit)ルンゲ・クッタ法ともいう. さらに, すべての対角成分が同一の値である場合, 単純対角陰的(singly diagonally implicit)ルンゲ・クッタ(SDIRK)法という.

DIRK法の場合, s個のn元非線形方程式を逐次解くことにより解を求めることができるようになる. n元非線形方程式をニュートン法で解くときa_iiを含む係数の線形連立方程式を解く必要があるが, SDIRK法であればLU分解を1回で済ませて結果を繰り返し使うことができるようになる.

DIRK法は次のように表される.

  
  ki = f(yn + hΣaijkj + aiiki) (Σ は j = 1, 2, ... , i-1) (i = 1, 2, ... , s)
  yn+1 = yn + hΣbiki (Σ は i = 1, 2, ... , s)

 
工夫を進めて次のように変えたのがローゼンブロック法(Rosenbrock method)である.

  
  (I - hγiiJ)ki = f(yn + hΣαijkj) + hJΣγijkj  (Σは j = 1, 2, ... , i-1) (i = 1, 2, ... , s)
  yn+1 = yn + hΣbiki  (Σ は i = 1, 2, ... , s)

 
α_ij, γ_ij, b_iは係数, J = f'(y_n) である. DIRK法を線形化したものとみることができるので線形陰的(linearly implicit)ルンゲ・クッタ法とよばれる. 各ステップではk_iを未知数とし(I – hγ_iiJ)を係数とする線形連立方程式を解くことで計算を進められるようになる.

これに属する種々の公式が提案されており, 公式によりA安定あるいはA(α)安定である.

補外法

補外法をスティフな方程式に適用するためにGBSアルゴリズムの線形陰的な(ローゼンブロックタイプの)拡張が研究されている.

  
  (I - hjJ)(y1 - y0) = hjf(t0, y0)  [線形陰的オイラー法]
  (I - hjJ)(yi+1 - yi) = -(I + hjJ)(yi - yi-1) + 2hjf(ti, yi)  [線形陰的中点則]

  Tj1 = (1/2)(ynj-1 + ynj+1)

 
Jはヤコビ行列の近似で, J = 0 とおくとGBSアルゴリズムに一致する.

線形陰的中点則や線形陰的オイラー法を使った補外法もいくつか実装されており, αが90°に近いA(α)安定である.

微分代数方程式

次のように常微分方程式と代数方程式が連立した方程式を微分代数方程式(Differential algebraic equations (DAE))という.

  
  x' = f(t, x, z)
  0 = g(t, x, z)

 
代数方程式は制約条件を表し, この制約のもとで常微分方程式を解くとみることができる. 代数制約式g()を解析的に微分と消去を繰り返し z’ = h(t, x, z) のような形に変換すれば, すべての未知数について陽的な連立常微分方程式が得られる(方程式が特異でなければ可能である). この変換に必要な微分回数を微分代数方程式の指数という.

単振り子の例

下図のように, 長さl[m]の糸に重さm[kg]のおもりがついた振り子を考える. 糸と鉛直線の角度がθ[rad]で重力加速度をg[m/s²]とする.

この振り子の運動はθを使って次の2階常微分方程式で表される.

  
  θ'' = -(g/l)sin(θ), θ(0) = θ0

 
θ₀で手を離すと振り子は周期運動を開始し, その周期Tは次のようになることがわかっている.

  
  T = 4sqrt(l/g)K(sin(θ0/2))

 
ただし, Kは第1種完全楕円積分である.

この問題は1階連立形に変換して解くことができて, 初期値を θ(0) = π/4 とし, Derkfを使って t = 0.1 ごとにプロットすると次のようになった. なお, 周期Tごとの値を赤で示した.

θ₀ = π/4, l = 1 として2周期分の振り子の動きを示すと次のようになる.

さて, 同じ問題を x₁ = l・sin(θ), x₂ = -l・cos(θ) として直交座標系で表すと次の方程式になる. zは未定乗数である.

  
  mx1'' = -zx1
  mx2'' = -zx2 - mg
  x12 + x22 = l2

 
この方程式はθを使って変数変換しzを消去すれば上のθを使った方程式に戻すことができるが, ここではこのまま計算を行うことにする.

q₁ = x₁, q₂ = x₂, v₁ = x₁‘, v₂ = x₂‘, λ = z/m と変数変換して1階連立形に書き直すと次のようになる.

  
  q1' = v1
  q2' = v2
  v1' = -λq1
  v2' = -λq2 - g
  0 = q12 + q22 - l2

 
これは微分代数方程式の形になっている. 最後の式が代数制約式であり, 位置の制約条件を表す(物理的にはおもりの軌道が半径lの円周上にあることを示す). これを1回微分すると次のようになる.

  
  q1q1' + q2q2' =
  q1v1 + q2v2 = 0

 
これは速度の制約条件を表す. もう1回微分すると次のようになる.

  
  q1v1' + q2v2' + v12 + v22 =
  -λ - gq2 + v12 + v22 = 0

 
これは加速度の制約条件を表す. これをもう1回微分するとλ’の式が得られるのでこの方程式の指数は3である.

微分代数方程式の解法 (1)

微分代数方程式は指数が高いほど解きにくいので, 制約式の解析的微分を繰り返し指数を下げて解くのが1つの方法である.

単振り子の制約式を2回微分して得られたλを元の連立方程式に代入すると次の4元連立常微分方程式(指数0)が得られる.

  
  q1' = v1
  q2' = v2
  v1' = (gq2 - v12 - v22)q1
  v2' = (gq2 - v12 - v22)q2 - g

 
この連立常微分方程式は普通に解くことができて, 初期値を q₁(0) = l・sin(π/4), q₂(0) = -l・cos(π/4) とし, Derkfを使って計算してみるとθを使った方程式と同じ解が得られた.

微分代数方程式の制約式を微分して得られた常微分方程式は不変性を持つ常微分方程式といわれ, 制約条件が陽には示されていないが連立常微分方程式の中に含まれている. このような方程式には安定性の問題があり, 計算を進めていったときに誤差の累積により制約条件が満たされなくなっていくことがある.

途中で誤差が蓄積してきたと仮定してq₂の初期値を -l*cos(π/4) + 0.05 としてDerkfで計算してみる(本来このような初期値は許されないが). 下図のように制約条件から外れた軌跡を描くようになるのがわかる. x-y平面での振り子の軌跡を表す. 赤線は制約条件を満たす本来の軌跡である.

これを安定化する方法の一例として, X’ = f(t, x) の代わりに次の方程式を解く方法がある. h(x)は不変方程式で, h(x) = 0 のときに制約条件を満たす. 従って, 制約条件を満たしていればこの方程式は元の方程式と同じ解を持つ.

  
  x' = f(t, x) - γF(x)h(x)

 
詳細は省略するが, 今考えている単振り子の問題の場合, 具体的には次のようにする(文献[4]参照).

  
  h(x) = ( q12 + q22 - l2 )
         ( q1v1 + q2v2    )

  F(x) = D(HD)-1
  DT = H = ( 2q1  2q2  0   0  )
           ( v1   v2   q1  q2 )

 
γ = 10 として上図と同様に計算してみると次のようになった.

制約条件を満たす方向に安定化しているのがわかる.

このように, 微分代数方程式は, 問題の定式化にいくつかのやり方がある場合には指数の小さなものを選び, 指数が大きければ解析的に微分を繰り返すことにより指数を下げることにより常微分方程式に変換して解くのがよい. なお, 問題によっては制約条件を満たすように安定化を図る方がよい.

微分代数方程式の解法 (2)

微分代数方程式を直接離散化すると非常にスティフな方程式になる. しかし, RADAU5(Radau IIA法を使用)やDASSL(BDFを使用)のように微分代数方程式を直接離散化して解くことができるプログラムもできている. 指数を下げるなどの再定式化には手間がかかるので, これを使って解けるのであればその方がよい. しかし, 多くのプログラムは指数1にしか使えないなど, この場合でも指数が小さくなるように定式化しておくのがよいことは変わらない.

XLPackには微分代数方程式に使えるいくつかのプログラムが収録されている. 計算例は実用ルーチンのベンチマーク編を参照せよ.