常微分方程式の数値解析法

常微分方程式の数値法は、常微分方程式(ODE)の解の数値近似値を求めるために使用される方法です。数値法の使用は「数値積分」とも呼ばれますが、この用語は積分の計算を指す場合もあります。

多くの微分方程式は正確に解くことができません。しかし、工学などの実用的目的においては、解の数値近似値で十分な場合が多くあります。ここで学習するアルゴリズムは、そのような近似値の計算に使用できます。別の方法としては、微積分の手法を使用して解の 級数展開を得る方法があります。

常微分方程式は、物理学、化学、生物学、経済学など、多くの科学分野で登場します。^[1]さらに、数値偏微分方程式のいくつかの方法は、偏微分方程式を常微分方程式に変換し、それを解く必要があります。

一次微分方程式は、次のような初期値問題（IVP）である。^[2]

ここで は関数 であり、初期条件は与えられたベクトルです。1 次とは、方程式に yの 1 次導関数のみが現れ、それより高次の導関数は存在しないことを意味します。${\displaystyle f}$${\displaystyle f:[t_{0},\infty )\times \mathbb {R} ^{d}\to \mathbb {R} ^{d}}$${\displaystyle y_{0}\in \mathbb {R} ^{d}}$

高次システムへの一般性を失うことなく、高次 ODE は追加の変数を導入することで、より大きな 1 次方程式のシステムに変換できるため、1 次微分方程式に限定します。たとえば、2 次方程式y '' = − yは、 y ′ = zとz ′ = − yという 2 つの 1 次方程式に書き直すことができます。

このセクションでは、境界値問題に対する数値的手法について説明し、境界値問題(BVP) には異なるツールセットが必要であることを指摘します。BVP では、解yの値または成分を複数の点で定義します。このため、BVP を解くには異なる手法を使用する必要があります。たとえば、シューティング法(およびその変種) や、有限差分法、^[3] ガラーキン法、^[4]または選点法などのグローバル法は、このクラスの問題に適しています。

ピカール・リンデレーフの定理は、f がリプシッツ連続である場合に、唯一の解が存在することを述べています。

1次の初期値問題を解く数値法は、多くの場合、^[5] 線型多段階法、またはルンゲ・クッタ法という2つの大きなカテゴリのいずれかに分類されます。さらに、方法を明示的なものおよび暗黙的なものに分けることで分類することができます。たとえば、暗黙的な線型多段階法には、アダムス・モールトン法、後方微分法(BDF)が含まれ、暗黙的なルンゲ・クッタ法^[6]には、対角暗黙ルンゲ・クッタ法 (DIRK) ^{[7] [8]}単対角暗黙ルンゲ・クッタ法 (SDIRK) [ ^9]およびガウス・ラダウ法^{[10] (}ガウス求積法^[11]に基づく) 数値法が含まれます。線型多段階法ファミリーの明示的な例としては、アダムス・バッシュフォース法があり、大まかな経験則では、硬い微分方程式には暗黙的なスキームの使用が必要であり、硬くない問題は明示的なスキームでより効率的に解くことができるとされています。

いわゆる一般線形法（GLM）は、上記の2つの大きなクラスの方法を一般化したものである。^[12]

曲線上の任意の点から、曲線の接線に沿って短い距離を移動することで、曲線上の近くの点の近似値を見つけることができます。

微分方程式（１ ）から始めて、導関数y ′を有限差分近似 で置き換える。

これを整理すると次の式が得られる。

