硬い方程式

計算数学において、スティッフ方程式は初期値問題である。

${\displaystyle {\dot {u}}=f(u)\,,\qquad u(0)=u_{0}\,,\qquad t\in [0,T]\,,}$

どこ${\displaystyle f:{\mathbb {R} }^{d}\rightarrow {\mathbb {R} }^{d}}$効率的な数値積分には専用の暗黙的時間ステップ法が必要となる。硬い方程式の最も単純な数学的特徴付けは、必要条件である。

${\displaystyle {\frac {T}{d}}{\big (}{\mathrm {div} }_{u}\,f{\big )}(u)\ll -1\,.}$

以来${\displaystyle {\big (}{\mathrm {div} }_{u}\,f{\big )}(u)={\mathrm {trace} }({\mathrm {grad} }_{u}\,f)(u)={\mathrm {trace} }\,f'(u)}$、 どこ${\displaystyle f'(u)\in {\mathbb {R} }^{d\times d}}$は、​${\displaystyle f}$その時点で${\displaystyle u}$上記の基準は容易に評価でき、剛性を定量化します。この基準は、以下の非線形剛性方程式について導出、説明、および図示されています。

定数係数を持つ 線形システムの場合${\displaystyle {\dot {u}}=Au}$発散は一定であり、剛性は時間スケールに関連するグローバルな特性となる。${\displaystyle T}$。

非線形システムの場合、剛性は通常、解の軌跡に沿って空間的にも時間的にも変化する。${\displaystyle u(t)}$ここで、基準は局所的な剛性を定量化する。実際の計算では、剛性方程式は必ず適応法を用いて解かれる。^{[ 1 ] [ 2 ]}

硬い微分方程式に関する文献は豊富にあるが、概念の厳密な定義を試みるよりも、直感的な説明やヒューリスティックの方がはるかに一般的である。HairerとWanner ^{[ 3 ]}は、最も明白な特徴を簡潔に説明している。

硬直方程式とは、明示的な解法が適用できない問題のことである。

これは、明示的な積分法では極めて小さな時間ステップを使用せざるを得ないという観察結果を指している。${\displaystyle h}$数値安定性を維持するため、このような方法は競争力を持ちません。各ステップは安価ですが、ステップの総数は${\displaystyle N=T/h}$非常に大きくなり、${\displaystyle [0,T]}$効果的に時間稼ぎをする。

対照的に、硬い方程式に対する陰解法では、各ステップでコストのかかる「代数」方程式の解法が必要となる。ステップごとの追加作業は、優れた安定性によって相殺され、大きな時間ステップの使用が可能となる。安定性に関する制約が厳しくなければ、全体の計算量は管理可能な範囲に収まり、効率を損なうことなく要求される精度を達成できる。一部の硬い方程式では、その効率は最良の陽解法よりも数桁高い場合がある。

剛性現象は、数値解析の他の分野におけるよく知られた性能問題と類似している。例えば、最適化においては、勾配法にも同様の制限がある。凸関数を最小化するために最急降下法を用いる場合${\displaystyle F(x)}$反復につながる

${\displaystyle x^{k+1}=x^{k}-h\,{\mathrm {grad} }_{x}F(x^{k})\,,}$

どこ${\displaystyle h}$は線探索のステップサイズです。これは単に微分方程式に対する明示的なオイラー法です。${\displaystyle {\dot {x}}=-{\mathrm {grad} }_{x}F(x^{k})}$急勾配の問題は、選択の制約により収束が遅いことで知られています。${\displaystyle h}$一般的な解決策は、勾配法をニュートン法のような方法に置き換えることであり、これは本質的に微分方程式に対して「陰解法」を用いることに相当する。

別の例としては、次の形式の非線形方程式を解く場合が挙げられる。${\displaystyle x=g(x)}$単純な不動点反復、${\displaystyle x^{k+1}=g(x^{k})}$、リプシッツ定数が${\displaystyle L[g]}$は大きいが、収縮の場合に限る。しかし、この方程式は、より弱いが一般的な条件、例えばマップが${\displaystyle gI}$は単調であり、対数リプシッツ定数は${\displaystyle g}$満たす${\displaystyle M[g]<1}$固定点反復法は発散するため、ニュートン法のような反復法を用いる必要がある。

3つ目の例は、（離散化された）ポアソン方程式の反復解法である。${\displaystyle -\Delta u=f}$ここで、ヤコビ法（行列演算を回避しようとする方法）は、拡散方程式の対応する線法離散化に対して、明示的なオイラー法（擬似時間）を使用することと同等である。${\displaystyle u_{t}=\Delta u+f}$定常解を求めます。拡散方程式は、スティッフ方程式の典型的な例です。明示的再帰の時間ステップは、CFL条件によって厳しく制限され、空間のメッシュ幅の2乗に比例する必要があります。これにより、ヤコビ反復の収束が非常に遅くなります。以前と同様に、解決策は、前処理法やマルチグリッド法など、何らかの形式の代数方程式を解くことです。

このように、計算効率を高めるために専用の数値解法（通常はかなりの複雑さを伴う）が必要となる数学的問題は数多く存在する。常微分方程式の初期値問題における具体的な例として、硬直方程式が挙げられる。しかしながら、適切に設計された適応型硬直方程式ソルバーを用いることで、ほとんどの硬直方程式の効率的な数値積分は、現在では日常的な作業となっている。

スティッフ方程式の最初の言及は、1952年のカーティスとヒルシュフェルダーの論文に見られる^{[ 4 ]}。著者らは、以下のスカラーモデル方程式について論じている。

${\displaystyle {\dot {x}}={\frac {1}{a(t,x)}}{\big (}xg(t){\big )}}$

そして、剛性を以下のように特徴づける。

もし${\displaystyle \Delta t}$望ましい解像度は${\displaystyle t}$または数値積分で使用される区間の場合、方程式は「硬い」です。

$\displaystyle \left|{\frac {a(t,x)}{\Delta t}}\right|\ll 1}$

そして${\displaystyle g}$行儀が良い。

この特徴付けは時間スケールに関連しています${\displaystyle \Delta t}$過渡現象の減衰率については、係数によって決まります。${\displaystyle a(t,x)}$これは負の値であると想定される。${\displaystyle a(t,x)}$時間と空間で変化するため、それに応じて剛性も変化する。プロセロとロビンソンは、類似しているがより有益なモデル問題を導入した [ ^{5 ] 。}${\displaystyle a(t,x)=1/\lambda <0}$定数と方程式の研究

