● 講義資料

2013年度・後期・数理解析・計算機数学３・第８回 1

● 講義資料

▼ 講義予定

• 常微分方程式の初期値問題の数値解法

– 前進オイラー法・後退オイラー法とその誤差評価 – テイラー展開法

– 修正オイラー法と数値的不安定性

● 前回の講義のまとめ

★ 前進オイラー法

• 最も単純な離散変数法は,前進オイラー法

xn+1=xn+hf(tn, xn)

である. （単に「オイラー法」と言ったときには,この前進オイラー法を指す.）

• 微分方程式を離散変数法で数値的に解く場合には,微分 dx

dt を何らかの形で離散化する必要 がある. そのため,以下のように差分化を行う.

dx

dt(t) = lim

h→0

x(t+h)−x(t) h をt=tⁿ で考えれば,h >0 が小さいと仮定して,

dx

dt(tⁿ) = lim

h→0

x(tⁿ+h)−x(tⁿ)

h ∼x(tn+1)−x(tⁿ) h

と書き直すことができる. （これを前進差分と呼ぶ.）これを, 微分方程式x^′ =f(t, x)に適 用すれば,

x(tⁿ+1)−x(tⁿ)

h ∼f(tⁿ, xⁿ) x(tn+1)∼x(tⁿ) +hf(tⁿ, x(tⁿ)) となる. したがって,前進オイラー法の計算スキーム

xⁿ+1=xⁿ+hf(tⁿ, xⁿ) に相当する式を得ることができる.

• 前進オイラー法のスキームの右辺は,直前のxⁿの値があれば直接計算可能である. その意味 で陽的解法(explicit method)となっている.

• 前進オイラー法のスキームは,直前のxⁿ の値のみで,次のxn+1 を計算可能である. その意 味で１段階法となっている.

Nov. 20, 2013, Version: 1.0 [email protected]

2013年度・後期・数理解析・計算機数学３・第８回 2

★ 後退オイラー法

• 前進オイラー法での差分化を行った点をt=tⁿではなくt=tn+1 での値と考える. すなわち dx

dt(tn+1) = lim

h→0

x(tⁿ+h)−x(tⁿ)

h ∼ x(tn+1)−x(tⁿ) h と書き直すことができる. これを,微分方程式x^′ =f(t, x)に適用すれば,

x(tⁿ+1)−x(tⁿ)

h ∼f(tn+1, xn+1) x(tn+1)∼x(tn) +hf(tn+1, x(tn+1))

となる計算スキームを得ることができる.（これを後退差分と呼ぶ.）このようにして得られ た計算スキーム

xn+1=xⁿ+hf(tn+1, xn+1) を後退オイラー法(backward Euler method)と呼ぶ.

• 後退オイラー法のスキームは,直前のxⁿ の値のみで,次のxn+1 を計算可能である. その意 味で１段階法となっている.

• 後退オイラー法のスキームの右辺は,直前のxⁿの値だけからは直接計算できない. その意味 で陰的解法(implicit method)となっている.

この詳細な意味は以下の通りである.

– 微分方程式x^′(t) =x(t)に後退オイラー法を適用すると, xⁿ+1=xⁿ+hxⁿ+1

を得る. ここで, xn から xn+1 を得るためには, この式を代数的に xn+1 について解 いて,

xn+1= 1 1−hxⁿ として計算する必要がある.

– 微分方程式x^′(t) = sin(x(t))に後退オイラー法を適用すると, xn+1=xn+hsin(xn+1)

を得る. ここで, xⁿ からxn+1 を得るために, この式を代数的にxn+1 について解くこ とはできない. したがって,逐次近似

y0=xⁿ,

y^k+1=xⁿ+hsin(y^k)

を用いて,その極限としてxn+1 を求める必要がある. （なお, hが十分に小さければ, この逐次近似が収束することを証明できる.）

Nov. 20, 2013, Version: 1.0 [email protected]

2013年度・後期・数理解析・計算機数学３・第８回 3

★ オイラー法の誤差

• 今後,考える微分方程式の初期値問題は,

x^′(t) =f(x(t))

