常微分方程式の数値解法

常微分方程式の数値解法 (じょうびぶんほうていしきのすうちかいほう、) は、数値解析において常微分方程式を数値的に解く技術の総称である森正武. 数値解析 第2版. 共立出版.。

数値解法の必要性 これまで様々な自然現象 (物理現象など) を記述するために多くの常微分方程式が作られ、多くの数学者たちがその解法を探求してきたが、フックス型微分方程式時弘哲治、工学における特殊関数、共立出版。坂井秀隆. (2015). 常微分方程式. 東京大学出版会.などを除いて、式変形による計算だけで厳密に解ける常微分方程式は多くない。そのため多くの研究者たちが常微分方程式を数値的に解く技術について研究をしてきた。最も標準的な手法はルンゲ・クッタ法であり遠藤理平 (2018) ルンゲ・クッタで行こう!―物理シミュレーションを基礎から学ぶButcher, J. C. (1996). A history of Runge-Kutta methods. Applied numerical mathematics, 20(3), 247-260.、MATLABにはode45として搭載されている。しかしこれは万能なソルバーとは言えない。例えばパンルヴェ方程式岡本和夫. (2009). パンルヴェ方程式. 岩波書店.野海正俊. (2000). パンルヴェ方程式-対称性からの入門. すうがくの風景 4. 朝倉書店.岡本和夫. (1985). パンルヴェ方程式序説. 上智大学数学講究録, 19.やリッカチ方程式リッカチのひ・み・つ : 解ける微分方程式の理由を探る. 井ノ口順一著. 日本評論社, 2010.9.などは非線形性によって精度の良い計算ができず、数値実験結果だけを見ていると間違った結論 (幻影解) にたどり着く危険がある。、

線型多段法 リープ・フロッグ法 オイラー法の進化系 後退オイラー法 (en) :en:semi-implicit Euler method :en:Euler–Maruyama method (確率微分方程式に特化した解法) :en:exponential integrator (行列指数関数を使う解法)Hochbruck, M., & Ostermann, A. (2010). Exponential integrators. :en:Acta Numerica, 19, 209-286.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. 狙い撃ち法 :en:Bulirsch–Stoer algorithmMonroe, J. L. (2002). Extrapolation and the Bulirsch-Stoer algorithm. Physical Review E, 65(6), 066116.Kirpekar, S. (2003). Implementation of the Bulirsch Stoer extrapolation method. Department of Mechanical Engineering, UC Berkeley/California. シンプレクティック数値積分法Hairer, E., Lubich, C., & Wanner, G. (2006). Geometric numerical integration: structure-preserving algorithms for ordinary differential equations. :en:Springer Science & Business Media.吉田春夫. (1995). シンプレクティック数値解法 (古典力学の輝き--未解決問題と新しい発見). 数理科学, 33(6), p37-46.Symplectic integrators: An introduction, American Journal of Physics 73, 938 (2005); https://doi.org/10.1119/1.2034523 Denis Donnelly.Y. B. Suris, Hamiltonian Runge-Kutta type methods and their variational formulation (1990) Matematicheskoe modelirovanie, 2(4), 78-87.Iserles, A., & Quispel, G. R. W. (2016). Why geometric integration?. arXiv preprint arXiv:1602.07755. テイラー級数を用いる方法平山弘, 小宮聖司, & 佐藤創太郎. (2002). Taylor 級数法による常微分方程式の解法. 日本応用数理学会論文誌, 12(1), 1-8.平山弘. (2013). Taylor 展開法による常微分方程式の高次並列計算. 研究報告ハイパフォーマンスコンピューティング (HPC), 2013(3), 1-6.平山弘, & 佐藤創太郎. (2002). 遅延微分方程式の級数による解法 (Computer Algebra: Algorithms, Implementations and Applications).Hirayama, H. (2002). Solution of ordinary differential equations by Taylor series method. JSIAM, 12, 1-8.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階常微分方程式、すなわち （またはその開集合）で定義されたベクトル値連続関数 により定まる次の方程式

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

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

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

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

