硬い方程式

数学において硬い方程式（英: 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$ を満足する方程式

$y'=f(t,y)$

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

上記の方程式に対してルンゲ＝クッタ法の安定性は 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-安定である。

参考文献

• 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 .
• 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