${\displaystyle y(t+h)\approx y(t)+hy'(t)}$

そして（１）を使うと次のようになる。

この式は通常、次のように適用される。ステップサイズhを選択し、シーケンスを構築する。を正確な解の数値推定値 で表す。( 3 )に基づいて、これらの推定値を次の再帰スキーム で計算する。${\displaystyle t_{0},t_{1}=t_{0}+h,t_{2}=t_{0}+2h,...}$${\displaystyle y_{n}}$${\displaystyle y(t_{n})}$

これはオイラー法（または、後述の後退オイラー法とは対照的に前進オイラー法）です。この方法は、1768 年にこれを説明した レオンハルト オイラーにちなんで名付けられました。

オイラー法は明示的方法の一例です。これは、新しい値y _{n +1が、} y _nのような既知のものに基づいて定義されることを意味します。

（２）の代わりに近似値を用いると、

後退オイラー法は次のようになります。

後退オイラー法は暗黙的な方法であり、 y _{n +1 を}求めるために方程式を解く必要があります。これを実現するには、固定小数点反復法またはニュートン・ラプソン法（の修正版）がよく使用されます。

この方程式を解くには、明示的な方法よりも時間がかかります。使用する方法を選択する際には、このコストを考慮する必要があります。( 6 )のような暗黙的な方法の利点は、通常、硬い方程式を解くのにより安定しており、より大きなステップサイズhを使用できること です。

指数積分器は、最近多くの開発が行われている積分器の大きなクラスを表します。^[13]指数積分器の 歴史は少なくとも1960年代に遡ります。

（ １ ）の代わりに、微分方程式が次のいずれかの形であると仮定する。

あるいは、背景状態について局所的に線形化されて、線形項と非線形項が生成されます。 ${\displaystyle -Ay}$${\displaystyle {\mathcal {N}}(y)}$

指数積分器は、( 7 )を乗じて、その結果を時間間隔にわたって正確に積分することによって構築されます。 ${\textstyle e^{At}}$${\displaystyle [t_{n},t_{n+1}=t_{n}+h]}$

${\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次指数積分器は、全区間にわたって一定に保つことによって実現できます。 ${\displaystyle {\mathcal {N}}(y(t_{n}+\tau ))}$

オイラー法は、多くの場合、十分な精度がありません。より正確に言えば、オイラー法は 1 次のみです (次数の概念については後述します)。このため、数学者はより高次の方法を求めました。

1 つの可能性は、 y _{n +1 を}決定するために、以前に計算された値y _nだけでなく、より多くの過去の値に依存させて解を求めることです。これにより、いわゆるマルチステップ法が得られます。おそらく最も単純なのは、2 次で (大まかに言えば) 2 つの時間値に依存するリープフロッグ法です。

実用的な多段階法のほとんどすべては線形多段階法のファミリーに属し、その形式は次の通りである。 $${\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}}}$$

もう 1 つの可能性は、区間 内のより多くの点を使用することです。これは、カール・ルンゲとマルティン・クッタにちなんで名付けられたルンゲ・クッタ法のファミリーにつながります。彼らの 4 次法の 1 つは特に人気があります。 ${\displaystyle [t_{n},t_{n+1}]}$

ODE を解くためのこれらの方法の 1 つを適切に実装するには、時間ステップの式以上のものが必要になります。

常に同じステップ サイズを使用するのは非効率的であることが多いため、可変ステップ サイズ メソッドが開発されました。通常、ステップ サイズは、ステップあたりの (ローカル) エラーが許容レベルを下回るように選択されます。つまり、メソッドでは、ローカル エラーの推定値で あるエラー インジケータも計算する必要があります。

この考え方の拡張は、異なる順序の異なる方法の間で動的に選択することです（これは可変順序法と呼ばれます）。リチャードソン外挿法^[14]に基づく方法、例えばブリルシュ・ストアアルゴリズム^{[15] [16]}は、異なる順序のさまざまな方法を構築するためによく使用されます。

その他の望ましい機能は次のとおりです。

• 高密度出力: t ₀、t ₁、t ₂、...の点だけでなく、積分区間全体に対する安価な数値近似値。
• イベントの位置: たとえば、特定の関数が消滅する時間を見つけます。これには通常、ルート検索アルゴリズムを使用する必要があります。
• 並列コンピューティングのサポート。
• 時間に関して積分する場合、時間の可逆性

多くの方法は、ここで説明したフレームワークに該当しません。代替方法のいくつかの種類は次のとおりです。

• 多重微分法は関数fだけでなくその導関数も使用します。このクラスには、エルミート・オブレシュコフ法やフェルバーグ法のほか、パーカー・ソチャッキ法^[17]やビシュコフ・シェルバコフ法など、解yのテイラー級数の係数を再帰的に計算する手法が含まれます。
• 2 階常微分方程式の Nyström 法。すべての高階常微分方程式は、(1) の形式の 1 階常微分方程式に変換できると述べました。これは確かに正しいのですが、最善の方法ではないかもしれません。特に、Nyström 法は2 階方程式を直接処理します。
• 幾何学的積分法^{[18] [19]}は、特別なクラスの常微分方程式（例えば、ハミルトン方程式の解のためのシンプレクティック積分器）のために特別に設計されている。これらの積分法は、数値解がこれらのクラスの基礎となる構造や幾何学を尊重するように配慮している。
• 量子化状態システム法は、状態量子化の考え方に基づいた ODE 積分法のファミリーです。不連続性が頻繁に発生するスパース システムをシミュレートする場合に効果的です。

一部の IVP では、非常に高い時間分解能や非常に長い時間間隔での積分が必要となるため、従来のシリアル タイム ステップ法ではリアルタイムで実行することが計算上不可能になります (数値天気予報、プラズマ モデリング、分子動力学における IVP など)。これらの問題に対応するために、並列コンピューティングを使用してシミュレーションの実行時間を短縮する Parallel-in-Time (PinT) 法が開発されました。

