コンテンツにスキップ

英文维基 | 中文维基 | 日文维基 | 草榴社区

「常微分方程式の数値解法」の版間の差分

出典: フリー百科事典『ウィキペディア(Wikipedia)』
削除された内容 追加された内容
重複する誘導を除去、typo修正
(同じ利用者による、間の1版が非表示)
1行目: 1行目:
{{Differential equations}}
{{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>。
'''常微分方程式の数値解法''' (じょうびぶんほうていしきのすうちかいほう、{{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}}
{{div col|rules=yes}}
*[[線型多段法]]<ref name="Yamamoto1"/><ref name="mori"/>
*[[線型多段法]]<ref name="Yamamoto1"/><ref name="mori"/>
20行目: 19行目:
{{div col end}}
{{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=C言語による数値計算プログラミング演習 10. 常微分方程式の初期値問題 |accessdate=2021-02-02}}</ref>(ただし Press らは著書において計算速度の観点から古典的ルンゲ=クッタ法に否定的に言及している<ref>Press et al., pp. 708-709, 712.</ref>)。

{{See also |ルンゲ=クッタ法}}

=== 多段法 ===
一段法は <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) 常微分方程式の数値解法, [[岩波書店]].
* 三井斌友 (2003) 常微分方程式の数値解法, [[岩波書店]].
* 三井斌友 et. al. (2004). [[微分方程式]]による[[計算科学]]入門. [[共立出版]].
* 三井斌友 et. al. (2004). [[微分方程式]]による[[計算科学]]入門. [[共立出版]].
* {{Cite book |和書 |author=齊藤宣一 |title=数値解析 (共立講座 数学探検 17) |publisher=共立出版 |date=2017 |isbn=9784320992740}}

===洋書===
===洋書===
{{Div col|rules=yes}}
{{Div col|rules=yes}}
62行目: 134行目:
* Shampine, L. F. (2018). Numerical solution of ordinary differential equations. Routledge.
* 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]].
* 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}}
{{Div col end}}
====微分代数方程式の数値解法====
====微分代数方程式の数値解法====

2021年2月3日 (水) 09:25時点における版

常微分方程式の数値解法 (じょうびぶんほうていしきのすうちかいほう、: Numerical methods for ODEs) は、数値解析において常微分方程式を数値的に解く技術の総称である[1][2]

数値解法の必要性

これまで様々な自然現象 (物理現象など) を記述するために多くの常微分方程式が作られ、多くの数学者たちがその解法を探求してきたが、フックス型微分方程式[3][4]などを除いて、手計算だけで厳密に解ける常微分方程式は多くない。そのため多くの研究者たちが常微分方程式を数値的に解く技術について研究をしてきた[1][2]。最も標準的な手法はルンゲ・クッタ法であり[1][2][5][A 1]MATLABにはode45として搭載されている。しかしこれは万能なソルバーとは言えない。例えばパンルヴェ方程式[6][7][8]リッカチ方程式[9]などは非線形性によって精度の良い計算ができず、数値実験結果だけを見ていると間違った結論 (幻影解) にたどり着く危険がある。そのため[要出典]

などの新しい解法に関する研究が進められている。

初期値問題

個の変数 に関する1階常微分方程式、すなわち (またはその開集合)で定義されたベクトル値連続関数 により定まる次の方程式

について考える。その初期値問題 (initial-value problem) とは、初期条件 を満たす関数 を求めることである[12]。関数 が第2引数についてリプシッツ連続であるとき、すなわち定数 が存在し

ノルム)を満たすとき、その初期値問題には解が一意に存在する[13]ピカール・リンデレフの定理)。本節では常微分方程式の初期値問題を数値的に解くことについて議論する。

初期値問題の例

例として、電気回路の研究から導かれたファン・デル・ポール振動子について考える。その運動方程式は

であり、これは2階の常微分方程式であるものの、 とおくと

という上述の形に帰着できる[14]

一段法

区間 の厳密解を とする。一段法[15]として知られるクラスの数値解法では、離散化した時刻 () での厳密解 の近似値

という漸化式によって定める[16]。関数 の選択が数値積分スキームを選択することに対応する。離散化した時刻の差分 を刻み幅あるいはステップサイズと呼ぶ[17]。なお、ここでは時間ステップ は一定としたが、これを動的に決定する適応刻み英語版という手法もある[18][19]

厳密解 から差分商

を導入するとき、数値積分スキーム の局所離散化誤差[20] (local discretization error) は

