硬い方程式

数学において、硬い方程式とは、ステップ サイズを極端に小さくしない限り、方程式を解く特定の数値方法が数値的に不安定になる微分方程式です。硬さの正確な定義を定式化することは困難であることが証明されていますが、主な考え方は、方程式に解の急激な変化につながる可能性のあるいくつかの項が含まれているということです。

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

初期値問題を考える

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

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

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

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

たとえば、このシステムを短い区間で扱う場合、数値積分に問題はありません。しかし、区間が非常に大きい場合 (たとえば 10 ¹¹ )、多くの標準コードでは正しく積分できません。 ${\displaystyle t\in [0,40]}$

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

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

ここで は任意の定数、 は特定の積分です。ここで ${\displaystyle \kappa_{t}}$${\displaystyle \mathbf {g} (x)}$

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

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

最も速い過渡応答と最も遅い過渡応答です。剛性比を次のように定義します^[2]${\displaystyle \kappa _{t}e^{{\overline {\lambda }}x}\mathbf {c} _{t}}$${\displaystyle \kappa _{t}e^{{\underline {\lambda }}x}\mathbf {c} _{t}}$

このセクションでは、剛性の現象のさまざまな側面について考察します。「現象」は「特性」よりも適切な言葉でしょう。なぜなら、後者はむしろ剛性を正確な数学的用語で定義できることを意味します。しかし、線形定数係数システムの限定されたクラスであっても、満足のいく方法で定義することはできないことが判明しています。また、剛性の概念を要約するために作成できる (そしてほとんど作成されている) 定性的な記述をいくつか見ていき、剛性の「定義」としておそらく最も満足のいくものを述べます。

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

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

硬い問題の多くの例で示される他の特性もありますが、それぞれに反例があるため、これらの特性は硬いことの適切な定義にはなりません。それでも、これらの特性に基づく定義は一部の著者によって一般的に使用されており、硬いことの存在に関する良い手がかりとなります。Lambert は、前述の理由から、これらを定義ではなく「ステートメント」と呼んでいます。これらのいくつかは次のとおりです。

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

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

厳密な解に（ は実部が負の複素数）の 形式の項が含まれる場合、標準的な数値手法を微分方程式の解の近似に適用すると、重大な問題が発生する可能性があります。${\displaystyle e^{\lambda t}}$${\displaystyle \lambda }$

. . .

急速に減衰する過渡解を含む問題は、バネや減衰システムの研究、制御システムの解析、化学反応速度論の問題など、さまざまな分野で自然に発生します。これらはすべて、大きなバネ定数（物理的剛性）を持つバネや質量システムの運動を解析する際に使用されるため、硬い（数学的剛性）微分方程式系と呼ばれる一連の問題の例です。 [ ^5]

例えば、初期値問題

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

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

これはかなり大きい。システム（１０）は確かにステートメント１と３を満たす。ここではバネ定数は大きく、減衰定数はさらに大きい。^[6]（「大きい」は明確に定義された用語ではないが、上記の量が大きいほど、剛性の影響はより顕著になる。）（１０）の正確な解は ${\displaystyle k}$${\displaystyle c}$

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

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

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

上記のオイラー法について考えてみましょう。テスト方程式に適用される明示的なオイラー法は、 ${\displaystyle y'=k\cdot y}$

${\displaystyle 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}.}$

したがって、となります。この方法の絶対安定領域は となり、これは右側に描かれた円です。オイラー法は A 安定ではありません。 ${\displaystyle y_{n}=(1+hk)^{n}\cdot y_{0}}$${\displaystyle \phi (z)=1+z}$${\displaystyle {\bigl \{}z\in \mathbb {C} \,{\big |}\,|1+z|<1{\bigr \}}}$

動機付けの例は でした。ステップ サイズを取ったときのzの値は であり、これは安定領域外です。実際、数値結果はゼロに収束しません。ただし、ステップ サイズ では となり、これは安定領域内にあり、数値結果は、かなりゆっくりでは ありますが、ゼロに収束します。${\displaystyle k=-15}$${\displaystyle h={\tfrac {1}{4}}}$${\displaystyle z=-15\times {\tfrac {1}{4}}=-3.75}$${\displaystyle h={\tfrac {1}{8}}}$${\displaystyle z=-1.875}$

台形法を検討する

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

