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{\bigl (}t,y(t){\bigr )}} 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}
高階プロセス オイラー法は一階常微分方程式を積分しますが、任意の階数の常微分方程式は 一階常微分方程式の連立方程式として表すことができます 。 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時 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]
高階の例 この 3 次例では、次の情報が与えられていると仮定します。
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'''を分離して次の方程式を得ることができます。
f ( t , y , y ′ , y ″ ) = y ‴ = sin t + ( cos t ) y + t 2 y ′ − 4 t y ″ {\displaystyle f\left(t,y,y',y''\right)=y'''=\sin {t}+(\cos {t})y+t^{2}y'-4ty''}
これを使用すると、 の解を得ることができます 。 また、 の解を使用すると 、 の解を得ることができます 。 必要なだけ同じ式を使用してこのプロセスを続けることで、希望するものを見つけることができます 。 y → 1 {\displaystyle {\vec {y}}_{1}} y → 1 = ( y 1 y 1 ′ y 1 ″ ) = ( y 0 + h ⋅ y 0 ′ y 0 ′ + h ⋅ y 0 ″ y 0 ″ + h ⋅ f ( t 0 , y 0 , y 0 ′ , y 0 ″ ) ) = ( 2 + 0.5 ⋅ − 1 − 1 + 0.5 ⋅ 3 3 + 0.5 ⋅ ( sin 0 + ( cos 0 ) ⋅ 2 + 0 2 ⋅ ( − 1 ) − 4 ⋅ 0 ⋅ 3 ) ) = ( 1.5 0.5 4 ) {\displaystyle {\vec {y}}_{1}={\begin{pmatrix}y_{1}\\y_{1}'\\y_{1}''\end{pmatrix}}={\begin{pmatrix}y_{0}+h\cdot y'_{0}\\y'_{0}+h\cdot y''_{0}\\y''_{0}+h\cdot f\left(t_{0},y_{0},y'_{0},y''_{0}\right)\end{pmatrix}}={\begin{pmatrix}2+0.5\cdot -1\\-1+0.5\cdot 3\\3+0.5\cdot \left(\sin {0}+(\cos {0})\cdot 2+0^{2}\cdot (-1)-4\cdot 0\cdot 3\right)\end{pmatrix}}={\begin{pmatrix}1.5\\0.5\\4\end{pmatrix}}} y → 1 {\displaystyle {\vec {y}}_{1}} y → 2 {\displaystyle {\vec {y}}_{2}} y → 2 = ( y 2 y 2 ′ y 2 ″ ) = ( y 1 + h ⋅ y 1 ′ y 1 ′ + h ⋅ y 1 ″ y 1 ″ + h ⋅ f ( t 1 , y 1 , y 1 ′ , y 1 ″ ) ) = ( 1.5 + 0.5 ⋅ 0.5 0.5 + 0.5 ⋅ 4 4 + 0.5 ⋅ ( sin 0.5 + ( cos 0.5 ) ⋅ 1.5 + 0.5 2 ⋅ 0.5 − 4 ⋅ 0.5 ⋅ 4 ) ) = ( 1.75 2.5 0.9604... ) {\displaystyle {\vec {y}}_{2}={\begin{pmatrix}y_{2}\\y_{2}'\\y_{2}''\end{pmatrix}}={\begin{pmatrix}y_{1}+h\cdot y'_{1}\\y'_{1}+h\cdot y''_{1}\\y''_{1}+h\cdot f\left(t_{1},y_{1},y'_{1},y''_{1}\right)\end{pmatrix}}={\begin{pmatrix}1.5+0.5\cdot 0.5\\0.5+0.5\cdot 4\\4+0.5\cdot \left(\sin {0.5}+(\cos {0.5})\cdot 1.5+0.5^{2}\cdot 0.5-4\cdot 0.5\cdot 4\right)\end{pmatrix}}={\begin{pmatrix}1.75\\2.5\\0.9604...\end{pmatrix}}} 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)}
以下では、テイラー展開を使用してオイラー法で発生する誤差を解析します。テイラー展開を拡張すると、 ルンゲ・クッタ法が 生成されます。
（３） これに密接に関連する導出法は、導関数を
前向き 差分 式に置き換えることである。
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}}} 微分方程式の 。これもまたオイラー法となる。 [8] y ′ = f ( t , y ) {\displaystyle y'=f(t,y)}
同様の計算により、 中点法 と 後方オイラー法 が導き出されます。
（４） 最後に、微分方程式をからまで積分し、微積分学の基本定理を適用すると次 の 式 が 得られる。 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{\bigl (}t,y(t){\bigr )}\,\mathrm {d} t.} ここで、積分を左側の 長方形法 （長方形は 1 つだけ）で近似します。
∫ 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{\bigl (}t,y(t){\bigr )}\,\mathrm {d} t\approx hf{\bigl (}t_{0},y(t_{0}){\bigr )}.} 両方の方程式を組み合わせると、再びオイラー法が見つかります。 [9]
この考え方を継続することで、さまざまな 線形多段階法 に到達できます。
ローカル切り捨てエラー オイラー法の局所 的打ち切り誤差は 、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})+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})+O\left(h^{3}\right).} この結果は、 3次導関数が有界である場合に有効である。 [10] y {\displaystyle y}
これは、 が小さい場合 、局所的な打ち切り誤差が にほぼ比例することを示しています 。そのため、オイラー法は、局所的な打ち切り誤差がステップサイズのより高いべき乗に比例する ルンゲ・クッタ法 や 線形マルチステップ法 などの高次手法よりも精度が低くなります。 h {\displaystyle h} h 2 {\displaystyle h^{2}}
局所的打ち切り誤差の若干異なる定式化は、 テイラーの定理 の剰余項のラグランジュ形式を用いることで得られる 。が 連続な2次微分を持つ場合、 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 ).} [11] 上記の誤差の式において、未知の厳密解の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}}{\bigl (}t_{0},y(t_{0}){\bigr )}+{\frac {\partial f}{\partial y}}{\bigl (}t_{0},y(t_{0}){\bigr )}\,f{\bigl (}t_{0},y(t_{0}){\bigr )}.}
グローバル切り捨てエラー グローバル 打ち切り誤差 とは、初期時刻からその時刻に到達するまでに何ステップ必要であっても、そのステップを経た後の、 固定時刻 における誤差のことである。グローバル打ち切り誤差は、各ステップで発生するローカル打ち切り誤差の累積効果である。 [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 |y(t_{i})-y_{i}|\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}}{\bigl (}|{\frac {d}{dy}}{\bigl [}f(t,y){\bigr ]}|{\bigr )}} t {\displaystyle t} M = max ( | d 2 d t 2 [ y ( t ) ] | ) {\textstyle M={\text{max}}{\bigl (}|{\frac {d^{2}}{dt^{2}}}{\bigl [}y(t){\bigr ]}|{\bigr )}} 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 L={\text{max}}{\bigl (}|{\frac {d}{dy}}{\bigl [}f(t,y){\bigr ]}|{\bigr )}=\max _{2\leq t\leq 3}{\bigl (}|{\frac {d}{dy}}{\bigl [}1+(t-y)^{2}{\bigr ]}|{\bigr )}=\max _{2\leq t\leq 3}{\bigl (}|2(t-y)|{\bigr )}=\max _{2\leq t\leq 3}{\bigl (}|2(t-[t+{\frac {1}{t-1}}])|{\bigr )}=\max _{2\leq t\leq 3}{\bigl (}|-{\frac {2}{t-1}}|{\bigr )}=2} 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 M={\text{max}}{\bigl (}|{\frac {d^{2}}{dt^{2}}}{\bigl [}y(t){\bigr ]}|{\bigr )}=\max _{2\leq t\leq 3}\left(|{\frac {d^{2}}{dt^{2}}}{\bigl [}t+{\frac {1}{1-t}}{\bigr ]}|\right)=\max _{2\leq t\leq 3}\left(|{\frac {2}{(-t+1)^{3}}}|\right)=2}
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 {\text{error bound}}={\frac {hM}{2L}}\left(e^{L(t_{i}-t_{0})}-1\right)={\frac {0.5\cdot 2}{2\cdot 2}}\left(e^{2(2.5-2)}-1\right)=0.42957} 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) ピンク色の円盤はオイラー法の安定領域を示しています。 オイラー法を線形方程式に適用すると、積が 領域外にある 場合、数値解は不安定になる。 y ′ = k y {\displaystyle y'=ky} h k {\displaystyle hk}
{ z ∈ C | | z + 1 | ≤ 1 } , {\displaystyle {\bigl \{}z\in \mathbf {C} \,{\big |}\,|z+1|\leq 1{\bigr \}},} 右に図示されているように、この領域は（線形） 安定領域 と呼ばれます。 [18] 例では、 となるので、 であればとなり 、これは安定領域外となり、したがって数値解は不安定になります。 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)} 。 これにより、ルンゲ・クッタ法 のファミリーが生まれます 。
もう 1 つの可能性は、2 段階の Adams-Bashforth 法で示されているように、より多くの過去の値を使用することです。
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}).} この手法は線形多段階法 のファミリーにつながる 。メモリ使用量を最小限に抑えるために圧縮センシングの技術を用いた他の改良法もある [21]。
大衆文化において 映画 『Hidden Figures』 では、 キャサリン・ジョンソンは宇宙飛行士 ジョン・グレン の地球軌道からの大気圏再突入を計算する際にオイラー法を利用している 。 [22]
参照
注記 ^ ブッチャー 2003, p. 45; ヘアラー、ノーセット、ワナー 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 \xi \in [t_{0},t_{0}+h]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/21216efa1ba4d1d6558385ff51b38a37f3ffc42f)
![{\textstyle L={\text{max}}{\bigl (}|{\frac {d}{dy}}{\bigl [}f(t,y){\bigr ]}|{\bigr )}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8795e9535eb3a8b2c7553a589ed946574b5d09f3)
![{\textstyle M={\text{max}}{\bigl (}|{\frac {d^{2}}{dt^{2}}}{\bigl [}y(t){\bigr ]}|{\bigr )}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ad0a378ffc1cbdbd1593d11a81adc87ef810f5f2)
![{\displaystyle L={\text{max}}{\bigl (}|{\frac {d}{dy}}{\bigl [}f(t,y){\bigr ]}|{\bigr )}=\max _{2\leq t\leq 3}{\bigl (}|{\frac {d}{dy}}{\bigl [}1+(ty)^{2}{\bigr ]}|{\bigr )}=\max _{2\leq t\leq 3}{\bigl (}|2(ty)|{\bigr )}=\max _{2\leq t\leq 3}{\bigl (}|2(t-[t+{\frac {1}{t-1}}])|{\bigr )}=\max _{2\leq t\leq 3}{\bigl (}|-{\frac {2}{t-1}}|{\bigr )}=2}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f92322351c8d87a6f5ba41e0335b86f502c795bf)
![{\displaystyle M={\text{max}}{\bigl (}|{\frac {d^{2}}{dt^{2}}}{\bigl [}y(t){\bigr ]}|{\bigr )}=\max _{2\leq t\leq 3}\left(|{\frac {d^{2}}{dt^{2}}}{\bigl [}t+{\frac {1}{1-t}}{\bigr ]}|\right)=\max _{2\leq t\leq 3}\left(|{\frac {2}{(-t+1)^{3}}}|\right)=2}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d13641f653285747f8f2895d42bbbce688895190)
ウィキメディア・コモンズのオイラー法関連メディア
