Approach to finding numerical solutions of ordinary differential equations
(図1) オイラー法の図解。未知の曲線は青で、その多角形近似は赤で示されています。 数学 および 計算科学 において 、 オイラー法 （ 前進オイラー法 とも呼ばれる）は、与えられた 初期値を持つ 常微分方程式 （ODE）を解くための 一次数値 計算手法である。これは 常微分方程式の数値積分 における 最も基本的な 明示的手法 であり、最も単純な ルンゲ・クッタ法である。オイラー法は、 レオンハルト・オイラー が著書 『積分計算の原理 』 （1768～1770年出版）で初めて提唱したこと にちなんで名付けられた。 [1]
オイラー法は一次近似法であり、局所誤差（ステップごとの誤差）はステップサイズの2乗に比例し、大域誤差（特定の時点における誤差）はステップサイズに比例します。オイラー法は、 予測子修正子法 など、より複雑な手法を構築するための基礎としてよく用いられます。
幾何学的記述
目的と機能の理由 与えられた点から始まり、与えられた微分方程式を満たす未知の曲線の形状を計算する問題を考えてみましょう。ここで、微分方程式とは、曲線上の任意の点の位置が計算されれば、その点における曲線の 接線 の 傾き を計算できる式と考えることができます
考え方としては、曲線は最初は未知ですが、その始点 （図1参照）は既知です。そして、微分方程式から、曲線の傾き（ ）を 計算し、接線を求めることができます。 A 0 , {\displaystyle A_{0},} A 0 {\displaystyle A_{0}}
その接線に沿って点まで小さなステップを踏みます。 この小さなステップでは傾きはあまり変化しないため、 曲線に近くなります。が まだ曲線上にあると仮定すると、 上記の点の場合と同じ推論が適用できます。いくつかのステップを踏むと、 折れ線曲線 （ ）が計算されます。一般に、この曲線は元の未知の曲線から大きく逸脱することはなく、ステップサイズが十分に小さく、計算間隔が有限であれば、2つの曲線間の誤差は小さく抑えることができます。 [2] A 1 . {\displaystyle A_{1}.} A 1 {\displaystyle A_{1}} A 1 {\displaystyle A_{1}} A 0 {\displaystyle A_{0}} A 0 , A 1 , A 2 , A 3 , … {\displaystyle A_{0},A_{1},A_{2},A_{3},\dots }
一次過程 と の値が与えられている場合 、 の導関数 は と の与えられた関数であり 、 と 表記されます 。 を設定することで処理を開始します。次に、 t軸に沿ったすべてのステップのサイズの 値を選択し、 （または同等の ）を設定します。ここで、オイラー法を使用して と から を求めます 。 [ 3 ] t 0 {\displaystyle t_{0}} y ( t 0 ) {\displaystyle y(t_{0})} y {\displaystyle y} t {\displaystyle t} y {\displaystyle y} y ′ ( t ) = f ( t , y ( t ) ) {\displaystyle y'(t)=f\left(t,y(t)\right)} y 0 = y ( t 0 ) {\displaystyle y_{0}=y(t_{0})} h {\displaystyle h} t n = t 0 + n h {\displaystyle t_{n}=t_{0}+nh} t n + 1 = t n + h {\displaystyle t_{n+1}=t_{n}+h} y n + 1 {\displaystyle y_{n+1}} y n {\displaystyle y_{n}} t n {\displaystyle t_{n}}
y n + 1 = y n + h f ( t n , y n ) . {\displaystyle y_{n+1}=y_{n}+hf(t_{n},y_{n}).}
の値は、 時刻 における解の近似値 、すなわち です 。オイラー法は 明示的 です。すなわち、 の解は に対して の明示的な関数です 。 y n {\displaystyle y_{n}} t n {\displaystyle t_{n}} y n ≈ y ( t n ) {\displaystyle y_{n}\approx y(t_{n})} y n + 1 {\displaystyle y_{n+1}} y i {\displaystyle y_{i}} i ≤ n {\displaystyle i\leq n}
高階過程 オイラー法は1階常微分方程式を積分しますが、任意の階数の常微分方程式は1階常微分方程式の系として表すことができます。 次のように定義される 階数の常微分方程式が与えられたとき N {\displaystyle N} N {\displaystyle N}
y ( N + 1 ) ( t ) = f ( t , y ( t ) , y ′ ( t ) , … , y ( N ) ( t ) ) , {\displaystyle y^{(N+1)}(t)=f\left(t,y(t),y'(t),\ldots ,y^{(N)}(t)\right),}
、、 およびと 同様に 、希望する時刻における常微分方程式の解の近似値に達するまで、次の式を実行します。 h {\displaystyle h} t 0 {\displaystyle t_{0}} y 0 , y 0 ′ , … , y 0 ( N ) {\displaystyle y_{0},y'_{0},\dots ,y_{0}^{(N)}}
y → i + 1 = ( y i + 1 y i + 1 ′ ⋮ y i + 1 ( N − 1 ) y i + 1 ( N ) ) = ( y i + h ⋅ y i ′ y i ′ + h ⋅ y i ″ ⋮ y i ( N − 1 ) + h ⋅ y i ( N ) y i ( N ) + h ⋅ f ( t i , y i , y i ′ , … , y i ( N ) ) ) {\displaystyle {\vec {y}}_{i+1}={\begin{pmatrix}y_{i+1}\\y'_{i+1}\\\vdots \\y_{i+1}^{(N-1)}\\y_{i+1}^{(N)}\end{pmatrix}}={\begin{pmatrix}y_{i}+h\cdot y'_{i}\\y'_{i}+h\cdot y''_{i}\\\vdots \\y_{i}^{(N-1)}+h\cdot y_{i}^{(N)}\\y_{i}^{(N)}+h\cdot f\left(t_{i},y_{i},y'_{i},\ldots ,y_{i}^{(N)}\right)\end{pmatrix}}}
これらの一次システムはオイラー法、あるいは実際には一次システムのための他の任意の手法で扱うことができる。 [4]
一次例 初期値問題が与えられた場合
y ′ = y , y ( 0 ) = 1 , {\displaystyle y'=y,\quad y(0)=1,}
オイラー法を用いて近似したい 。 [5] y ( 4 ) {\displaystyle y(4)}
ステップサイズを1（ h = 1 ） （図2）方程式の数値積分の図解 青はオイラー法、緑は 中点法 、赤は厳密解、 ステップサイズは y ′ = y , y ( 0 ) = 1. {\displaystyle y'=y,y(0)=1.} y = e t . {\displaystyle y=e^{t}.} h = 1.0. {\displaystyle h=1.0.} オイラー法は
y n + 1 = y n + h f ( t n , y n ) . {\displaystyle y_{n+1}=y_{n}+hf(t_{n},y_{n}).}
まず を計算しなければなりません 。この単純な微分方程式では、関数 は で定義されます 。 f ( t 0 , y 0 ) {\displaystyle f(t_{0},y_{0})} f {\displaystyle f} f ( t , y ) = y {\displaystyle f(t,y)=y}
f ( t 0 , y 0 ) = f ( 0 , 1 ) = 1. {\displaystyle f(t_{0},y_{0})=f(0,1)=1.}
上記の手順により、点 における解曲線に接する直線の傾きを求めました 。傾きは の変化を の変化で割ったもの 、つまり として定義されることを思い出してください 。 ( 0 , 1 ) {\displaystyle (0,1)} y {\displaystyle y} t {\displaystyle t} Δ y Δ t {\textstyle {\frac {\Delta y}{\Delta t}}}
次のステップでは、上記の値にステップ サイズを掛けます 。ここではステップ サイズは 1 とします。 h {\displaystyle h}
h ⋅ f ( y 0 ) = 1 ⋅ 1 = 1. {\displaystyle h\cdot f(y_{0})=1\cdot 1=1.}
ステップサイズは の変化なので 、ステップサイズと接線の傾きを掛け合わせると、 値の変化が得られます。この値は初期 値に加算され、計算に使用する次の値が得られます。 t {\displaystyle t} y {\displaystyle y} y {\displaystyle y}
y 0 + h f ( y 0 ) = y 1 = 1 + 1 ⋅ 1 = 2. {\displaystyle y_{0}+hf(y_{0})=y_{1}=1+1\cdot 1=2.}
上記の手順を繰り返して 、、 およびを見つけます 。 y 2 {\displaystyle y_{2}} y 3 {\displaystyle y_{3}} y 4 {\displaystyle y_{4}}
y 2 = y 1 + h f ( y 1 ) = 2 + 1 ⋅ 2 = 4 , y 3 = y 2 + h f ( y 2 ) = 4 + 1 ⋅ 4 = 8 , y 4 = y 3 + h f ( y 3 ) = 8 + 1 ⋅ 8 = 16. {\displaystyle {\begin{aligned}y_{2}&=y_{1}+hf(y_{1})=2+1\cdot 2=4,\\y_{3}&=y_{2}+hf(y_{2})=4+1\cdot 4=8,\\y_{4}&=y_{3}+hf(y_{3})=8+1\cdot 8=16.\end{aligned}}}
このアルゴリズムは反復的な性質を持つため、エラーを避けるために、以下に示すように計算をチャート形式で整理すると役立ちます。
n {\displaystyle n} y n {\displaystyle y_{n}} t n {\displaystyle t_{n}} f ( t n , y n ) {\displaystyle f(t_{n},y_{n})} h {\displaystyle h} Δ y {\displaystyle \Delta y} y n + 1 {\displaystyle y_{n+1}} 0 1 0 1 1 1 2 1 2 1 2 1 2 4 2 4 2 4 1 4 8 3 8 3 8 1 8 16
この計算の結果は です 。微分方程式の正確な解は なので 、 となります 。オイラー法の近似値は、特に大きな値のステップサイズ のために、この特定のケースではそれほど正確ではありませんでしたが 、図が示すように、その動作は定性的に正しいです。 y 4 = 16 {\displaystyle y_{4}=16} y ( t ) = e t {\displaystyle y(t)=e^{t}} y ( 4 ) = e 4 ≈ 54.598 {\displaystyle y(4)=e^{4}\approx 54.598} h {\displaystyle h}
他のステップサイズの使用 （図3）同じ図 h = 0.25. {\displaystyle h=0.25.} 導入部で述べたように、オイラー法はステップサイズ が小さいほど精度が高くなります。下の表は、ステップサイズを変えた場合の結果を示しています。一番上の行は前のセクションの例に対応し、2行目は図に示されています。 h {\displaystyle h}
ステップサイズ オイラー法の結果 誤差 1 16.00 38.60 0.25 35.53 19.07 0.1 45.26 0 9.34 0.05 49.56 0 5.04 0.025 51.98 0 2.62 0.0125 53.26 0 1.34
表の最後の列に記録されている誤差は、における正確な解 とオイラー近似との差です。表の下部では、ステップサイズは前の行のステップサイズの半分であり、誤差も前の行の誤差の約半分です。これは、少なくともステップサイズがかなり小さい場合、誤差がステップサイズにほぼ比例することを示唆しています。これは一般に当てはまり、他の方程式にも当てはまります。詳細については、「全体的打ち切り誤差」のセクションを参照してください t = 4 {\displaystyle t=4}
図にも示されている 中点法 などの他の手法は、より好ましい挙動を示します。中点法の全体誤差は、ステップサイズの 2乗 にほぼ比例します。このため、オイラー法は一次法、中点法は二次法と呼ばれます。
上記の表から、小数点以下3桁まで正確な答えを得るために必要なステップサイズは約0.00001であることが推測できます。これは40万ステップ必要であることを意味します。この大きなステップ数は計算コストを増大させます。そのため、特に高い精度が求められる場合は、 ルンゲ・クッタ法 や 線形多段階法 などの高次手法が用いられます。 [6]
高次の例 この三次の例では、以下の情報が与えられていると仮定します
y ‴ + 4 t y ″ − t 2 y ′ − ( cos t ) y = sin t t 0 = 0 y 0 = y ( t 0 ) = 2 y 0 ′ = y ′ ( t 0 ) = − 1 y 0 ″ = y ″ ( t 0 ) = 3 h = 0.5 {\displaystyle {\begin{aligned}&y'''+4ty''-t^{2}y'-(\cos {t})y=\sin {t}\\&t_{0}=0\\&y_{0}=y(t_{0})=2\\&y'_{0}=y'(t_{0})=-1\\&y''_{0}=y''(t_{0})=3\\&h=0.5\end{aligned}}}
これから y ''' を分離して次の方程式を得ることができます。
y ‴ = f ( t , y , y ′ , y ″ ) = sin t + ( cos t ) y + t 2 y ′ − 4 t y ″ {\displaystyle y'''=f{\left(t,y,y',y''\right)}=\sin {t}+\left(\cos {t}\right)y+t^{2}y'-4ty''}
これを使用すると、 の解を得ることができます 。 また、 の解を使用すると 、 の解を得ることができます 。 必要なだけ同じ式を使用してこのプロセスを続けることで、希望するものを見つけることができます 。 y → 1 {\displaystyle {\vec {y}}_{1}} y → 1 = ( y 1 y 1 ′ y 1 ″ ) = ( y 0 y 0 ′ y 0 ″ ) + h ( y 0 ′ y 0 ″ f ( t 0 , y 0 , y 0 ′ , y 0 ″ ) ) = ( 2 − 1 3 ) + 0.5 ⋅ ( − 1 3 sin 0 + ( cos 0 ) ⋅ 2 + 0 2 ⋅ ( − 1 ) − 4 ⋅ 0 ⋅ 3 ) = ( 1.5 0.5 4 ) {\displaystyle {\begin{aligned}{\vec {y}}_{1}&={\begin{pmatrix}y_{1}\\y_{1}'\\y_{1}''\end{pmatrix}}={\begin{pmatrix}y_{0}\\y'_{0}\\y''_{0}\end{pmatrix}}+h{\begin{pmatrix}y'_{0}\\y''_{0}\\f{\left(t_{0},y_{0},y'_{0},y''_{0}\right)}\end{pmatrix}}\\[1ex]&={\begin{pmatrix}2\\-1\\3\end{pmatrix}}+0.5\cdot {\begin{pmatrix}-1\\3\\\sin {0}+(\cos {0})\cdot 2+0^{2}\cdot (-1)-4\cdot 0\cdot 3\end{pmatrix}}\\[1ex]&={\begin{pmatrix}1.5\\0.5\\4\end{pmatrix}}\end{aligned}}} y → 1 {\displaystyle {\vec {y}}_{1}} y → 2 {\displaystyle {\vec {y}}_{2}} y → 2 = ( y 2 y 2 ′ y 2 ″ ) = ( y 1 y 1 ′ y 1 ″ ) + h ( y 1 ′ y 1 ″ f ( t 1 , y 1 , y 1 ′ , y 1 ″ ) ) = ( 1.5 0.5 4 ) + 0.5 ⋅ ( 0.5 4 sin 0.5 + ( cos 0.5 ) ⋅ 1.5 + 0.5 2 ⋅ 0.5 − 4 ⋅ 0.5 ⋅ 4 ) = ( 1.75 2.5 0.9604... ) {\displaystyle {\begin{aligned}{\vec {y}}_{2}&={\begin{pmatrix}y_{2}\\y_{2}'\\y_{2}''\end{pmatrix}}={\begin{pmatrix}y_{1}\\y'_{1}\\y''_{1}\end{pmatrix}}+h{\begin{pmatrix}y'_{1}\\y''_{1}\\f{\left(t_{1},y_{1},y'_{1},y''_{1}\right)}\end{pmatrix}}\\[1ex]&={\begin{pmatrix}1.5\\0.5\\4\end{pmatrix}}+0.5\cdot {\begin{pmatrix}0.5\\4\\\sin {0.5}+(\cos {0.5})\cdot 1.5+0.5^{2}\cdot 0.5-4\cdot 0.5\cdot 4\end{pmatrix}}\\[1ex]&={\begin{pmatrix}1.75\\2.5\\0.9604...\end{pmatrix}}\end{aligned}}} y → i {\displaystyle {\vec {y}}_{i}}
導出 オイラー法はいくつかの方法で導出できます。
まず、上記の幾何学的記述があります もう一つの可能性は、 の周りの 関数の テイラー展開 を考えることです 。 y {\displaystyle y} t 0 {\displaystyle t_{0}} y ( t 0 + h ) = y ( t 0 ) + h y ′ ( t 0 ) + 1 2 h 2 y ″ ( t 0 ) + O ( h 3 ) . {\displaystyle y(t_{0}+h)=y(t_{0})+hy'(t_{0})+{\tfrac {1}{2}}h^{2}y''(t_{0})+O\left(h^{3}\right).}
微分方程式は と述べている 。これをテイラー展開に代入し、二次項と高次項を無視すると、オイラー法が得られる。 [7] y ′ = f ( t , y ) {\displaystyle y'=f(t,y)}
以下では、テイラー展開を使用してオイラー法で発生する誤差を解析しますが、これを拡張して ルンゲ・クッタ法 を作成することもできます。 これに密接に関連する導出法として 、微分方程式の導関数を
前向き 差分 公式に代入する方法が あります。これもまた、オイラー法となります。 [8] 同様の計算から、 中点法 と 後向きオイラー法 が得られます。 y ′ ( t 0 ) ≈ y ( t 0 + h ) − y ( t 0 ) h {\displaystyle y'(t_{0})\approx {\frac {y(t_{0}+h)-y(t_{0})}{h}}} y ′ = f ( t , y ) {\displaystyle y'=f(t,y)} 最後に、微分方程式を から まで積分し 、 微積分の基本定理を 適用すると次の式 が得られます。 ここで、積分を左辺の 長方形法 （長方形は1つだけ）
で近似します 。両方の式を組み合わせると、再びオイラー法が見つかります。 [9] t 0 {\displaystyle t_{0}} t 0 + h {\displaystyle t_{0}+h} y ( t 0 + h ) − y ( t 0 ) = ∫ t 0 t 0 + h f ( t , y ( t ) ) d t . {\displaystyle y(t_{0}+h)-y(t_{0})=\int _{t_{0}}^{t_{0}+h}f\left(t,y(t)\right)\,\mathrm {d} t.} ∫ t 0 t 0 + h f ( t , y ( t ) ) d t ≈ h f ( t 0 , y ( t 0 ) ) . {\displaystyle \int _{t_{0}}^{t_{0}+h}f\left(t,y(t)\right)\,\mathrm {d} t\approx hf\left(t_{0},y(t_{0})\right).} この考え方を継続することで、さまざまな 線形多段階法 に到達できます。
局所的打ち切り誤差 オイラー法の局所 的打ち切り誤差 は、1ステップで生じる誤差です。これは、1ステップ後の数値解 と時刻における正確な解との差です 。数値解は次のように与えられます y 1 {\displaystyle y_{1}} t 1 = t 0 + h {\displaystyle t_{1}=t_{0}+h}
y 1 = y 0 + h f ( t 0 , y 0 ) . {\displaystyle y_{1}=y_{0}+hf(t_{0},y_{0}).}
正確な解を得るには、上記の導出セクションで説明したテイラー展開を使用します。
y ( t 0 + h ) = y ( t 0 ) + h y ′ ( t 0 ) + 1 2 h 2 y ″ ( t 0 ) + O ( h 3 ) . {\displaystyle y(t_{0}+h)=y(t_{0})+hy'(t_{0})+{\tfrac {1}{2}}h^{2}y''(t_{0})+{\mathcal {O}}{\left(h^{3}\right)}.}
オイラー法によって導入される局所的切り捨て誤差 (LTE) は、次の方程式の差で表されます。
L T E = y ( t 0 + h ) − y 1 = 1 2 h 2 y ″ ( t 0 ) + O ( h 3 ) . {\displaystyle \mathrm {LTE} =y(t_{0}+h)-y_{1}={\tfrac {1}{2}}h^{2}y''(t_{0})+{\mathcal {O}}{\left(h^{3}\right)}.}
この結果は、 3次導関数が有界である場合に有効である。 [10] y {\displaystyle y}
これは、 が小さい場合 、局所的な打ち切り誤差が にほぼ比例することを示しています 。そのため、オイラー法は、局所的な打ち切り誤差がステップサイズのより高いべき乗に比例する ルンゲ・クッタ法 や 線形マルチステップ法 などの高次手法よりも精度が低くなります。 h {\displaystyle h} h 2 {\displaystyle h^{2}}
局所的打ち切り誤差の若干異なる定式化は、テイラーの定理 の剰余項にラグランジュ形式を用いることで得られる 。が連続な2次微分を持つ場合、次 式を満たすような が 存在する [11]。 y {\displaystyle y} ξ ∈ [ t 0 , t 0 + h ] {\displaystyle \xi \in [t_{0},t_{0}+h]}
L T E = y ( t 0 + h ) − y 1 = 1 2 h 2 y ″ ( ξ ) . {\displaystyle \mathrm {LTE} =y(t_{0}+h)-y_{1}={\tfrac {1}{2}}h^{2}y''(\xi ).}
上記の誤差の式において、未知の厳密解の2階微分は、 微分方程式の右辺を含む式に置き換えることができる。実際、この式から [ 12]が導かれる。 y {\displaystyle y} y ′ = f ( t , y ) {\displaystyle y'=f(t,y)}
y ″ ( t 0 ) = ∂ f ∂ t ( t 0 , y ( t 0 ) ) + ∂ f ∂ y ( t 0 , y ( t 0 ) ) f ( t 0 , y ( t 0 ) ) . {\displaystyle y''(t_{0})={\frac {\partial f}{\partial t}}\left(t_{0},y(t_{0})\right)+{\frac {\partial f}{\partial y}}\left(t_{0},y(t_{0})\right)\,f\left(t_{0},y(t_{0})\right).}
全体的切り捨て誤差 全体 的切り捨て誤差 とは、初期時刻からその時刻に到達するまでに何ステップ必要になったかに関わらず、 固定時刻における誤差です。全体的切り捨て誤差は、各ステップで発生する局所的切り捨て誤差の累積的な影響です。 [13] ステップ数は と簡単に決定でき 、これは に比例し 、各ステップで発生する誤差は に比例します （前のセクションを参照）。したがって、全体的切り捨て誤差は に比例すると予想されます 。 [14] t i {\displaystyle t_{i}} t i − t 0 h {\textstyle {\frac {t_{i}-t_{0}}{h}}} 1 h {\textstyle {\frac {1}{h}}} h 2 {\displaystyle h^{2}} h {\displaystyle h}
この直感的な推論は明確に説明できる。解が 有界な2次導関数を持ち、その2番目の引数において リプシッツ連続で ある場合 、大域的打ち切り誤差（ と表記 ）は y {\displaystyle y} f {\displaystyle f} | y ( t i ) − y i | {\displaystyle |y(t_{i})-y_{i}|}
| y ( t i ) − y i | ≤ h M 2 L ( e L ( t i − t 0 ) − 1 ) {\displaystyle \left|y(t_{i})-y_{i}\right|\leq {\frac {hM}{2L}}\left(e^{L(t_{i}-t_{0})}-1\right)}
ここで、は 与えられた区間における の二階微分の上界であり、 は のリプシッツ定数である 。 [15] あるいはもっと簡単に言えば、 のとき 、値 （ は 定数として扱われる）である。対照的に、 の関数は変数 のみを含む正確な解である 。 M {\displaystyle M} y {\displaystyle y} L {\displaystyle L} f {\displaystyle f} y ′ ( t ) = f ( t , y ) {\displaystyle y'(t)=f(t,y)} L = max ( | d d y [ f ( t , y ) ] | ) {\textstyle L={\text{max}}\left(|{\frac {d}{dy}}\left[f(t,y)\right]|\right)} t {\displaystyle t} M = max ( | d 2 d t 2 [ y ( t ) ] | ) {\textstyle M=\max \left(\left|{\frac {d^{2}}{dt^{2}}}\left[y(t)\right]\right|\right)} y ( t ) {\displaystyle y(t)} t {\displaystyle t}
この境界値の正確な形は実用上あまり重要ではない。なぜなら、ほとんどの場合、この境界値はオイラー法によって生じる実際の誤差を大幅に過大評価するからである。 [16] 重要なのは、この境界値が、大域的打ち切り誤差が（近似的に）に比例することを示していることである 。このため、オイラー法は一次法であると言われている。 [17] h {\displaystyle h}
例 微分方程式と厳密解 が あり、 のときの と を 求めたいとします 。
したがって、 t =2.5、 h =0.5 における誤差境界を求めることができます y ′ = 1 + ( t − y ) 2 {\displaystyle y'=1+(t-y)^{2}} y = t + 1 t − 1 {\displaystyle y=t+{\frac {1}{t-1}}} M {\displaystyle M} L {\displaystyle L} 2 ≤ t ≤ 3 {\displaystyle 2\leq t\leq 3} L = max | d d y f ( t , y ) | = max 2 ≤ t ≤ 3 | d d y [ 1 + ( t − y ) 2 ] | = max 2 ≤ t ≤ 3 | 2 ( t − y ) | = max 2 ≤ t ≤ 3 | 2 ( t − [ t + 1 t − 1 ] ) | = max 2 ≤ t ≤ 3 | − 2 t − 1 | = 2 {\displaystyle {\begin{aligned}L&=\max \left|{\frac {d}{dy}}f(t,y)\right|=\max _{2\leq t\leq 3}\left|{\frac {d}{dy}}\left[1+\left(t-y\right)^{2}\right]\right|\\[1ex]&=\max _{2\leq t\leq 3}\left|2\left(t-y\right)\right|=\max _{2\leq t\leq 3}\left|2\left(t-\left[t+{\frac {1}{t-1}}\right]\right)\right|\\[1ex]&=\max _{2\leq t\leq 3}\left|-{\frac {2}{t-1}}\right|=2\end{aligned}}} M = max | d 2 d t 2 [ y ( t ) ] | = max 2 ≤ t ≤ 3 | d 2 d t 2 ( t + 1 1 − t ) | = max 2 ≤ t ≤ 3 | 2 ( − t + 1 ) 3 | = 2 {\displaystyle {\begin{aligned}M&=\max \left|{\frac {d^{2}}{dt^{2}}}\left[y(t)\right]\right|\\&=\max _{2\leq t\leq 3}\left|{\frac {d^{2}}{dt^{2}}}\left(t+{\frac {1}{1-t}}\right)\right|\\&=\max _{2\leq t\leq 3}\left|{\frac {2}{\left(-t+1\right)^{3}}}\right|=2\end{aligned}}}
error bound = h M 2 L ( e L ( t i − t 0 ) − 1 ) = 0.5 ⋅ 2 2 ⋅ 2 ( e 2 ( 2.5 − 2 ) − 1 ) = 0.42957 {\displaystyle {\begin{aligned}{\text{error bound}}&={\frac {hM}{2L}}\left(e^{L(t_{i}-t_{0})}-1\right)\\[1ex]&={\frac {0.5\cdot 2}{2\cdot 2}}\left(e^{2(2.5-2)}-1\right)=0.42957\end{aligned}}} t 0 は 2 に等しい
ことに注意してください 。これは、 における t の下限値だからです 。 2 ≤ t ≤ 3 {\displaystyle 2\leq t\leq 3}
数値安定性 （図4）ステップサイズ（青い四角）と（赤い丸） を用いてオイラー法で計算された 解 。黒い曲線は厳密解を示しています y ′ = − 2.3 y {\displaystyle y'=-2.3y} h = 1 {\displaystyle h=1} h = 0.7 {\displaystyle h=0.7} オイラー法は数値的に 不安定に なる場合があり、特に 硬い方程式 の場合に不安定になります。つまり、正確な解が ではない方程式の数値解が非常に大きくなります。これは、線形方程式 を使用して説明できます。 正確な解は であり 、 のときにゼロに減衰します 。ただし、ステップ サイズ でこの方程式にオイラー法を適用すると 、数値解は質的に間違っています。つまり、振動して増大します (図を参照)。これが不安定であるということです。より小さなステップ サイズ (たとえば ) を使用すると、 数値解はゼロに減衰します。 y ′ = − 2.3 y , y ( 0 ) = 1. {\displaystyle y'=-2.3y,\qquad y(0)=1.} y ( t ) = e − 2.3 t {\displaystyle y(t)=e^{-2.3t}} t → ∞ {\displaystyle t\to \infty } h = 1 {\displaystyle h=1} h = 0.7 {\displaystyle h=0.7}
(図 5) ピンク色の円盤はオイラー法の安定領域を示しています。 線形方程式 にオイラー法を適用すると、積が 右図の 領域外にある 場合、数値解は不安定になります。この領域は（線形） 安定領域 と呼ばれます。 [18] この例では、 と なるので、 となると となり、 は安定領域外にあるため、数値解は不安定になります。 y ′ = k y {\displaystyle y'=ky} h k {\displaystyle hk} { z ∈ C | | z + 1 | ≤ 1 } , {\displaystyle \left\{z\in \mathbf {C} \,{\big |}\,|z+1|\leq 1\right\},} k = − 2.3 {\displaystyle k=-2.3} h = 1 {\displaystyle h=1} h k = − 2.3 {\displaystyle hk=-2.3}
この制限と、誤差の収束が遅いことから、 オイラー法は数値積分の単純な例を除いて、あまり用いられません [ 要出典 ] 。物理システムのモデルには、減衰の速い要素（つまり、大きな負の指数関数を持つ要素）を表す項が含まれることがよくあります。これらの要素が全体の解において重要でない場合でも、それらが引き起こす不安定性により、オイラー法を用いる場合は非常に小さなタイムステップが必要となる場合があります。 h {\displaystyle h}
丸め誤差 オイラー法の ステップでは、 丸め誤差 はおおよそ の大きさです。 ここで、は マシンイプシロン です 。丸め誤差が独立した確率変数であると仮定すると、期待される合計丸め誤差は に比例します 。 [19] したがって、ステップサイズが非常に小さい値の場合、切り捨て誤差は小さくなりますが、丸め誤差の影響は大きくなる可能性があります。 オイラー法の式で 補正和を使用すれば、丸め誤差の影響のほとんどは簡単に回避できます。 [20] n {\displaystyle n} ε y n {\displaystyle \varepsilon y_{n}} ε {\displaystyle \varepsilon } ε h {\textstyle {\frac {\varepsilon }{\sqrt {h}}}}
変更と拡張 上述の安定性の問題を解決するオイラー法の単純な修正法として、 後退オイラー法 があります。 これは、関数が ステップの開始点ではなく終了点で評価される点で、（標準的な、あるいは前進的な）オイラー法と異なります。後退オイラー法は 暗黙法 であり、つまり、後退オイラー法の式は 両辺に を持つため、後退オイラー法を適用する際には方程式を解く必要があります。そのため、実装コストは高くなります。 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}).} f {\displaystyle f} y n + 1 {\displaystyle y_{n+1}}
安定性を高めるためにオイラー法に他の修正を加えたものが、 指数オイラー法 または 半暗黙的オイラー法です 。
より複雑な手法を用いることで、より高い次数（およびより高い精度）を実現できます。一つの可能性として、関数評価の回数を増やすことが挙げられます。これは、 本稿で既に述べた 中点法によって例示されます。これは 、ルンゲ・クッタ法 の族へと繋がります 。 y n + 1 = y n + h f ( t n + 1 2 h , y n + 1 2 h f ( t n , y n ) ) {\displaystyle y_{n+1}=y_{n}+hf\left(t_{n}+{\tfrac {1}{2}}h,y_{n}+{\tfrac {1}{2}}hf(t_{n},y_{n})\right)}
もう一つの可能性は、2段階のアダムズ・バッシュフォース法で示されるように、より多くの過去の値を使用することです。
これは 線形多段階法 のファミリーにつながります 。メモリ使用量を最小限に抑えるために圧縮センシングの技術を使用する他の修正法もあります [21] y n + 1 = y n + 3 2 h f ( t n , y n ) − 1 2 h f ( t n − 1 , y n − 1 ) . {\displaystyle y_{n+1}=y_{n}+{\tfrac {3}{2}}hf(t_{n},y_{n})-{\tfrac {1}{2}}hf(t_{n-1},y_{n-1}).}
大衆文化において 映画 『Hidden Figures』 の中で、 キャサリン・ジョンソンは宇宙飛行士 ジョン・グレン の地球軌道からの大気圏再突入を計算する際にオイラー法を用いています 。 [22]
参照
注釈 ^ Butcher 2003, p. 45; Hairer, Nørsett & Wanner 1993, p. 35 ^ アトキンソン 1989, p. 342; ブッチャー 2003, p. 60 ^ ブッチャー 2003, p. 45; ヘアラー、ノーセット、ワナー 1993, p. 36 ^ ブッチャー 2003, p. 3; ヘアラー、ノーセット、ワナー 1993, p. 2 ^ アトキンソン 1989, p. 344 も参照 ^ ハイラー、ノーセット、ワナー、1993、p. 40 ^ アトキンソン 1989, p. 342; ヘアラー、ノーセット、ワナー 1993, p. 36 ^ アトキンソン 1989、342ページ ^ アトキンソン 1989、343ページ ^ ブッチャー 2003、60ページ ^ アトキンソン 1989、342ページ ^ Stoer & Bulirsch 2002、p. 474 ^ アトキンソン 1989、344ページ ^ ブッチャー 2003、49ページ ^ アトキンソン、1989、p. 346;ラコバ 2012、方程式 (1.16) ^ イゼルレス 1996、7ページ ^ ブッチャー 2003、63ページ ^ ブッチャー 2003, p. 70; イザールズ 1996, p. 57 ^ ブッチャー 2003、74～75ページ ^ ブッチャー 2003、75～78ページ ^ Unni, MP; Chandra, MG; Kumar, AA (2017年3月). 「圧縮センシングを用いた微分方程式の数値解法におけるメモリ削減」. 2017 IEEE 第13回国際信号処理とその応用コロキウム (CSPA) . pp. 79– 84. doi :10.1109/CSPA.2017.8064928. ISBN 978-1-5090-1184-1 . S2CID 13082456. ^ カーン、アミナ（2017年1月9日）「アメリカ人を宇宙に送り込んだ『Hidden Figures』の数学者に会う」 ロサンゼルス・タイムズ 。 2017年 2月12日 閲覧
参考文献
外部リンク ウィキブック 微積分学には、 「オイラー法」 に関するページがあります
![{\displaystyle {\begin{aligned}{\vec {y}}_{1}&={\begin{pmatrix}y_{1}\\y_{1}'\\y_{1}''\end{pmatrix}}={\begin{pmatrix}y_{0}\\y'_{0}\\y''_{0}\end{pmatrix}}+h{\begin{pmatrix}y'_{0}\\y''_{0}\\f{\left(t_{0},y_{0},y'_{0},y''_{0}\right)}\end{pmatrix}}\\[1ex]&={\begin{pmatrix}2\\-1\\3\end{pmatrix}}+0.5\cdot {\begin{pmatrix}-1\\3\\sin {0}+(\cos {0})\cdot 2+0^{2}\cdot (-1)-4\cdot 0\cdot 3\end{pmatrix}}\\[1ex]&={\begin{pmatrix}1.5\\0.5\\4\end{pmatrix}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/12d39d9f7b06f2891cd612caf12d43f1895a4c97)
![{\displaystyle {\begin{aligned}{\vec {y}}_{2}&={\begin{pmatrix}y_{2}\\y_{2}'\\y_{2}''\end{pmatrix}}={\begin{pmatrix}y_{1}\\y'_{1}\\y''_{1}\end{pmatrix}}+h{\begin{pmatrix}y'_{1}\\y''_{1}\\f{\left(t_{1},y_{1},y'_{1},y''_{1}\right)}\end{pmatrix}}\\[1ex]&={\begin{pmatrix}1.5\\0.5\\4\end{pmatrix}}+0.5\cdot {\begin{pmatrix}0.5\\4\\sin {0.5}+(\cos {0.5})\cdot 1.5+0.5^{2}\cdot 0.5-4\cdot 0.5\cdot 4\end{pmatrix}}\\[1ex]&={\begin{pmatrix}1.75\\2.5\\0.9604...\end{pmatrix}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b2977cf7737e8f5d2bf33689bc668c308ce85f14)
![{\displaystyle \xi \in [t_{0},t_{0}+h]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/21216efa1ba4d1d6558385ff51b38a37f3ffc42f)
![{\textstyle L={\text{max}}\left(|{\frac {d}{dy}}\left[f(t,y)\right]|\right)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/918c33596aaf1730984ce5857ee16c113dbd1bb4)
![{\textstyle M=\max \left(\left|{\frac {d^{2}}{dt^{2}}}\left[y(t)\right]\right|\right)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0d7390889a132e632a6f49b7a3e07a4638d63817)
![{\displaystyle {\begin{aligned}L&=\max \left|{\frac {d}{dy}}f(t,y)\right|=\max _{2\leq t\leq 3}\left|{\frac {d}{dy}}\left[1+\left(ty\right)^{2}\right]\right|\\[1ex]&=\max _{2\leq t\leq 3}\left|2\left(ty\right)\right|=\max _{2\leq t\leq 3}\left|2\left(t-\left[t+{\frac {1}{t-1}}\right]\right)\right|\\[1ex]&=\max _{2\leq t\leq 3}\left|-{\frac {2}{t-1}}\right|=2\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/06b0eee0130d84eded626517c689bb00244b5b06)
![{\displaystyle {\begin{aligned}M&=\max \left|{\frac {d^{2}}{dt^{2}}}\left[y(t)\right]\right|\\&=\max _{2\leq t\leq 3}\left|{\frac {d^{2}}{dt^{2}}}\left(t+{\frac {1}{1-t}}\right)\right|\\&=\max _{2\leq t\leq 3}\left|{\frac {2}{\left(-t+1\right)^{3}}}\right|=2\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7f40aadeb1ff24778f87c6b1b8bf32cbf55f8bf9)
![{\displaystyle {\begin{aligned}{\text{誤差限界}}&={\frac {hM}{2L}}\left(e^{L(t_{i}-t_{0})}-1\right)\\[1ex]&={\frac {0.5\cdot 2}{2\cdot 2}}\left(e^{2(2.5-2)}-1\right)=0.42957\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/65ddc6a993d9eee4e821c97a268269ca8890e1b4)
ウィキメディア・コモンズのオイラー法に関連するメディア
