硬い方程式

数学・数値解析において硬い方程式（英: stiff equation）は、常微分方程式の数値解法・偏微分方程式の数値解法において、刻み幅を極めて小さくしない限り、数値的不安定になる微分方程式である。硬さを的確に定義するのが困難であると判明したが、方程式に解の急激な変化を起こせる項が含まれていることは確かである。

導入の例

硬い方程式に対する異なる数値解の安定性

下記の初期値問題を考える。

$y'=-15y,\quad t\geq 0,\quad y(0)=1.$

この問題は直接に解くことができ、厳密解 (水色の曲線) が次の公式で与えられる。

$y(t)=e^{-15t}.$

公式によって、$\lim _{t\to \infty }y(t)=0$ も明らかである。

同じ振舞いを持つ数値解を求めよう。様々な数値的方法を用いて得られる数値解は右側の画像に表示される。

1. 刻み幅 h = 1/4 のオイラー法 に対応する解（赤色の曲線）は激しく振動し、素早くグラフの範囲を超える。
2. 刻み幅 h = 1/8 のオイラー法 に対応する解（緑色の曲線）も振動するが、グラフの範囲にいる。
3. 刻み幅 h = 1/8 の台形公式に対応する解（青色の曲線、凡例ではAdams-Moulton法と名付けられたが同じ方法である）は振動せず、期待通りに0に減衰していく。

よって、オイラー法は上記の硬い方程式に対し数値的不安定である。一方、台形公式は数値的安定である。

他の例として、もっとも有名な硬い方程式の一つは、Robertsonの化学反応を支配する方程式系である。

${\dot {x}}=-0.04x+10^{4}y\cdot z,$

${\dot {y}}=0.04x-10^{4}y\cdot z-3\cdot 10^{7}y^{2},$

${\dot {z}}=3\cdot 10^{7}y^{2}.$

[0, 40] のような短い区間では、上記の方程式系を数値的に積分することに問題はない。しかし区間が極めて大きい場合（例えば 10¹¹）、多数のコードは方程式系を正しく積分することができなくなる。

常微分方程式

→「ルンゲ＝クッタ法 § 安定性」、「線型多段法」、および「常微分方程式の数値解法」も参照

上述の例の示すように、硬い常微分方程式の近似解を計算するとき、数値的に安定な方法を使うべきである。常微分方程式における数値的安定性に複数の定義が存在する。特に、線型方程式に対する安定性と非線型方程式に対する安定性を分けて考える必要がある。

硬さの比例

線形常微分方程式系の硬さは簡単に測ることができる。一般的な線型方程式系

$y'=Ay,\quad y(t_{0})=y_{0}$

を考える。上記方程式に対する 硬さの比例 (stiffness ratio) は、行列 A の最大固有値（の絶対値）を最小固有値（の絶対値）で割った商である。つまり、A の固有値を λ₁ ≥ λ₂ ≥ ... ≥ λ_n としたとき、方程式に対する硬さの比例を ‖ λ₁ ‖ / ‖ λ_n ‖ と定義する。

典型的な硬さの比例は、10¹⁷ あたりである。極端な場合に、その数は 10³¹ にもなる。

非線型方程式の場合は、代わりに関数のヤコビ行列の固有値を使って比例を同じ公式で計算する。

線型安定性

線型常微分方程式に対する安定性は 線型安定性 (linear stability)、あるいは 絶対安定性 (absolute stability) という。線型テスト方程式

$y'=\lambda y,\quad t\geq 0,\quad y(0)=1,\quad \lambda \in \mathbb {C}$

を考える。この方程式は簡単に解くことができ、厳密解は y(t) = e^λt である。Re λ < 0 が成立するとき、y の t → ∞ の極限も 0 である。理想的に、近似解にもそのような振舞いを期待できる。しかし刻み幅 h が一定のとき、すべての方法に対する近似解がそのような振舞いを持つわけではない。それを区別するのが線型安定性である。

