オイラー法のソースを表示

'''オイラー法'''（オイラーほう、{{lang-en-short|Euler method}}）とは、[[常微分方程式の数値解法]]の一つである。この方法は、数学的に理解しやすく、[[プログラム (コンピュータ)|プログラム]]的にも簡単なので、[[数値解析]]の初歩的な学習問題としてよく取りあげられる。

== 定義と公式の導出 ==
常微分方程式とその[[初期値問題]]を次のように定める。
: <math> y' = f(t,y), \quad y(t_0) = y_0. </math>
時間の刻み幅を {{mvar|h}} とする。{{mvar|T{{sub|F}}}} 時間後の数値解を求めるために、まずは時間を離散化し、{{math|1=''t{{sub|n}}'' = ''t''{{sub|0}} + ''nh''}} とすると、オイラー法は次の公式で定義される。
:<math> y_{n+1} = y_n + hf(t_n,y_n), \quad t_n \le T_F. </math>
ここで、{{mvar|y{{sub|n}}}} は {{mvar|t{{sub|n}}}} での数値解である。この公式を導出するために、解の存在性と[[滑らかさ (数学)|滑らかさ]]は{{仮リンク|ピカール・リンデレフの定理|en|Picard–Lindelöf theorem}}より保証されると想定する（特に、{{math|''f''(''t'', ''y'')}} は[[リプシッツ連続]]である）。上記の初期値問題の厳密解（ときには[[解析解]]）を {{mvar|y}} にし、{{math|''y''(''t'' + ''h'')}} の[[テイラー展開]]を考える：
:<math> y(t+h) = y(t) + y'(t)h + \frac{1}{2}y''(t)h^2 + O(h^3). </math>
ここで {{math|''y&prime;''(''t'')}} を微分方程式により {{math|''f''(''t'', ''y'')}} に変換すると上記式が
:<math> y(t+h) = y(t) + f(t,y)h + O(h^2) </math>
となる。{{math|''O''(''h''{{sup|2}})}} の項を切り捨てて、{{mvar|t}} を {{mvar|t{{sub|n}}}} に、{{math|''y''(''t{{sub|n}}'')}}（厳密解）を {{mvar|y{{sub|n}}}}（数値解）に置き換えるとオイラー法の公式である。他に微分の定義から公式を導出する方法も存在する。
<!--
常微分方程式系の場合、{{mvar|y}} を <math>\vec{y} = (y_1,\ldots,y_k)</math> に置き換えて方法はそのまま適用できる。
-->

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

== 収束性と安定性 == 
[[Image:Stability region for Euler method.svg|thumb|複素平面でzがピンク色の円板内部領域はオイラー法の絶対安定性領域である。]]

数値解析における[[収束性]]は、おおよそ刻み幅 {{mvar|h}} を十分に小さくすると、方法の局所[[誤差]]（の絶対値）も小さくなることを意味する。時間 {{mvar|t{{sub|n}}}} での局所誤差を <math>e_{n,h} = y(t_n) - y_n</math> とする。数学的に、収束性は
:<math> \lim_{h \to 0} \|e_{n,h}\| = 0 </math>
を意味する。原則として、収束しない、また収束性を証明できない方法は絶対に使ってはいけない{{要出典|date=2022年3月}}。そのため、オイラー法の収束性を示す必要がある。

{{math|''y''(''t{{sub|n}}'' + ''h'')}} のテイラー展開からオイラー法の公式を引いて両辺の[[絶対値]]を取ると
:<math>\|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\|</math>
となる。解の滑らかさの仮設よりリプシッツ連続を用いて、不等式
:<math> \|f(t_n,y(t_n)) - f(t_n,y_n)\| \le \lambda \|e_{n,h}\| </math>
を得る。ここでの <math>\lambda</math> は <math> f </math> のリプシッツ定数である。[[三角不等式]]より上記の両式を合わせて、
:<math> \|e_{n+1,h}\| \le (1+h\lambda) \|e_{n,h}\| + Ch^2 </math>
という[[漸化式]]になる。{{mvar|C}} は定数であり、{{math|''h''{{sup|2}}}} の係数の絶対値と考えても大差はない。<math> e_{0,h} = y_0 - y(t_0) = 0 </math> を使ってこの漸化式を解くと[[上界 (数学)|上界]]
:<math> \|e_{n,h}\| \le \frac{C_{\mathrm{max}} h}{\lambda}\left((1 + h\lambda)^{n} - 1\right) </math>
がある（[[数学的帰納法|帰納法]]による証明も可能である）。ここで、{{math|''C''{{sub|max}}}} も定数である。固定された時間 <math> t_n = t_0 + nh</math> での局所誤差の上界はゆえに
:<math> \|e_{n,h}\|\le \frac{C_{\mathrm{max}} h}{\lambda}(e^{(t_n-t_0)\lambda}- 1) </math>
（ここで不等式 <math> 1+x \le e^x </math>を使った）。上記式から {{mvar|h}} が {{math|0}} の極限で局所誤差も {{math|0}} に収束する。すなわち、オイラー法は収束である。 そのうえ、{{mvar|e{{sup|x}}}} のテイラー展開を用いて、<math>\|e_{n,h}\| = O(h^2)</math> であることも明らかになる。したがって、オイラー法は1次方法となる。

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

線形微分方程式<math>y' = k y</math>をオイラー法により求める場合
<math>z=hk</math>を複素平面にとったとき図の円より外側の領域で数値解が不安定となる。
:<math> \{ z \in \mathbf{C} \mid |z+1| \le 1 \}, </math>