により定義される[21][22]。整合性のために の極限で局所離散化誤差(厳密にはその上限)は 0 に収束することが要求される[23][24]。さらに、この極限で局所離散化誤差が を満足するとき、この積分スキーム 次精度であるという[25][26]

一方、 を固定して とするとき、大域離散化誤差[27] (global discretization error)

が 0 に収束するならば、その積分スキームは収束するという[28][29] 次精度 () の1段法を十分に滑らかな関数 に適用するとき、そのスキームは収束し大域離散化誤差は

のように振る舞うことが保証されている、すなわち大域離散化誤差のオーダーは局所離散化誤差のオーダーに等しい[30][31]。この結果はすべての整合的な一段法が で漸近的に安定であることを意味するものの、ただし現実的に可能な で1段法が安定であることは必ずしも保証されない(#安定性節を参照)[32]

オイラー法、ホイン法、古典的ルンゲ=クッタ法 (RK4) の相対誤差の比較。初期値 のもとでの常微分方程式 の数値解の での値をステップサイズの関数として対数プロットした。各手法がそれぞれ1次精度、2次精度、4次精度であることに対応して、傾き 1, 2, 4 で誤差が減少している[33]

積分スキーム としては以下のものが知られている。

  • オイラー法[34]: とするもの。これは1次精度のスキームである。
  • ホイン法英語版[35]: とするもの。これは2次精度のスキームであり、関数 の評価を2回必要とする。
  • 古典的ルンゲ=クッタ法: これは4次精度のスキームであり、関数 の評価を4回必要とする[36]

このうち古典的ルンゲ=クッタ法は適用可能範囲の広さ[37]やプログラミングの容易さのために広く用いられている[38](ただし Press らは著書において計算速度の観点から古典的ルンゲ=クッタ法に否定的に言及している[39])。

多段法

一段法は の値を だけから定めるものであったが、より多くのステップでの値 , ..., を使う積分スキームは多段法[15] (multistep method) と呼ばれる[40]。この場合、最初の , ..., はこのスキームでは定めることができず、1段法などの他の方法を用いる必要がある[41]。多段法としてはアダムス・バッシュフォース法や、それを応用する予測子修正子法などがある[42](これはどちらも , ..., の線型関数として定まるため、線型多段法と呼ばれる[43])。

安定性

数値積分スキームの安定性はしばしばそれを初期値問題

は複素数の定数)に適用することで定量化される[44]。問題のスキームを一定の刻み幅 で適用することで得られる数列 について、それが極限 で 0 に収束するような がなす複素平面上の領域を絶対安定領域[45]と呼ぶ[46]。もとの初期値問題の厳密解は であり、これは のとき で 0 に収束する[47]。そこである積分スキームの絶対安定領域が左半平面を含むとき、そのスキームはA-安定 (A-stable) であるという[48][49]

関数 が大きな値を取り積分スキームによっては極端に時間刻み幅 を小さくする必要がある場合、その常微分方程式は硬い方程式と呼ばれる[50]。この場合には十分に安定なスキームを用いる必要が生じる[49]。これはしばしば 陰関数的に 等から定める陰解法によって実現される[51]

幾何学的解法

解析対象となる微分方程式が何らかの特別な性質を持つとき、汎用の数値積分法ではなく、その性質を尊重するように構成された数値解法を用いることがあり、そのような手法は構造保存型解法または幾何学的数値解法と呼ばれる[52][53]。例えば古典力学ハミルトン力学)において時間発展はシンプレクティック写像であり、エネルギーハミルトニアン)等の保存量が存在する[53]。この場合にはシンプレクティック数値積分法というハミルトン系にのみ適用可能な数値解法が存在し、良好な性質を持つことが知られている[54]

境界値問題

常微分方程式

の解 で, に対する境界条件

を満足するものを求める問題を境界値問題 (boundary-value problem) と呼ぶ[55]。初期値問題と異なり、境界値問題では複数の解が存在すること、あるいは解が存在しないことがあり得る。また、スツルム=リウヴィル型微分方程式のように常微分方程式にパラメータ が含まれ、境界値問題に解が存在するようにパラメータ を同時に定める問題もある[56]

これらの問題を数値的に解く最も単純な方法が狙い撃ち法 (shooting method) である[57]。一方の端点(例えば )において初期条件 を適当に定めて微分方程式をもう一方の端点 まで解き、境界条件を満たすように未定の初期条件(および固有値)を適切に選ぶ[58][59]。これはニュートン法などの関数の根を求めるアルゴリズム(これはしばしば逐次反復を伴う)を常微分方程式の数値解法と組み合わせることを意味する[60][61]。ただし初期条件によっては区間 全体で定義された解が存在しないことがあり、そのために改良された手法がある[62]。あるいは有限要素法 (finite element method) などの手法も用いられる[63]

