数値計算をやっていると「スティフな問題」というものに出くわすことがある。
ではスティフは何なのか。
手元にある本でラルストン・ラビノヴィッツ著 “電子計算機のための数値解析の理論と応用”から引用すると
その物理的システムにおいて，意味のある異なった時間の尺度（または時定数）をもった微分方程式で表されるプロセスが存在する。これに対して通常の数値解法は不適当である，なぜならば，正確な解のそれぞれ対応する成分が徐々に消えてしまった後であっても，時定数の大きいものが残っているので最も速いプロセスを考えなければならないからである。
とある。日本語訳が分かりづらいが、桁が異なる時定数を持つ複数のプロセスを連立方程式系で解く場合に起きる。 雑な言い方をすると、時間刻みを細かくしないと不安定になって(発散したりして)計算できない場合をスティフと言う。
一般的に陽的な解法はスティフな問題が苦手で、陰的な解法を選ぶ。一般的に陰的な解法は計算コストが高いためスティフな問題は計算時間がかかる。
似たような記述がMATLABのヘルプドキュメントにもある: (ODE ソルバーの選択 | 基本的なソルバー選択)
スティッフとは明確な定義が難しい用語ですが、一般的には、問題中のいずれかの部分でスケールに差があるときにスティッフになります。たとえば、時間スケールに大幅な違いがあると変化する解の要素が ODE に 2 つある場合、その方程式はスティッフな可能性があります。
カーネギー・メロン大学のLecture Note An Introduction to Physically Based Modeling に “Energy Functions and Stiffness“というchapterがある。この中からも引用しておく:
Note that when we stretch the energy function, we change the time constants on the true solutions to the differential equations. In the symmetric case, the solution is $x = x_0 e^{−t}$, $y = y_0e^{−t}$. In the stretched case, the solution is $x = x_0 e^{−\varepsilon t}$, $y = y_0e^{−t}$. The time constants for the symmetric case are the same, but in the stretched case, the time constants have a ratio of $\varepsilon$. This is the usual root of “stiffness” in differential equations. As we saw, the maximum step size is set by the process with the fastest time constant, and this can make simulation of the longer time constant unbearably slow.
In order to deal with the difficulties of stiffness, we have two main choices. We can either use a method based on a more accurate local model of the energy landscape, or we can try to formulate a different but related set of differential equations capable of solving our original problem. Reformulating the problem to eliminate the stiffness is preferred when possible, but very often there is no choice besides simulating the stiff equation.
機械訳すると
エネルギー関数を引き伸ばすと、微分方程式の真の解の時定数が変化することに注意しよう。 対称の場合、解は$x = x_0e^{-t}$, $y = y_0e^{-t}$となる。 引き伸ばした場合、解は$x = x_0e^{-\varepsilon t}$, $y = y_0 e^{-t}$となる。 対称の場合、時定数は同じだが、引き伸ばした場合は時定数の比が$\varepsilon$になる。 通常、これが微分方程式における「硬さ」の根源だ。 これまで見てきたように、最大ステップサイズは時定数が最も速いプロセスで設定されるため、時定数が長い方のシミュレーションが耐えられないほど遅くなることがある。
スティフネスの難しさに対処するためには、主に2つの選択肢がある。 エネルギーランドスケープのより正確なローカルモデルに基づく方法を使うか、元の問題を解くことができる、異なるが関連する微分方程式のセットを定式化することを試みるかである。 可能であれば、剛性を排除するために問題を再定式化することが望ましいが、硬い方程式をシミュレーションする以外に選択肢がない場合が非常に多い。
“電子計算機のための数値解析の理論と応用”に戻ると、この本の上記の記述の前では一つの方程式での問題も触れている。これがWikipediaで書かれているような例で、刻み幅の変更で計算結果が大きく変わってしまうような不安定な例だ。