一つの方法による時刻 t_n での近似解を y_n とする。複素数平面上の集合

$D=\{h\lambda \in \mathbb {C} \mid \lim _{n\to \infty }y_{n}=0\}$

は方法に対する 線型安定性領域 (linear stability domain)、あるいは 絶対安定性領域 (region of absolute stability) という。この集合はすなわち、与えられた方法による近似解が期待通りの振舞いを持つすべての hλ （からなる集合）である。特に、ルンゲ＝クッタ法に対する線型安定性領域は以下の形で与えられる。

$D=\{z\in \mathbb {C} \mid |r(z)|<1\}\quad (z=h\lambda ).$

ここで、r(z) は等式 y_n = (r(z))ⁿ を成立させる関数であり、時々方法に対する 安定性関数 という。例えば、オイラー法に対応する関数は r(z) = 1 + z である。

一般的に、方法に対する安定性領域（の面積）が大きいほど、その方法はより安定である。よってもっとも安定な方法に対する安定性領域は左複素数平面すべてを含めるべきである。そのような方法を A-安定 (A-stable) という。A-安定な方法は（すくなくとも線型）硬い方程式の場合でも、刻み幅 h を精度のみの考慮で選択することができ、よって硬い方程式を解くために適切な方法だと考えられる。しかし、優れる安定性を持つ方法を実装するには通常高い計算コストが所要される。そのため、実践では常にA-安定な方法を使うわけではなく、方程式の性質、精度の要件や計算コストの制限などの条件を共に考えてから適切な方法を選ぶのが必要となる。

非線型安定性

上述の安定性理論に考察されたのは線型方程式のみである。その理論は時折り非線型方程式にも適用できるが、決して正しいわけではない。非線形方程式の研究を完全に一般化するのが困難であるように、すべての方程式に対する安定性を考察するのもほぼ不可能である。現在非線形方程式に対する安定性はほとんど単調性条件 $\langle f(t,y)-f(t,z),y-z\rangle \leq 0$ を満足する方程式

$y'=f(t,y)$

のみを考える。ただし、$\langle \cdot ,\cdot \rangle$ は標準内積である。この発想は、ダールキストによるものである。また、ルンゲ＝クッタ法と線型多段法に対する安定性の定義は異なる。なぜならば、線型多段法は時刻毎に多数の成分からベクトルを記憶する必要があり、偏差を測るには標準内積と異なる内積を定義しなければならないからである（比べて、ルンゲ＝クッタ法は時刻毎に単一の実数を記憶し、よって標準内積が使われる）。

上記の方程式に対してルンゲ＝クッタ法の安定性は B-安定性 (B-stability) という。方程式にルンゲ＝クッタ法を適用するときに、異なる初期値 y₀ と ˆy₀ に対し不等式

$\lVert y_{1}-{\hat {y}}_{1}\rVert \leq \lVert y_{0}-{\hat {y}}_{0}\rVert$

が成立すれば、その方法をB-安定と呼ぶ。ここで、y₁ と ˆy₁ は時刻 t₁ でのそれぞれの近似解である。B-安定な方法は必ずA-安定である。

さらに、ルンゲ＝クッタ法の係数が b_i ≥ 0 かつ行列

$M=(b_{i}a_{ij}-b_{j}a_{ji}-b_{i}b_{j})_{ij}$

が半正定値であるという条件を満足するとき、その方法を 代数的安定 (algebraically stable) という。代数的安定な方法は必ずB-安定である。

線型多段法の安定性は G-安定性 (G-stability) という。G-安定性は、B-安定性と同じアイディアを持つが、上述通り標準内積が通用しないので同じように定義することができない。

k 段線型多段法の一般形式は次の公式で与えられる。

$\sum _{i=0}^{k}\alpha _{i}y_{n-i}=h\sum _{i=0}^{k}\beta _{i}f(t_{n-i},y_{n-i}).$