解の存在検証

高精度に解く技術が追求されている一方で、「計算機で解の存在を検証する」という研究もおこなわれている[64]。このような研究が必要となるのは、近似解が求まったとしてもそれが幻影解である危険性があるからである。偏微分方程式ではすでに幻影解が報告されているので[64][B 1][B 2]常微分方程式でも警戒が必要である。偏微分方程式の時と同様に関数解析学的な手法[64][65][66][67]も考えられるが、関数解析学に頼らない手法 (例えば狙い撃ち法スペクトル法アフィン演算など) に基づく研究が主流であり[64]、欧米などの海外[B 3][B 4][B 5][B 6][B 7][B 8][B 9][B 10][B 11]のみならず日本国内でも研究されている[B 12][B 13][B 14][B 15]。また、爆発解 (: Blow-up solution) に特化した精度保証付き解法も探求されている[B 16][B 17]

解の存在検証・計算機援用証明が行われた方程式

関連ソフトウェア・ライブラリ

出典

  1. ^ a b c d 山本哲朗『数値解析入門』(増訂版)サイエンス社〈サイエンスライブラリ 現代数学への入門 14〉、2003年6月。ISBN 4-7819-1038-6 
  2. ^ a b c d 森正武. 数値解析 第2版. 共立出版.
  3. ^ 時弘哲治、工学における特殊関数共立出版
  4. ^ 坂井秀隆. (2015). 常微分方程式. 東京大学出版会.
  5. ^ 遠藤理平 (2018) ルンゲ・クッタで行こう!―物理シミュレーションを基礎から学ぶ
  6. ^ 岡本和夫. (2009). パンルヴェ方程式. 岩波書店.
  7. ^ 野海正俊. (2000). パンルヴェ方程式-対称性からの入門. すうがくの風景 4. 朝倉書店.
  8. ^ 岡本和夫. (1985). パンルヴェ方程式序説. 上智大学数学講究録, 19.
  9. ^ リッカチのひ・み・つ : 解ける微分方程式の理由を探る. 井ノ口順一著. 日本評論社, 2010.9.
  10. ^ Hairer, E., Lubich, C., & Wanner, G. (2006). Geometric numerical integration: structure-preserving algorithms for ordinary differential equations. en:Springer Science & Business Media.
  11. ^ 吉田春夫. (1995). シンプレクティック数値解法 (古典力学の輝き--未解決問題と新しい発見). 数理科学, 33(6), p37-46.
  12. ^ Stoer & Bulirsch, p. 465.
  13. ^ Stoer & Bulirsch, p. 467.
  14. ^ Iserles, pp. 106-107.
  15. ^ a b 加古富志雄. “数値解析”. 2021年2月2日閲覧。
  16. ^ Stoer & Bulirsch, p. 473.
  17. ^ 齊藤, pp. 97-98.
  18. ^ Press et al., p. 709.
  19. ^ Jun Makino. “5.6 刻み幅調節と埋め込み型公式”. 2021年2月2日閲覧。
  20. ^ 齊藤, pp. 99, 102.
  21. ^ Stoer & Bulirsch, pp. 473-474.
  22. ^ Hackbusch, p. 71.
  23. ^ Stoer & Bulirsch, p. 474.
  24. ^ Hackbusch, p. 72.
  25. ^ 齊藤, p. 103.
  26. ^ Stoer & Bulirsch, p. 474.
  27. ^ 齊藤, p. 102.
  28. ^ Stoer & Bulirsch, p. 477.
  29. ^ Hackbusch, p. 72.
  30. ^ Stoer & Bulirsch, pp. 477-479.
  31. ^ Hackbusch, pp. 73-74.
  32. ^ Hackbusch, p. 74.
  33. ^ 齊藤, p. 107.
  34. ^ 齊藤, p. 97-98, 103.
  35. ^ 齊藤, p. 106.
  36. ^ Stoer & Bulirsch, pp. 475-476.
  37. ^ Press et al., p. 709.
  38. ^ C言語による数値計算プログラミング演習 10. 常微分方程式の初期値問題”. 2021年2月2日閲覧。
  39. ^ Press et al., pp. 708-709, 712.
  40. ^ Stoer & Bulirsch, p. 492.
  41. ^ Stoer & Bulirsch, p. 492.
  42. ^ Press et al., pp. 747-750.
  43. ^ Stoer & Bulirsch, p. 508.
  44. ^ Iserles, pp. 56-57.
  45. ^ Jun Makino (1998年). “2 数値解の安定性”. 2021年2月1日閲覧。
  46. ^ Iserles, pp. 56-57.
  47. ^ Iserles, p. 56.
  48. ^ Iserles, p. 59.
  49. ^ a b Jun Makino (1998年). “11.2 A-安定性”. 2021年2月1日閲覧。
  50. ^ Iserles, pp. 53-56.
  51. ^ Press et al., pp. 735-736.
  52. ^ 齊藤, p. 120.
  53. ^ a b 宮武勇登. “保存則に即した数値計算手法”. 2021年2月2日閲覧。
  54. ^ 吉田春夫. “可変時間ステップによるシンプレクティック数値解法”. 2021年2月2日閲覧。
  55. ^ Stoer & Bulirsch, pp. 539-540.
  56. ^ Stoer & Bulirsch, p. 541.
  57. ^ Stoer & Bulirsch, pp. 542-548.
  58. ^ Press et al., pp. 753-755, 757-759.
  59. ^ Hairer et al. (1993), p. 105.
  60. ^ Hairer et al. (1993), p. 106.
  61. ^ Press et al., pp. 754-755.
  62. ^ Hairer et al. (1993), pp. 106-107.
  63. ^ Iserles, pp. 171-176.
  64. ^ a b c d 精度保証付き数値計算の基礎』大石進一 編著、コロナ社、2018年。
  65. ^ M. Nakao, M. Plum, Y. Watanabe (2019) Numerical Verification Methods and Computer-Assisted Proofs for Partial Differential Equations (Springer Series in Computational Mathematics).
  66. ^ 中尾充宏, & 山本野人. (1998). 精度保証付き数値計算 チュートリアル: 応用数理最前線.
  67. ^ 中尾充宏, & 渡部善隆. (2011). 実例で学ぶ精度保証付き数値計算, サイエンス社.
  68. ^ 大石進一、非線形解析入門、コロナ社。

