常微分方程式の数値解法

テンプレート:Differential equations 常微分方程式の数値解法 (じょうびぶんほうていしきのすうちかいほう、テンプレート:Lang-en-short) は、数値解析において常微分方程式を数値的に解く技術の総称である。

数値解法の必要性

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

• 線型多段法
• リープ・フロッグ法
• オイラー法の進化系
  • 後退オイラー法 (en)
  • en:semi-implicit Euler method
  • en:Euler–Maruyama method (確率微分方程式に特化した解法)
• en:exponential integrator (行列指数関数を使う解法)
• 狙い撃ち法
• en:Bulirsch–Stoer algorithm
• シンプレクティック数値積分法
• テイラー級数を用いる方法

テンプレート:Div col end などの新しい解法に関する研究が進められている。

初期値問題

$N$ 個の変数 $y=(y_1,\cdots ,y_N)\in {ℝ}^N$ に関する1階常微分方程式、すなわち $[t_0,b]\times {ℝ}^N$ （またはその開集合）で定義されたベクトル値連続関数 $f(t,y)=(f_1(t,y),\cdots ,f_N(t,y))$ により定まる次の方程式 テンプレート:Indent について考える。その初期値問題 (initial-value problem) とは、初期条件 $y(t_0)=y_0$ を満たす関数 $y(t)$ を求めることである。関数 $f$ が第2引数についてリプシッツ連続であるとき、すなわち定数 $L$ が存在し テンプレート:Indent （$\Vert \cdot \Vert$ は ${ℝ}^N$ のノルム）を満たすとき、その初期値問題には解が一意に存在する（ピカール・リンデレフの定理）。本節では常微分方程式の初期値問題を数値的に解くことについて議論する。

初期値問題の例

例として、電気回路の研究から導かれたファン・デル・ポール振動子について考える。その運動方程式は テンプレート:Indent であり、これは2階の常微分方程式であるものの、$v=dx/dt$ とおくと テンプレート:Indent という上述の形に帰着できる。

一段法

区間 $[t_0,t]$ の厳密解を $y(t)$ とする。一段法として知られるクラスの数値解法では、離散化した時刻 $t_i=t_0+ih$ ($h=(t-t_0)/n$) での厳密解 $y(t_i)$ の近似値 ${\eta }_i$ を テンプレート:Indent という漸化式によって定める。関数 $\mathrm{\Phi }$ の選択が数値積分スキームを選択することに対応する。離散化した時刻の差分 $h=t_{i+1}-t_i$ を刻み幅あるいはステップサイズと呼ぶ。なお、ここでは時間ステップ $h=t_{i+1}-t_i$ は一定としたが、これを動的に決定するテンプレート:仮リンクという手法もある。

厳密解 $y(t)$ から差分商 テンプレート:Indent を導入するとき、数値積分スキーム $\mathrm{\Phi }$ の局所離散化誤差 (local discretization error) は テンプレート:Indent により定義される。整合性のために $h\to 0$ の極限で局所離散化誤差（厳密にはその上限）は 0 に収束することが要求される。さらに、この極限で局所離散化誤差が $\mathrm{\Delta }(t,\eta ;h)-\mathrm{\Phi }(t,\eta ;h)=O(h^p)$ を満足するとき、この積分スキーム $\mathrm{\Phi }$ は $p$ 次精度であるという。

一方、$t$ を固定して $n\to \mathrm{\infty }$ とするとき、大域離散化誤差 (global discretization error) テンプレート:Indent が 0 に収束するならば、その積分スキームは収束するという。$p$ 次精度 ($p>0$) の1段法を十分に滑らかな関数 $f$ に適用するとき、そのスキームは収束し大域離散化誤差は テンプレート:Indent のように振る舞うことが保証されている、すなわち大域離散化誤差のオーダーは局所離散化誤差のオーダーに等しい。この結果はすべての整合的な一段法が $h\to 0$ で漸近的に安定であることを意味するものの、ただし現実的に可能な $h$ で1段法が安定であることは必ずしも保証されない（#安定性節を参照）。

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

積分スキーム $\mathrm{\Phi }$ としては以下のものが知られている。

• オイラー法: $\mathrm{\Phi }(t,\eta ;h)=f(t,\eta )$ とするもの。これは1次精度のスキームである。
• テンプレート:仮リンク: $\mathrm{\Phi }(t,\eta ;h)=[f(t,\eta )+f(t+h,\eta +hf(t,\eta ))]/2$ とするもの。これは2次精度のスキームであり、関数 $f$ の評価を2回必要とする。
• 古典的ルンゲ＝クッタ法: これは4次精度のスキームであり、関数 $f$ の評価を4回必要とする。

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

多段法

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

安定性

数値積分スキームの安定性はしばしばそれを初期値問題 テンプレート:Indent （$\lambda$ は複素数の定数）に適用することで定量化される。問題のスキームを一定の刻み幅 $h$ で適用することで得られる数列 $\{y_n{\}}_{n=0,1,\cdots }$ について、それが極限 $n\to \mathrm{\infty }$ で 0 に収束するような $h\lambda$ がなす複素平面上の領域を絶対安定領域と呼ぶ。もとの初期値問題の厳密解は $y(t)=e^{\lambda t}$ であり、これは $\mathrm{R}\mathrm{e}\{\lambda \}<0$ のとき $t\to \mathrm{\infty }$ で 0 に収束する。そこである積分スキームの絶対安定領域が左半平面を含むとき、そのスキームはA-安定 (A-stable) であるという。