という上述の形に帰着できるIserles, pp. 106-107.。

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

という漸化式によって定めるStoer & Bulirsch, p. 473.。関数 の選択が数値積分スキームを選択することに対応する。離散化した時刻の差分 を刻み幅あるいはステップサイズと呼ぶ齊藤, pp. 97-98.。なお、ここでは時間ステップ は一定としたが、これを動的に決定するという手法もあるPress et al., p. 709.。

厳密解 から差分商

を導入するとき、数値積分スキーム の局所離散化誤差齊藤, pp. 99, 102. (local discretization error) は

により定義されるStoer & Bulirsch, pp. 473-474.Hackbusch, p. 71.。整合性のために の極限で局所離散化誤差（厳密にはその上限）は 0 に収束することが要求されるStoer & Bulirsch, p. 474.Hackbusch, p. 72.。さらに、この極限で局所離散化誤差が

を満足するとき、この積分スキーム は 次精度であるという齊藤, p. 103.。

一方、 を固定して とするとき、大域離散化誤差齊藤, p. 102. (global discretization error)

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

のように振る舞うことが保証されている、すなわち大域離散化誤差のオーダーは局所離散化誤差のオーダーに等しいStoer & Bulirsch, pp. 477-479.Hackbusch, pp. 73-74.。この結果はすべての整合的な一段法が で漸近的に安定であることを意味するものの、ただし現実的に可能な で1段法が安定であることは必ずしも保証されない（#安定性節を参照）Hackbusch, p. 74.。

積分スキーム としては以下のものが知られている。 オイラー法齊藤, p. 97-98, 103.: とするもの。これは1次精度のスキームである。 齊藤, p. 106.: とするもの。これは2次精度のスキームであり、関数 の評価を2回必要とする。 古典的ルンゲ＝クッタ法: これは4次精度のスキームであり、関数 の評価を4回必要とするStoer & Bulirsch, pp. 475-476.。 このうち古典的ルンゲ＝クッタ法は適用可能範囲の広さやプログラミングの容易さのために広く用いられている（ただし Press らは著書『Numerical Recipes』において計算速度の観点から古典的ルンゲ＝クッタ法に否定的に言及しているPress et al., pp. 708-709, 712.）。

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

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

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

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

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

境界値問題 常微分方程式

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

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

これらの問題を数値的に解く最も単純な方法が狙い撃ち法 (shooting method) であるStoer & Bulirsch, pp. 542-548.。一方の端点（例えば ）において初期条件 を適当に定めて微分方程式をもう一方の端点 まで解き、境界条件を満たすように未定の初期条件（および固有値）を適切に選ぶPress et al., pp. 753-755, 757-759.Hairer et al. (1993), p. 105.。これはニュートン法などの関数の根を求めるアルゴリズム（これはしばしば逐次反復を伴う）を常微分方程式の数値解法と組み合わせることを意味するHairer et al. (1993), p. 106.Press et al., pp. 754-755.。ただし初期条件によっては区間 全体で定義された解が存在しないことがあり、そのために改良された手法があるHairer et al. (1993), pp. 106-107.。あるいは有限要素法 (finite element method) などの手法も用いられるIserles, pp. 171-176.。

解の存在検証