近似解法に関する論文

  1. ^ Butcher, J. C. (1996). A history of Runge-Kutta methods. Applied numerical mathematics, 20(3), 247-260.
  2. ^ Hochbruck, M., & Ostermann, A. (2010). Exponential integrators. en:Acta Numerica, 19, 209-286.
  3. ^ Al-Mohy, A. H., & Higham, N. J. (2011). Computing the action of the matrix exponential, with an application to exponential integrators. en:SIAM journal on scientific computing, 33(2), 488-511.
  4. ^ Monroe, J. L. (2002). Extrapolation and the Bulirsch-Stoer algorithm. Physical Review E, 65(6), 066116.
  5. ^ Kirpekar, S. (2003). Implementation of the Bulirsch Stoer extrapolation method. Department of Mechanical Engineering, UC Berkeley/California.
  6. ^ Symplectic integrators: An introduction, American Journal of Physics 73, 938 (2005); https://doi.org/10.1119/1.2034523 Denis Donnelly.
  7. ^ Y. B. Suris, Hamiltonian Runge-Kutta type methods and their variational formulation (1990) Matematicheskoe modelirovanie, 2(4), 78-87.
  8. ^ Iserles, A., & Quispel, G. R. W. (2016). Why geometric integration?. arXiv preprint arXiv:1602.07755.
  9. ^ 平山弘, 小宮聖司, & 佐藤創太郎. (2002). Taylor 級数法による常微分方程式の解法. 日本応用数理学会論文誌, 12(1), 1-8.
  10. ^ 平山弘. (2013). Taylor 展開法による常微分方程式の高次並列計算. 研究報告ハイパフォーマンスコンピューティング (HPC), 2013(3), 1-6.
  11. ^ 平山弘, & 佐藤創太郎. (2002). 遅延微分方程式の級数による解法 (Computer Algebra: Algorithms, Implementations and Applications).
  12. ^ Hirayama, H. (2002). Solution of ordinary differential equations by Taylor series method. JSIAM, 12, 1-8.
  13. ^ Hirayama, H. (2015). Performance of a Higher-Order Numerical Method for Solving Ordinary Differential Equations by Taylor Series. In Integral Methods in Science and Engineering (pp. 321-328). Birkhäuser, Cham.

