微分方程式の数値積分の例y ′ = y 、 y （ 0 ） = 1. {\displaystyle y'=y,y(0)=1.} 青：オイラー法
赤：正確な解：
y = e t y=e^{t}} 。
h = 1.0 {\displaystyle h=1.0} 同じイラストh = 0.25。 {\displaystyle h=0.25.} h → 0 {\displaystyle h\to 0} 常微分方程式の数値解法とは、 常微分方程式 （ODE）の解の数値 近似値を求めるために用いられる手法である。この手法は「数値積分 」とも呼ばれるが、この用語は積分 の計算を指す場合もある。
多くの微分方程式は厳密に解くことができません。しかし、 工学などの実用的な目的においては、 解の数値近似値で十分な場合がよくあります。ここで研究したアルゴリズムは 、そのような近似値を計算するのに使用できます。別の方法としては、微積分学 の手法を用いて解の級数展開 を得る方法があります。
常微分方程式は、物理学 、化学 、生物学 、経済学 など、多くの科学分野で出現します。[ 1 ] 数値偏微分方程式の一部の手法では、 偏微分方程式を 常微分方程式に変換し、それを解く必要があります。
問題 1階微分方程式は、次の形式の初期値問題 （IVP）である。[ 2 ]
どこf {\displaystyle f} f : [ t 0 、 ∞ ） × R d → R d {\displaystyle f:[t_{0},\infty )\times \mathbb {R} ^{d}\to \mathbb {R} ^{d}} y 0 ∈ R d \displaystyle y_{0}\in \mathbb {R} ^{d}} 次とは、 方程式にy の1階微分のみが現れ、それ以上の高階微分は存在しないことを意味します。
高階システムへの一般性を失うことなく、ここでは1階 微分方程式に限定します。なぜなら、高階常微分方程式は、追加の変数を導入することで、より大きな1階方程式のシステムに変換できるからです。例えば、2階方程式y ′′ = − y は y ′ = z z ′ = − y に書き換えることができます。
このセクションでは、初期値問題 (IVP) の数値解法について説明し、境界値問題 y の値、つまり成分を複数の点で定義します。そのため、BVP を解くには異なる方法を使用する必要があります。たとえば、シューティング法 (およびその変種) や、有限差分法 [ 3 ] ガラーキン法 [ 4 ] 、コロケーション法 などの大域的手法が、この種の問題に適しています。
ピカール・リンデレーフの定理は、 fが リプシッツ連続 である場合、一意の解が存在することを述べている。
方法 1 階 IVP を解くための数値解法は、多くの場合、大きく 2 つのカテゴリに分類されます。[ 5 ] 線形多段階法 、またはRunge–Kutta 法です 。さらに、法を明示的法と暗黙的法に分けることで、より細分化できます。たとえば、暗黙的線形多段階法には Adams–Moulton 法 や後退微分法 (BDF)が含まれますが、暗黙的 Runge–Kutta 法 [ 6 ] [ 7 ] [ 8 ] [ 9 ] [ 10 ] (ガウス求積法 [ 11 ] 線形多段階法の明示的な例としては 、Adams–Bashforth 法が あり、下対角のButcher タブローを持つ Runge–Kutta 法は すべて明示的 です。大まかな経験則として、硬い 微分方程式には陰解法を用いる必要があり、一方、硬くない問題は陽解法を用いることでより効率的に解くことができる。
いわゆる一般線形法 （GLM）は、上記の2つの大きなクラスの方法を一般化したものである。[ 12 ]
オイラー法 曲線上の任意の点から、曲線に接する 直線に沿って少し移動することで、曲線上の近くの点の近似値を求めることができます。
微分方程式（1 ）から始めて、導関数y ′ を有限差分 近似で置き換えます。
これを整理すると、次の式が得られる。 y （ t + h ） ≈ y （ t ） + h y ′ （ t ） {\displaystyle y(t+h)\approx y(t)+hy'(t)} １
この式は通常、次のように適用されます。ステップサイズh を選択し、数列を構築します。t 0 、 t 1 = t 0 + h 、 t 2 = t 0 + 2 h 、 … {\displaystyle t_{0},t_{1}=t_{0}+h,t_{2}=t_{0}+2h,\dots } y n \displaystyle y_n}} y （ t n ） {\displaystyle y(t_{n})} ３ 再帰的 スキームによってこれらの推定値を計算する。
これはオイラー法 前進オイラー法 後退オイラー法 とは対照的である）である。この方法は、1768年にこれを考案したレオンハルト・オイラーにちなんで名付けられた。
オイラー法は明示的な y n +1が、 y n
後退オイラー法 ( 2
後退オイラー法 が得られます。
後退オイラー法は陰 解法であり、 y n +1固定点反復法や ニュートン・ラフソン法 （またはその改良版）がよく用いられる。
この方程式を解くには、明示的な方法よりも時間がかかります。使用する方法を選択する際には、このコストを考慮に入れる必要があります。( 6 )のような暗黙的な方法の利点は、通常、硬い方程式を 解くのに安定していることであり、つまり、より大きなステップサイズh を使用できます。
一次指数積分法 指数積分器は、近年多くの発展を遂げた積分器の大きなクラスを説明する。[ 13 ]
（ １
あるいは、背景状態の周りで局所的に線形化されて線形項が生成されている。− A y {\displaystyle -Ay} N （ y ） {\displaystyle {\mathcal {N}}(y)}
指数積分器は、( 7 e A t {\textstyle e^{At}} [ t n 、 t n + 1 ] {\displaystyle [t_{n},t_{n+1}]} t n + 1 = t n + h {\displaystyle t_{n+1}=t_{n}{+}h} y n + 1 = e − A h y n + ∫ 0 h e − （ h − τ ） A N （ y （ t n + τ ） ） d τ 。 {\displaystyle y_{n+1}=e^{-Ah}y_{n}+\int _{0}^{h}e^{-(h-\tau )A}{\mathcal {N}}{\left(y\left(t_{n}+\tau \right)\right)}\,d\tau .}
1次指数積分器は、保持することによって実現できます。N （ y （ t n + τ ） ） {\displaystyle {\mathcal {N}}(y(t_{n}+\tau ))}
一般化 オイラー法はしばしば十分な精度が得られない。より正確に言えば、オイラー法の次数は1である（次数 の概念については後述する）。このため、数学者たちはより高次の方法を模索するようになった。
一つの可能性として、以前に計算された値y n を使ってy n +1多段階法 である。おそらく最も単純なのは、2次精度で（大まかに言えば）2つの時間値に依存するリープフロッグ法だろう。
実用的な多段階法のほとんどすべては、線形多段階法 のファミリーに属し、その形式は次のようになる。 α k y n + k + α k − 1 y n + k − 1 + ⋯ + α 0 y n = h [ β k f （ t n + k 、 y n + k ） + β k − 1 f （ t n + k − 1 、 y n + k − 1 ） + ⋯ + β 0 f （ t n 、 y n ） ] 。 {\displaystyle {\begin{aligned}&{}\alpha _{k}y_{n+k}+\alpha _{k-1}y_{n+k-1}+\cdots +\alpha _{0}y_{n}\\&{}\quad =h\left[\beta _{k}f(t_{n+k},y_{n+k})+\beta _{k-1}f(t_{n+k-1},y_{n+k-1})+\cdots +\beta _{0}f(t_{n},y_{n})\right].\end{aligned}}}
別の可能性としては、区間内のより多くの点を使用することです。[ t n 、 t n + 1 ] {\displaystyle [t_{n},t_{n+1}]} これは、カール・ルンゲ とマルティン・クッタ にちなんで名付けられたルンゲ・クッタ法と呼ばれる手法 群につながります。彼らの4次法の一つは特に広く用いられています。
高度な機能 常微分方程式を解くためのこれらの方法のいずれかを適切に実装するには、時間ステップ公式以上のものが必要となる。
常に同じステップサイズを使用するのは非効率的な場合が多いため、可変ステップサイズ法 が開発されてきた。通常、ステップサイズは、各ステップにおける（局所的な）誤差が許容レベルを下回るように選択される。つまり、これらの方法では、局所誤差の推定値である誤差指標も計算する必要がある。
このアイデアの拡張として、異なる次数を持つさまざまな方法を動的に選択する方法があります（これは可変次数法と呼ばれます）。 リチャードソン外挿法 [ 14 ] ブリルシュ・ストーアアルゴリズム [ 15 ] [ 16 ]
その他の望ましい機能としては、以下のようなものがあります。
高密度出力：積分区間全体に対する安価な数値近似値であり、 t 0 、t 1 、t 2 、...の点だけでなく、積分区間全体に対する数値近似値である。イベント位置特定 ：例えば、特定の関数が消滅する時刻を特定すること。これは通常、根探索アルゴリズム の使用を必要とする。並列コンピューティング のサポート。時間に関して積分する場合、時間反転性
代替方法 ここで議論されている枠組みには当てはまらない方法も数多く存在する。代替方法の例としては、以下のようなものがある。
分析 数値解析 とは、数値計算手法の設計だけでなく、その解析も含む。この解析における3つの中心的な概念は以下のとおりである。
収束性 ：その方法が解を近似するかどうか、次数 ：解をどれだけよく近似しているか、そして安定性 [ 22 ]
収束 数値解がステップサイズh が 0 に近づくにつれて厳密解に近づく場合、数値解法は収束する と言われます。より正確には、リプシッツ 関数f を持つすべての ODE (1)とすべてのt * > 0に対して、
リム h → 0 + 最大 n = 0 、 1 、 … 、 ⌊ t * / h ⌋ ‖ y n 、 h − y （ t n ） ‖ = 0. {\displaystyle \lim _{h\to 0^{+}}\max _{n=0,1,\dots ,\lfloor t^{*}/h\rfloor }\left\|y_{n,h}-y(t_{n})\right\|=0.}
上記の方法はすべて収束する。
安定性と剛性 一部の微分方程式では、オイラー法、明示的ルンゲ・クッタ法 、多段階法 （例えば、アダムス・バッシュフォース法）などの標準的な方法を適用すると、解に不安定性が生じますが、他の方法では安定した解が得られる場合があります。方程式におけるこの「困難な挙動」（必ずしもそれ自体が複雑であるとは限りません）は、剛性 と呼ばれ、多くの場合、基礎となる問題に異なる時間スケールが存在することによって引き起こされます。[ 23 ] 衝撃振動子 のような機械システムにおける衝突は、通常、物体の運動時間よりもはるかに短い時間スケールで発生します。この不一致により、状態パラメータの曲線に非常に「急激な変化」が生じます。
硬い問題は、化学反応速度論 、制御理論 、固体力学 、天気予報 、生物学 、プラズマ物理学 、電子工学 など、あらゆる分野で普遍的に見られます。硬さを克服する一つの方法は、微分方程式の概念を微分包含 の概念に拡張することです。微分包含は、滑らかでない状態を許容し、モデル化します。[ 24 ] [ 25 ]
2階1次元境界値問題の数値解法 境界値問題（BVP）は通常、元のBVPを離散化して得られる近似的に等価な行列問題を解くことによって数値的に解かれます。[ 28 ] 有限差分法 と呼ばれています。[ 3 ] 有限差分係数 を構築します。例えば、1階導関数の2次中心差分 近似は次のように与えられます。
u ′ （ x 私 ） = u 私 + 1 − u 私 − 1 2 h + O （ h 2 ） 、 {\displaystyle u'(x_{i})={\frac {u_{i+1}-u_{i-1}}{2h}}+{\mathcal {O}}(h^{2}),}
また、2階微分に対する2次中心差分 は次のように表される。
u 」 （ x 私 ） = u 私 + 1 − 2 u 私 + u 私 − 1 h 2 + O （ h 2 ） 。 {\displaystyle u''(x_{i})={\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}+{\mathcal {O}}(h^{2}).}
これらの式では、h = x 私 − x 私 − 1 {\displaystyle h=x_{i}-x_{i-1}} は、離散化された領域における隣接するx 値間の距離です。次に、標準的な行列法 で解くことができる線形システムを構築します。たとえば、解くべき方程式が次のようになっているとします。
d 2 u d x 2 = u 、 u （ 0 ） = 0 、 u （ 1 ） = 1. {\displaystyle {\begin{aligned}&{\frac {d^{2}u}{dx^{2}}}=u,\\[1ex]&u(0)=0,\\&u(1)=1.\end{aligned}}}
次のステップは、問題を離散化し、次のような線形微分近似を使用することです。
u 私 」 = u 私 + 1 − 2 u 私 + u 私 − 1 h 2 {\displaystyle u''_{i}={\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}}
そして、得られた連立一次方程式を解きます。すると、次のような方程式が得られます。
u 私 + 1 − 2 u 私 + u 私 − 1 h 2 − u 私 = 0 、 ∀ 私 = 1 、 2 、 3 、 … 、 n − 1 。 {\displaystyle {\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}-u_{i}=0,\quad \forall i={1,2,3,\dots ,n-1}.}
一見すると、この方程式系は変数で乗算されていない項が含まれていないという事実に関連する困難を抱えているように見えるが、実際にはこれは誤りである。i = 1 および n − 1 には境界 値を 含む項が存在する。u （ 0 ） = u 0 {\displaystyle u(0)=u_{0}} u （ 1 ） = u n {\displaystyle u(1)=u_{n}} 線形方程式系 が得られます。
注記 ↑ Chicone, C. (2006). 応用を伴う常微分方程式（第34巻）。Springer Science & Business Media。 ↑ ブラディ（2006年 、533～655ページ ） 1 2 LeVeque, RJ (2007). 常微分方程式および偏微分方程式に対する有限差分法：定常状態および時間依存問題（第98巻）。SIAM。 ↑ Slimane Adjerid および Mahboub Baccouch (2010) Galerkin 法。Scholarpedia、5(10):10056。 ↑ Griffiths, DF、& Higham, DJ (2010). 常微分方程式の数値解法：初期値問題。Springer Science & Business Media。 ↑ ハイラー、ノーセット、 ワナー (1993 、pp. 204–215) ↑ Alexander, R. (1977). 硬い常微分方程式に対する対角陰解法ルンゲ・クッタ法. SIAM Journal on Numerical Analysis, 14(6), 1006-1021. ↑ Cash, JR (1979). 誤差推定を伴う対角的に陰的なルンゲ・クッタ公式. IMA Journal of Applied Mathematics, 24(3), 293-301. ↑ Ferracina, L., & Spijker, MN (2008). 単対角陰解法ルンゲ・クッタ法の強い安定性。応用数値数学、58(11)、1675-1686。 ↑ Everhart, E. (1985). ガウス・ラドー間隔を用いた効率的な積分器。国際天文学連合コロキウム（第83巻、185-202頁）。ケンブリッジ大学出版局。 ↑ Weisstein, Eric W.「ガウス求積法」MathWorld（Wolfram Web Resource）より。https ://mathworld.wolfram.com/GaussianQuadrature.html ↑ Butcher, JC (1987). 常微分方程式の数値解析：ルンゲ・クッタ法と一般線形法。Wiley-Interscience。 ↑ Hochbruck & Ostermann (2010 , pp. 209–286) これは指数積分器に関する現代的で包括的なレビュー論文である。 ↑ Brezinski, C., & Zaglia, MR (2013). Extrapolation methods: theory and practice. Elsevier. ↑ Monroe, JL (2002). 外挿とBulirsch-Stoerアルゴリズム。Physical Review E、65(6)、066116。 ↑ Kirpekar, S. (2003). Bulirsch Stoer 外挿法の実装。カリフォルニア大学バークレー校機械工学科。 ↑ Nurminskii, EA、Buryi, AA (2011)。グラフィックスプロセッサを用いた常微分方程式系の解法のためのパーカー・ソチャッキ法。数値解析と応用、4(3)、223。 ↑ Hairer, E., Lubich, C., & Wanner, G. (2006). 幾何学的数値積分：常微分方程式のための構造保存アルゴリズム（第31巻）。Springer Science & Business Media。 ↑ Hairer, E., Lubich, C., & Wanner, G. (2003). Störmer–Verlet法による幾何学的数値積分。Acta Numerica、12、399-450。 ↑ Nievergelt, Jürg (1964). "並列法による常微分方程式の積分" . Communications of the ACM . 7 (12): 731–733 . doi : 10.1145/355588.365137 . S2CID 6361754 . ↑ "Parallel-in-Time.org" . Parallel-in-Time.org . 2023年 11月15日 取得 . ↑ Higham, NJ (2002). 数値アルゴリズムの精度と安定性 (Vol. 80). SIAM. ↑ Miranker, A. (2001). Numerical Methods for Stiff Equations and Singular Perturbation Problems: and singular perturbation problems (Vol. 5). Springer Science & Business Media. ↑ Markus Kunze; Tassilo Kupper (2001). "非平滑力学系：概要". Bernold Fiedler (編)『 力学系のエルゴード理論、解析、および効率的なシミュレーション 』Springer Science & Business Media、p . 431。ISBN 978-3-540-41290-8 。↑ Thao Dang (2011). 「ハイブリッドシステムのモデルベーステスト」。Justyna Zander 、 Ina Schieferdecker、Pieter J. Mosterman (編) 『組み込みシステムのためのモデルベーステスト』 CRC Press、p. 411。ISBN 978-1-4398-1845-9 。↑ Brezinski, C., & Wuytack, L. (2012). 数値解析：20世紀の歴史的発展。Elsevier。 ↑ Butcher, JC (1996). ルンゲ・クッタ法の歴史。応用数値数学、20(3)、247-260。 ↑ Ascher, UM、Mattheij, RM、Russell, RD (1995)。常微分方程式の境界値問題の数値解法。応用数理学会。
参考文献 ブラディ、ブライアン（2006）。数値解析入門 。ニュージャージー州アッパーサドルリバー：ピアソン・プレンティスホール。ISBN 978-0-13-013054-9 。 JC Butcher著 、『常微分方程式の数値解法 』 、ISBN 0-471-96758-0 Hairer, E.; Nørsett, SP; Wanner, G. (1993).常微分方程式の解法 I. 非剛性問題 . Springer Series in Computational Mathematics. Vol. 8 (第 2 版). Springer-Verlag, Berlin. ISBN 3-540-56670-8 MR 1227985 . エルンスト・ヘアラー、ゲルハルト・ワナー著『常微分方程式の解法 II：硬い問題と微分代数問題』 第2版、シュプリンガー・フェルラーク、ベルリン、1996年。ISBN 3-540-60452-9 （この2巻からなるモノグラフは、当該分野のあらゆる側面を体系的に網羅している。 ） Hochbruck, Marlis ; Ostermann, Alexander (2010 年 5 月). "指数積分器". Acta Numerica . 19 : 209– 286. Bibcode : 2010AcNum..19..209H . CiteSeerX 10.1.1.187.6794 . doi : 10.1017/S0962492910000048 . S2CID 4841957 . アリエ・イザーレス著『微分方程式の数値解析入門』 ケンブリッジ大学出版局、1996年。ISBN 0-521-55376-8 （ハードカバー）、ISBN 0-521-55655-4 （ペーパーバック版） （数学を専攻する上級学部生および大学院生を対象とした教科書で、数値偏微分方程式 についても解説している。） ジョン・デンホルム・ランバート著、『常微分方程式系の数値解法』、 ジョン・ワイリー・アンド・サンズ社、チチェスター、1991年。ISBN 0-471-92990-5 （教科書。 イゼルレスの著書よりやや難易度が高い。）
外部リンク Joseph W. Rudmin、「天体力学へのパーカー・ソチャッキ法の応用」、 1998年、2016年5月16日にポルトガル語ウェブアーカイブに アーカイブされました 。 Dominique Tournès、L'intégration approchée des équations différentielles ordinaires (1671–1914) 、パリ大学博士課程 7 - Denis Diderot、1996 年。ヴィルヌーヴ・ダスク : Presses universitaires du Septentrion、1997 年、468 ページ。 (ODE 数値解析の歴史に関する広範なオンライン資料。ODE 数値解析の歴史に関する英語資料については、たとえば、彼が引用した Chabert と Goldstine による紙の本を参照してください。) Pchelintsev, AN (2020). 「カオスシステムの解を構築するための正確な数値的方法とアルゴリズム」. Journal of Applied Nonlinear Dynamics . 9 (2): 207–221 . arXiv : 2011.10664 . doi : 10.5890/JAND.2020.06.004 . S2CID 225853788 . GitHub 上の kv（厳密な常微分方程式ソルバーを備えたC++ ライブラリ）INTLAB （ MATLAB / GNU Octave で作成された、厳密な常微分方程式ソルバーを含むライブラリ）
![{\displaystyle [t_{n},t_{n+1}]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ea7558633c567b46e0e6a00c0284e61646613660)
![{\displaystyle {\begin{aligned}&{}\alpha _{k}y_{n+k}+\alpha _{k-1}y_{n+k-1}+\cdots +\alpha _{0}y_{n}\\&{}\quad =h\left[\beta _{k}f(t_{n+k},y_{n+k})+\beta _{k-1}f(t_{n+k-1},y_{n+k-1})+\cdots +\beta _{0}f(t_{n},y_{n})\right].\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/47bb8bad1854dddf30061661932c07d46850ceff)
![{\displaystyle {\begin{aligned}&{\frac {d^{2}u}{dx^{2}}}=u,\\[1ex]&u(0)=0,\\&u(1)=1.\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/adb905f356d1c605b4d1250a44f291fb04effa89)