関数 $f(t,y)$ が大きな値を取り積分スキームによっては極端に時間刻み幅 $h$ を小さくする必要がある場合、その常微分方程式は硬い方程式と呼ばれる。この場合には十分に安定なスキームを用いる必要が生じる。これはしばしば $y_{i+1}$ を陰関数的に $y_i$ 等から定める陰解法によって実現される。

幾何学的解法

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

境界値問題

常微分方程式 テンプレート:Indent の解 $y(x)$ で, $a\ne b$ に対する境界条件 テンプレート:Indent を満足するものを求める問題を境界値問題 (boundary-value problem) と呼ぶ。初期値問題と異なり、境界値問題では複数の解が存在すること、あるいは解が存在しないことがあり得る。また、スツルム＝リウヴィル型微分方程式のように常微分方程式にパラメータ $\lambda$ が含まれ、境界値問題に解が存在するようにパラメータ $\lambda$ を同時に定める問題もある。

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

解の存在検証

テンプレート:Seealso 高精度に解く技術が追求されている一方で、「計算機で解の存在を検証する」という研究もおこなわれている。このような研究が必要となるのは、近似解が求まったとしてもそれが幻影解である危険性があるからである。偏微分方程式ではすでに幻影解が報告されているので、常微分方程式でも警戒が必要である。偏微分方程式の時と同様に関数解析学的な手法も考えられるが、関数解析学に頼らない手法 (例えば狙い撃ち法、スペクトル法やアフィン演算など) に基づく研究が主流であり、欧米などの海外のみならず日本国内でも研究されている。また、爆発解 (テンプレート:Lang-en-short) に特化した精度保証付き解法も探求されている。

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

• ローレンツ方程式
• レスラー方程式
• ダフィング方程式

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

• COSY INFINITY
• INTLAB: MATLAB/GNU Octaveを使用する区間演算ライブラリであり、version 11 からODEソルバー (AWA toolbox と Taylor model toolbox の二つ) が搭載された。
• kv C++による区間演算ライブラリであり、ODEソルバーを搭載している。
  • テンプレート:GitHub
• 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日)。
• E.ハイラー、S.P.ネルセット、G.バンナー, 三井斌友(訳)：「常微分方程式の数値解法 I　(基礎編）」、丸善出版、ISBN 978-4-621-06282-1（2007年12月）。
• E.ハイラー、G.バンナー, 三井斌友(訳)：「常微分方程式の数値解法 II　(発展編）」、丸善出版、ISBN 978-4-621-06317-0（2008年8月）。
• テンプレート:Cite book
• 神永正博:「Pythonと実例で学ぶ微分方程式：はりの方程式から感染症の数理モデルまで」、コロナ社、ISBN 978-4-339-06123-9 (2021年10月22日)。※ 応用とPythonのライブラリの利用に主眼がある。

洋書

• 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, テンプレート:ISBN2. (日本語版が丸善出版から発売されている、三井斌友が翻訳を担当)
• 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, テンプレート:ISBN2.
• John D. Lambert, Numerical Methods for Ordinary Differential Systems, John Wiley & Sons, Chichester, 1991. テンプレート:ISBN2.
• 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.
• テンプレート:Cite book
• テンプレート:Cite book
• テンプレート:Cite book
• John C. Butcher: "B-Series : Algebraic Analysis of Numerical Methods", Springer（SSCM, volume 55), ISBN 978-3030709556 (April, 2021).
• Hairer Ernst, Lubich Christian, Wanner Gerhard: "Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations" (2nd Ed.), Springer, ISBN 978-3-540-30663-4 (2006).
• Kees Vuik, Fred Vermolen,Martin van Gijzen and Thea Vuik: "Numerical Methods for Ordinary Differential Equations", ISBN 978-94-6366-665-7 (2023年1月8日）

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

• 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.

外部リンク

• Verified ODE (IVP) Solver
• 常微分方程式で記述される非線形力学系に関する精度保証付き計算法の研究
• 硬い常微分方程式が与える連続力学系の周期解の数値計算
• 宮武勇登：「今更聞けない数値計算アルゴリズム（常微分⽅程式 編）」

解説記事

近似解法

• 常微分方程式の初期値問題の数値解法
• 常微分方程式の数値解法
• symplectic integrator in nLab
• Symplectic数値積分法
• Symplectic Integratorの紹介

精度保証

• テンプレート:PDFlink
• テンプレート:PDFlink (ファン・デル・ポール振動子を扱った実験例が掲載されている)

研究集会

• 常微分方程式の数値解法とその周辺 2019
• 常微分方程式の数値解法とその周辺 2018

ソフトウェア

• INTLAB
• Short demonstration of the AWA toolbox
• Short demonstration of the Taylor model toolbox