精度保証・計算機援用証明に関する論文

  1. ^ Breuer, B., Plum, M., & McKenna, P. J. (2001). "Inclusions and existence proofs for solutions of a nonlinear boundary value problem by spectral numerical methods." In Topics in Numerical Analysis (pp. 61–77). Springer, Vienna.
  2. ^ Gidas, B., Ni, W. M., & Nirenberg, L. (1979). "Symmetry and related properties via the maximum principle." Communications in Mathematical Physics, 68(3), 209–243.
  3. ^ a b c Lohner,R.J.,Enclosing the Solution of Ordinary lnitial and Boundary Value Problems, Computer arithmetic:Scientific Computation and Programming Languages,Kaucher,E.,Kulisch,U., Ullrich,Ch.(eds.), B.G.Teubner,Stuttgart (1987), 255−286.
  4. ^ Rihm, R. (1994). Interval methods for initial value problems in ODEs. Topics in Validated Computations, 173-207.
  5. ^ Hungria, A., Lessard, J. P., & Mireles-James, J. D. (2014). Radii polynomial approach for analytic solutions of differential equations: Theory, examples, and comparisons. Math. Comp.
  6. ^ Nedialkov, N. S., Jackson, K. R., & Pryce, J. D. (2001). An effective high-order interval method for validating existence and uniqueness of the solution of an IVP for an ODE. Reliable Computing, 7(6), 449-465.
  7. ^ Corliss, G. F. (1989). Survey of interval algorithms for ordinary differential equations. Applied Mathematics and Computation, 31, 112-120.
  8. ^ Nedialkov, N. S. (2000). Computing rigorous bounds on the solution of an initial value problem for an ordinary differential equation (Ph.D. thesis). University of Toronto.
  9. ^ Eijgenraam, P. (1981). The solution of initial value problems using interval arithmetic: formulation and analysis of an algorithm. MC Tracts.
  10. ^ Nedialkov, N. S., & Jackson, K. R. (1999). An interval Hermite-Obreschkoff method for computing rigorous bounds on the solution of an initial value problem for an ordinary differential equation. Reliable Computing, 5(3), 289-310.
  11. ^ Nedialkov, N. S., Jackson, K. R., & Corliss, G. F. (1999). Validated solutions of initial value problems for ordinary differential equations. Applied Mathematics and Computation, 105(1), 21-68.
  12. ^ Berz, M., & Makino, K. (1998). Verified integration of ODEs and flows using differential algebraic methods on high-order Taylor models. Reliable computing, 4(4), 361-369.
  13. ^ 柏木啓一郎, & 柏木雅英. (2011). 平均値形式とアフィン演算を用いた常微分方程式の精度保証法. 日本応用数理学会論文誌, 21(1), 37-58.
  14. ^ Kashiwagi, M. (1995). Numerical Validation for Ordinary Differential Equations using Power Series Arithmetic. In Numerical Analysis Of Ordinary Differential Equations And Its Applications (pp. 213-218).
  15. ^ 相馬隆郎, & 大石進一. (2003). 精度保証付き数値計算法を用いた常微分方程式の全解探索アルゴリズム. 電子情報通信学会論文誌 A, 86(6), 663-673.
  16. ^ Takayasu, A., Matsue, K., Sasaki, T., Tanaka, K., Mizuguchi, M., & Oishi, S. I. (2017). Numerical validation of blow-up solutions of ordinary differential equations. en:Journal of Computational and Applied Mathematics, 314, 10-29.
  17. ^ Matsue, K., & Takayasu, A. (2019). Rigorous numerics of blow-up solutions for ODEs with exponential nonlinearity. arXiv preprint arXiv:1902.01842.
  18. ^ Hassard, B., Zhang, J., Hastings, S. P., & Troy, W. C. (1994). A computer proof that the Lorenz equations have “chaotic” solutions. Applied Mathematics Letters, 7(1), 79-83.
  19. ^ Mischaikow, K., & Mrozek, M. (1995). Chaos in the Lorenz equations: a computer-assisted proof. en:Bulletin of the American Mathematical Society, 32(1), 66-72.
  20. ^ Mischaikow, K., & Mrozek, M. (1998). Chaos in the Lorenz equations: A computer assisted proof. Part II: Details. en:Mathematics of Computation, 67(223), 1023-1046.
  21. ^ Mischaikow, K., Mrozek, M., & Szymczak, A. (2001). Chaos in the lorenz equations: A computer assisted proof part iii: Classical parameter values. Journal of Differential Equations, 169(1), 17-56.
  22. ^ Galias, Z., & Zgliczyński, P. (1998). Computer assisted proof of chaos in the Lorenz equations. Physica D: Nonlinear Phenomena, 115(3-4), 165-188.
  23. ^ Tucker, W. (1999). The Lorenz attractor exists. Comptes Rendus de l'Académie des Sciences-Series I-Mathematics, 328(12), 1197-1202.
  24. ^ Zgliczynski, P. (1997). Computer assisted proof of chaos in the Rössler equations and in the Hénon map. Nonlinearity, 10(1), 243.
  25. ^ 大石進一. (1993). 非線形常微分方程式の周期解の数値的存在検証と近似解の精度保証. 電子情報通信学会技術研究報告. CAS, 回路とシステム, 93(102), 91-96.
  26. ^ 大石進一. (1993). 非線形常微分方程式の周期解の数値的存在検証と近似解の精度保証 (常微分方程式系の数値解析とその周辺). 京都大学数理解析研究所講究録.
  27. ^ Makino, K., & Berz, M. (2006). Cosy infinity version 9. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 558(1), 346-350.
  28. ^ Berz, M., Makino, K., Shamseddine, K., Hoffstätter, G. H., & Wan, W. (1996). 32. COSY INFINITY and Its Applications in Nonlinear Dynamics.
  29. ^ S.M. Rump: INTLAB - INTerval LABoratory. In Tibor Csendes, editor, Developments in Reliable Computing, pages 77-104. Kluwer Academic Publishers, Dordrecht, 1999.
  30. ^ Overview of kv – a C++ library for verified numerical computation, Masahide Kashiwagi, SCAN 2018.

