硬い方程式

数学において、スティッフ方程式とは、微分方程式を解くための特定の数値解法が、刻み幅を極めて小さくしない限り数値的に不安定となる方程式のことである。スティッフネスの正確な定義を定式化することは困難であることが判明しているが、その主要な考え方は、方程式に解の急激な変化につながる可能性のある項が含まれているというものである。

微分方程式を数値的に積分する場合、解曲線が大きく変化する領域では必要なステップサイズが比較的小さく、解曲線が直線になって傾きがほぼ 0 の線に近づく領域では必要なステップサイズが比較的大きくなることが予想されます。しかし、問題によってはそうならないことがあります。数値法で微分方程式に信頼性の高い解を与えるには、解曲線が非常に滑らかな領域ではステップサイズが許容できないほど小さい必要がある場合があります。この現象は剛性として知られています。同じ解を持つ 2 つの異なる問題があり、一方が剛性でもう一方が剛性である場合があります。この現象は両方の問題で同じであるため、厳密な解の特性であるはずがなく、微分方程式自体の特性でなければなりません。そのため、このようなシステムは剛性システムとして知られています。

動機付けの例

初期値問題を考える

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

1

正確な解（水色で表示）は

y(t)=e^{-15t},\quad y(t)\to 0{\text{ as }}t\to \infty .

2

同じ動作を示す数値解を求めます。

図 (右) は、方程式に適用されたさまざまな数値積分器の数値的問題を示しています。

ステップサイズがのオイラー法は大きく振動し、すぐにグラフの範囲から外れます (赤で表示)。 $h={\tfrac {1}{4}}$
ステップサイズを半分にしたオイラー法では、グラフの境界内で解が生成されますが、ゼロの周りで振動します (緑色で表示)。 $h={\tfrac {1}{8}}$

台形法（すなわち2段階アダムス・モールトン法）は次のように表される。

y_{n+1}=y_{n}+{\tfrac {1}{2}}h{\bigl (}f(t_{n},y_{n})+f(t_{n+1},y_{n+1}){\bigr )},

3

ここで、オイラー法の代わりにこの方法を適用すると、はるかに良い結果（青）が得られます。数値結果は、厳密解と同様に、単調にゼロまで減少します。

y'=f(t,y)

別の例

硬い常微分方程式（ODE）の最も顕著な例の一つは、ロバートソンの化学反応を記述するシステムである。 ^[1]

{\begin{aligned}\mathrm {A} \xrightarrow {0.04} &\mathrm {B} \\\mathrm {B} +\mathrm {B} \xrightarrow {3\cdot 10^{7}} &\mathrm {C} +\mathrm {B} \\\mathrm {B} +\mathrm {C} \xrightarrow {1\cdot 10^{4}} &\mathrm {A} +\mathrm {C} \end{aligned}}

${\begin{aligned}{\dot {x}}&=-0.04x+10^{4}y\cdot z\\{\dot {y}}&={\hphantom {-}}0.04x-10^{4}y\cdot z-3\cdot 10^{7}y^{2}\\{\dot {z}}&=3\cdot 10^{7}y^{2}\end{aligned}}$

4

例えば、この系を短い区間で扱う場合、数値積分は問題ありません。しかし、区間が非常に大きい場合（例えば10の¹¹乗）、多くの標準コードでは正しく積分できません。 $t\in [0,40]$

剛性比

線形定数係数不均質系を考える

\mathbf {y} '=\mathbf {A} \mathbf {y} +\mathbf {f} (x),

5

ここで、およびは、固有値（異なると仮定）とそれに対応する固有ベクトルを持つ定数対角化可能な行列である。（ 5 ）の一般解は以下の形をとる。 $\mathbf {y} ,\mathbf {f} \in \mathbb {R} ^{n}$ $\mathbf {A}$ $n\times n$ $\lambda _{t}\in \mathbb {C} ,t=1,2,\ldots ,n$ $\mathbf {c} _{t}\in \mathbb {C} ^{n},t=1,2,\ldots ,n$

\mathbf {y} (x)=\sum _{t=1}^{n}\kappa _{t}e^{\lambda _{t}x}\mathbf {c} _{t}+\mathbf {g} (x),