${\displaystyle {\dot {x}}=\lambda {\big (}xg(t){\big )}+{\dot {g}}(t)\,,\qquad x(0)=x_{0}\neq g(0)\,,\qquad t\in [0,T].}$

解は減衰する過渡現象から構成される$\displaystyle {\big (}x_{0}-g(0){\big )}\,{\mathrm {e} }^{t\lambda }}$特定の解決策とともに${\displaystyle g(t)}$これは有界かつ滑らかであると想定されている。数学的解${\displaystyle x(t)}$は明らかに有界である。明示的オイラー法と暗黙的オイラー法を比較することで、数値解法が数学的解の特性を再現する条件を調べることができる。カーティス・ヒルシュフェルダー基準については、${\displaystyle \Delta t=T}$そして${\displaystyle \Delta t=h}$条件に応じて${\displaystyle T\lambda \ll -1}$そして${\displaystyle h\lambda \ll -1}$、 それぞれ。

微分方程式に対する明示的なオイラー法${\displaystyle {\dot {u}}=f(u)}$定義される

${\displaystyle u_{n+1}=u_{n}+hf(u_{n})\,,}$

暗黙のオイラー法は再帰によって定義される

${\displaystyle u_{n+1}=u_{n}+hf(u_{n+1})\,.}$

ここ${\displaystyle h}$は（一定の）時間ステップであり、$u_n$正確な解を近似する${\displaystyle u(t_{n})}$その時${\displaystyle t_{n}=nh}$明示的な方法で${\displaystyle u_{n+1}}$ベクトル場の評価によって直接計算される${\displaystyle f}$しかし、陰解法では、「代数」方程式を解かなければならない。${\displaystyle u_{n+1}-hf(u_{n+1})=u_{n}}$のために${\displaystyle u_{n+1}}$これにより、再帰処理の安定性が向上します。

プロセロ・ロビンソン問題の場合、明示的なオイラー法は数値解（積分形式）を生成する。

${\displaystyle x_{n}=\left(1+h\lambda \right)^{n}x_{0}+\sum _{k=0}^{n-1}\left(1+h\lambda \right)^{n-k-1}h{\big (}{\dot {g}}(t_{k})-\lambda g(t_{k}){\big )}\,.}$

障害は明らかです。${\displaystyle 1+h\lambda }$安定性条件を課さない限り、指数関数的成長（数値的不安定性）が生じる。${\displaystyle |1+h\lambda |\leq 1}$。 以来${\displaystyle \lambda }$これは現実であり、否定的である。これは、${\displaystyle -2\leq h\lambda \leq 0}$したがって、明示的なオイラー法では、規則性に適応したステップサイズを使用することはできません。${\displaystyle g(t)}$統合を完了するための手順の総数は、${\displaystyle N=T/h\sim T|\lambda |}$、非常に大きくなります。効率は、積で定量化される剛性に反比例します。${\displaystyle T|\lambda |\gg 1}$。

その結果、完全に滑らかな数学的解を持つスカラー問題を解くのに、明示的なオイラー法では「永遠に」かかる可能性がある。${\displaystyle \sim g(t)}$（高速な）過渡現象が減衰すると、初期値が${\displaystyle x(0)=g(0)}$したがって、ステップサイズの制限は過渡現象の有無とは関係ありません。これは、任意のステップサイズが${\displaystyle h}$そのため${\displaystyle |1+h\lambda |>1}$必然的に数値的不安定性を引き起こす。

代わりに陰解法オイラー法を用いると、数値解は次のようになる。

${\displaystyle x_{n}=\left(1-h\lambda \right)^{-n}x_{0}+\sum _{k=0}^{n-1}\left(1-h\lambda \right)^{k-n}h{\big (}{\dot {g}}(t_{k+1})-\lambda g(t_{k+1}){\big )}\,.}$

ここでの安定性要件は${\displaystyle |1-h\lambda |^{-1}\leq 1}$これはすべての人にとって満たされる${\displaystyle h>0}$いつ${\displaystyle {\mathrm {Re} }\,\lambda <0}$どんなに大きくても${\displaystyle h|\lambda |}$そうです。実際には、過渡現象を解決する必要があるという小さな問題があります（これは最初は短いステップサイズを必要とするかもしれません）が、それが減衰したら、「大きな」ステップサイズを使用できます。${\displaystyle h}$規則性に適応している${\displaystyle g(t)}$大きさに関係なく${\displaystyle |h\lambda |}$上記の2つの方法のうち、硬い方程式を効率的に扱うことができるのは、陰解法オイラー法のみです。

例： プロセロ・ロビンソン検定問題を考えてみよう。

${\displaystyle {\dot {x}}\,=\,\lambda (x-\sin \omega t)+\omega \cos \omega t\,;\qquad x(0)=x_{0}\,,}$

