微分方程式の数値積分の図解y ′ = y 、 y ( 0 ) = 1. {\displaystyle y'=y,y(0)=1.} 青:オイラー法
赤: 正確な解: 。
y = e t {\textstyle 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 ] さらに、数値偏微分方程式のいくつかの手法では、 偏微分方程式を 常微分方程式に変換し、それを解く必要があります。
問題 一次微分方程式は、初期値問題 (IVP)であり、[ 2 ]
y ′ ( t ) = f ( t 、 y ( t ) ) 、 y ( t 0 ) = y 0 、 {\displaystyle y'(t)=f(t,y(t)),\qquad y(t_{0})=y_{0},} 1
ここでは関数 であり、初期条件は与えられたベクトルです。1 次とは、 方程式に y の1次導関数のみが現れ、それ以上の高次の導関数は現れないことを意味します。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}}
高階微分方程式系への一般性を失うことなく、ここでは1階 微分方程式に限定します。これは、高階常微分方程式は追加の変数を導入することで、より大きな1階方程式系に変換できるためです。例えば、2階方程式y " = − yは、 y ′ = z とz ′ = − y という2つの1階方程式に書き直すことができます。
この節では、境界値問題 (IVP)の数値解析手法について説明し、境界値問題(BVP)では異なるツールセットが必要であることを指摘します。BVPでは、解y の成分、つまり値(複数の点)を定義します。そのため、BVPを解くには異なる手法を用いる必要があります。例えば、シューティング法(およびその派生法)や、 有限差分法 [ 3 ] 、ガラーキン法 [ 4 ] 、コロケーション法 などの大域的手法は、この種の問題に適しています。
ピカール・リンデレフの定理は、 fが リプシッツ連続 である場合、唯一の解が存在することを述べています。
方法 1 階初期値問題 (IVP) を解く数値法は、多くの場合、[ 5 ] 線型多段階法 、またはルンゲ・クッタ法 という 2 つの大きなカテゴリのいずれかに分類されます。さらに、陽的手法と陰的手法に分けることができます。たとえば、陰的線型多段階法には、 アダムス・モールトン法 、後方微分法 (BDF)が含まれ、陰的ルンゲ・クッタ法 [ 6 ] には、対角陰的ルンゲ・クッタ法 (DIRK)、[ 7 ] [ 8 ] 単対角陰的ルンゲ・クッタ法 (SDIRK)、[ 9 ] およびガウス・ラダウ法[ 10 ] (ガウス求積法 [ 11 ] に基づく) 数値法が含まれます。線型多段階法ファミリーの明示的な例としては、 アダムス・バッシュフォース法 があり、下対角のブッチャー テーブルを持つルンゲ・クッタ法はすべて 陽的 です。大まかな経験則では、硬い 微分方程式には暗黙的なスキームの使用が必要ですが、硬くない問題は明示的なスキームを使用してより効率的に解決できます。
いわゆる一般線形法 (GLM)は、上記の2つの大きなクラスの手法を一般化したものである。[ 12 ]
オイラー法 曲線上の任意の点から、曲線の接線 に沿って短い距離を移動することで、曲線上の近くの点の近似値を見つけることができます。
微分方程式(1 )から始めて、導関数y ′を有限差分 近似 で置き換える。
y ′ ( t ) ≈ y ( t + h ) − y ( t ) h 、 {\displaystyle y'(t)\approx {\frac {y(t+h)-y(t)}{h}},} 2
これを整理すると次の式が得られ 、(1 )を用いると次のようになる。 y ( t + h ) ≈ y ( t ) + h y ′ ( t ) {\displaystyle y(t+h)\approx y(t)+hy'(t)}
y ( t + h ) ≈ y ( t ) + h f ( t 、 y ( t ) ) 。 {\displaystyle y(t+h)\approx y(t)+hf(t,y(t)).} 3
この式は通常、次のように適用される。ステップサイズh を選択し、数列 を構築する。を厳密解の数値推定値 とする。( 3 ) に基づき、これらの推定値を以下の再帰的 手法 で計算する。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})}
y n + 1 = y n + h f ( t n 、 y n ) 。 {\displaystyle y_{n+1}=y_{n}+hf(t_{n},y_{n}).} 4
これはオイラー法 (または後述の後退オイラー法 とは対照的に前進オイラー法 )です。この法は、1768年にこの法を考案したレオンハルト・オイラー にちなんで名付けられました。
オイラー法は明示的 手法の一例です。これは、新しい値y n +1が、 y n のような既知の値に基づいて定義されることを意味します。
後退オイラー法 ( 2 )の代わりに近似値を用いると、
y ′ ( t ) ≈ y ( t ) − y ( t − h ) h 、 {\displaystyle y'(t)\approx {\frac {y(t)-y(th)}{h}},} 5
後退オイラー法は 次のようになります。
y n + 1 = y n + h f ( t n + 1 、 y n + 1 ) 。 {\displaystyle y_{n+1}=y_{n}+hf(t_{n+1},y_{n+1}).} 6
後退オイラー法は暗黙的な 解法であり、 y n +1 を求めるには方程式を解く必要があります。この解を求めるには、固定小数点反復法や ニュートン・ラプソン法 (の修正版)がよく用いられます。
この方程式を解くには、陽解法よりも時間がかかります。使用する手法を選択する際には、このコストを考慮する必要があります。( 6 )のような陰解法の利点は、通常、硬い方程式を 解く際により安定しており、より大きなステップサイズh を使用できることです。
一次指数積分法 指数積分器は、近年大きく発展した積分器の大きなクラスを表します。[ 13 ] その起源は少なくとも1960年代に遡ります。
( 1 )の代わりに、微分方程式が以下のいずれかの形式であると仮定する。
y ′ ( t ) = − あ y + 北 ( y ) 、 {\displaystyle y'(t)=-A\,y+{\mathcal {N}}(y),} 7
あるいは、背景状態について局所的に線形化されて、線形項と非線形項が生成されます。 − あ y {\displaystyle -Ay} 北 ( y ) {\displaystyle {\mathcal {N}}(y)}
指数積分器は、( 7 )を で乗じて、その結果を時間間隔にわたって正確に積分することによって構築されます。 ここで、 この積分方程式は正確ですが、積分を定義しません。 e あ 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 − あ h y n + ∫ 0 h e − ( h − τ ) あ 北 ( 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次指数積分器は、全区間にわたって一定に保つことで実現できます。 北 ( y ( t n + τ ) ) {\displaystyle {\mathcal {N}}(y(t_{n}+\tau ))}
y n + 1 = e − あ h y n + あ − 1 ( 1 − e − あ h ) 北 ( y ( t n ) ) 。 {\displaystyle y_{n+1}=e^{-Ah}y_{n}+A^{-1}\left(1-e^{-Ah}\right){\mathcal {N}}{\left(y(t_{n})\right)}\ .} 8
一般化 オイラー法はしばしば十分な精度を示さない。より正確に言えば、オイラー法は1階しか持たない(階数 の概念については後述)。そのため、数学者たちはより高階な解法を模索するようになった。
一つの可能性は、 y n +1 を決定する際に、以前に計算された値y n だけでなく、より多くの過去の値にも依存するようにすることです。これはいわゆるマルチステップ法 となります。おそらく最も単純なのは、 2次のリープフロッグ法 で、(大まかに言えば)2つの時間値に依存します。
ほとんどすべての実用的な多段階法は線形多段階法 のファミリーに属し、その形式は次の通りである。 α け y n + け + α け − 1 y n + け − 1 + ⋯ + α 0 y n = h [ β け f ( t n + け 、 y n + け ) + β け − 1 f ( t n + け − 1 、 y n + け − 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}]}
高度な機能 ODE を解くためのこれらの方法の 1 つを適切に実装するには、時間ステップの式以上のものが必要になります。
常に同じステップサイズを用いるのは非効率な場合が多いため、可変ステップサイズ法 が開発されています。通常、ステップサイズは、ステップあたりの(局所的な)誤差が許容レベルを下回るように選択されます。つまり、これらの手法では、局所的な誤差の推定値 である誤差指標も計算する必要があります。
この考え方の拡張として、異なる順序の異なる手法を動的に選択する方法があります(これは可変順序法と呼ばれます)。 リチャードソン外挿法 [ 14 ] に基づく手法、例えばブリルシュ・ストアーアルゴリズム [ 15 ] [ 16 ] は、異なる順序の様々な手法を構築するためによく用いられます。
その他の望ましい機能は次のとおりです。
稠密な出力: t 0 、t 1 、t 2 、 ...の点だけでなく、積分区間全体に対する簡便な数値近似値。イベントの位置 :例えば、特定の関数が消滅する時刻を見つけること。これには通常、根を求めるアルゴリズム の使用が必要となる。並列コンピューティング のサポート。時間に関して積分する場合、時間可逆性
代替方法 ここで議論されている枠組みに当てはまらない方法も数多くあります。代替方法には以下のようなものがあります。
多重微分法は 、関数f だけでなくその導関数も用いる。このクラスには、エルミート・オブレシュコフ法 やフェルバーグ法の ほか、パーカー・ソチャッキ法 [ 17 ] やビシュコフ・シェルバコフ法といった、解yの テイラー級数 の係数を再帰的に計算する手法が含まれる。二階常微分方程式に対するニストローム法 。すべての高階常微分方程式は(1)の形の一階常微分方程式に変換できると述べました。これは確かに正しいのですが、必ずしも最善の方法とは限りません。特に、ニストローム法は 二階方程式に直接適用できます。幾何学的積分法 [ 18 ] [ 19 ] は、特別なクラスの常微分方程式(例えば、ハミルトン方程式 の解を求めるシンプレクティック積分法 )のために特別に設計されている。これらの積分法は、数値解がこれらのクラスの基礎となる構造や幾何学を尊重するように配慮している。量子化状態システム法は 、状態量子化の考え方に基づく常微分方程式の積分法の一種です。頻繁に不連続点を持つ疎なシステムをシミュレーションする際に効果的です。
時間並列法 一部のIVPでは、非常に高い時間分解能や非常に長い時間間隔での積分が必要となるため、従来のシリアルタイムステップ法ではリアルタイム実行が計算上不可能となります(数値気象予報、プラズマモデリング、分子動力学におけるIVPの例)。これらの問題に対処するため、並列計算を用いてシミュレーション実行時間を短縮する Parallel-in-Time (PinT)法が開発されました。
初期のPinT法(最も古いものは1960年代に提案された)[ 20 ] は、必要とする並列コンピューティングアーキテクチャがまだ広く利用できていなかったため、研究者に当初は見過ごされていました。利用可能なコンピューティングパワーが増えるにつれて、2000年代初頭にParareal が開発されたことで関心が再燃しました。PinTアルゴリズムは、さまざまなIVPを解くのに適しています。エクサスケールコンピューティングの出現により、PinTアルゴリズムはますます研究の注目を集めており、世界で最も強力な スーパーコンピュータ を利用できるように開発されています。2023年現在で最も人気のある方法には、Parareal、PFASST、ParaDiag、MGRITなどがあります。[ 21 ]
分析 数値解析 とは、数値手法の設計だけでなく、その解析も意味します。この解析における中心となる概念は以下の3つです。
収束性 :方法が解を近似するかどうか次数 :解をどれだけ正確に近似しているか、そして安定性 :誤差が抑制されるかどうか。 [ 22 ]
収束 数値解がステップサイズhが0に近づくにつれて正確な解に近づくとき、数値法は 収束する と言われる。より正確には、リプシッツ 関数f とt * > 0を持つすべての常微分方程式(1)に対して 、
lim h → 0 + max 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.}
上記の方法はすべて収束的です。
一貫性と秩序 数値計算法が
y n + k = Ψ ( t n + k ; y n , y n + 1 , … , y n + k − 1 ; h ) . {\displaystyle y_{n+k}=\Psi (t_{n+k};y_{n},y_{n+1},\dots ,y_{n+k-1};h).\,}
手法の局所的(打ち切り)誤差 とは、手法の1つのステップで発生する誤差のことです。つまり、それ以前のステップで誤差が発生しなかったと仮定した場合の手法によって得られる結果と、正確な解との差です。
δ n + k h = Ψ ( t n + k ; y ( t n ) , y ( t n + 1 ) , … , y ( t n + k − 1 ) ; h ) − y ( t n + k ) . {\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}).}
この方法は、次の 場合、 無矛盾 であると言える。 この方法は、次の場合、次数 を持つ。したがって、次数が0より大きい場合、その方法は無矛盾である。上で紹介した(前進)オイラー法(4)と後退オイラー法(6)はどちらも次数が1であるため、無矛盾である。実際に使用されているほとんどの方法は、より高い次数を達成している。無矛盾は収束の必要条件だが、十分条件ではない。方法が収束するためには、無矛盾かつ零点安定性の 両方を満たしていなければならない。 lim h → 0 δ n + k h h = 0. {\displaystyle \lim _{h\to 0}{\frac {\delta _{n+k}^{h}}{h}}=0.} p {\displaystyle p} δ n + k h = O ( h p + 1 ) as h → 0. {\displaystyle \delta _{n+k}^{h}=O(h^{p+1})\quad {\text{as }}h\to 0.}
関連する概念として、大域誤差(打ち切り誤差) があります。これは、固定された時刻 に到達するために必要なすべてのステップで発生する誤差です。明示的には、時刻 における大域誤差は となります。階のワンステップ法の大域誤差は です。特に、このような方法は収束します。ただし、この記述はマルチステップ法では必ずしも当てはまりません。 t {\displaystyle t} t {\displaystyle t} y N − y ( t ) {\displaystyle y_{N}-y(t)} N = ( t − t 0 ) / h {\displaystyle N=(t-t_{0})/h} p {\displaystyle p} O ( h p ) {\displaystyle O(h^{p})}
安定性と剛性 一部の微分方程式では、オイラー法、陽的ルンゲ・クッタ法 、多段階法 (例えばアダムス・バッシュフォース法)などの標準的な手法を適用すると、解が不安定になるが、他の手法では安定した解が得られる場合がある。方程式におけるこの「困難な挙動」(それ自体は必ずしも複雑ではない)は剛性 として記述され、多くの場合、基礎にある問題に異なる時間スケールが存在することによって引き起こされる。[ 23 ] 例えば、衝突振動子 のような機械システムにおける衝突は、通常、物体の運動の時間よりもはるかに短い時間スケールで発生する。この不一致により、状態パラメータの曲線に非常に「急激な変化」が生じる。
硬い問題は、化学反応速度論 、制御理論 、固体力学 、気象予報 、生物学 、プラズマ物理学 、電子工学 など、あらゆる分野で広く見られる。この硬直性を克服する一つの方法は、微分方程式の概念を微分包含 の概念に拡張することである。これは、非滑らかさを許容し、モデル化する。[ 24 ] [ 25 ]
歴史 以下はこの分野におけるいくつかの重要な進展のタイムラインです。 [ 26 ] [ 27 ]
2次1次元境界値問題の数値解 境界値問題(BVP)は通常、元のBVPを離散化することで得られる近似的に等価な行列問題を解くことで数値的に解かれる。[ 28 ] 1次元のBVPを数値的に解くために最も一般的に用いられる方法は、有限差分法 と呼ばれる。[ 3 ] この方法は、点の値の線形結合を利用して、関数の導関数を記述する有限差分係数 を構築する。例えば、1次導関数の2次中心差分 近似は次のように表される。
u ′ ( x i ) = u i + 1 − u i − 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 i ) = u i + 1 − 2 u i + u i − 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}).}
これらの式の両方において、は離散化領域上の隣接するx 値間の距離です。次に、標準的な行列法 で解ける線形連立方程式を構築します。例えば、解くべき方程式が次の式であるとします。 h = x i − x i − 1 {\displaystyle h=x_{i}-x_{i-1}}
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 i ″ = u i + 1 − 2 u i + u i − 1 h 2 {\displaystyle u''_{i}={\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}}
得られた連立一次方程式を解きます。すると、次のような方程式が得られます。
u i + 1 − 2 u i + u i − 1 h 2 − u i = 0 , ∀ i = 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 には境界値とを含む項があり、これら2つの値は既知であるため、この方程式に代入するだけで、結果として非自明な解を持つ非同次線形連立方程式が 得られます。 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ページ)^ a b LeVeque, RJ (2007). 常微分方程式と偏微分方程式の有限差分法:定常問題と時間依存問題(第98巻). SIAM. ^ Slimane AdjeridとMahboub Baccouch(2010)「ガラーキン法」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巻, pp. 185–202). ケンブリッジ大学出版局. ^ Weisstein, Eric W. 「ガウス求積法」 MathWorld(Wolfram Webリソース)より。https ://mathworld.wolfram.com/GaussianQuadrature.html ^ Butcher, JC (1987). 常微分方程式の数値解析:ルンゲ・クッタ法と一般線形法. Wiley-Interscience. ^ Hochbruck & Ostermann (2010 , pp. 209–286) これは指数積分器に関する最新の広範なレビュー論文である。^ Brezinski, C., Zaglia, MR (2013). 外挿法:理論と実践. エルゼビア. ^ 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). 数値アルゴリズムの精度と安定性 (第80巻). SIAM. ^ Miranker, A. (2001). スティフ方程式と特異摂動問題のための数値解析法:特異摂動問題(第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世紀における歴史的発展. エルゼビア. ^ Butcher, JC (1996). ルンゲ・クッタ法の歴史. 応用数値数学, 20(3), 247-260. ^ Ascher, UM, Mattheij, RM, & Russell, RD (1995). 常微分方程式の境界値問題の数値解法. Society for Industrial and Applied Mathematics.
参考文献
外部リンク