図の領域は線形安定領域と呼ばれる。<ref>{{harvnb|Butcher|2003|p=70}}; {{harvnb|Iserles|1996|p=57}}</ref>
例えば <math>k = -2.3</math>の場合では、時間の刻み幅<math>h=1</math>では<math>z = hk = -2.3</math>となる。

よってこの <math> z </math> では安定領域より外側の領域のためオイラー法の数値解は不安定となる。
<!--
残念ながらオイラー法は、絶対安定性領域 (region of absolute stability) が <math>\{z \in \mathbb{C}\mid |1+z| < 1 \}</math> という左複素数平面の円板であって、A-安定な方法ではない。そのせいで、実装するときに方程式の硬さを予め計算しなければならない。硬い方程式の場合では、他の方法を試そう。
-->

== 例 ==
簡単な例として、次の1階常微分方程式を考えよう：
:<math> y'= 3y + 2, \quad y(0) =y_0. </math>
この方程式に対するのオイラー法は
:<math> y_{n+1} = y_n + h(3y_n +2)</math>
という漸化式になる。この漸化式は厳密解を求めることができて、初期値を用いて、
:<math> y_n=\left(y_0+\frac{2}{3}\right)(1+3h)^n-\frac{2}{3}</math>
となる。この例では、微分方程式の厳密解を直接に求めることができる。解いて厳密解は
:<math> y=\left(y_0+\frac{2}{3}\right)e^{3x}-\frac{2}{3}</math>
となる。{{math|1=''x'' = ''nh''}} のときの誤差は
:<math>y_n-y(x)=\left(y_0+\frac{2}{3}\right)\left((1+3h)^{n}-e^{3nh}\right)</math>
となる。{{mvar|e{{sup|x}}}} のテイラー展開と[[二項定理]]を用いて等式
:<math>e^{3nh} - (1+3h)^n = O(h^2)</math>
は明らかになる。故に、この例でオイラー法の局所誤差は {{math|''O''(''h''{{sup|2}})}} であり、1次方法となる。すなわち、
:<math>y_n-y(x)=\left(y_0+\frac{2}{3}\right)O(h^2)</math>
したがって、この表示から {{mvar|h}} が {{math|0}} の極限で局所誤差が {{math|0}} に収束することが分かる。

== 後退オイラー法 ==
前に述べたように、オイラー法は[[数値的安定性|数値的不安定]]である。硬い微分方程式の解を計算するときに刻み幅を非常に小さくしないと近似解は厳密解に収束しない。そのような方程式の近似解を求めるために、数値的安定な方法が必要となる。安定性を持つ方法の中では、後退オイラー法（{{Lang-en-short|backward Euler method}}）が一番シンプルな方法である。{{math|''y''(''t'' + ''h'')}} の代わりに、{{math|''y''(''t'')}} の {{math|''t'' + ''h''}} に対するテイラー展開を考える
:<math>y(t) = y(t+h) - y'(t+h)h + \frac{1}{2}y''(t+h)h^2 + O(h^3)</math>
前と同じ置き換えをすると、次の公式となる。
:<math> y_{n+1} = y_n + hf(t_{n+1},y_{n+1}) + O(h^2) </math>
<math> O(h^2)</math> の項を切り捨てて、それが後退オイラー法の公式となる。

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

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

== 拡張 ==
オイラー法は１次方法である。１次方法は極端に精度が低いので、実践には向かない。そのため、もっと次数の高い方法が必要となる。オイラー法と後退オイラー法の公式をもとにして平均値を取ると次の公式となる。
:<math> y_{n+1} = y_n + \frac{1}{2}h(f(t_n,y_n)+f(t_{n+1},y_{n+1})) </math>
これは微分方程式での[[台形公式]]である。[[有限差分法]]を用いて[[偏微分方程式]]に適用するときはクランク・ニコルソン法とも呼ばれる。この方法が2次方法でA-安定であることを証明できる<ref>{{harvnb|Iserles|2008|loc = Section 1.3}}</ref>。

台形公式は見ての通り陰公式である。それを陽公式にも変換できる。{{math|''y{{sub|n+1}}''}}をもう一度オイラー法で近似すると次の公式となる。
:<math> 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)) </math>
これは{{仮リンク|ホイン法|en|Heun's method}}（または改良オイラー法、{{Lang-en-short|improved Euler method}}）とよばれる陽公式である。台形公式と同じ2次方法であるけど、A-安定な方法ではない。修正オイラー法は、もっともシンプルな2段[[ルンゲ＝クッタ法]]である。

しかし、2次方法はオイラー法より精度が遥かに高いが、実践で使うにはまだ精度が足りない。現在よく使われているのが高次のルンゲ＝クッタ法である（MATLABのコマンドode23は3段3次ルンゲ＝クッタ法で、ode45は6段5次のルンゲ＝クッタ法である<ref>{{Cite web |url=http://blogs.mathworks.com/cleve/2014/05/26/ordinary-differential-equation-solvers-ode23-and-ode45/ |title=Ordinary Differential Equation Solvers ODE23 and ODE45 |author=Cleve Moler |accessdate=2016-12-16}}</ref>）。

== 脚注 ==
{{Reflist}}

== 参考文献 ==
* {{Citation|title=A First Course in the Numerical Analysis of Differential Equations (Second Edition)|year=2008|last1=Iserles|first1=Arieh|publisher=[[Cambridge University Press]]|isbn=978-0-521-73490-5}}.

== 外部リンク ==
{{Commonscat|Euler method}}

{{DEFAULTSORT:おいらあほう}}
[[Category:解析学]]
[[Category:数値解析]]
[[Category:数値微分方程式]]
[[Category:数学に関する記事]]
[[Category:レオンハルト・オイラー]]