硬い方程式のソースを表示

数学・[[数値解析]]において硬い方程式（{{Lang-en-short|stiff equation}}）は、[[常微分方程式の数値解法]]・[[偏微分方程式の数値解法]]において、刻み幅を極めて小さくしない限り、[[数値的安定性|数値的不安定]]になる[[微分方程式]]である。硬さを的確に定義するのが困難であると判明したが、方程式に解の急激な変化を起こせる項が含まれていることは確かである。
<!-- When integrating a differential equation numerically, one would expect the requisite step size to be relatively small in a region where the [[solution curve]] displays much variation and to be relatively large where the solution curve straightens out to approach a line with slope nearly zero.  For some problems this is not the case.  Sometimes the step size is forced down to an unacceptably small level in a region where the solution curve is very smooth.  The phenomenon being exhibited here is known as '''stiffness'''.  In some cases we may have two different problems with the same solution, yet problem one is ''not'' stiff and problem two ''is'' stiff.  Clearly the phenomenon cannot be a property of the exact solution, since this is the same for both problems, and must  be a property of the differential system itself.  It is thus appropriate to speak of ''stiff systems''.

微分方程式の解を数値的に計算するとき、解曲線が急激に変化する領域に対し必要な刻み幅が比較的に小さく、緩やかに変化する領域に対し必要な刻み幅が比較的に大きいであると期待できる。しかし、ある問題に対してそれは正しくない。これらの問題に対し、たとえ解曲線がとても滑らかな領域の中でも時々刻み幅を許容できないほど小さくしなければならない。このような現象は、硬さとして知られている。また、異なる問題が同じ解を持つ場合に、一つの問題が硬く、もう一つの問題が硬くないこともある。そのため、硬さは明らかに方程式の解の性質ではなく、方程式自体の性質である。したがって、硬い方程式について語るのが適切である。 -->

== 導入の例 ==
[[Image:StiffEquationNumericalSolvers.svg|thumb|500px|right|硬い方程式に対する異なる数値解の安定性]]

下記の[[初期値問題]]を考える。
: <math>y'=-15y,\quad t \ge 0, \quad y(0)=1.</math>
この問題は直接に解くことができ、厳密解 (水色の曲線) が次の公式で与えられる。
: <math>y(t)=e^{-15t}.</math>
公式によって、<math> \lim_{t \to \infty} y(t) = 0 </math> も明らかである。

同じ振舞いを持つ数値解を求めよう。様々な数値的方法を用いて得られる数値解は右側の画像に表示される。 

{{ordered list
|1= 刻み幅 {{math|''h'' {{=}} 1/4}} の[[オイラー法]] に対応する解（赤色の曲線）は激しく振動し、素早くグラフの範囲を超える。
|2= 刻み幅 {{math|''h'' {{=}} 1/8}} の[[オイラー法]] に対応する解（緑色の曲線）も振動するが、グラフの範囲にいる。
|3= 刻み幅 {{math|''h'' {{=}} 1/8}} の[[台形公式]]に対応する解（青色の曲線、凡例ではAdams-Moulton法と名付けられたが同じ方法である）は振動せず、期待通りに0に減衰していく。
}}

よって、オイラー法は上記の硬い方程式に対し数値的不安定である。一方、台形公式は数値的安定である。

他の例として、もっとも有名な硬い方程式の一つは、Robertsonの[[化学反応]]を支配する方程式系である。

: <math> \dot x = -0.04 x + 10^4 y \cdot z,</math>
: <math> \dot y = 0.04 x - 10^4 y \cdot z - 3\cdot 10^7 y^2, </math>
: <math> \dot z = 3\cdot 10^7 y^2.</math>

{{math|[0, 40]}} のような短い区間では、上記の方程式系を数値的に積分することに問題はない。しかし区間が極めて大きい場合（例えば {{math|10{{sup|11}}}}）、多数のコードは方程式系を正しく積分することができなくなる。