正確な解${\displaystyle x(t)=x_{0}\,{\mathrm {e} }^{\lambda t}+\sin \omega t}$. 取る${\displaystyle x_{0}=0}$均一解は厳密解には存在しない${\displaystyle x(t)=\sin \omega t}$しかし、指数項は、点を通過する局所解に現れる。${\displaystyle x_{n}}$数値計算法によって生成された。${\displaystyle \omega =\pi }$問題は解決しました${\displaystyle [0,1]}$使用${\displaystyle N=5}$そして${\displaystyle 10}$それぞれステップ (${\displaystyle h=0.2}$そして${\displaystyle h=0.1}$この実験は、${\displaystyle h\lambda }$これは数値手法の安定性/減衰に影響を与えます。${\displaystyle \lambda =-20}$陽解法と陰解法のオイラー法には同じ計算設定を使用する。

結果を図1に示す。${\displaystyle T\lambda =-20}$局所解が急速に収縮していることから明らかなように、中程度の強い指数関数的減衰が見られます。この設定は視覚化のために意図的に簡略化されています。実際の硬い方程式では、はるかに大きな値が一般的です。それでも、${\displaystyle \lambda }$は、明示的オイラー法のステップサイズ制限を示すのに十分な大きさで負の値である。${\displaystyle N=5}$計算された解は振動し、初期値が厳密解から取られているにもかかわらず、厳密解から乖離する。同じ演習を繰り返すと、${\displaystyle T\lambda =-2000}$例えば、図は、剛性の増加に伴って2つの方法の差がどのように大きくなるかを示しています。

より現実的な問題では、不安定性はより劇的になることが多く、解の「爆発」と呼ばれます。ここで、最大安定ステップサイズ（${\displaystyle h=0.1}$）はわずかな係数でしか上回られていない。一方、陰解法オイラー法は安定性を損なうことなく進行し、より大きなステップサイズでも正確な結果を得ることができる。

実際の計算では、硬い方程式は適応法を用いて解かれます。これらの方法は、計算プロセスを適切に制御するためにステップサイズを自動的に選択し、通常は精度と安定性の両方を制御します。したがって、適応型陽解法を用いる場合、適応機能が安定する最大ステップサイズを超えないように対策を講じるため、解が不安定になることはありません。その代わりに、この方法は（非常に大きな）効率低下という代償を伴います。

分析の本質は、問題パラメータ、${\displaystyle \lambda }$これはベクトル場の特徴であり、解の過渡的挙動に関連しており、時間スケール、すなわちステップサイズと結びついている。${\displaystyle h}$または積分範囲${\displaystyle T}$直接。${\displaystyle h\lambda }$または${\displaystyle T\lambda }$これらは無次元であり、スケーリング不変であり、負の値をとる。一般的な剛性概念はこれらの特性を反映し、非線形微分方程式における剛性も扱えるものでなければならない。

スカラープロセロ・ロビンソン問題は、線形方程式系に拡張することができ、

${\displaystyle {\dot {x}}=A{\big (}x-g(t){\big )}+{\dot {g}}(t)\,,\quad A\in {\mathbb {R} }^{d\times d}}$

非線形システムと同様に、

${\displaystyle {\dot {x}}=f{\big (}x-g(t){\big )}+{\dot {g}}(t)\,,\quad f:{\mathbb {R} }^{d}\rightarrow {\mathbb {R} }^{d}}$

ただし、${\displaystyle f(0)=0}$一般的な剛性概念では、ベクトル場が「大きくて負」であるとはどういう意味かを明確にする必要がある。上記の解析では、滑らかな関数も示されている。${\displaystyle g}$せいぜい小さな役割しか果たさない。したがって、${\displaystyle u}$のために${\displaystyle x-g(t)}$問題には3種類あります。

${\displaystyle {\begin{aligned}{\dot {u}}\,&=\,\lambda u\\{\dot {u}}\,&=\,Au\\{\dot {u}}\,&=\,f(u)\,.\end{aligned}}}$

ここでは、マトリックスの剛性を特性評価する必要がある。${\displaystyle A}$非線形マップの場合も同様${\displaystyle f}$上記の最初の式は、よく知られているダールキスト線形テスト方程式であり、主に時間ステップ法の安定領域を決定するために使用されます。安定領域${\displaystyle S}$メソッドの集合は、すべての${\displaystyle h\lambda \in \mathbb {C} }$数値解法が有界解を生成するようにする。剛性に関するこれまでの議論は、実質的にダールクイストのテスト方程式に基づいている。

2番目の線形方程式系を考察するためには、固有値を利用するのが一般的である。${\displaystyle \lambda _{k}}$行列の${\displaystyle A}$結局、数値解法が有界（安定）解を生成するためには、${\displaystyle {\dot {u}}=Au}$すべての固有値が以下を満たすことが必要かつ十分である${\displaystyle h\lambda _{k}\in S}$前述の議論は、一度に 1 つの固有値を考慮することで継続されます。線形問題は、すべての固有値の実部が負であり、かつ少なくとも 1 つの固有値が次の条件を満たす場合に、スティッフとなります。${\displaystyle T\lambda _{k}\ll -1}$; 後者の固有値は、すべての明示的方法には境界のある安定領域があり、${\displaystyle h\lambda _{k}\in S}$ステップサイズは、この要件を満たすために十分に小さく選択する必要があります（図2を参照）。

しかし、剛性の異なる特徴付けとして、「剛性比」という形で定義されるものがよく見られます。

${\displaystyle {\frac {\max _{k}|{\mathrm {Re} }\,\lambda _{k}|}{\min _{k}|{\mathrm {Re} }\,\lambda _{k}|}}\,,}$

と仮定すると${\displaystyle \lambda _{k}\in {\mathbb {C} }^{-}}$^{[ 6 ]}剛性システムはしばしば大きな剛性比を持つが、ByrneとHindmarshは、これは必要条件でも十分条件でもないと強調している。${\displaystyle {\dot {u}}=Au}$硬い方程式である。^{[ 7 ]}

まず、時間スケールとの比較はできません${\displaystyle T}$しかし、それは「時定数」の本質的な比較に過ぎない。${\displaystyle 1/\lambda _{k}}$ここで、剛性比が大きいということは、単に「速い」過渡成分と「遅い」過渡成分が存在することを意味するだけであり、これらはしばしば大きく異なる時間スケールと呼ばれます。第二に、上述のように、スカラー剛性方程式があります。これらは剛性比を持ちます。${\displaystyle 1}$第三に、特異行列を持つ問題では剛性比が破綻する。${\displaystyle A}$これは、一部のアプリケーションでよく見られる現象です。（例えば、化学反応速度論では、特異点は質量保存則に対応します。）

したがって、剛性比は誤解を招くものであり、剛性方程式の最も重要な側面を特徴づけることができない。

同様に、ヤコビ行列が${\displaystyle f'(u)}$剛性比が大きい。^{[ 8 ]} 当然ながら、アルテミエフとアヴェリナが指摘しているように、これも失敗する。^{[ 9 ]}

例えば、非線形自律常微分方程式系${\displaystyle f(0)=0}$ヤコビ行列の固有値の実部が [sic] の場合、スティッフと呼ぶことができる。${\displaystyle f'(0)}$上記の条件を満たす。しかし、硬い常微分方程式の例としてよく用いられる有名なファン・デル・ポール方程式は、この定義を満たさない。非線形システムの場合、ヤコビ行列の固有値のみで剛性を決定することは不可能である。

この議論も部分的には妥当ではあるものの、ファンデルポール方程式のゼロ解は不安定な平衡点であるため、的外れである。したがって、剛性は原点付近では発生せず、後述するようにリミットサイクルに沿ってのみ発生する。ファンデルポール方程式は、このリミットサイクルにおいて詳細に解析される。そのため、正確な剛性特性評価は局所的である必要があり、実際の解付近における（小さな）摂動の挙動を捉える必要がある。

剛性比の代わりに、現代的なアプローチは機能的な${\displaystyle s_{2}[\cdot ]}$剛性指標と呼ばれる基準${\displaystyle s_{2}[A]\!\cdot \!T\ll -1}$そして${\displaystyle s_{2}[f]\!\cdot \!T\ll -1}$それぞれ、異なるベクトル場特性を時間スケールに関連付けることで剛性を定量化する。${\displaystyle T}$この特徴づけは、追加の観察に基づいている。したがって、DekkerとVerwer ^{[ 10 ]}は次のように述べている（原文強調）。

剛性の本質は、計算対象となる解がゆっくりと変化する一方で、急速に減衰する摂動が存在するという点にある。

シャンパイン^{[ 11 ]}も同様の見解を示しているが、概念的にはより具体的で、剛性を散逸と関連付け、実際には不可逆性とは何かを説明している。

後者の状態を説明するのに私たちが好む方法は、[解]が[時間の]逆方向において非常に不安定である、というものです。

逆不安定性は対数ノルムを用いて容易に特徴づけることができる。したがって、${\displaystyle {\dot {u}}=Au}$、微分不等式を用いて解のノルムを制限できる。${\displaystyle \|u\|_{2}^{2}=u^{*}u}$ユークリッドノルム、場合によっては${\displaystyle u\in {\mathbb {C} }^{d}}$。 それから

${\displaystyle m_{2}[A]\cdot \|u\|_{2}\,\leq \,{\mathrm {D} }_{t}\,\|u\|_{2}\,\leq \,M_{2}[A]\cdot \|u\|_{2}\,,}$

どこ${\displaystyle {\mathrm {D} }_{t}}$は時間微分を表し、${\displaystyle m_{2}[A]}$そして${\displaystyle M_{2}[A]}$は、下側および上側の対数ノルムです。${\displaystyle A}$それぞれ、これらは対称二次形式の極値であり、

${\displaystyle m_{2}[A]\,\leq \,{\frac {u^{*}{\mathrm {He} }(A)u}{u^{*}u}}\,\leq \,M_{2}[A]\,,}$

どこ${\displaystyle {\mathrm {He} }(A)=(A+A^{*})/2}$エルミート部分を表す${\displaystyle A}$その停留点は、対称固有値問題から求められる。

${\displaystyle {\mathrm {He} }(A)v=\mu v\,,}$

だれの${\displaystyle d}$実固有値${\displaystyle \mu _{1}\geq \mu _{2}\geq \dots \geq \mu _{d}}$これらは対数値と呼ばれます${\displaystyle A}$固有値に関して${\displaystyle \lambda _{k}}$の${\displaystyle A}$対数値は以下を満たす

${\displaystyle m_{2}[A]\,\equiv \,\mu _{d}\,\leq \,\min _{k}{\mathrm {Re} }\,\lambda _{k}\,;\qquad \max _{k}{\mathrm {Re} }\,\lambda _{k}\,\leq \,\mu _{1}\,\equiv \,M_{2}[A]\,,}$

トレースIDと共に

${\displaystyle {\mathrm {Re} }\,{\mathrm {trace} }\,A\,=\,\sum _{k}{\mathrm {Re} }\,\lambda _{k}\,=\,\sum _{k}\mu _{k}\,=\,{\mathrm {trace} }\,{\mathrm {He} }(A)\,.}$

以来${\displaystyle u(t)={\mathrm {e} }^{tA}u_{0}}$上記の微分不等式から次のことが導かれる。

${\displaystyle {\mathrm {e} }^{tm_{2}[A]}\,\leq \,\|{\mathrm {e} }^{tA}\|_{2}\,\leq \,{\mathrm {e} }^{tM_{2}[A]}\,\qquad t\geq 0\,,}$

システムの最大成長率と最大減衰率を制限する。したがって、長さの時間間隔にわたって${\displaystyle T>0}$最大可能な前方時間成長率は${\displaystyle {\mathrm {e} }^{TM_{2}[A]}}$最大逆時間成長率は${\displaystyle {\mathrm {e} }^{-Tm_{2}[A]}}$。

Söderlind ら^{[ 12 ]}は、Shampine の観察を定量化しました。したがって、硬いシステムでは、

${\displaystyle {\mathrm {e} }^{-Tm_{2}[A]}\,\gg \,{\mathrm {e} }^{TM_{2}[A]}\,,}$

これから次のことが導かれる${\displaystyle {\mathrm {e} }^{T{\big (}m_{2}[A]+M_{2}[A]{\big )}}\ll 1}$つまり、${\displaystyle T{\big (}m_{2}[A]+M_{2}[A]{\big )}\ll -1}$剛性指標は^{[ 13 ]}で定義されている。

${\displaystyle s_{2}[A]\,=\,{\frac {m_{2}[A]+M_{2}[A]}{2}}\,,}$

硬い方程式は、次の条件によって特徴付けられる。

${\displaystyle s_{2}[A]\cdot T\ll -1\,.}$

含める理由${\displaystyle M_{2}[A]}$剛性指標の定義において、非線形システムではしばしば次のようなことが起こる可能性がある。${\displaystyle M_{2}[A]}$は正であり、剛性を相殺します。これは、剛性指標を奇パリティを持つように定義することによって考慮されます。つまり、${\displaystyle s_{2}[-A]=-s_{2}[A]}$。

シャンパインの基準は、${\displaystyle m_{2}[A]}$は大きく負の値であり、つまり、の下限は${\displaystyle \|{\mathrm {e} }^{tA}\|_{2}}$は極めて小さい。言い換えれば、流れはほぼ半群であり、${\displaystyle {\mathrm {e} }^{tA}}$これは「特異点に近い」状態である。リプシッツ系では、（理論的には）逆時間でも微分方程式を解くことができるが、スティッフ系ではこれは事実上不可能である。逆時間スティッフ方程式は、極めて条件の悪い「逆問題」に対応する。

非線形写像に対数ノルムを用いると^{[ 14 ]、} 同じ議論が非線形ベクトル場にも適用できる。${\displaystyle f}$しかし、${\displaystyle s_{2}[f]}$はグローバルに定義され、代わりに剛性指標はローカルに定義される。${\displaystyle s_{2}[f'(u)]}$軌道に沿って。

剛性指標に加えて、Söderlindら^{[ 15 ]}は局所参照時間スケールを導入している。${\displaystyle \Delta t(u)}$定義される

${\displaystyle \Delta t(u)\,=\,1/\max {\big (}1/T,-s_{2}[f'(u)]{\big )}\,.}$

剛性は、方法に依存しない条件によって特徴付けられる。${\displaystyle T/\Delta t(u)\gg 1}$（局所）剛性係数と呼ばれる。剛性係数が大きいことは剛性の必要条件である。しかし、方法依存の基準も評価できる。数値計算では、${\displaystyle \Delta t(u)}$実際のステップサイズと比較できます${\displaystyle h}$局所的なステップサイズ剛性係数${\displaystyle h/\Delta t(u)}$積分方法の選択と精度要件によります。陰解法のみがステップサイズを使用できます。${\displaystyle h/\Delta t(u)\gg 1}$しかし、積分過程の一部では、ステップサイズ剛性係数が中程度にとどまる場合もある。

暗黙的時間ステップ法の重要な側面は、次の形式の方程式を解く必要があることです。

${\displaystyle u=\gamma hf(u)+\psi }$

のために${\displaystyle u}$あらゆる段階で。${\displaystyle \gamma >0}$は中程度の大きさの方法特性定数であり、${\displaystyle \psi }$は既知のベクトルです。この方程式が（一意の）解を持つかどうかという疑問が生じます。${\displaystyle hf}$「大きい」とは、${\displaystyle hs_{2}[f']\ll -1}$これは、計算が困難であることを示しています。方程式を次のように書き換えます。

${\displaystyle {\big (}\gamma hf-I{\big )}(u)=-\psi \,,}$

一様単調性定理が適用される：写像が${\displaystyle \gamma hf-I}$は（負の）単調である、つまり、${\displaystyle M_{2}[\gamma hf-I]<0}$。 以来

${\displaystyle M_{2}[\gamma hf-I]<0\quad \Leftrightarrow \quad M_{2}[\gamma hf]<1\,,}$

この条件は、厳しい計算では通常、かなりの余裕をもって満たされます。${\displaystyle M_{2}[f]<0}$全員にとって満足のいくものです${\displaystyle h>0}$しかし、もし${\displaystyle M_{2}[f]>0}$そのため、わずかなステップサイズの制限が必要になる場合があります。シャンパイン基準および剛性指標では、次のことが成り立ちます。

${\displaystyle {\frac {m_{2}[\gamma hf]+M_{2}[\gamma hf]}{2}}\,=\,\gamma h\!\cdot \!s_{2}[f]\ll -1\,,}$

どこ${\displaystyle m_{2}[\gamma hf]\ll -1}$、そして${\displaystyle M_{2}[\gamma hf]}$正の場合は中程度である。したがって、この方程式は通常、硬い計算ではニュートン反復を必要とするが（${\displaystyle \gamma hf}$ is not a contraction), existence and uniqueness are usually guaranteed also when the step size stiffness factor is large.

Because it is relatively expensive to compute the stiffness indicator by solving the symmetric eigenvalue problem ${\displaystyle \,{\mathrm {He} }{\big (}f'(u){\big )}v=\mu v\,}$ for ${\displaystyle \mu _{1}}$ and ${\displaystyle \mu _{d}}$, an inexpensive alternative is of interest. Noting that the stiffness indicator is the arithmetic average of the largest and the smallest logarithmic value, an option is to replace this average by the average of all logarithmic values. With the aforementioned trace identity (see above), we define the scaled (arithmetic average) real trace functional^[16]

${\displaystyle \tau [A]\,=\,{\frac {1}{d}}\,{\mathrm {Re} }\,{\mathrm {trace} }(A)\,=\,{\frac {1}{d}}\,\sum _{k=1}^{d}{\mathrm {Re} }\,a_{kk}\,=\,{\frac {1}{d}}\,\sum _{k=1}^{d}{\mathrm {Re} }\,\lambda _{k}\,=\,{\frac {1}{d}}\,\sum _{k=1}^{d}\mu _{k}\,,}$

noting that the computation of ${\displaystyle \tau [A]}$ is both inexpensive and trivial: no eigenvalue computations are required. The scaled trace (as well as the divergence) is a linear functional, hence ${\displaystyle \tau [-A]=-\tau [A]}$, sharing the odd parity of ${\displaystyle s_{2}[A]}$. Moreover, for a nonlinear map,

${\displaystyle \tau [f'(u)]\,=\,{\frac {1}{d}}\,{\mathrm {trace} }\,f'(u)\,=\,{\frac {1}{d}}\,{\mathrm {trace} }({\mathrm {grad} }_{u}\,f)(u)\,=\,{\frac {1}{d}}\,({\mathrm {div} }_{u}\,f)(u)\,.}$

The functional ${\displaystyle \,\tau [f'(u)]\,}$ serves as an alternative stiffness indicator. Stiffness can therefore also be characterized by the divergence condition

${\displaystyle {\frac {T}{d}}\,{\big (}{\mathrm {div} }_{u}\,f{\big )}(u)\,\ll \,-1\,,}$

which is readily evaluated for most problems.

For a linear system ${\displaystyle {\dot {u}}=Au}$, it holds that

${\displaystyle {\mathrm {D} }_{t}\,{\mathrm {det} }\,{\mathrm {e} }^{tA}\,=\,{\mathrm {trace} }\,A\cdot {\mathrm {det} }\,{\mathrm {e} }^{tA}\,.}$

Defining a scaled (geometric average) absolute determinant ${\displaystyle \,\delta [A]=|{\mathrm {det} }\,A|^{1/d}\,}$ we have ${\displaystyle \,\delta [I]=1\,}$ and ${\displaystyle \,\delta [\alpha A]=|\alpha |\!\cdot \!\delta [A]\,}$, avoiding the standard determinant's direct dependence on the dimension ${\displaystyle d}$. (This is of special importance when the differential equation is derived from a method-of-lines discretization of a partial differential equation.) Here ${\displaystyle \,\delta [A]\,}$ is the linear "size" of the phase space volume spanned by the column vectors of ${\displaystyle A}$, such as in ${\displaystyle \,\delta [2I]=2\,}$, independent of ${\displaystyle d}$. The differential equation for the determinant of the flow can now be rewritten

${\displaystyle {\mathrm {D} }_{t}\,\delta {\big [}{\mathrm {e} }^{tA}{\big ]}\,=\,\tau [A]\cdot \delta {\big [}{\mathrm {e} }^{tA}{\big ]}\,.}$

Hence

${\displaystyle \delta {\big [}{\mathrm {e} }^{tA}{\big ]}\,=\,{\mathrm {e} }^{t\,\tau [A]}\,,}$

expressing that an initial (linear) phase volume changes by a factor ${\displaystyle \,\delta {\big [}{\mathrm {e} }^{tA}{\big ]}\,}$ over a time ${\displaystyle t}$. If ${\displaystyle \,\tau [A]<0\,}$, phase volume is compressed by the flow ${\displaystyle \,{\mathrm {e} }^{tA}\,}$ for ${\displaystyle t\geq 0}$.

It follows that ${\displaystyle \,\tau [A]\leq 0\,}$ is a necessary condition for stability. (Compare the sufficient criterion ${\displaystyle \,M_{2}[A]\leq 0\,}$ for the stability of the zero solution.) Shampine's characterization of stiffness, in terms of strongly unstable solutions in reverse time, can be written as ${\displaystyle \,T\!\cdot \tau [-A]\gg 1\,}$, which, due to the odd parity of the trace, is equivalent to the stiffness condition ${\displaystyle T\!\cdot \tau [A]\ll -1}$. These observations, concepts and constructions provide the key to characterizing and quantifying stiffness, with the divergence criterion ${\displaystyle \,T\!\cdot \tau [f']\ll -1\,}$ being the simplest.

A large number of nontrivial test problems have been collected in the Bari Test Set for Initial Value Problem Solvers.^[17] This contains detailed examples of problems, solvers, and results achieved under specified conditions by well-known codes.

ここでは、剛性指標（簡略化されたもの）を用いて剛性を正しく識別・特徴付ける方法を示すために、3つの剛性の高い非線形問題を取り上げます。これには、ファン・デル・ポール方程式やオレゴネーター方程式に見られるような、（局所的に不安定な）転換点における複雑で急速に変化する挙動を区別する能力が含まれます。さらに、剛性指標は、剛性が顕著な部分区間と、問題が非剛性である部分区間を区別することもできます。最後に、場合によっては、どの従属変数が剛性に最も寄与しているかを特定することも可能です。

適応時間ステップが精度要件にどのように依存し、解に沿ってどのように変化するかといった運用基準も、ステップサイズ剛性係数を定量化しながら説明されています。これらを総合すると、剛性指標から詳細な洞察が得られ、その一部は実際の計算が始まる前に得られます。

問題1.（火炎伝播） スカラー火炎伝播モデル^{[ 18 ]}について考える。${\displaystyle t\in [0,1]}$、

${\displaystyle \varepsilon \,{\dot {u}}\,=\,2\,(u^{2}-u^{3})\,;\qquad u(0)=\varepsilon \ll 1\,.}$

右辺の2つの項は、炎内部の燃焼には炎の表面を通して供給される酸素が必要であるという事実を反映しています。${\displaystyle u}$表面積を増やす${\displaystyle \sim u^{2}}$最終的には、より大きな体積での燃焼を維持するには小さすぎるものになる。${\displaystyle \sim u^{3}}$この簡略化されたモデルは、硬い方程式に対するソフトウェアの時間ステップ適応性を検証するためのテスト問題としてよく使用され、非硬い領域と硬い領域の間の急激な遷移が正しく解決されているかどうかを確認します。

のために${\displaystyle u\in (0,1)}$右辺は正であり、解は${\displaystyle u(t)}$は単調増加です。${\displaystyle u(0)=\varepsilon }$不安定平衡点から離れていく${\displaystyle u=0}$ベクトル場（剛性指標）の発散は単純に${\displaystyle f'(u)=2\,(2u-3u^{2})/\varepsilon }$。 以来

${\displaystyle f'(\epsilon )\approx 4\,;\qquad f'(1)=-\,{\frac {2}{\varepsilon }}\ll -1\,,}$

解は最初は非剛性です。${\displaystyle u\ll 1}$。 で${\displaystyle t\approx 0.5}$いつ${\displaystyle u\sim 1/2}$不安定性は破滅的なものとなり、第二の安定平衡状態への極めて急速な移行が起こっている。${\displaystyle u=1}$。

いつ${\displaystyle u>2/3}$安定性が回復し、解はスティッフ領域に入ります（図3参照）。ここでは、${\displaystyle \varepsilon =10^{-3}}$方程式は、小さいほどより厳格になる。${\displaystyle \,\varepsilon \,}$が選ばれ、${\displaystyle \,u=1\,}$解決策が移動する。

このデモンストレーションでは、3次ルンゲ・クッタ法と2次誤差推定器、および3次A安定陰的ルンゲ・クッタ法と2次誤差推定器を比較します。陰的方法は非剛性領域では利点がありませんが、剛性領域ではより大きなステップを使用できます。2つの方法の差は、ステップが小さいほど大きくなります。${\displaystyle \,\varepsilon }$。

問題2. (ファンデルポール方程式) ファンデルポール方程式は、通常、異なるスケーリングを用いたいくつかの代替形式で表されます。これは、解が極限サイクルに近づく非線形問題です。ここでは、時間をスケーリングしてシステムを次のように記述します。

${\displaystyle {\begin{aligned}{\dot {x}}&=2\kappa \cdot y\\{\dot {y}}&=2\kappa ^{2}\cdot (1-x^{2})\,y-2\kappa \cdot x\,,\end{aligned}}}$

初期条件付き${\displaystyle x(0)=2}$、${\displaystyle y(0)=0}$リミットサイクル上で選択された。時間スケールは、リミットサイクルの周期が${\displaystyle \mathrm {O} (1)}$パラメータにほとんど依存しない${\displaystyle \kappa }$問題は次のように解決できます。${\displaystyle [0,1]}$これは、全周期にわたる現象である。この問題は、電気回路における非線形振動の研究に端を発する。

大きな値の場合${\displaystyle \kappa }$（図4を参照）${\displaystyle \kappa =200}$）、リミットサイクルは、安定性の壊滅的な喪失を示すほぼ不連続な遷移によって接続された2つの安定な分岐から構成されます。トレースの同一性により、剛性指標は${\displaystyle s_{2}[f'(u)]}$スケールされたトレースと一致する${\displaystyle \tau [f'(u)]}$によって与えられた

${\displaystyle s_{2}[f'(u)]\,=\,\tau [f'(u)]\,=\,\kappa ^{2}\,(1-x^{2})\,,}$

どこ${\displaystyle u}$ベクトル${\displaystyle (x,y)^{\mathrm {T} }}$剛性指標は、ファンデルポール系の2番目の式の中に、一見して分かりやすい形で隠されている。

スケーリングされたトレースには定数項があり、${\displaystyle \kappa ^{2}}$これはシステムの線形部分からの寄与を表します。この項は正であるため、不安定化します。平衡${\displaystyle x=y=0}$不安定です。${\displaystyle \tau [f'(0)]=\kappa ^{2}>0}$スケーリングされたトレースには負の項も含まれる。${\displaystyle -\kappa ^{2}x^{2}}$剛性に対する非線形寄与を表す。スケーリングされたトレースは${\displaystyle y}$後者の変数は、剛性に特に影響を与えない。

解が硬い領域に入るのは、${\displaystyle \tau [f'(u)]\ll -1}$。 以来${\displaystyle \tau [f'(u)]=\kappa ^{2}\,(1-x^{2})}$したがって、剛性は${\displaystyle \,|x|>1}$すなわち、2 つの安定枝に沿って、剛性は${\displaystyle \,\kappa ^{2}}$、そして、任意に大きくなる可能性があります。特に、安定ブランチに沿って問題を解くために明示的な方法を使用する場合、計算量は次のようになります。${\displaystyle \mathrm {O} (\kappa ^{2})}$対照的に、努力は${\displaystyle \mathrm {O} (1)}$硬いソルバーの場合。^{[ 19 ]}

なお、これらの情報はすべて事前に入手可能であり、問​​題の数値解によって確認されています。

問題3.（オレゴネーター方程式）化学反応速度論は、硬い微分方程式の宝庫である。オレゴネーター方程式は、3つの反応を含む自己触媒反応をモデル化する。適切な初期条件の下では、周期の極限サイクルが存在する。${\displaystyle T\approx 305}$方程式は次のとおりです。

${\displaystyle {\begin{aligned}{\dot {x}}\,&=\,s\cdot (x-xy+y-q\,x^{2})\\{\dot {y}}\,&=\,(z-y-xy)/s\\{\dot {z}}\,&=\,w\cdot (x-z)\,.\end{aligned}}}$

テスト問題では、以下のパラメータを使用します。${\displaystyle s=77.27}$、${\displaystyle q=8.375\cdot 10^{-6}}$、${\displaystyle w=0.161}$初期条件とともに${\displaystyle x(0)=3.2,\,y(0)=1.45,\,z(0)=2.5}$[ ^{20 ]​}

システムの次元は${\displaystyle d=3}$剛性指標${\displaystyle s_{2}[f']}$そしてスケーリングされたトレース（発散）${\displaystyle \tau [f'(u)}$両者は一致しないが、その差は小さい。さらに重要なのは、この方程式は剛性のより複雑な例であるため、解析においてパラメータをできるだけ長く保持することで予備的な洞察を得る必要があるということである。そこで、この利点を提供するスケーリングされた発散を選択する。したがって、次のことがわかる。

${\displaystyle T\cdot \tau [f'(u)]\,=\,\alpha +\beta \,x+\gamma \,y,}$

どこ

${\displaystyle {\begin{aligned}\alpha \,&=\,T\cdot (s^{2}-sw-1)/(3s)\,\approx \,7.838\cdot 10^{3}\\\beta \,&=\,-T\cdot (2qs^{2}+1)/(3s)\,\approx \,-1.447\\\gamma \,&=\,-T\cdot s/3\,\approx \,-7.856\cdot 10^{3}.\end{aligned}}}$

ここで定数項${\displaystyle \alpha }$は発散への線形寄与です。これは正の値なので不安定化します。他の 2 つの項は非線形寄与を表し、${\displaystyle x}$そして${\displaystyle y}$それらの係数${\displaystyle \beta }$そして${\displaystyle \gamma }$負の値であり、剛性に寄与します。しかし、発散は${\displaystyle z}$, the latter variable does not contribute to stiffness other than to the constant term ${\displaystyle \alpha }$. The largest coefficient is ${\displaystyle \gamma }$, indicating that stiffness will be particularly strong when ${\displaystyle y}$ is large and positive. Whether ${\displaystyle x}$ will contribute depends on what magnitude it will reach. In Figure 5, the actual data are collected during one full period ${\displaystyle [0,305]}$ plotted on the normalized interval ${\displaystyle [0,1]}$.

During a brief interval when ${\displaystyle x}$ is large, the scaled divergence is of order ${\displaystyle T\cdot \tau [f'(u)]\approx -10^{5}}$ (just barely visible as a notch in the upper left corner of the high resolution scaled trace graph), making the problem stiff there. But the worse is yet to come; when ${\displaystyle y}$ becomes large stiffness becomes severe, with scaled divergence at ${\displaystyle T\cdot \tau [f'(u)]\approx -10^{7}}$.

When an adaptive third order A-stable Runge-Kutta method is used, the step sizes exceed the reference time scale by three to four orders of magnitude. This corresponds to the efficiency gain of the stiff solver over the potential use of an explicit method.

As noted above, stiff equations require implicit time stepping methods, either Runge-Kutta methods or linear multistep methods. There are also alternatives based on extrapolation techniques. But not all implicit methods are suitable. A good method must have a stability region ${\displaystyle S}$ covering all or most of the negative half-plane ${\displaystyle {\mathbb {C} }^{-}}$.

The stability region is determined by applying the method to the Dahlquist test equation^[21]${\displaystyle \,{\dot {u}}=\lambda u}$, and the stability region ${\displaystyle S}$ is the set of ${\displaystyle h\lambda \in {\mathbb {C} }}$ for which the method produces bounded solutions. A method with ${\displaystyle S\supset {\mathbb {C} }^{-}}$, i.e., where ${\displaystyle S}$ contains the entire left half-plane, is called A-stable.

There are many A-stable Runge-Kutta methods of high order to choose from, but unfortunately, for A-stable linear multistep methods the convergence order is limited to ${\displaystyle p=2}$. Therefore one usually has to settle for less. Among multistep methods the best choice for stiff equations are the backward differentiation (BDF) methods, to which the implicit Euler method belongs. There are BDF methods of convergence orders ${\displaystyle p=1:6}$, but only orders ${\displaystyle p\leq 2}$ are A-stable. While the stability region remains large for higher orders, it eventually deteriorates, and only methods up to order ${\displaystyle p=5}$ are regularly used.

Software based on BDF-type methods (or similar) are Matlab's ode15s, and C or Fortran codes such as CVODE, LSODE, MEBDF, DASSL. The software is elaborate, robust and of high complexity, uses variable order and variable time step adaptivity, and usually offers many options for problems with special structure, or for special handling of Jacobian matrices. Accuracy criteria can be specified both for absolute and relative error estimates. Some codes include options for nonstiff equations, or for differential-algebraic equations.

For implicit Runge-Kutta methods applied to the test equation, the differential equation is replaced by a recursion ${\displaystyle u_{n+1}=R(h\lambda )\,u_{n}}$ここで、有理関数${\displaystyle R}$は安定性関数と呼ばれます。

${\displaystyle {\mathrm {Re} }\,z\leq 0\,\Rightarrow \,|R(z)|\leq 1\,,}$

この方法はA安定です。高次のA安定法はありますが、実装例は少ないです。A安定ルンゲ・クッタ法の利点は、

${\displaystyle M_{2}[hA]\leq 0\,\Rightarrow \,\|R(hA)\|_{2}\leq 1\,,}$

つまり、もし${\displaystyle A}$が負定値である場合、この方法は縮小再帰を生成します。（これはフォン・ノイマンの不等式の変形です。）残念ながら、これは追加の条件の下でのみ非線形問題に適用できます。B安定ルンゲ・クッタ法はA安定法のサブセットであり、次のようになります。${\displaystyle \,M_{2}[hf]\leq 0\,}$非線形問題の場合、この方法は縮小再帰を生成します。^{[ 22 ] [ 23 ] [ 24 ]} L安定性などの他の特別な安定性要件も一般的です。^{[ 25 ]}これは、追加の要件を持つA安定法の別のサブセットです。${\displaystyle \,R(\infty )=0\,}$固有モードの減衰を改善するために${\displaystyle \,{\mathrm {Re} }\,h\lambda \ll -1}$。

暗黙的ルンゲ・クッタ法は計算複雑度が高いため、ソフトウェア開発時には特別な効率要件も考慮する必要がある。硬い問題に対する効率的なルンゲ・クッタ法ソフトウェアとしては、次数 10 の Matlab ode23 が挙げられる。${\displaystyle 2}$そして${\displaystyle 3}$また、 B-およびL-安定を実装したFortranコードRADAU5 ^{[ 26 ]も存在する。}${\displaystyle 5^{\mathrm {th} }}$ラダウ IIa 法の次数。硬い方程式に使用される対角陰解法ルンゲ・クッタ法の一般的な概説については、ケネディとカーペンターを参照してください。^{[ 27 ]}

いずれの場合も、硬方程式の数値解法には専用の専門ソフトウェアが必要です。これは日常的な作業ではありますが、成功にはコードの適切な設定を綿密に検討し、解決すべき問題を正しく理解することがしばしば求められます。

1. 剛性は完全に理解されているのか？文献では、剛性の正確な定義は存在しないと示唆されることがある。これは、剛性比が広く言及または使用されているためと思われるが、その明らかな欠点はずっと以前から認識されている。実際の離散化方法や精度要件などの運用基準が問題解決のために持ち込まれることもある。しかし、これは基本的に単純な問題を過度に複雑化させている。したがって、実務家にとって、非剛性方程式と剛性方程式の区別があることは明らかであり、この区別を数学的に記述できる必要がある。剛性指標は、この区別を特徴付け、定量化するためのシンプルで必要な基準を提供する。今日では、剛性は複雑ではあるが完全に理解された現象であり、優れた効率的で信頼性の高い専用ソフトウェアが利用可能であると言っても差し支えない。

2. 剛性インジケータ。剛性インジケータ${\displaystyle s_{2}[f']}$スケーリングされた発散よりも頑健である${\displaystyle \tau [f']}$前者は両極端の対数値のみを使用するのに対し、後者は安定性に同様の影響を与えないにもかかわらず、中間の対数値も使用します。スケーリングされた発散の利点は、上記の例で示されているように、固有値の計算を回避し、ダイナミクスの事前の解析的理解をサポートすることです。剛性は通常桁違いであるため、おおよその定量化で十分です。

3. 散逸系と保存系。上記で説明した剛性の概念は、散逸性、つまり減衰またはエネルギー損失の尺度です。いくつかの文献では、特に双曲型偏微分方程式の線法離散化から導出された場合には、剛性のある保存系も存在すると示唆しています。しかし、そのようなシステムは、（潜在的に）非常に振動的なシステムという別のカテゴリに属します。固有値が虚軸の遠くに位置する場合、大きなステップを使用する陰解法を選択することで、対応する高周波を抑制できます。ただし、サンプリング定理に従って、高周波現象はもはや解像されず、波形が歪む可能性があります。散逸問題と保存問題の区別は重要であり、計算アプローチが異なります。剛性は放物型問題の不可逆性と密接に関連していますが、双曲型問題には固有の減衰がなく、エネルギー保存などの不変量があります。同様に、分離可能なハミルトン系は発散のないベクトル場を持つため、「面積保存性」（位相体積の保存）を持つ。この構造のため、剛性指標はゼロとなる。

4. 語源。「stiff」という用語は、カーティスとヒルシュフェルダーによって導入されました。ヒルシュフェルダー^{[ 28 ]}によると、この用語が選ばれたのは、最初の例がサーボシステムに関連しており、サーボと制御システムの間に「密結合」があったためです。同様の効果は、高ゲイン負帰還制御システムでも見られ、ゲインが高いほど剛性が高まります。もう1つの関連性として、2次方程式で表される（機械的な）剛性の概念が挙げられます。

${\displaystyle m{\ddot {x}}+c{\dot {x}}+kx=F\,,}$

どこ${\displaystyle m}$質量を表す、${\displaystyle c}$減衰係数、${\displaystyle k}$ばね定数、そして${\displaystyle F}$外部から加えられた力。その考え方は、「硬いバネ」（大きな）を含む方程式です。${\displaystyle k}$数学的な意味での硬直を引き起こす。${\displaystyle m=1}$方程式は次のように書き換えられます。

${\displaystyle {\begin{aligned}{\dot {x}}&=y\\{\dot {y}}&=-kx-cy+F\,,\end{aligned}}}$

これは次の形式のシステムです${\displaystyle {\dot {u}}=Au+g}$の固有値は${\displaystyle A}$は

${\displaystyle \lambda _{1,2}\,=\,{\frac {-c\pm {\sqrt {c^{2}-4k}}}{2}}\,;}$

これらは大きいと主張されています${\displaystyle k}$は大きい。臨界減衰の場合、パラメータは以下を満たす必要がある。${\displaystyle c=2{\sqrt {mk}}}$と${\displaystyle m=1}$につながる

${\displaystyle \lambda _{1,2}\,=\,-c/2\,.}$

興味深いことに、ベクトル場の剛性指標（ここではスケーリングされた発散と同一）を計算すると、次の式が得られる。

${\displaystyle s_{2}[A]\,=\,\tau [A]\,=\,-c/2\,,}$

これは、減衰定数が${\displaystyle c}$は大きく、ばね定数とは無関係である。したがって、剛性の唯一の原因は、エネルギーを散逸させるダンパーである。 「剛性」という用語は誤称であり、すぐに別の意味で新しい文脈で定着したという結論を避けることは難しい。機械工学では、大きなばね定数は通常、臨界減衰に近い値を得るために大きな減衰定数と組み合わされるが、この混乱は理解できる。

• 後退微分公式は、特に硬い微分方程式の解法に用いられる陰解法の一種である。
• 状態番号
• 微分包含は、微分方程式の概念を拡張したもので、不連続性を許容し、剛性の問題を回避する手段の一つとなる。
• 明示的メソッドと暗黙的メソッド

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