6

ここで、は任意の定数、は特定の積分である。ここで、 $\kappa _{t}$ $\mathbf {g} (x)$

\operatorname {Re} (\lambda _{t})<0,\qquad t=1,2,\ldots ,n,

7

これは、各項がとして成り立つことを意味し、解はとして漸近的に近づきます。項はが実数の場合には単調に減少し、が複素数の場合には正弦的に減少します。 $e^{\lambda _{t}x}\mathbf {c} _{t}\to 0$ $x\to \infty$ $\mathbf {y} (x)$ $\mathbf {g} (x)$ $x\to \infty$ $e^{\lambda _{t}x}\mathbf {c} _{t}$ $\lambda _{t}$ $\lambda _{t}$

を時間として解釈する（物理問題ではよくあるように）ことを過渡解、定常解と呼ぶ。が大きい場合、対応する項はが増加するにつれて急速に減少するため、高速過渡と呼ばれる。が小さい場合、対応する項はが緩やかに減少するため、低速過渡と呼ばれる。をで定義する。 $x$ ${\textstyle \sum _{t=1}^{n}\kappa _{t}e^{\lambda _{t}x}\mathbf {c} _{t}}$ $\mathbf {g} (x)$ $\left|\operatorname {Re} (\lambda _{t})\right|$ $\kappa _{t}e^{\lambda _{t}x}\mathbf {c} _{t}$ $x$ $\left|\operatorname {Re} (\lambda _{t})\right|$ $\kappa _{t}e^{\lambda _{t}x}\mathbf {c} _{t}$ ${\overline {\lambda }},{\underline {\lambda }}\in \{\lambda _{t},t=1,2,\ldots ,n\}$

\left|\operatorname {Re} ({\overline {\lambda }})\right|\geq \left|\operatorname {Re} (\lambda _{t})\right|\geq \left|\operatorname {Re} ({\underline {\lambda }})\right|,\qquad t=1,2,\ldots ,n

8

最も速い過渡応答と最も遅い過渡応答である。ここで剛性比を^[2]で定義する。 $\kappa _{t}e^{{\overline {\lambda }}x}\mathbf {c} _{t}$ $\kappa _{t}e^{{\underline {\lambda }}x}\mathbf {c} _{t}$

{\frac {\left|\operatorname {Re} ({\overline {\lambda }})\right|}{\left|\operatorname {Re} ({\underline {\lambda }})\right|}}.

9

剛性の特性評価

この節では、剛性という現象の様々な側面について考察します。「現象」という言葉は「性質」という言葉よりも適切でしょう。なぜなら、後者はむしろ剛性を厳密な数学的用語で定義できることを意味するからです。しかし、線形定数係数系という限定されたクラスにおいてさえ、満足のいく方法で定義することは不可能であることが判明しています。また、剛性の概念を要約するために提示できる（そしてほとんど提示されてきた）いくつかの定性的な記述についても考察し、それらの中で剛性の「定義」としておそらく最も適切なものを述べます。

JD Lambert は剛性を次のように定義しています。

絶対安定の有限領域を持つ数値法を、任意の初期条件を持つシステムに適用し、特定の積分区間で、その区間での正確な解の滑らかさに比べて過度に小さいステップ長を使用せざるを得ない場合、その区間ではシステムは硬いと言われます。

スティフな問題の多くの例に見られる他の特性もありますが、それぞれに反例が存在するため、これらの特性はスティフネスの適切な定義とはなり得ません。とはいえ、これらの特性に基づく定義は一部の著者によって一般的に用いられており、スティフネスの存在を示す優れた手がかりとなります。ランバートは前述の理由から、これらを定義ではなく「ステートメント」と呼んでいます。そのいくつかを以下に示します。

線形定数係数システムは、その固有値すべてが負の実部を持ち、剛性比が大きい場合、剛性です。
剛性は、精度の要件ではなく、安定性の要件によって歩幅が制限されるときに発生します。
硬直性は、溶液のいくつかの成分が他の成分よりもはるかに急速に崩壊するときに発生します。^[3]

語源