の形をしていると仮定する. もし,x^′(t) =f(t, x(t))の形をしている場合には,X(t) = (t, x(t)) とおき,X(t)に関する微分方程式を考えることにより,x^′(t) =f(x(t))に帰着できるので,こ の形の微分方程式だけを考えても一般性を失わないことに注意する必要がある.

• 関数xのテイラー展開と,前進オイラー法の計算スキームを書くと,

x(tⁿ+1) =x(tⁿ) +hx^′(tⁿ) +h²

2 x^′′(tⁿ) +O(h³),

=x(tⁿ) +hf(x(tⁿ)) +h²

2 f(x^′(tⁿ))f(x(tⁿ)) +O(h³), xn+1=xⁿ+hf(xⁿ)

x(tn+1)−xn+1=x(tⁿ)−xⁿ+h(f(x(tⁿ))−f(xⁿ)) +h²

2 f^′(x(tⁿ))f(x(tⁿ)) +O(h³) と計算できる. ここで,f がリプシッツ条件

|f(p)−f(q)| ≤L|p−q| をみたし,f はなめらかであると仮定すると,

|x(tn+1)−xn+1| ≤(1 +Lh)|x(tⁿ)−xⁿ|+M h² が成り立つ.

• 後退オイラー法の場合には,

x(tⁿ) =x(tn+1)−hx^′(tn+1) +h²

2 x^′′(tn+1) +O(h³),

=x(tⁿ+1)−hf(x(tⁿ+1)) +h²

2 f(x^′(tⁿ+1))f(x(tⁿ+1)) +O(h³), xⁿ+1=xⁿ+hf(xⁿ+1)

x(tⁿ+1)−xⁿ+1=x(tⁿ)−xⁿ+h(f(x(tⁿ+1))−f(xⁿ+1)) +h²

2 f^′(x(tⁿ+1))f(x(tⁿ+1)) +O(h³) より,

|x(tn+1)−xn+1| ≤(1−Lh)⁻¹|x(tⁿ)−xⁿ|+M h² が成り立つ.

• 前進オイラー法,後退オイラー法のいずれの場合にも,

ǫⁿ+1=Aǫⁿ+M h², ǫⁿ=|x(tⁿ)−xⁿ|

の形の不等式を得ることができた. この評価を局所離散化誤差と呼び,ある時刻t=tⁿ にお いて,数値解の値 xⁿ と厳密解の値 x(tⁿ)が一致しているとき, １ステップ進んだ際の誤差が O(h²)であることを示している.

Nov. 20, 2013, Version: 1.0 [email protected]

2013年度・後期・数理解析・計算機数学３・第８回 4

• この局所離散化誤差の不等式からt=T =N h においての誤差

|x(T)−x^N| ≤Che^{T L}, 前進オイラー法,

|x(T)−x^N| ≤Che^{T L}, 後退オイラー法,

を得ることができる. すなわち,オイラー法は誤差がh に関して１次の公式であることがわ かる.

● 講義資料

▼ 修正オイラー法による数値計算

0 0.5 1 1.5 2 2.5 3

0 0.2 0.4 0.6 0.8 1

x’ = -x (Modified Euler, h=0.01)

x^′(t) =x(t),x(0) = 1

-1 -0.5 0 0.5 1 1.5

0 0.5 1 1.5 2 2.5 3 3.5 4 x’ = -x (Modified Euler, h=0.01)

-1 -0.5 0 0.5 1 1.5

0 1 2 3 4 5 6 7

x’ = -x (Modified Euler, h=0.001)

x^′(t) =x(t),x(0) = 1,h= 0.01 x^′(t) =−x(t),x(0) = 1,h= 0.001

● 実習内容

1. （★）テイラー展開法を用いて,６回目資料の微分方程式1, 2を解きなさい.

2. （★★）修正オイラー法を用いて,６回目資料の微分方程式1 (λ=−1) を解きなさい.

3. （★★）Fibinacchi 数列xn+2=xn+1+xⁿ,x0= 1,x1= (1−√

5)/2 の一般項の近似値を, この漸化式を使って計算しなさい. その結果が,本来の挙動（一般項をexplicitに計算した値 の挙動）と異ることを確認しなさい.

Nov. 20, 2013, Version: 1.0 [email protected]