ここで、α_i と β_i は定数であり、ベクトル α = (α₀, ..., α_k) と β = (β₀, ..., β_k) は方法の 生成ペア (generating pair) という。 この方法に対応する one-leg法 (one-leg method) は次の公式で与えられる。

$\sum _{i=0}^{k}\alpha _{i}y_{n-i}=h{\Biggl (}\sum _{i=0}^{k}\beta _{i}{\Biggr )}\,f{\Biggl (}x_{n}-\theta h,{\biggl (}\sum _{i=0}^{k}\beta _{i}{\biggr )}^{-1}\sum _{i=0}^{k}\beta _{i}y_{n-i}{\Biggr )}$

ただし、

$\theta ={\frac {\sum _{i=0}^{k}i\beta _{i}}{\sum _{i=0}^{k}\beta _{i}}}$

である。線型多段法は対応するone-leg法と同じ線型安定性を持つため、同じ非線形安定性を持つことも期待できる。しかし上記の方程式に対する安定性を分析するには、線型多段法より対応するone-leg法を用いる方が遥かに簡単である。よって以下の定義はone-leg法に対するものである。

与えられた k 次正定値対称行列 G に対応する $\mathbb {R} ^{kn}$ 上の内積を以下のように定義できる：

$\langle U,V\rangle _{G}=\sum _{i=1}^{k}\sum _{j=1}^{k}g_{ij}\langle U_{i},V_{j}\rangle ,\quad U,V\in \mathbb {R} ^{kn}$

ここで、$U=(U_{1},\ldots ,U_{k})^{\mathrm {T} },\;U_{i}\in \mathbb {R} ^{n}$である。

one-leg法に対し、行列

$M=\alpha \beta ^{\mathrm {T} }+\beta \alpha ^{\mathrm {T} }-{\begin{bmatrix}G&0\\0&0\end{bmatrix}}+{\begin{bmatrix}0&0\\0&G\end{bmatrix}}$

を半正定値にする G が存在するとき、その方法をG-安定という。この定義は初見でルンゲ＝クッタ法の安定性とは全く違うものと思われるかもしれないが、本質的には同じものである。なぜならば、上記の定義を以て、不等式

$\lVert y_{1}-{\hat {y}}_{1}\rVert _{G}\leq \lVert y_{0}-{\hat {y}}_{0}\rVert _{G}$

を証明できるからである。

また、G-安定性もB-安定性のようにA-安定性より強い条件に見えるかもしれないが、実際にA-安定性とは同値である。すなわち、線型多段法がA-安定であることは対応するone-leg法がG-安定であることの必要十分条件となる。

参考文献

• Burden, Richard L.; Faires, J. Douglas (1993), Numerical Analysis (5th ed.), Boston: Prindle, Weber and Schmidt, ISBN 0-534-93219-3 .
• Butcher, John C. (2008), Numerical Methods for Ordinary Differential Equations, New York: John Wiley & Sons, ISBN 978-0-470-72335-7 .
• Dahlquist, Germund (1976), Error analysis for a class of methods for stiff non-linear initial value problems, Lecture Notes in Mathematics, 506, pp. 60-72, doi:10.1007/BFb0080115 .
• Dahlquist, Germund (1978), “G-stability is equivalent to A-stability”, BIT 18 (4): 384-401 .
• Hairer, Ernst; Wanner, Gerhard (1996), Solving ordinary differential equations II: Stiff and differential-algebraic problems (second ed.), Berlin: Springer-Verlag, ISBN 978-3-540-60452-5 .
• Iserles, Arieh (2008), A First Course in the Numerical Analysis of Differential Equations (Second Edition), Cambridge University Press, ISBN 978-0-521-73490-5 .
• Iserles, Arieh; Nørsett, Syvert (1991), Order Stars, Chapman & Hall, ISBN 978-0-412-35260-7 .

外部リンク

• An Introduction to Physically Based Modeling: Energy Functions and Stiffness
• Stiff systems Lawrence F. Shampine and Skip Thompson Scholarpedia, 2(3):2855. doi:10.4249/scholarpedia.2855