「剛性」という用語の起源は明確に解明されていない。ジョセフ・オークランド・ヒルシュフェルダーによれば、「剛性」という用語が用いられるのは、このようなシステムが駆動部と被駆動部との間の密接な結合に対応するためである^[4]。リチャード・L・バーデンとJ・ダグラス・フェアーズによれば、

厳密解に（は負の実部を持つ複素数）の形式の項が含まれる場合、標準的な数値手法を使用して微分方程式の解を近似すると、重大な問題が発生する可能性があります。 $e^{\lambda t}$ $\lambda$
. . .
急速に減衰する過渡解を含む問題は、バネや減衰システムの研究、制御システムの解析、化学反応速度論の問題など、幅広い応用分野で自然に発生します。これらはすべて、大きなバネ定数（物理的剛性）を持つバネや質量システムの運動解析に応用されるため、「硬い（数学的剛性）」微分方程式系と呼ばれる一連の問題の例です。 [ 5 ^]

例えば、初期値問題

m{\ddot {x}}+c{\dot {x}}+kx=0,\qquad x(0)=x_{0},\qquad {\dot {x}}(0)=0,

10

、、は（５）の形で表すことができ、および $m=1$ $c=1001$ $k=1000$ $n=2$

\mathbf {A} ={\begin{pmatrix}0&1\\-1000&-1001\end{pmatrix}},\qquad \mathbf {f} (t)={\begin{pmatrix}0\\0\end{pmatrix}},\qquad \mathbf {x} (0)={\begin{pmatrix}x_{0}\\0\end{pmatrix}},

11

は固有値を持つ。両方の固有値は負の実部を持ち、剛性比は ${\overline {\lambda }}=-1000,{\underline {\lambda }}=-1$

{\frac {\left|-1000\right|}{\left|-1\right|}}=1000,

12

これはかなり大きい。するとシステム（10）は確かに記述1と3を満たす。ここではバネ定数は大きく、減衰定数はさらに大きい。^[6]（「大きい」という言葉は明確に定義されていないが、上記の量が大きいほど剛性の影響はより顕著になる。）（10）の正確な解は $k$ $c$

x(t)=x_{0}\left(-{\frac {1}{999}}e^{-1000t}+{\frac {1000}{999}}e^{-t}\right)\approx x_{0}e^{-t}.

13

式13は単純な指数関数と非常によく似た挙動を示すが、係数が小さくても項が存在するため、数値計算はステップサイズに非常に敏感になる。( 10 )の安定積分には、解曲線の滑らかな部分まで非常に小さなステップサイズが必要であるため、精度に必要な誤差よりもはるかに小さい誤差となる。したがって、この系はステートメント2とランバートの定義も満たしている。 $x_{0}e^{-t}$ $e^{-1000t}$

A安定性

硬い問題に対する数値法の挙動は、という初期条件を条件とするテスト方程式にこれらの方法を適用することで分析できます。この方程式の解はです。この解は、のときに 0 に近づきます。数値法がこの挙動 (固定ステップサイズに対して) を示す場合、その方法は A 安定であるといわれます。^[7] L 安定 (以下を参照) な数値法は、ステップサイズが無限大になると 1 ステップで解が 0 に近づくという、より強い特性があります。A 安定な方法では、例で説明したような不安定性の問題は発生しません。 $y'=ky$ $y(0)=1$ $k\in \mathbb {C}$ $y(t)=e^{kt}$ $t\to \infty$ $\operatorname {Re} (k)<0.$

ルンゲ・クッタ法

ルンゲ・クッタ法をテスト方程式に適用すると、の形をとり、帰納法によってとなる。この関数は安定関数と呼ばれる。したがって、がである条件はと等価である。これが、の集合である絶対安定領域（単に安定領域と呼ばれることもある）の定義の根拠となる。この方法は、絶対安定領域にの集合、すなわち左半平面が含まれる場合、A安定である。 $y'=k\cdot y$ $y_{n+1}=\phi (hk)\cdot y_{n}$ $y_{n}={\bigl (}\phi (hk){\bigr )}^{n}\cdot y_{0}$ $\phi$ $y_{n}\to 0$ $n\to \infty$ $|\phi (hk)|<1$ ${\bigl \{}z\in \mathbb {C} \,{\big |}\,|\phi (z)|<1{\bigr \}}$ ${\bigl \{}z\in \mathbb {C} \,{\big |}\,\operatorname {Re} (z)<0{\bigr \}}$

