オイラー法

提供: testwiki
ナビゲーションに移動 検索に移動

オイラー法(オイラーほう、テンプレート:Lang-en-short)とは、常微分方程式の数値解法の一つである。この方法は、数学的に理解しやすく、プログラム的にも簡単なので、数値解析の初歩的な学習問題としてよく取りあげられる。

定義と公式の導出

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

y=f(t,y),y(t0)=y0.

時間の刻み幅を テンプレート:Mvar とする。テンプレート:Mvar 時間後の数値解を求めるために、まずは時間を離散化し、テンプレート:Math とすると、オイラー法は次の公式で定義される。

yn+1=yn+hf(tn,yn),tnTF.

ここで、テンプレート:Mvarテンプレート:Mvar での数値解である。この公式を導出するために、解の存在性と滑らかさテンプレート:仮リンクより保証されると想定する(特に、テンプレート:Mathリプシッツ連続である)。上記の初期値問題の厳密解(ときには解析解)を テンプレート:Mvar にし、テンプレート:Mathテイラー展開を考える:

y(t+h)=y(t)+y(t)h+12y(t)h2+O(h3).

ここで テンプレート:Math を微分方程式により テンプレート:Math に変換すると上記式が

y(t+h)=y(t)+f(t,y)h+O(h2)

となる。テンプレート:Math の項を切り捨てて、テンプレート:Mvarテンプレート:Mvar に、テンプレート:Math(厳密解)を テンプレート:Mvar(数値解)に置き換えるとオイラー法の公式である。他に微分の定義から公式を導出する方法も存在する。

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

収束性と安定性

複素平面でzがピンク色の円板内部領域はオイラー法の絶対安定性領域である。

数値解析における収束性は、おおよそ刻み幅 テンプレート:Mvar を十分に小さくすると、方法の局所誤差(の絶対値)も小さくなることを意味する。時間 テンプレート:Mvar での局所誤差を en,h=y(tn)yn とする。数学的に、収束性は

limh0en,h=0

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

テンプレート:Math のテイラー展開からオイラー法の公式を引いて両辺の絶対値を取ると

en+1,h=en+h+h(f(tn,y(tn))f(tn,yn))+O(h2)

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

f(tn,y(tn))f(tn,yn)λen,h

を得る。ここでの λf のリプシッツ定数である。三角不等式より上記の両式を合わせて、

en+1,h(1+hλ)en,h+Ch2

という漸化式になる。テンプレート:Mvar は定数であり、テンプレート:Math の係数の絶対値と考えても大差はない。e0,h=y0y(t0)=0 を使ってこの漸化式を解くと上界

en,hCmaxhλ((1+hλ)n1)

がある(帰納法による証明も可能である)。ここで、テンプレート:Math も定数である。固定された時間 tn=t0+nh での局所誤差の上界はゆえに

en,hCmaxhλ(e(tnt0)λ1)

(ここで不等式 1+xexを使った)。上記式から テンプレート:Mvarテンプレート:Math の極限で局所誤差も テンプレート:Math に収束する。すなわち、オイラー法は収束である。 そのうえ、テンプレート:Mvar のテイラー展開を用いて、en,h=O(h2) であることも明らかになる。したがって、オイラー法は1次方法となる。

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

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

{z𝐂|z+1|1},

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

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

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

y=3y+2,y(0)=y0.

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

yn+1=yn+h(3yn+2)

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

yn=(y0+23)(1+3h)n23

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

y=(y0+23)e3x23

となる。テンプレート:Math のときの誤差は

yny(x)=(y0+23)((1+3h)ne3nh)

となる。テンプレート:Mvar のテイラー展開と二項定理を用いて等式

e3nh(1+3h)n=O(h2)

は明らかになる。故に、この例でオイラー法の局所誤差は テンプレート:Math であり、1次方法となる。すなわち、

yny(x)=(y0+23)O(h2)

したがって、この表示から テンプレート:Mvarテンプレート:Math の極限で局所誤差が テンプレート:Math に収束することが分かる。

後退オイラー法

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

y(t)=y(t+h)y(t+h)h+12y(t+h)h2+O(h3)

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

yn+1=yn+hf(tn+1,yn+1)+O(h2)

O(h2) の項を切り捨てて、それが後退オイラー法の公式となる。

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

とはいえ、この方法はA-安定(テンプレート:Lang-en-short)である。よって近似解は刻み幅に関係なく収束する。後退オイラー法も普通は使われない。

拡張

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

yn+1=yn+12h(f(tn,yn)+f(tn+1,yn+1))

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

台形公式は見ての通り陰公式である。それを陽公式にも変換できる。テンプレート:Mathをもう一度オイラー法で近似すると次の公式となる。

yn+1=yn+12h(f(tn,yn)+f(tn+1,yn+hf(tn,yn))

これはテンプレート:仮リンク(または改良オイラー法、テンプレート:Lang-en-short)とよばれる陽公式である。台形公式と同じ2次方法であるけど、A-安定な方法ではない。修正オイラー法は、もっともシンプルな2段ルンゲ=クッタ法である。

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

脚注

テンプレート:Reflist

参考文献

外部リンク

テンプレート:Commonscat