== 常微分方程式 == 
{{See also|1=ルンゲ＝クッタ法#安定性|2=線型多段法|3=常微分方程式の数値解法}}
上述の例の示すように、硬い[[常微分方程式]]の近似解を計算するとき、数値的に安定な方法を使うべきである。常微分方程式における数値的安定性に複数の定義が存在する。特に、線型方程式に対する安定性と非線型方程式に対する安定性を分けて考える必要がある。

=== 硬さの比例 ===
線形常微分方程式系の硬さは簡単に測ることができる。一般的な線型方程式系
: <math> y' = Ay, \quad y(t_0) = y_0</math> 
を考える。上記方程式に対する '''硬さの比例'''<!--比？--> (stiffness ratio) は、[[行列]] {{mvar|A}} の最大[[固有値]]（の[[絶対値]]）を最小固有値（の絶対値）で割った商である。つまり、{{mvar|A}} の固有値を {{math|''λ{{sub|1}}'' ≥ ''λ{{sub|2}}'' ≥ ... ≥ ''λ{{sub|n}}''}} としたとき、方程式に対する硬さの比例を {{math| {{norm|''λ{{sub|1}}''}} / {{norm|''λ{{sub|n}}''}}}} と定義する。

典型的な硬さの比例は、{{math|10{{sup|17}}}} あたりである。極端な場合に、その数は {{math|10{{sup|31}}}} にもなる{{sfn|Iserles|2008|p=56}}。

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

=== 線型安定性 ===
[[線型微分方程式|線型常微分方程式]]に対する安定性は '''線型安定性''' (linear stability)、あるいは '''絶対安定性''' (absolute stability) という。線型テスト方程式
: <math> y' = \lambda y, \quad t \ge 0, \quad y(0)=1, \quad \lambda \in \mathbb{C} </math>
を考える。この方程式は簡単に解くことができ、厳密解は {{math|''y''(''t'') {{=}} ''e{{sup|λt}}''}} である。{{math|Re ''λ'' < 0}} が成立するとき、{{mvar|y}} の {{math|''t'' → ∞}} の[[極限]]も {{math|0}} である。理想的に、近似解にもそのような振舞いを期待できる。しかし刻み幅 {{mvar|h}} が一定のとき、すべての方法に対する近似解がそのような振舞いを持つわけではない。それを区別するのが線型安定性である。

一つの方法による時刻 {{mvar|t{{sub|n}}}} での近似解を {{mvar|y{{sub|n}}}} とする。[[複素数平面]]上の[[集合]]
: <math> D = \{h\lambda \in \mathbb{C} \mid \lim_{n\to \infty} y_n = 0\} </math>
は方法に対する '''線型安定性領域''' (linear stability domain)、あるいは '''絶対安定性領域''' (region of absolute stability) という{{sfn|Iserles|2008|pp=56–57}}{{sfn|Burden|Faires|1993|p=344}}。この集合はすなわち、与えられた方法による近似解が期待通りの振舞いを持つすべての {{math|''hλ''}} （からなる集合）である。特に、[[ルンゲ＝クッタ法]]に対する線型安定性領域は以下の形で与えられる。
: <math> D = \{z \in \mathbb{C} \mid |r(z)| < 1\} \quad (z = h\lambda) .</math>
ここで、{{math|''r''(''z'')}} は等式 {{math|''y{{sub|n}}'' {{=}} (''r''(''z''))''{{sup|n}}''}} を成立させる関数であり、時々方法に対する '''安定性関数''' という。例えば、オイラー法に対応する関数は {{math|''r''(''z'') {{=}} 1 + ''z''}} である。

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

<!-- また、ルンゲ＝クッタ法にはA-安定性より強い '''L-安定性'''（L-stability）という概念が存在する。線型多段法にもA-安定性とは異なる '''零点安定性'''（zero-stability）が存在する。これらの安定性の定義や意味については、それぞれの記事を参照。-->
=== 非線型安定性 ===

上述の安定性理論に考察されたのは線型方程式のみである。その理論は時折り非線型方程式にも適用できるが、決して正しいわけではない。非線形方程式の研究を完全に一般化するのが困難であるように、すべての方程式に対する安定性を考察するのもほぼ不可能である。現在非線形方程式に対する安定性はほとんど単調性条件 <math>\langle f(t,y) - f(t,z), y - z\rangle \le 0</math> を満足する方程式
: <math> y' = f(t,y) </math>
のみを考える。ただし、<math>\langle \cdot, \cdot \rangle</math> は標準内積である。この発想は、ダールキストによるものである{{sfn|Dahlquist|1976}}。また、[[ルンゲ＝クッタ法]]と[[線型多段法]]に対する安定性の定義は異なる。なぜならば、線型多段法は時刻毎に多数の成分からベクトルを記憶する必要があり、[[偏差]]を測るには[[標準内積]]と異なる[[内積]]を定義しなければならないからである（比べて、ルンゲ＝クッタ法は時刻毎に単一の実数を記憶し、よって標準内積が使われる）<ref name="B2008361">{{harvnb|Butcher|2008|p=361}}.</ref>。

上記の方程式に対して[[ルンゲ＝クッタ法]]の安定性は '''B-安定性''' (B-stability) という。方程式にルンゲ＝クッタ法を適用するときに、異なる初期値 {{math|''y''{{sub|0}}}} と {{math|{{hat|''y''}}{{sub|0}}}} に対し不等式
: <math> \lVert y_1 - \hat{y}_1 \rVert \le \lVert y_0 - \hat{y}_0 \rVert </math>
が成立すれば、その方法をB-安定と呼ぶ<ref name="HW96181">{{harvnb|Hairer|Wanner|1996|p=181}}.</ref>。ここで、{{math|''y''{{sub|1}}}} と {{math|{{hat|''y''}}{{sub|1}}}} は時刻 {{math|''t''{{sub|1}}}} でのそれぞれの近似解である。B-安定な方法は必ずA-安定である<ref name="HW96181" />。

さらに、[[ルンゲ＝クッタ法]]の係数が {{math|''b{{sub|i}}'' ≥ 0}} かつ[[行列]] 
: <math> M = (b_ia_{ij} - b_ja_{ji} - b_ib_j)_{ij} </math>
が[[半正定値]]であるという条件を満足するとき、その方法を '''代数的安定''' (algebraically stable) という。代数的安定な方法は必ずB-安定である{{sfn|Hairer|Wanner|1996|p=182}}。

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

{{mvar|k}} 段線型多段法の一般形式は次の公式で与えられる。
: <math> \sum_{i=0}^k \alpha_i y_{n-i} = h\sum_{i=0}^k \beta_i f(t_{n-i},y_{n-i}). </math>
ここで、{{mvar|α{{sub|i}}}} と {{mvar|β{{sub|i}}}} は定数であり、ベクトル {{math|''α'' {{=}} (''α{{sub|0}}'', ..., ''α{{sub|k}}'')}} と {{math|''β'' {{=}} (''β{{sub|0}}'', ..., ''β{{sub|k}}'')}} は方法の '''生成ペア''' (generating pair) という<ref group="注">[[線型多段法]]は生成ペアで一意に定まる。また、これらの係数で定義される多項式 <math> \rho(\omega) = \sum_{i=0}^k \alpha_i\omega^{k-i}</math> と <math> \sigma(\omega) = \sum_{i=0}^k \beta_i\omega^{k-i} </math> も[[線型多段法]]の分析に重要である。</ref>。 この方法に対応する '''one-leg法''' (one-leg method) は次の公式で与えられる{{sfn|Butcher|2008|p=360}}。
: <math> \sum_{i=0}^k \alpha_i y_{n-i} = h\Biggl(\sum_{i=0}^k \beta_i\Biggr)\,
f\Biggl(x_{n} - \theta h,\biggl(\sum_{i=0}^k \beta_i\biggr)^{-1}\sum_{i=0}^k \beta_i y_{n-i}\Biggr) </math>
ただし、
: <math> \theta = \frac{\sum_{i=0}^k i\beta_i}{\sum_{i=0}^k \beta_i} </math>
である。[[線型多段法]]は対応するone-leg法と同じ線型安定性を持つため、同じ非線形安定性を持つことも期待できる。しかし上記の方程式に対する安定性を分析するには、線型多段法より対応するone-leg法を用いる方が遥かに簡単である<ref name="B2008361"/>。よって以下の定義はone-leg法に対するものである。

与えられた {{mvar|k}} 次正定値[[対称行列]] {{mvar|G}} に対応する <math> \mathbb{R}^{kn} </math> 上の内積を以下のように定義できる：
: <math> \langle U, V \rangle_G = \sum_{i=1}^k \sum_{j=1}^k g_{ij}\langle U_i, V_j \rangle, \quad U, V \in \mathbb{R}^{kn} </math>
ここで、<math> U = (U_1,\ldots, U_k)^\mathrm{T}, \; U_i \in \mathbb{R}^n </math>である。

one-leg法に対し、行列
: <math> M = \alpha \beta^\mathrm{T} + \beta \alpha^\mathrm{T} - \begin{bmatrix} G & 0 \\ 0 & 0 \end{bmatrix}
+ \begin{bmatrix} 0 & 0 \\ 0 & G \end{bmatrix} </math>
を半正定値にする {{mvar|G}} が存在するとき、その方法をG-安定という{{sfn|Butcher|2008|p=364}}。この定義は初見で[[ルンゲ＝クッタ法]]の安定性とは全く違うものと思われるかもしれないが、本質的には同じものである。なぜならば、上記の定義を以て、不等式
: <math>\lVert y_1 - \hat{y}_1 \rVert_G \le \lVert y_0 - \hat{y}_0 \rVert_G </math>
を証明できるからである<ref>{{harvnb|Butcher|2008|pp=362-364}}</ref>。

また、G-安定性もB-安定性のようにA-安定性より強い条件に見えるかもしれないが、実際にA-安定性とは[[同値]]である{{sfn|Dahlquist|1978}}。すなわち、[[線型多段法]]がA-安定であることは対応するone-leg法がG-安定であることの必要十分条件となる。

== 脚注 ==
===注===
{{reflist|group="注"}}
===出典===
{{reflist|30em}}

== 参考文献 ==
* {{citation | first1 = Richard L. | last1 = Burden | first2 = J. Douglas | last2 = Faires | year = 1993 | isbn = 0-534-93219-3 | title = Numerical Analysis | edition = 5th | publisher = [[Prindle, Weber and Schmidt]] | location = Boston }}.
* {{Citation | last1=Butcher | first1=John C. | author1-link=John C. Butcher | title=Numerical Methods for Ordinary Differential Equations | publisher=[[John Wiley & Sons]] | location=New York | isbn=978-0-470-72335-7 | year=2008}}.
* {{citation | first1 = Germund | last1 = Dahlquist | year = 1976 | title = Error analysis for a class of methods for stiff non-linear initial value problems | series = Lecture Notes in Mathematics | volume = 506 | pages = 60-72 | doi = 10.1007/BFb0080115 }}.
* {{citation | first1 = Germund | last1 = Dahlquist | year = 1978 | title = G-stability is equivalent to A-stability |journal = BIT | volume = 18 | issue = 4 | pages = 384-401 | }}.
<!--* {{citation | first1 = Germund | last1 = Dahlquist | author1-link = Germund Dahlquist | year = 1963 | title = A special stability problem for linear multistep methods | journal = BIT | volume = 3 | issue = 1 | pages = 27–43 | doi = 10.1007/BF01963532 }}.
* {{citation | first1 = David | last1 = Eberly| year = 2008 | title = Stability analysis for systems of differential equations  }}.[http://www.geometrictools.com/Documentation/StabilityAnalysis.pdf]
* {{citation | first1 = B. L. | last1 = Ehle | year = 1969 | title = On Padé approximations to the exponential function and A-stable methods for the numerical solution of initial value problems | publisher = Report 2010, [[University of Waterloo]] }}.
* {{citation | first1 = C. W. | last1 = Gear | year = 1971 | title = Numerical Initial-Value Problems in Ordinary Differential Equations | publisher = [[Prentice Hall]] | location = Englewood Cliffs }}.
* {{citation | first1 = C. W. | last1 = Gear | year = 1981 | title = Numerical solution of ordinary differential equations: Is there anything left to do? | journal = [[Society for Industrial and Applied Mathematics|SIAM Review]] | volume = 23 | issue = 1 | pages = 10–24 }}. -->
* {{citation | first1 = Ernst | last1 = Hairer | first2 = Gerhard | last2 = Wanner | year = 1996 | title = Solving ordinary differential equations II: Stiff and differential-algebraic problems | edition = second | publisher = [[Springer-Verlag]] | location = Berlin | isbn = 978-3-540-60452-5 }}.
<!-- * {{citation | first1 = J. O. | last1 = Hirshfelder | year = 1963 | title = Applied Mathematics as used in Theoretical Chemistry | journal = [[American Mathematical Society|American Mathematical Society Symposium]] | pages = 367–376 }}. -->
* {{Citation|title=A First Course in the Numerical Analysis of Differential Equations (Second Edition)|year=2008|last1=Iserles|first1=Arieh|publisher=Cambridge University Press|isbn=978-0-521-73490-5}}.
* {{citation | first1 = Arieh | last1 = Iserles | first2 = Syvert | last2 = Nørsett | year = 1991 | title = Order Stars | publisher = [[Chapman & Hall]] | isbn = 978-0-412-35260-7 }}.
<!-- * {{citation | first1 = Erwin | last1 = Kreyszig | year = 1972 | isbn = 0-471-50728-8 | title = Advanced Engineering Mathematics | edition = 3rd | publisher = [[John Wiley & Sons|Wiley]] | location = New York }}.
* {{citation | first1 = J. D. | last1 = Lambert | year = 1977 | title = The initial value problem for ordinary differential equations | journal = The State of the Art in Numerical Analysis | editor = D. Jacobs | publisher = [[Academic Press]] | location = New York | pages = 451–501 }}.
* {{citation | first1 = J. D. | last1 = Lambert | year = 1992 | isbn = 978-0-471-92990-1 | title = Numerical Methods for Ordinary Differential Systems | publisher = [[John Wiley & Sons|Wiley]] | location = New York }}.
*{{citation | first1 = John  | last1 = Mathews |first2= Kurtis|last2= Fink| year = 1992 | title = Numerical methods using MATLAB}}.
*{{Cite book | last1=Press | first1=WH | last2=Teukolsky | first2=SA | last3=Vetterling | first3=WT | last4=Flannery | first4=BP | year=2007 | title=Numerical Recipes: The Art of Scientific Computing | edition=3rd | publisher=[[Cambridge University Press]] |  publication-place=New York | isbn=978-0-521-88068-8 | chapter=Section 17.5. Stiff Sets of Equations | chapter-url=http://apps.nrbook.com/empanel/index.html#pg=931}}
* {{citation | first1 = L. F. | last1 = Shampine | first2 = C. W. | last2 = Gear | year = 1979 | title = A user's view of solving stiff ordinary differential equations | journal = [[Society for Industrial and Applied Mathematics|SIAM Review]] | volume = 21 | issue = 1 | pages = 1–17 }}.
* {{citation | first1 = Gerhard | last1 = Wanner | first2 = Ernst | last2 = Hairer | first3 = Syvert | last3 = Nørsett | year = 1978 | title = Order stars and stability theory | journal = BIT | volume = 18 | issue = 4 | pages = 475–489 | doi = 10.1007/BF01932026 }}.
*Stability of Runge-Kutta Methods [http://homepages.cwi.nl/~jason/Classes/numwisk/ch10.pdf ] -->

== 外部リンク ==
* [https://www.cs.cmu.edu/~baraff/pbm/energons.pdf An Introduction to Physically Based Modeling: Energy Functions and Stiffness]
* [http://www.scholarpedia.org/article/Stiff_systems Stiff systems] [[Lawrence F. Shampine]] and [[Skip Thompson]] [[ Scholarpedia]], 2(3):2855.	[[doi:10.4249/scholarpedia.2855]]



{{DEFAULTSORT:かたいほうていしき}}
[[Category:解析学]]
[[Category:数値解析]]
[[Category:数値微分方程式]]
[[Category:数学に関する記事]]