初期のPinT法（最も古いものは1960年代に提案された）^[20]は、必要な並列コンピューティングアーキテクチャがまだ広く利用できなかったため、研究者に当初は見過ごされていました。利用可能なコンピューティングパワーが増えると、2000年代初頭にPararealが開発され、関心が再燃しました。Pararealは、さまざまなIVPを解くのに適した、柔軟で使いやすいPinTアルゴリズムです。エクサスケールコンピューティングの登場により、PinTアルゴリズムはますます研究の注目を集めており、世界で最も強力なスーパーコンピューターを活用できるように開発されています。2023年現在、最も人気のある方法には、Parareal、PFASST、ParaDiag、MGRITなどがあります。^[21]

数値解析は数値手法の設計だけでなく、その解析も含まれます。この解析における中心となる概念は次の 3 つです。

• 収束: 方法が解を近似するかどうか、
• 次数：解をどれだけ正確に近似しているか、そして
• 安定性：誤差が抑えられているかどうか。 ^[22]

数値解がステップサイズhが0に近づくにつれて正確な解に近づく場合、数値法は収束すると言われる。より正確には、リプシッツ関数fとt ^*  > 0 のすべての常微分方程式(1)に対して、

${\displaystyle \lim _{h\to 0^{+}}\max _{n=0,1,\dots ,\lfloor t^{*}/h\rfloor }\left\|y_{n,h}-y(t_{n})\right\|=0.}$

上記の方法はすべて収束的です。

数値計算法が

${\displaystyle y_{n+k}=\Psi (t_{n+k};y_{n},y_{n+1},\dots ,y_{n+k-1};h).\,}$

メソッドのローカル(切り捨て) エラーは、メソッドの 1 つのステップで発生するエラーです。つまり、前のステップでエラーが発生しなかったと仮定してメソッドによって返される結果と、正確な解との差です。

${\displaystyle \delta _{n+k}^{h}=\Psi \left(t_{n+k};y(t_{n}),y(t_{n+1}),\dots ,y(t_{n+k-1});h\right)-y(t_{n+k}).}$

この方法は、

${\displaystyle \lim _{h\to 0}{\frac {\delta _{n+k}^{h}}{h}}=0.}$

この方法は 、 ${\displaystyle p}$

${\displaystyle \delta _{n+k}^{h}=O(h^{p+1})\quad {\mbox{as }}h\to 0.}$

したがって、ある方法は、次数が 0 より大きい場合に矛盾がないと言えます。上で紹介した (前進) オイラー法 (4) と後退オイラー法 (6) は、どちらも次数が 1 であるため矛盾がありません。実際に使用されている方法のほとんどは、より高い次数を達成しています。矛盾がないことは収束に必要な条件ですが^[要出典]、十分ではありません。つまり、ある方法が収束するためには、矛盾がなく、かつゼロ安定性も備えていなければなりません。

関連する概念は、グローバル（打ち切り）誤差であり、これは、固定された時間 に到達するために必要なすべてのステップで発生する誤差です。明示的には、時間 におけるグローバル誤差は で、ここでです。 a次の 1 ステップ法のグローバル誤差は です。特に、このような方法は収束します。この記述は、マルチステップ法では必ずしも当てはまりません。 ${\displaystyle t}$${\displaystyle t}$${\displaystyle y_{N}-y(t)}$${\displaystyle N=(t-t_{0})/h}$${\displaystyle p}$${\displaystyle O(h^{p})}$

一部の微分方程式では、オイラー法、陽的ルンゲ・クッタ法、多段階法（アダムス・バッシュフォース法など）などの標準的な方法を適用すると、解が不安定になるが、他の方法では安定した解が得られる可能性がある。方程式（それ自体は必ずしも複雑ではない）におけるこの「困難な動作」は剛性として説明され、多くの場合、基礎となる問題に異なる時間スケールが存在することによって引き起こされる。^[23]たとえば、衝撃振動子のような機械システムにおける衝突は、通常、物体の運動の時間よりもはるかに小さい時間スケールで発生する。この不一致により、状態パラメータの曲線に非常に「急な曲がり」が生じる。

硬い問題は、化学反応速度論、制御理論、固体力学、天気予報、生物学、プラズマ物理学、電子工学のいたるところで見られます。硬直性を克服する一つの方法は、微分方程式の概念を微分包含の概念に拡張し、非滑らかさを許容してモデル化することです。^{[24] [25]}

以下はこの分野におけるいくつかの重要な進展のタイムラインです。 ^{[26] [27]}