参考文献

和書

洋書

  • Mitsui, T., & Shinohara, Y. (1995). Numerical analysis of ordinary differential equations and its applications. en:World Scientific.
  • Iserles, A. (2009). A first course in the numerical analysis of differential equations. Cambridge University Press.
  • Hairer, Ernst; Nørsett, Syvert Paul; Wanner, Gerhard (1993), Solving ordinary differential equations I: Nonstiff problems, Berlin, New York: Springer-Verlag, ISBN 978-3-540-56670-0. (日本語版が丸善出版から発売されている、三井斌友が翻訳を担当)
  • Wanner, G. & Hairer, E. (1996), Solving ordinary differential equations II: Stiff and differential-algebraic problems (2nd ed.). Springer Berlin Heidelberg. (日本語版が丸善出版から発売されている、三井斌友が翻訳を担当)
  • Butcher, John C. (2008), Numerical Methods for Ordinary Differential Equations, New York: John Wiley & Sons, ISBN 978-0-470-72335-7.
  • John D. Lambert, Numerical Methods for Ordinary Differential Systems, John Wiley & Sons, Chichester, 1991. ISBN 0-471-92990-5.
  • Deuflhard, P., & Bornemann, F. (2012). Scientific computing with ordinary differential equations. en:Springer Science & Business Media.
  • 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.
  • Press, William H.; Teukolsky, Saul A.; Vetterling, William T.; Flannery, Brian P. (1992). Numerical Recipes in C: The Art of Scientific Computing (2nd ed.). Cambridge University Press. doi:10.2277/0521431085. ISBN 978-0-521-43108-8 
  • Stoer, Josef; Bulirsch, R. (2002). Introduction to Numerical Analysis. Springer. doi:10.1007/978-0-387-21738-3. ISBN 978-0-387-21738-3 
  • Hackbusch, Wolfgang (2014). The Concept of Stability in Numerical Mathematics. Springer. doi:10.1007/978-3-642-39386-0. ISBN 978-3-642-39386-0 

微分代数方程式の数値解法

  • Brenan, K. E., Campbell, S. L., & Petzold, L. R. (1996). Numerical solution of initial-value problems in differential-algebraic equations. SIAM.
  • Hairer, E., Lubich, C., & Roche, M. (2006). The numerical solution of differential-algebraic systems by Runge-Kutta methods. Springer.
  • Kunkel, P., & Mehrmann, V. (2006). Differential-algebraic equations: analysis and numerical solution. European Mathematical Society.
  • Marz, R. (1992). Numerical methods for differential algebraic equations. en:Acta Numerica, 1, 141-198.

遅延微分方程式の数値解法

  • Bellen, A., & Zennaro, M. (2013). Numerical methods for delay differential equations. Oxford University Press.
  • Zennaro, M. (1995). Delay differential equations: theory and numerics. Theory and numerics of ordinary and partial differential equations, 291-333.

外部リンク

解説記事

近似解法

精度保証

研究集会

ソフトウェア