例: オイラー法

上記のオイラー法について考えてみましょう。テスト方程式に適用される陽的オイラー法はです。したがって、となります。この方法の絶対安定領域はであり、これは右図に示す円板に相当します。オイラー法はA安定ではありません。 $y'=k\cdot y$ $y_{n+1}=y_{n}+h\cdot f(t_{n},y_{n})=y_{n}+h\cdot (ky_{n})=y_{n}+h\cdot k\cdot y_{n}=(1+h\cdot k)y_{n}.$ $y_{n}=(1+hk)^{n}\cdot y_{0}$ $\phi (z)=1+z$ ${\bigl \{}z\in \mathbb {C} \,{\big |}\,|1+z|<1{\bigr \}}$

例としてはがありました。ステップサイズをとした場合のzの値はとなり、これは安定領域外にあります。実際、数値結果はゼロに収束しません。しかし、ステップサイズをとした場合、となり、これは安定領域内にあり、数値結果はゼロに収束しますが、収束速度はかなり遅いです。 $k=-15$ $h={\tfrac {1}{4}}$ $z=-15\times {\tfrac {1}{4}}=-3.75$ $h={\tfrac {1}{8}}$ $z=-1.875$

例: 台形法

台形法をテスト方程式に適用する場合を考えてみると、はとなる。を解くと、が得られる。したがって、安定関数はとなり、絶対安定領域はとなる。この領域には左半平面が含まれるため、台形法は A 安定である。実際、安定領域は左半平面と同じであるため、の数値解がゼロに収束するのは、厳密な解がゼロに収束する場合に限る。ただし、台形法は完璧な動作をするわけではない。すべての減衰成分を減衰させるが、のときはとなるため、急速に減衰する成分はごくわずかにしか減衰しない。このことからL 安定性の概念が生まれた。つまり、ある方法が L 安定とは、A 安定でのときはとなる。台形法は A 安定だが L 安定ではない。暗黙のオイラー法は L 安定法の一例である。^[8] $y_{n+1}=y_{n}+{\tfrac {1}{2}}h\cdot {\bigl (}f(t_{n},y_{n})+f(t_{n+1},y_{n+1}){\bigr )},$ $y'=k\cdot y$ $y_{n+1}=y_{n}+{\tfrac {1}{2}}h\cdot \left(ky_{n}+ky_{n+1}\right).$ $y_{n+1}$ $y_{n+1}={\frac {1+{\frac {1}{2}}hk}{1-{\frac {1}{2}}hk}}\cdot y_{n}.$ $\phi (z)={\frac {1+{\frac {1}{2}}z}{1-{\frac {1}{2}}z}}$ $\left\{z\in \mathbb {C} \ \left|\ \left|{\frac {1+{\frac {1}{2}}z}{1-{\frac {1}{2}}z}}\right|<1\right.\right\}.$ $y'=k\cdot y$ $\phi (z)\to 1$ $z\to -\infty$ $|\phi (z)|\to 0$ $z\to \infty$

一般理論

係数とを持つルンゲ・クッタ法の安定関数は、で与えられます。ここではすべて1のベクトルを表します。これは有理関数（ある多項式を別の多項式で割る）です。 $\mathbf {A}$ $\mathbf {b}$ $\phi (z)={\frac {\det \left(\mathbf {I} -z\mathbf {A} +z\mathbf {e} \mathbf {b} ^{\mathsf {T}}\right)}{\det \left(\mathbf {I} -z\mathbf {A} \right)}},$ $\mathbf {e}$

陽的ルンゲ・クッタ法は、厳密に下三角の係数行列を持つため、その安定性関数は多項式となる。したがって、陽的ルンゲ・クッタ法はA安定ではない。 $\mathbf {A}$