• 1768年 -レオンハルト・オイラーが自身の方法を発表。
• 1824年 -オーギュスタン＝ルイ・コーシーがオイラー法の収束を証明。この証明では、コーシーは暗黙のオイラー法を使用している。
• 1855年 -フランシス・バッシュフォースが書いた手紙の中で、ジョン・カウチ・アダムスの多段階法が初めて言及された。
• 1895年 -カール・ルンゲが最初のルンゲ・クッタ法を発表。
• 1901年 -マルティン・クッタが人気の高い4次ルンゲ・クッタ法を説明。
• 1910年 -ルイス・フライ・リチャードソンが外挿法、リチャードソン外挿法を発表。
• 1952年 - チャールズ・F・カーティスとジョセフ・オークランド・ヒルシュフェルダーが「スティッフ方程式」という用語を作った。
• 1963年 -ゲルムンド・ダールキストが積分法のA安定性を導入。

境界値問題（BVP）は通常、元のBVPを離散化して得られるほぼ同等の行列問題を解くことによって数値的に解かれます。^[28] 1次元でBVPを数値的に解くために最も一般的に使用される方法は、有限差分法と呼ばれます。^[3] この方法は、点の値の線形結合を利用して、関数の導関数を記述する有限差分係数を構築します。たとえば、1次導関数の2次中心差分近似は次のように表されます。

${\displaystyle {\frac {u_{i+1}-u_{i-1}}{2h}}=u'(x_{i})+{\mathcal {O}}(h^{2}),}$

2次導関数の2次中心差分は次のように表される。

${\displaystyle {\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}=u''(x_{i})+{\mathcal {O}}(h^{2}).}$

これらの式の両方において、は離散化された領域上の隣接するx値間の距離です。次に、標準的な行列法で解くことができる線形システムを構築します。たとえば、解くべき方程式が次のとおりであるとします。 ${\displaystyle h=x_{i}-x_{i-1}}$

${\displaystyle {\begin{aligned}&{}{\frac {d^{2}u}{dx^{2}}}-u=0,\\&{}u(0)=0,\\&{}u(1)=1.\end{aligned}}}$

次のステップは問題を離散化し、次のような線形微分近似を使用することです。

${\displaystyle u''_{i}={\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}}$

得られた線形方程式を解きます。次のような方程式が得られます。

${\displaystyle {\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}-u_{i}=0,\quad \forall i={1,2,3,...,n-1}.}$

一見すると、この方程式系は、変数が乗算されない項を含まないという事実に関連する難しさがあるように見えますが、実際にはそうではありません。i = 1 および n − 1 には境界値を含む項があり、これら 2 つの値は既知であるため、この方程式にそれらを代入するだけで、結果として、非自明な解を持つ非同次線形方程式系が得られます。 ${\displaystyle u(0)=u_{0}}$${\displaystyle u(1)=u_{n}}$

• クーラント・フリードリヒス・レヴィ条件
• エネルギードリフト
• 一般線形法
• 数値解析トピックの一覧#常微分方程式の数値解析法
• 可逆参照システム伝播アルゴリズム
• Modelica言語とOpenModelicaソフトウェア

• ブレイディ、ブライアン (2006)。『数値解析入門』アッパーサドルリバー、ニュージャージー州：ピアソン・プレンティス・ホール。ISBN 978-0-13-013054-9。
• JC Butcher著『常微分方程式の数値解析法』ISBN 0-471-96758-0 
• Ernst Hairer、Syvert Paul Nørsett、Gerhard Wanner、『常微分方程式の解法 I: ノンスティッフ問題』第2 版、Springer Verlag、ベルリン、1993 年。ISBN 3-540-56670-8。 
• Ernst Hairer および Gerhard Wanner、「常微分方程式の解法 II: スティフ問題と微分代数問題」、第 2 版、Springer Verlag、ベルリン、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.
• Arieh Iserles、『微分方程式の数値解析入門』、ケンブリッジ大学出版局、1996 年。ISBN 0-521-55376-8 (ハードカバー)、ISBN 0-521-55655-4 (ペーパーバック)。 (数学の上級学部生および大学院生を対象とした教科書で、数値偏微分方程式についても説明しています。)  
  
• John Denholm Lambert 著『Numerical Methods for Ordinary Differential Systems』、 John Wiley & Sons、Chichester、1991 年。ISBN 0-471-92990-5 。 (教科書。Iserles 著より少し難易度が高い。) 
  

• 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). 「カオスシステムの解を構築するための正確な数値手法とアルゴリズム」.応用非線形力学ジャーナル. 9 (2): 207–221. arXiv : 2011.10664 . doi :10.5890/JAND.2020.06.004. S2CID  225853788.
• GitHubの kv (厳密な ODE ソルバーを備えたC++ライブラリ)
• INTLAB (厳密な ODE ソルバーを含むMATLAB / GNU Octave製のライブラリ)