これをテスト方程式に適用すると、 ${\displaystyle y'=k\cdot y}$

${\displaystyle y_{n+1}=y_{n}+{\tfrac {1}{2}}h\cdot \left(ky_{n}+ky_{n+1}\right).}$

利回り の解決${\displaystyle y_{n+1}}$

${\displaystyle y_{n+1}={\frac {1+{\frac {1}{2}}hk}{1-{\frac {1}{2}}hk}}\cdot y_{n}.}$

したがって、安定性関数は

${\displaystyle \phi (z)={\frac {1+{\frac {1}{2}}z}{1-{\frac {1}{2}}z}}}$

そして絶対安定領域は

${\displaystyle \left\{z\in \mathbb {C} \ \left|\ \left|{\frac {1+{\frac {1}{2}}z}{1-{\frac {1}{2}}z}}\right|<1\right.\right\}.}$

この領域には左半平面が含まれるため、台形法は A 安定です。実際、安定領域は左半平面と同一であり、したがって の数値解がゼロに収束するのは、厳密な解がゼロに収束する場合のみです。ただし、台形法は完全な動作をしません。台形法はすべての減衰成分を減衰しますが、急速に減衰する成分は として非常に緩やかに減衰するだけです。これが L 安定性の概念につながりました。つまり、ある方法が L 安定であるのは、それが A 安定であり、としてである場合です。台形法は A 安定ですが、L 安定ではありません。暗黙のオイラー法はL 安定法の例です。^[8]${\displaystyle y'=k\cdot y}$${\displaystyle \phi (z)\to 1}$${\displaystyle z\to -\infty }$${\displaystyle |\phi (z)|\to 0}$${\displaystyle z\to \infty }$

係数を持つルンゲ・クッタ法の安定関数は次のように与えられる 。${\displaystyle \mathbf {A} }$${\displaystyle \mathbf {b} }$

${\displaystyle \phi (z)={\frac {\det \left(\mathbf {I} -z\mathbf {A} +z\mathbf {e} \mathbf {b} ^{\mathsf {T}}\right)}{\det(\mathbf {I} -z\mathbf {A} )}},}$

ここで、すべて 1 のベクトルを表します。これは有理関数(1 つの多項式を別の多項式で割ったもの) です。 ${\displaystyle \mathbf {e} }$

明示的ルンゲ・クッタ法は厳密に下三角の係数行列を持つため、その安定性関数は多項式になります。したがって、明示的ルンゲ・クッタ法は A 安定ではありません。 ${\displaystyle \mathbf {A} }$

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

線形多段階法は次のような形式をとる。

${\displaystyle 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).}$

これをテスト方程式に適用すると、

${\displaystyle y_{n+1}=\sum _{i=0}^{s}a_{i}y_{n-i}+hk\sum _{j=-1}^{s}b_{j}y_{n-j},}$

これを簡略化すると

${\displaystyle \left(1-b_{-1}z\right)y_{n+1}-\sum _{j=0}^{s}\left(a_{j}+b_{j}z\right)y_{n-j}=0}$

ここで、これは線形再帰関係である。再帰関係のすべての解が のときにゼロに収束する場合、この方法は A 安定である。特性多項式は ${\displaystyle z=hk}$${\displaystyle \{y_{n}\}}$${\displaystyle \operatorname {Re} (z)<0}$

${\displaystyle \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}.}$

のすべての解が単位円内にある 場合、の与えられた値に対してすべての解はゼロに収束します。${\displaystyle z}$${\displaystyle w}$${\displaystyle \Phi (z,w)=0}$

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

2段階アダムス・バッシュフォース法の絶対安定領域を決定しよう

${\displaystyle 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).}$

特性多項式は

${\displaystyle \Phi (w,z)=w^{2}-\left(1+{\tfrac {3}{2}}z\right)w+{\tfrac {1}{2}}z=0}$

根がある

${\displaystyle w={\tfrac {1}{2}}\left(1+{\tfrac {3}{2}}z\pm {\sqrt {1+z+{\tfrac {9}{4}}z^{2}}}\right),}$

したがって絶対安定領域は

${\displaystyle \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\}.}$

この領域は右側に示されています。この領域には左半平面全体は含まれません（実際には 間の実軸のみが含まれます）ので、Adams-Bashforth 法は A 安定ではありません。 ${\displaystyle -1\leq z\leq 0}$

明示的マルチステップ法は、明示的ルンゲ・クッタ法と同様に、A安定になることはありません。暗黙的マルチステップ法は、次数が最大2の場合にのみA安定になります。後者の結果は、第2ダルキスト障壁として知られており、硬い方程式に対する線形マルチステップ法の有用性を制限します。2次のA安定法の例は、前述の台形則であり、これも線形マルチステップ法と見なすことができます。^[10]

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

• バーデン、リチャード 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。
• Lambert, JD (1977)、D. Jacobs (編)、「常微分方程式の初期値問題」、数値解析の最新技術、ニューヨーク: Academic Press : 451–501。
• ランバート、JD（1992）、常微分システムの数値解析法、ニューヨーク：ワイリー、ISBN 978-0-471-92990-1。
• マシューズ、ジョン、フィンク、カーティス (1992)、MATLAB を使用した数値解析法。
• Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007)。「セクション 17.5. 方程式のスティフ セット」。数値レシピ: 科学計算の技法(第 3 版)。ニューヨーク: Cambridge University Press。ISBN 978-0-521-88068-8. 2011年8月11日時点のオリジナルよりアーカイブ。2011年8月17日閲覧。
• Shampine, LF; Gear, 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]

• 物理ベースモデリング入門: エネルギー関数と剛性
• スティフシステム Lawrence F. Shampine および Skip Thompson Scholarpedia、2(3):2855. doi:10.4249/scholarpedia.2855