暗黙的ルンゲ・クッタ法の安定性関数は、しばしば順序スターを用いて解析される。安定性関数を持つルンゲ・クッタ法の順序スターは、集合と定義される。ある方法がA安定であるためには、その安定性関数が左平面に極を持たず、かつその順序スターに純虚数が含まれない必要がある。^[9] $\phi$ ${\bigl \{}z\in \mathbb {C} \,{\big |}\,|\phi (z)|>|e^{z}|{\bigr \}}$

多段階法

線形多段階法は、という形をとります。これをテスト方程式に適用すると、となり、はと簡略化できます。これは線形再帰関係です。のとき、再帰関係のすべての解がゼロに収束する場合、この方法はA安定です。特性多項式はです。のすべての解が単位円内にある場合、の与えられた値に対してすべての解がゼロに収束します。 $y_{n+1}=\sum _{i=0}^{s}a_{i}y_{n-i}+h\sum _{j=-1}^{s}b_{j}f\left(t_{n-j},y_{n-j}\right).$ $y_{n+1}=\sum _{i=0}^{s}a_{i}y_{n-i}+hk\sum _{j=-1}^{s}b_{j}y_{n-j},$ $\left(1-b_{-1}z\right)y_{n+1}-\sum _{j=0}^{s}\left(a_{j}+b_{j}z\right)y_{n-j}=0$ $z=hk$ $\{y_{n}\}$ $\operatorname {Re} (z)<0$ $\Phi (z,w)=w^{s+1}-\sum _{i=0}^{s}a_{i}w^{s-i}-z\sum _{j=-1}^{s}b_{j}w^{s-j}.$ $z$ $w$ $\Phi (z,w)=0$

上記の形式の多段階法における絶対安定領域は、を満たすすべてのの集合である。また、この集合に左半平面が含まれる場合、多段階法はA安定であると言われる。 $z\in \mathbb {C}$ $w$ $\Phi (z,w)=0$ $|w|<1$

例: 2次アダムズ・バッシュフォース法

2段階アダムズ・バッシュフォース法の絶対安定領域を決定しましょう。特性多項式はであり、根を持つので、絶対安定領域はです。この領域は右側に示されています。この領域は左半平面全体を含んでいません（実際には間の実軸のみを含みます）。そのため、アダムズ・バッシュフォース法はA安定ではありません。 $y_{n+1}=y_{n}+h\left({\tfrac {3}{2}}f(t_{n},y_{n})-{\tfrac {1}{2}}f(t_{n-1},y_{n-1})\right).$ $\Phi (w,z)=w^{2}-\left(1+{\tfrac {3}{2}}z\right)w+{\tfrac {1}{2}}z=0$ $w={\tfrac {1}{2}}\left(1+{\tfrac {3}{2}}z\pm {\sqrt {1+z+{\tfrac {9}{4}}z^{2}}}\right),$ $\left\{z\in \mathbb {C} \ \left|\ \left|{\tfrac {1}{2}}\left(1+{\tfrac {3}{2}}z\pm {\sqrt {1+z+{\tfrac {9}{4}}z^{2}}}\right)\right|<1\right.\right\}.$ $-1\leq z\leq 0$

一般理論

陽的多段階法は、陽的ルンゲ・クッタ法と同様に、A安定にはなり得ません。陰的多段階法は、その次数が2以下の場合にのみA安定となります。後者の結果は第二ダールキスト障壁として知られており、硬い方程式に対する線形多段階法の有用性を制限します。2次のA安定法の例として、前述の台形則が挙げられますが、これも線形多段階法とみなすことができます。^[10]

参照

後退微分公式、特に硬い微分方程式の解法に使用される暗黙法の一種
条件数
微分包含は、不連続性を許容する微分方程式の概念の拡張であり、剛性の問題を回避するための方法でもある。
明示的メソッドと暗黙的メソッド

注記

