硬い方程式

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

導入の例

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

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

$y^{\prime }=-15y,t\ge 0,y(0)=1.$

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

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

公式によって、$\underset{t\to \mathrm{\infty }}{\lim }y(t)=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.$

テンプレート:Math のような短い区間では、上記の方程式系を数値的に積分することに問題はない。しかし区間が極めて大きい場合（例えば テンプレート:Math）、多数のコードは方程式系を正しく積分することができなくなる。

常微分方程式

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

硬さの比例

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

$y^{\prime }=Ay,y(t_0)=y_0$

を考える。上記方程式に対する 硬さの比例 (stiffness ratio) は、行列 テンプレート:Mvar の最大固有値（の絶対値）を最小固有値（の絶対値）で割った商である。つまり、テンプレート:Mvar の固有値を テンプレート:Math としたとき、方程式に対する硬さの比例を テンプレート:Math と定義する。

典型的な硬さの比例は、テンプレート:Math あたりである。極端な場合に、その数は テンプレート:Math にもなるテンプレート:Sfn。

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

線型安定性

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

$y^{\prime }=\lambda y,t\ge 0,y(0)=1,\lambda \in ℂ$

を考える。この方程式は簡単に解くことができ、厳密解は テンプレート:Math である。テンプレート:Math が成立するとき、テンプレート:Mvar の テンプレート:Math の極限も テンプレート:Math である。理想的に、近似解にもそのような振舞いを期待できる。しかし刻み幅 テンプレート:Mvar が一定のとき、すべての方法に対する近似解がそのような振舞いを持つわけではない。それを区別するのが線型安定性である。

一つの方法による時刻 テンプレート:Mvar での近似解を テンプレート:Mvar とする。複素数平面上の集合

$D=\{h\lambda \in ℂ\mid \underset{n\to \mathrm{\infty }}{\lim }{y}_n=0\}$

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

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

ここで、テンプレート:Math は等式 テンプレート:Math を成立させる関数であり、時々方法に対する 安定性関数 という。例えば、オイラー法に対応する関数は テンプレート:Math である。

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

非線型安定性

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

$y^{\prime }=f(t,y)$

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

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

$\Vert y_1-{\hat{y}}_1\Vert \le \Vert y_0-{\hat{y}}_0\Vert$

が成立すれば、その方法をB-安定と呼ぶ。ここで、テンプレート:Math と テンプレート:Math は時刻 テンプレート:Math でのそれぞれの近似解である。B-安定な方法は必ずA-安定である。

さらに、ルンゲ＝クッタ法の係数が テンプレート:Math かつ行列

$M=(b_i{a}_{ij}-b_j{a}_{ji}-b_i{b}_j{)}_{ij}$

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

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

テンプレート:Mvar 段線型多段法の一般形式は次の公式で与えられる。

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

ここで、テンプレート:Mvar と テンプレート:Mvar は定数であり、ベクトル テンプレート:Math と テンプレート:Math は方法の 生成ペア (generating pair) という。 この方法に対応する one-leg法 (one-leg method) は次の公式で与えられるテンプレート:Sfn。

$\sum _{i=0}^k{\alpha }_i{y}_{n-i}=h\left(\sum _{i=0}^k{\beta }_i\right)f(x_n-\theta h,(\sum _{i=0}^k{\beta }_i{)}^{-1}\sum _{i=0}^k{\beta }_i{y}_{n-i})$

ただし、

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

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

与えられた テンプレート:Mvar 次正定値対称行列 テンプレート:Mvar に対応する ${ℝ}^{kn}$ 上の内積を以下のように定義できる：

$⟨U,V{⟩}_G=\sum _{i=1}^k\sum _{j=1}^k{g}_{ij}⟨U_i,V_j⟩,U,V\in {ℝ}^{kn}$

ここで、$U=(U_1,\dots ,U_k{)}^{\mathrm{T}},U_i\in {ℝ}^n$である。

one-leg法に対し、行列

$M=\alpha {\beta }^{\mathrm{T}}+\beta {\alpha }^{\mathrm{T}}-\left[\begin{array}{cc}G& 0\\ 0& 0\end{array}\right]+\left[\begin{array}{cc}0& 0\\ 0& G\end{array}\right]$

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

$\Vert y_1-{\hat{y}}_1{\Vert }_G\le \Vert y_0-{\hat{y}}_0{\Vert }_G$

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

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

参考文献

• テンプレート:Citation.
• テンプレート:Citation.
• テンプレート:Citation.
• テンプレート:Citation.
• テンプレート:Citation.
• テンプレート:Citation.
• テンプレート:Citation.

外部リンク

• 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