高精度に解く技術が追求されている一方で、「計算機で解の存在を検証する」という研究もおこなわれている精度保証付き数値計算の基礎』大石進一 編著、コロナ社、2018年。。このような研究が必要となるのは、近似解が求まったとしてもそれが幻影解である危険性があるからである。偏微分方程式ではすでに幻影解が報告されているので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. , . 、常微分方程式でも警戒が必要である。偏微分方程式の時と同様に関数解析学的な手法中尾充宏, & 山本野人. (1998). 精度保証付き数値計算 チュートリアル: 応用数理最前線.中尾充宏, & 渡部善隆. (2011). 実例で学ぶ精度保証付き数値計算, サイエンス社.も考えられるが、関数解析学に頼らない手法 (例えば狙い撃ち法、スペクトル法やアフィン演算など) に基づく研究が主流であり、欧米などの海外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.Rihm, R. (1994). Interval methods for initial value problems in ODEs. Topics in Validated Computations, 173-207.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.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.Corliss, G. F. (1989). Survey of interval algorithms for ordinary differential equations. Applied Mathematics and Computation, 31, 112-120.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.Eijgenraam, P. (1981). The solution of initial value problems using interval arithmetic: formulation and analysis of an algorithm. MC Tracts.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.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.のみならず日本国内でも研究されている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.柏木啓一郎, & 柏木雅英. (2011). 平均値形式とアフィン演算を用いた常微分方程式の精度保証法. 日本応用数理学会論文誌, 21(1), 37-58.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).相馬隆郎, & 大石進一. (2003). 精度保証付き数値計算法を用いた常微分方程式の全解探索アルゴリズム. 電子情報通信学会論文誌 A, 86(6), 663-673.。また、爆発解 () に特化した精度保証付き解法も探求されている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.Matsue, K., & Takayasu, A. (2019). Rigorous numerics of blow-up solutions for ODEs with exponential nonlinearity. arXiv preprint arXiv:1902.01842.。 解の存在検証・計算機援用証明が行われた方程式 ローレンツ方程式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.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.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.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.Galias, Z., & Zgliczyński, P. (1998). Computer assisted proof of chaos in the Lorenz equations. Physica D: Nonlinear Phenomena, 115(3-4), 165-188.Tucker, W. (1999). The Lorenz attractor exists. Comptes Rendus de l'Académie des Sciences-Series I-Mathematics, 328(12), 1197-1202. レスラー方程式Zgliczynski, P. (1997). Computer assisted proof of chaos in the Rössler equations and in the Hénon map. Nonlinearity, 10(1), 243. ダフィング方程式大石進一、非線形解析入門、コロナ社。大石進一. (1993). 非線形常微分方程式の周期解の数値的存在検証と近似解の精度保証. 電子情報通信学会技術研究報告. CAS, 回路とシステム, 93(102), 91-96.大石進一. (1993). 非線形常微分方程式の周期解の数値的存在検証と近似解の精度保証 (常微分方程式系の数値解析とその周辺). 京都大学数理解析研究所講究録.

関連ソフトウェア・ライブラリ COSY INFINITYMakino, 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.Berz, M., Makino, K., Shamseddine, K., Hoffstätter, G. H., & Wan, W. (1996). 32. COSY INFINITY and Its Applications in Nonlinear Dynamics. INTLAB: MATLAB/GNU Octaveを使用する区間演算ライブラリであり、version 11 からODEソルバー (AWA toolbox と Taylor model toolbox の二つ) が搭載されたS.M. Rump: INTLAB - INTerval LABoratory. In Tibor Csendes, editor, Developments in Reliable Computing, pages 77-104. Kluwer Academic Publishers, Dordrecht, 1999.。 kv C++による区間演算ライブラリであり、ODEソルバーを搭載している。Overview of kv – a C++ library for verified numerical computation, Masahide Kashiwagi, SCAN 2018.

CAPD C++による区間演算ライブラリであり、力学系への応用のために開発された :en:Chebfun

出典

近似解法に関する論文

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

参考文献 和書 戸川隼人：「微分方程式の数値計算：有限要素法と差分法」、オーム社（1973年8月）。 一松信：「微分方程式と解法」、教育出版、ISBN 978-4-316-37661-5（1976年11月）。　 三井斌友：「常微分方程式の数値解法」, 岩波書店、ISBN 978-4-00-005453-9（2003年7月29日）。 三井斌友、小藤俊幸、斉藤善：「微分方程式による計算科学入門」、共立出版 、ISBN 978-4-320-01753-5 (2004年2月25日)。※ ハミルトン系に対するシンプレクティック法、遅延微分方程式、確率微分方程式が扱われている。 U.M.アッシャー、L.R.ペツォルド：「常微分方程式と微分代数方程式の数値解法」、培風館、ISBN 4-563-01125-8 (2006年7月14日)。