^ Robertson, HH (1966). 「反応速度方程式の解」.数値解析入門. Academic Press. pp. 178– 182.
^ ランバート（1992年、216～217ページ）
^ ランバート（1992年、217～220ページ）
^ ハーシュフェルダー（1963）
^ Burden & Faires (1993、p. 314)
^ クライスツィヒ (1972, pp. 62–68)
^ この定義はDahlquist（1963）によるものです。
^ L-安定性の定義はEhle（1969）による。
^ この定義はWanner、Hairer、Nørsett（1978）によるものです。また、IserlesとNørsett（1991）も参照してください。
^ Dahlquist（1963）を参照。

参考文献

バーデン、リチャード・L.; フェアーズ、J.ダグラス (1993)、『数値解析』（第5版）、ボストン：プリンドル、ウェーバー、シュミット、ISBN 0-534-93219-3。
Dahlquist, Germund (1963)、「線形多段階法の特殊な安定性問題」、BIT、3 (1): 27– 43、doi :10.1007/BF01963532、hdl : 10338.dmlcz/103497、S2CID 120241743。
エバリー、デイビッド（2008）、微分方程式系の安定性解析（PDF）。
Ehle, BL (1969), 指数関数のパデ近似と初期値問題の数値解法におけるA安定法について(PDF)、ウォータールー大学。
ギア、CW（1971）、常微分方程式の数値初期値問題、エングルウッドクリフス：プレンティスホール、Bibcode：1971nivp.book.....G。
ギア、CW（1981）、「常微分方程式の数値解：やるべきことは残っているか？」SIAMレビュー、23（1）：10-24、doi：10.1137/1023002。
ヘアラー、エルンスト、ワーナー、ゲルハルト（1996年）、常微分方程式の解法II：スティフ問題と微分代数問題（第2版）、ベルリン：シュプリンガー・フェアラーク、ISBN 978-3-540-60452-5。
ハーシュフェルダー, JO ( 1963)、「理論化学における応用数学」アメリカ数学会シンポジウム: 367–376。
イゼルレス、アリエ。 Nørsett、Syvert (1991)、Order Stars、Chapman & Hall、ISBN 978-0-412-35260-7。
クレイジグ、エルウィン（1972）、Advanced Engineering Mathematics（第3版）、ニューヨーク：Wiley、ISBN 0-471-50728-8。
ランバート、JD （1977）、D.ジェイコブス（編）、「常微分方程式の初期値問題」、数値解析の最新技術、ニューヨーク：アカデミックプレス：451-501。
ランバート、JD（1992）、常微分システムの数値解析法、ニューヨーク：ワイリー、ISBN 978-0-471-92990-1。
マシューズ、ジョン; フィンク、カーティス (1992)「MATLABを用いた数値解析法」。
Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). 「Section 17.5. Stiff Sets of Equations」. Numerical Recipes: The Art of Scientific Computing (第3版). ニューヨーク: Cambridge University Press. ISBN 978-0-521-88068-8. 2011年8月11日時点のオリジナルよりアーカイブ。2011年8月17日閲覧。
シャンピン, LF; ギア, CW (1979)「スティフ常微分方程式の解法に関するユーザーの視点」SIAM Review , 21 (1): 1– 17, doi :10.1137/1021001。
ワナー、ゲルハルト。ハイラー、エルンスト。 Nørsett、Syvert (1978)、「秩序星と安定性理論」、BIT、18 (4): 475–489、doi :10.1007/BF01932026、S2CID 8824105。
ルンゲ・クッタ法の安定性 [1]

外部リンク

物理ベースモデリング入門：エネルギー関数と剛性
スティフシステムローレンス・F・シャンパインとスキップ・トンプソンScholarpedia、2(3):2855. doi:10.4249/scholarpedia.2855

[1] Robertson, HH (1966). 「反応速度方程式の解」.数値解析入門. Academic Press. pp. 178– 182.

[2] ランバート（1992年、216～217ページ）

[3] ランバート（1992年、217～220ページ）

[4] ハーシュフェルダー（1963）

[5] Burden & Faires (1993、p. 314)

[6] クライスツィヒ (1972, pp. 62–68)

[7] この定義はDahlquist（1963）によるものです。

[8] L-安定性の定義はEhle（1969）による。

[9] この定義はWanner、Hairer、Nørsett（1978）によるものです。また、IserlesとNørsett（1991）も参照してください。

[10] Dahlquist（1963）を参照。