「常微分方程式の数値解法」の版間の差分
削除された内容 追加された内容
MathXplore (会話 | 投稿記録) m →洋書 |
重複する誘導を除去、typo修正 |
||
| (同じ利用者による、間の1版が非表示) | |||
1行目:
{{Differential equations}}
'''常微分方程式の数値解法''' (じょうびぶんほうていしきのすうちかいほう、{{lang-en-short|Numerical methods for ODEs}}) は、[[数値解析]]において[[常微分方程式]]を数値的に解く技術の総称である<ref name="Yamamoto1">{{Cite book |和書 |author=山本哲朗 |title=数値解析入門 |edition=増訂版 |date=2003-06 |publisher=[[サイエンス社]] |series=サイエンスライブラリ 現代数学への入門 14 |ISBN=4-7819-1038-6}}</ref><ref name="mori">[[森正武]]. 数値解析 第2版. [[共立出版]].</ref>。
{{seealso|常微分方程式|数値解析}}▼
==数値解法の必要性==
これまで様々な自然現象 (物理現象など) を記述するために多くの[[常微分方程式]]が作られ、多くの数学者たちがその解法を探求してきたが、[[フックス型微分方程式]]<ref name="toki">時弘哲治、工学における[[特殊関数]]、[[共立出版]]。</ref><ref name="sakai">坂井秀隆. (2015). [[常微分方程式]]. [[東京大学出版会]].</ref>などを除いて、手計算だけで厳密に解ける常微分方程式は多くない。そのため多くの研究者たちが常微分方程式を数値的に解く技術について研究をしてきた<ref name="Yamamoto1"/><ref name="mori"/>。最も標準的な手法は[[ルンゲ・クッタ法]]であり<ref name="Yamamoto1"/><ref name="mori"/><ref>遠藤理平 (2018) ルンゲ・クッタで行こう!―物理シミュレーションを基礎から学ぶ</ref><ref group="A">Butcher, J. C. (1996). A history of Runge-Kutta methods. Applied numerical mathematics, 20(3), 247-260.</ref>、[[MATLAB]]にはode45として搭載されている。しかしこれは万能なソルバーとは言えない。例えば[[パンルヴェ方程式]]<ref>岡本和夫. (2009). [[パンルヴェ方程式]]. [[岩波書店]].</ref><ref>野海正俊. (2000). パンルヴェ方程式-対称性からの入門. すうがくの風景 4. [[朝倉書店]].</ref><ref>岡本和夫. (1985). パンルヴェ方程式序説. [[上智大学]]数学講究録, 19.</ref>や[[リッカチ方程式]]<ref>リッカチのひ・み・つ : 解ける微分方程式の理由を探る. 井ノ口順一著. [[日本評論社]], 2010.9.</ref>などは非線形性によって精度の良い計算ができず、数値実験結果だけを見ていると間違った結論 ([[精度保証付き数値計算#幻影解|幻影解]]) にたどり着く危険がある。{{要出典|範囲=そのため|date=2021年2月}}、▼
{{seealso|硬い方程式|ルンゲ＝クッタ法のリスト}}▼
▲これまで様々な自然現象 (物理現象など) を記述するために多くの[[常微分方程式]]が作られ、多くの数学者たちがその解法を探求してきたが、[[フックス型微分方程式]]<ref name="toki">時弘哲治、工学における[[特殊関数]]、[[共立出版]]。</ref><ref name="sakai">坂井秀隆. (2015). [[常微分方程式]]. [[東京大学出版会]].</ref>などを除いて、手計算だけで厳密に解ける常微分方程式は多くない。そのため多くの研究者たちが常微分方程式を数値的に解く技術について研究をしてきた<ref name="Yamamoto1"/><ref name="mori"/>。最も標準的な手法は[[ルンゲ・クッタ法]]であり<ref name="Yamamoto1"/><ref name="mori"/><ref>遠藤理平 (2018) ルンゲ・クッタで行こう!―物理シミュレーションを基礎から学ぶ</ref><ref group="A">Butcher, J. C. (1996). A history of Runge-Kutta methods. Applied numerical mathematics, 20(3), 247-260.</ref>、[[MATLAB]]にはode45として搭載されている。しかしこれは万能なソルバーとは言えない。例えば[[パンルヴェ方程式]]<ref>岡本和夫. (2009). [[パンルヴェ方程式]]. [[岩波書店]].</ref><ref>野海正俊. (2000). パンルヴェ方程式-対称性からの入門. すうがくの風景 4. [[朝倉書店]].</ref><ref>岡本和夫. (1985). パンルヴェ方程式序説. [[上智大学]]数学講究録, 19.</ref>や[[リッカチ方程式]]<ref>リッカチのひ・み・つ : 解ける微分方程式の理由を探る. 井ノ口順一著. [[日本評論社]], 2010.9.</ref>などは非線形性によって精度の良い計算ができず、数値実験結果だけを見ていると間違った結論 ([[精度保証付き数値計算#幻影解|幻影解]]) にたどり着く危険がある。そのため、
{{div col|rules=yes}}
*[[線型多段法]]<ref name="Yamamoto1"/><ref name="mori"/>
20 ⟶ 19行目:
{{div col end}}
などの新しい解法に関する研究が進められている。
== 初期値問題==
<math>N</math> 個の変数 <math>y = ( y_1, \cdots, y_N ) \in \mathbb{R}^N</math> に関する1階[[常微分方程式]]、すなわち <math>[ t_0, b ] \times \mathbb{R}^N</math> （またはその開集合）で定義されたベクトル値連続関数 <math>f ( t, y ) = ( f_1 ( t, y ), \cdots, f_N ( t, y ) )</math> により定まる次の方程式
{{Indent|<math>\frac{ d y }{ d t } = f ( t, y )</math>}}
について考える。その[[初期値問題]] (initial-value problem) とは、初期条件 <math>y ( t_0 ) = y_0</math> を満たす関数 <math>y ( t )</math> を求めることである<ref>Stoer & Bulirsch, p. 465.</ref>。関数 <math>f</math> が第2引数について[[リプシッツ連続]]であるとき、すなわち定数 <math>L</math> が存在し
{{Indent|<math>\left\| f ( t, y_1 ) - f ( t, y_2 ) \right\| \leq L \| y_1 - y_2 \|</math>}}
（<math>\| \cdot \|</math> は <math>\mathbb{R}^N</math> の[[ノルム]]）を満たすとき、その初期値問題には解が一意に存在する<ref>Stoer & Bulirsch, p. 467.</ref>（[[ピカール・リンデレフの定理]]）。本節では常微分方程式の初期値問題を数値的に解くことについて議論する。
=== 初期値問題の例 ===
例として、[[電気回路]]の研究から導かれた[[ファン・デル・ポール振動子]]について考える。その運動方程式は
{{Indent|<math>\frac{ d^2 x }{ dt^2 }- \mu (1 - x^2) \frac{ dx }{ dt } + x= 0</math>}}
であり、これは2階の常微分方程式であるものの、<math>v = d x / d t</math> とおくと
{{Indent|<math>\frac{ d x }{ d t } = v , \ \ \frac{ d v }{ d t } = - x + \mu ( 1 - x^2 ) v</math>}}
という上述の形に帰着できる<ref>Iserles, pp. 106-107.</ref>。
=== 一段法 ===
区間 <math>[ t_0, t ]</math> の厳密解を <math>y ( t )</math> とする。一段法<ref name="加古">{{Cite web |author=加古富志雄 |url=https://www.ics.nara-wu.ac.jp/~kako/teaching/na/chap8.pdf |title=数値解析 |accessdate=2021-02-02}}</ref>として知られるクラスの数値解法では、離散化した時刻 <math>t_i = t_0 + i h</math> (<math>h = ( t - t_0 )/n</math>) での厳密解 <math>y ( t_i )</math> の近似値 <math>\eta_i</math> を
{{Indent|<math>\eta_{i + 1} = \eta_i + h \Phi ( t_i, \eta_i; h )</math>}}
という漸化式によって定める<ref>Stoer & Bulirsch, p. 473.</ref>。関数 <math>\Phi</math> の選択が数値積分スキームを選択することに対応する。離散化した時刻の差分 <math>h = t_{i+1} - t_i</math> を刻み幅あるいはステップサイズと呼ぶ<ref>齊藤, pp. 97-98.</ref>。なお、ここでは時間ステップ <math>h = t_{i+1} - t_i</math> は一定としたが、これを動的に決定する{{仮リンク|適応刻み|en|Adaptive step size|Adaptive step size}}という手法もある<ref>Press et al., p. 709.</ref><ref>{{Cite web |author=Jun Makino |url=http://jun-makino.sakura.ne.jp/kougi/system_suuri4_1998/all/node24.html |title=5.6 刻み幅調節と埋め込み型公式 |accessdate=2021-02-02}}</ref>。
厳密解 <math>y ( t )</math> から[[差分商]]
{{Indent|<math>\Delta ( t, \eta; h ) = \begin{cases} \frac{ y ( t + h ) - \eta }{ h } & h \neq 0 \\ f ( t, \eta ) & h = 0 \end{cases}</math>}}
を導入するとき、数値積分スキーム <math>\Phi</math> の局所離散化誤差<ref>齊藤, pp. 99, 102.</ref> (local discretization error) は
{{Indent|<math>\Delta ( t, \eta; h ) - \Phi ( t, \eta; h )</math>}}
により定義される<ref>Stoer & Bulirsch, pp. 473-474.</ref><ref>Hackbusch, p. 71.</ref>。整合性のために <math>h \to 0</math> の極限で局所離散化誤差（厳密にはその上限）は 0 に収束することが要求される<ref>Stoer & Bulirsch, p. 474.</ref><ref>Hackbusch, p. 72.</ref>。さらに、この極限で局所離散化誤差が
<math>\Delta ( t, \eta; h ) - \Phi ( t, \eta; h ) = O ( h^p )</math>
を満足するとき、この積分スキーム <math>\Phi</math> は <math>p</math> 次精度であるという<ref>齊藤, p. 103.</ref><ref>Stoer & Bulirsch, p. 474.</ref>。
一方、<math>t</math> を固定して <math>n \to \infty</math> とするとき、大域離散化誤差<ref>齊藤, p. 102.</ref> (global discretization error)
{{Indent|<math>e ( t; h_n ) = \eta ( t; h_n ) - y ( t ) , \ \ h_n = \frac{ t - t_0 }{ n }</math>}}
が 0 に収束するならば、その積分スキームは収束するという<ref>Stoer & Bulirsch, p. 477.</ref><ref>Hackbusch, p. 72.</ref>。<math>p</math> 次精度 (<math>p > 0</math>) の1段法を十分に滑らかな関数 <math>f</math> に適用するとき、そのスキームは収束し大域離散化誤差は
{{Indent|<math>e ( t; h_n ) = O ( h_n^p )</math>}}
のように振る舞うことが保証されている、すなわち大域離散化誤差のオーダーは局所離散化誤差のオーダーに等しい<ref>Stoer & Bulirsch, pp. 477-479.</ref><ref>Hackbusch, pp. 73-74.</ref>。この結果はすべての整合的な一段法が <math>h \to 0</math> で漸近的に[[数値的安定性|安定]]であることを意味するものの、ただし現実的に可能な <math>h</math> で1段法が安定であることは必ずしも保証されない（[[#安定性]]節を参照）<ref>Hackbusch, p. 74.</ref>。
[[File:Error of ODE Solvers (Japanese label).svg|thumb|320px|オイラー法、ホイン法、古典的ルンゲ＝クッタ法 (RK4) の相対誤差の比較。初期値 <math>y ( 0 ) = 0</math> のもとでの常微分方程式 <math>\frac{ d y }{ d t } = \cos ( y )</math> の数値解の <math>t = 1</math> での値をステップサイズの関数として対数プロットした。各手法がそれぞれ1次精度、2次精度、4次精度であることに対応して、傾き 1, 2, 4 で誤差が減少している<ref>齊藤, p. 107.</ref>。]]
積分スキーム <math>\Phi</math> としては以下のものが知られている。
* [[オイラー法]]<ref>齊藤, p. 97-98, 103.</ref>: <math>\Phi( t, \eta; h ) = f ( t, \eta )</math> とするもの。これは1次精度のスキームである。
* {{仮リンク|ホイン法|en|Heun's method}}<ref>齊藤, p. 106.</ref>: <math>\Phi ( t, \eta; h ) = [ f ( t, \eta ) + f ( t + h, \eta + h f ( t, \eta ) ) ] / 2</math> とするもの。これは2次精度のスキームであり、関数 <math>f</math> の評価を2回必要とする。
* [[ルンゲ＝クッタ法|古典的ルンゲ＝クッタ法]]: これは4次精度のスキームであり、関数 <math>f</math> の評価を4回必要とする<ref>Stoer & Bulirsch, pp. 475-476.</ref>。
このうち古典的ルンゲ＝クッタ法は適用可能範囲の広さ<ref>Press et al., p. 709.</ref>やプログラミングの容易さのために広く用いられている<ref>{{Cite web |url=http://www.ed.u-tokai.ac.jp/takakura/jp/takakura/book_SuuchiKeisan/RC10.pdf |title=Ｃ言語による数値計算プログラミング演習 10. 常微分方程式の初期値問題 |accessdate=2021-02-02}}</ref>（ただし Press らは著書において計算速度の観点から古典的ルンゲ＝クッタ法に否定的に言及している<ref>Press et al., pp. 708-709, 712.</ref>）。
=== 多段法 ===
一段法は <math>\eta_{i+1}</math> の値を <math>\eta_i</math> だけから定めるものであったが、より多くのステップでの値 <math>\eta_{i-1}</math>, ..., <math>\eta_{i-r+1}</math> を使う積分スキームは多段法<ref name="加古"/> (multistep method) と呼ばれる<ref>Stoer & Bulirsch, p. 492.</ref>。この場合、最初の <math>\eta_0</math>, ..., <math>\eta_{r-1}</math> はこのスキームでは定めることができず、1段法などの他の方法を用いる必要がある<ref>Stoer & Bulirsch, p. 492.</ref>。多段法としては[[アダムス・バッシュフォース法]]や、それを応用する[[予測子修正子法]]などがある<ref>Press et al., pp. 747-750.</ref>（これはどちらも <math>\eta_{i+1}</math> が <math>\eta_{i-r+1}</math>, ..., <math>\eta_{i}</math> の線型関数として定まるため、[[線型多段法]]と呼ばれる<ref>Stoer & Bulirsch, p. 508.</ref>）。
{{詳細記事|線型多段法}}
=== 安定性 ===
数値積分スキームの安定性はしばしばそれを初期値問題
{{Indent|<math>\frac{ d y }{ d t } = \lambda y , \ \ y ( 0 ) = 1</math>}}
（<math>\lambda</math> は複素数の定数）に適用することで定量化される<ref>Iserles, pp. 56-57.</ref>。問題のスキームを一定の刻み幅 <math>h</math> で適用することで得られる数列 <math>\{ y_n \}_{n = 0, 1, \cdots}</math> について、それが極限 <math>n \to \infty</math> で 0 に収束するような <math>h \lambda</math> がなす[[複素平面]]上の領域を絶対安定領域<ref>{{Cite web |author=Jun Makino |date=1998 |url=http://jun-makino.sakura.ne.jp/kougi/system_suuri4_1999/note7/node2.html |title=2 数値解の安定性 |accessdate=2021-02-01}}</ref>と呼ぶ<ref>Iserles, pp. 56-57.</ref>。もとの初期値問題の厳密解は <math>y ( t ) = e^{\lambda t}</math> であり、これは <math>\mathrm{Re} \{ \lambda \} < 0</math> のとき <math>t \to \infty</math> で 0 に収束する<ref>Iserles, p. 56.</ref>。そこである積分スキームの絶対安定領域が左半平面を含むとき、そのスキームはA-安定 (A-stable) であるという<ref>Iserles, p. 59.</ref><ref name="牧野-A安定性">{{Cite web |author=Jun Makino |date=1998 |url=http://jun-makino.sakura.ne.jp/kougi/system_suuri4_1998/all/node62.html |title=11.2 A-安定性 |accessdate=2021-02-01}}</ref>。
関数 <math>f ( t, y )</math> が大きな値を取り積分スキームによっては極端に時間刻み幅 <math>h</math> を小さくする必要がある場合、その常微分方程式は[[硬い方程式]]と呼ばれる<ref>Iserles, pp. 53-56.</ref>。この場合には十分に安定なスキームを用いる必要が生じる<ref name="牧野-A安定性"/>。これはしばしば <math>y_{i+1}</math> を[[陰関数]]的に <math>y_i</math> 等から定める[[陰解法]]によって実現される<ref>Press et al., pp. 735-736.</ref>。
{{詳細記事|硬い方程式}}
=== 幾何学的解法 ===
解析対象となる微分方程式が何らかの特別な性質を持つとき、汎用の数値積分法ではなく、その性質を尊重するように構成された数値解法を用いることがあり、そのような手法は構造保存型解法または幾何学的数値解法と呼ばれる<ref>齊藤, p. 120.</ref><ref name="宮武">{{Cite web |author=宮武勇登 |url=http://coop-math.ism.ac.jp/files/190/slide-Upscaling2016-miyatake.pdf |title=保存則に即した数値計算手法 |accessdate=2021-02-02}}</ref>。例えば[[古典力学]]（[[ハミルトン力学]]）において時間発展は[[シンプレクティック同相写像|シンプレクティック写像]]であり、[[エネルギー]]（[[ハミルトニアン]]）等の保存量が存在する<ref name="宮武"/>。この場合には[[シンプレクティック数値積分法]]というハミルトン系にのみ適用可能な数値解法が存在し、良好な性質を持つことが知られている<ref name="Yoshida1994">{{Cite web |url=http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/0889-06.pdf |title=可変時間ステップによるシンプレクティック数値解法 |author=吉田春夫 |accessdate=2021-02-02}}</ref>。
{{詳細記事|シンプレクティック数値積分法}}
== 境界値問題 ==
{{Indent|<math>\frac{ d y }{ d x } = f ( x, y )</math>}}
の解 <math>y ( x )</math> で, <math>a \neq b</math> に対する境界条件
{{Indent|<math>r ( y ( a ), y ( b ) ) = c</math>}}
を満足するものを求める問題を[[境界値問題]] (boundary-value problem) と呼ぶ<ref>Stoer & Bulirsch, pp. 539-540.</ref>。初期値問題と異なり、境界値問題では複数の解が存在すること、あるいは解が存在しないことがあり得る。また、[[スツルム＝リウヴィル型微分方程式]]のように常微分方程式にパラメータ <math>\lambda</math> が含まれ、境界値問題に解が存在するようにパラメータ <math>\lambda</math> を同時に定める問題もある<ref>Stoer & Bulirsch, p. 541.</ref>。
これらの問題を数値的に解く最も単純な方法が[[狙い撃ち法]] (shooting method) である<ref>Stoer & Bulirsch, pp. 542-548.</ref>。一方の端点（例えば <math>x = a</math>）において初期条件 <math>y ( a )</math> を適当に定めて微分方程式をもう一方の端点 <math>x = b</math> まで解き、境界条件を満たすように未定の初期条件（および固有値）を適切に選ぶ<ref>Press et al., pp. 753-755, 757-759.</ref><ref>Hairer et al. (1993), p. 105.</ref>。これは[[ニュートン法]]などの関数の根を求めるアルゴリズム（これはしばしば逐次反復を伴う）を常微分方程式の数値解法と組み合わせることを意味する<ref>Hairer et al. (1993), p. 106.</ref><ref>Press et al., pp. 754-755.</ref>。ただし初期条件によっては区間 <math>[ a, b ]</math> 全体で定義された解が存在しないことがあり、そのために改良された手法がある<ref>Hairer et al. (1993), pp. 106-107.</ref>。あるいは[[有限要素法]] (finite element method) などの手法も用いられる<ref>Iserles, pp. 171-176.</ref>。
==解の存在検証==
51 ⟶ 121行目:
* 三井斌友 (2003) 常微分方程式の数値解法, [[岩波書店]].
* 三井斌友 et. al. (2004). [[微分方程式]]による[[計算科学]]入門. [[共立出版]].
* {{Cite book |和書 |author=齊藤宣一 |title=数値解析 （共立講座 数学探検 17） |publisher=共立出版 |date=2017 |isbn=9784320992740}}
===洋書===
{{Div col|rules=yes}}
62 ⟶ 134行目:
* Shampine, L. F. (2018). Numerical solution of ordinary differential equations. Routledge.
* Dormand, John R. (1996), Numerical Methods for Differential Equations: A Computational Approach, Boca Raton: [[:en:CRC Press]].
* {{Cite book |last1=Press |first1=William H. |first2=Saul A. |last2=Teukolsky |first3=William T. |last3=Vetterling |first4=Brian P. |last4=Flannery |date=1992 |title=Numerical Recipes in C: The Art of Scientific Computing |edition=2nd |publisher=Cambridge University Press |doi=10.2277/0521431085 |isbn=978-0-521-43108-8}}
* {{Cite book |last1=Stoer |first1=Josef |last2=Bulirsch |first2=R. |title=Introduction to Numerical Analysis |publisher=Springer |date=2002 |isbn=978-0-387-21738-3 |doi=10.1007/978-0-387-21738-3}}
* {{Cite book |first=Wolfgang |last=Hackbusch |title=The Concept of Stability in Numerical Mathematics |date=2014 |publisher=Springer |isbn=978-3-642-39386-0 |doi=10.1007/978-3-642-39386-0}}
{{Div col end}}
====微分代数方程式の数値解法====
| |||