オイラー法

オイラー法（オイラーほう、英: Euler method）とは、常微分方程式の数値解法の一つである。この方法は、数学的に理解しやすく、プログラム的にも簡単なので、数値解析の初歩的な学習問題としてよく取りあげられる。

定義と公式の導出[編集]

常微分方程式とその初期値問題を次のように定める。

${\displaystyle y'=f(t,y),\quad y(t_{0})=y_{0}.}$

時間の刻み幅を h とする。T_F 時間後の数値解を求めるために、まずは時間を離散化し、t_n = t₀ + nh とすると、オイラー法は次の公式で定義される。

${\displaystyle y_{n+1}=y_{n}+hf(t_{n},y_{n}),\quad t_{n}\leq T_{F}.}$

ここで、y_n は t_n での数値解である。この公式を導出するために、解の存在性と滑らかさはピカール・リンデレフの定理より保証されると想定する（特に、f(t, y) はリプシッツ連続である）。上記の初期値問題の厳密解（ときには解析解）を y にし、y(t + h) のテイラー展開を考える：

${\displaystyle y(t+h)=y(t)+y'(t)h+{\frac {1}{2}}y''(t)h^{2}+O(h^{3}).}$

ここで y′(t) を微分方程式により f(t, y) に変換すると上記式が

${\displaystyle y(t+h)=y(t)+f(t,y)h+O(h^{2})}$

となる。O(h²) の項を切り捨てて、t を t_n に、y(t_n)（厳密解）を y_n（数値解）に置き換えるとオイラー法の公式である。他に微分の定義から公式を導出する方法も存在する。

オイラー法は陽公式である。すなわち、過去の値のみが未来の値の計算に必要である。

収束性と安定性[編集]

数値解析における収束性は、おおよそ刻み幅 h を十分に小さくすると、方法の局所誤差（の絶対値）も小さくなることを意味する。時間 t_n での局所誤差を ${\displaystyle e_{n,h}=y(t_{n})-y_{n}}$ とする。数学的に、収束性は

${\displaystyle \lim _{h\to 0}\|e_{n,h}\|=0}$

を意味する。原則として、収束しない、また収束性を証明できない方法は絶対に使ってはいけない^[要出典]。そのため、オイラー法の収束性を示す必要がある。

y(t_n + h) のテイラー展開からオイラー法の公式を引いて両辺の絶対値を取ると

${\displaystyle \|e_{n+1,h}\|=\left\|e_{n+h}+h{\Bigl (}f{\bigl (}t_{n},y(t_{n}){\bigr )}-f(t_{n},y_{n}){\Bigr )}+O(h^{2})\right\|}$

となる。解の滑らかさの仮設よりリプシッツ連続を用いて、不等式

${\displaystyle \|f(t_{n},y(t_{n}))-f(t_{n},y_{n})\|\leq \lambda \|e_{n,h}\|}$

を得る。ここでの ${\displaystyle \lambda }$ は ${\displaystyle f}$ のリプシッツ定数である。三角不等式より上記の両式を合わせて、

${\displaystyle \|e_{n+1,h}\|\leq (1+h\lambda )\|e_{n,h}\|+Ch^{2}}$

という漸化式になる。C は定数であり、h² の係数の絶対値と考えても大差はない。${\displaystyle e_{0,h}=y_{0}-y(t_{0})=0}$ を使ってこの漸化式を解くと上界

${\displaystyle \|e_{n,h}\|\leq {\frac {C_{\mathrm {max} }h}{\lambda }}\left((1+h\lambda )^{n}-1\right)}$

がある（帰納法による証明も可能である）。ここで、C_max も定数である。固定された時間 ${\displaystyle t_{n}=t_{0}+nh}$ での局所誤差の上界はゆえに

${\displaystyle \|e_{n,h}\|\leq {\frac {C_{\mathrm {max} }h}{\lambda }}(e^{(t_{n}-t_{0})\lambda }-1)}$

（ここで不等式 ${\displaystyle 1+x\leq e^{x}}$を使った）。上記式から h が 0 の極限で局所誤差も 0 に収束する。すなわち、オイラー法は収束である。 そのうえ、e^x のテイラー展開を用いて、${\displaystyle \|e_{n,h}\|=O(h^{2})}$ であることも明らかになる。したがって、オイラー法は1次方法となる。

収束性を示したことで、方法が使えるようになる。しかし、収束性が保証できるのは、h が十分小さい場合、近似解が厳密解に収束することのみである。一体 h をどれだけ小さくすれば正しい近似解を得られるのかは一切伝えていない。例えば h を 10⁻¹² 以下にしないと近似解が厳密解に近付かない場合、最低限でも 10¹² 時間での解を計算しなければならないので効率が大きく下がる。そのため、もし h に関係なく近似解が思わぬ行動を取れないことを示せるなら、h を自由に設定できてそのような心配はいらなくなる。上述の条件が満たされる方法は、おおよそ数値的に安定（正しくいうとA-安定）である。厳密な定義や他の安定性（L-安定、零点安定他）については、硬い方程式を参照。

線形微分方程式${\displaystyle y'=ky}$をオイラー法により求める場合 ${\displaystyle z=hk}$を複素平面にとったとき図の円より外側の領域で数値解が不安定となる。

${\displaystyle \{z\in \mathbf {C} \mid |z+1|\leq 1\},}$

図の領域は線形安定領域と呼ばれる。^[1] 例えば ${\displaystyle k=-2.3}$の場合では、時間の刻み幅${\displaystyle h=1}$では${\displaystyle z=hk=-2.3}$となる。

よってこの ${\displaystyle z}$ では安定領域より外側の領域のためオイラー法の数値解は不安定となる。

例[編集]

簡単な例として、次の1階常微分方程式を考えよう：

${\displaystyle y'=3y+2,\quad y(0)=y_{0}.}$

この方程式に対するのオイラー法は

${\displaystyle y_{n+1}=y_{n}+h(3y_{n}+2)}$

という漸化式になる。この漸化式は厳密解を求めることができて、初期値を用いて、

${\displaystyle y_{n}=\left(y_{0}+{\frac {2}{3}}\right)(1+3h)^{n}-{\frac {2}{3}}}$

となる。この例では、微分方程式の厳密解を直接に求めることができる。解いて厳密解は

${\displaystyle y=\left(y_{0}+{\frac {2}{3}}\right)e^{3x}-{\frac {2}{3}}}$

となる。x = nh のときの誤差は

${\displaystyle y_{n}-y(x)=\left(y_{0}+{\frac {2}{3}}\right)\left((1+3h)^{n}-e^{3nh}\right)}$

となる。e^x のテイラー展開と二項定理を用いて等式

${\displaystyle e^{3nh}-(1+3h)^{n}=O(h^{2})}$

は明らかになる。故に、この例でオイラー法の局所誤差は O(h²) であり、1次方法となる。すなわち、

${\displaystyle y_{n}-y(x)=\left(y_{0}+{\frac {2}{3}}\right)O(h^{2})}$

したがって、この表示から h が 0 の極限で局所誤差が 0 に収束することが分かる。

後退オイラー法[編集]

前に述べたように、オイラー法は数値的不安定である。硬い微分方程式の解を計算するときに刻み幅を非常に小さくしないと近似解は厳密解に収束しない。そのような方程式の近似解を求めるために、数値的安定な方法が必要となる。安定性を持つ方法の中では、後退オイラー法（英: backward Euler method）が一番シンプルな方法である。y(t + h) の代わりに、y(t) の t + h に対するテイラー展開を考える

${\displaystyle y(t)=y(t+h)-y'(t+h)h+{\frac {1}{2}}y''(t+h)h^{2}+O(h^{3})}$

前と同じ置き換えをすると、次の公式となる。

${\displaystyle y_{n+1}=y_{n}+hf(t_{n+1},y_{n+1})+O(h^{2})}$

${\displaystyle O(h^{2})}$ の項を切り捨てて、それが後退オイラー法の公式となる。

右側に ${\displaystyle f(t_{n+1},y_{n+1})}$ の項が存在するため、後退オイラー法は陰公式となる。すなわち、現在値を求めるために過去の数値だけではなく、現在や未来の数値を知らなければならない（後退オイラー法では未来値に依存しない）。結果として、一時刻ごとに方程式系を解かなければならない。ゆえに、特に非線型方程の場合、計算コストが大分上がる。

とはいえ、この方法はA-安定（英: A-stable）である。よって近似解は刻み幅に関係なく収束する。後退オイラー法も普通は使われない。

拡張[編集]

オイラー法は１次方法である。１次方法は極端に精度が低いので、実践には向かない。そのため、もっと次数の高い方法が必要となる。オイラー法と後退オイラー法の公式をもとにして平均値を取ると次の公式となる。

${\displaystyle y_{n+1}=y_{n}+{\frac {1}{2}}h(f(t_{n},y_{n})+f(t_{n+1},y_{n+1}))}$

これは微分方程式での台形公式である。有限差分法を用いて偏微分方程式に適用するときはクランク・ニコルソン法とも呼ばれる。この方法が2次方法でA-安定であることを証明できる^[2]。

台形公式は見ての通り陰公式である。それを陽公式にも変換できる。y_n+1をもう一度オイラー法で近似すると次の公式となる。

${\displaystyle y_{n+1}=y_{n}+{\frac {1}{2}}h(f(t_{n},y_{n})+f(t_{n+1},y_{n}+hf(t_{n},y_{n}))}$

これはホイン法（英語版）（または改良オイラー法、英: improved Euler method）とよばれる陽公式である。台形公式と同じ2次方法であるけど、A-安定な方法ではない。修正オイラー法は、もっともシンプルな2段ルンゲ＝クッタ法である。

しかし、2次方法はオイラー法より精度が遥かに高いが、実践で使うにはまだ精度が足りない。現在よく使われているのが高次のルンゲ＝クッタ法である（MATLABのコマンドode23は3段3次ルンゲ＝クッタ法で、ode45は6段5次のルンゲ＝クッタ法である^[3]）。

脚注[編集]

参考文献[編集]

• Iserles, Arieh (2008), A First Course in the Numerical Analysis of Differential Equations (Second Edition), Cambridge University Press, ISBN 978-0-521-73490-5 .

外部リンク[編集]

ウィキメディア・コモンズには、オイラー法に関連するカテゴリがあります。