● 前回の講義のまとめ

● 前回の講義のまとめ

常微分方程式の初期値問題

x(t) =f(t, x(t)), x(0) =x0

(1)

の数値解を構成する計算スキームを考える. 以下では, f(t, p)は滑らかであると仮定する.

【後退オイラー法】

1. 差分化を

dx

dt(tk)←→ xk−xk+1

−h と後退微分で置き換えると,

xk+1=xk+hf(tk+1, xk+1)

によって{xk}を与える計算スキームを得ることができる. これを後退オイラー法と呼ぶ.

2. 後退オイラー法に対しても,pに関するリプシッツ条件と ^∂f_∂t, ^∂f_∂p の有界性があれば,局所離散 化誤差ek と,大域的離散化誤差Ek =|x(tk)−xk|は

ek=O(h²), Ek =O(h) をみたす.

【テイラー展開法】

1. オイラー法よりも高次の近似を実現する計算スキームとして, xk+1=xk+hf(tk, xk) +h²

2!F2(tk, xk) +· · ·+h^p

p!Fp(tk, xk)

によって{xk} を与えるテイラー法を考えることができる. ここで, Fj(t, x) は, ^d_dt^j−1j f(t, x(t)) の x(t)などに, 微分方程式の条件x =f(t, x)を代入したものである.

2. テイラー法は,x(t)のテイラー展開

x(tk+1) =x(tk) +hx(tk) +h²

2!x(tk) +· · ·h^p

p!x^(p)(t) +O(h^p+1)

を利用した近似スキームである. （オイラー法は,テイラー法のp= 1に相当するスキームであ る.）したがって,f(t, p)に必要なだけの滑らかさがあれば,

ek =O(h^p+1), Ek =O(h^p) が成り立つ.

3. しかし,テイラー法を実際に適用するには f(t, p)の高次導関数を計算する必要があるため, 実 用上は用いられない.

【修正オイラー法】

1. オイラー法で用いた x(t)の近似である前進差分

x(tk)∼ x(tk+1)−x(tk) h だけでなく,後退差分

x(tk)∼ x(tk−1)−x(tk)

−h を組み合わせた中心差分（中点則）

x(tk)∼x(t_k+1)−x(t_k−1) 2h

を用いることを考える. これは,以下のように考えることができる.

x(t_k+1) =x(tk) +hx(tk) +h²

2 x(tk) +h³

3!x⁽³⁾(tk) +O(h⁴), x(tk−1) =x(tk)−hx(tk) +h²

2 x(tk)−h³

3!x⁽³⁾(tk) +O(h⁴), の上の式から下を引けば,

x(tk+1)−x(tk−1) = 2hx(tk) +h³

3 x⁽³⁾(tk) +O(h⁵) となる. このh³以上の項を無視すれば,

xk+1=xk−1+ 2hf(t_k, xk) というスキームを得る. これを修正オイラー法と呼ぶ.

2. 修正オイラー法は,初期値x0 だけではなく,x1 （出発値と呼ぶ）を必要とする. 通常はx1 を 求めるためにオイラー法を用いる.

3. 修正オイラー法に関しては,

ek =O(h³), Ek=O(h²) が成り立つ.

4. いま,最も安定と考えられる方程式

x(t) =−x(t), x(0) =x0

(2) に対して修正オイラー法を適用しよう. このとき,修正オイラー法のスキームは,３項間漸化式

xk+1=xk−1−2hxk (3)

と等しい. いま, (3) の一般項は, (3) の特性方程式

t²+ 2ht−1 = 0 の２つの解

η1=−h+ h²+ 1, η2=−h−

h²+ 1

ex08.tex,v 1.1 2007-05-27 16:19:24+09 naito Exp

を用いて,

xk =c1η^k₁+c2η^k₂ (4)

と書くことができる. いま, η2 <−1 より, c2 = 0 ならば, |xk| → ∞ となり, (1) の解の挙動 x(t)→0 とは大きく異なることになる.

ここで, 初期条件x0と出発値x1 を用いて, c1,c2 は c2=x1−x0

η2−η1 = x0−x1

2√ h²+ 1, c1=x0−c1

と書ける. したがって, x0=x1 であるかぎりは,任意のh >0 に対してc2= 0である.

5. 修正オイラー法を用いて(3) の数値解を構成すると, 真の解とは全く異なる振動しながら「発 散」する解があらわれる. このような現象を数値的不安定性と呼び,本来「安定」な数学的対象 を「離散化」または「数値化」したことにより発生する不安定性である. 最も安定な対象である と思われるx(t) =−x(t)でさえも不安定となってしまうため,修正オイラー法も実用上は用い られない.

● 講義資料

x(t) =λx(t), x(0) =x0, λ= 0. (5)

x(t) = (1−x(t)²), x(0) = 0. (6)

x(t) +ax(t) +bx(t) =f(x), x(0) =x0, x(0) =v0. (7) x(t) + sin(x(t)) = 0, x(0) =x0, x(0) =v0. (8) x(t) + (x(t)²−1)x(t) +x(t) = 0, x(0) =x0, x(0) =v0. (9)

★ 改良オイラー法

• x=x, x0= 1.0に対して, h= 0.5から順にhを 1/2倍したときの数値解と,真の解との誤差.

1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8

0 0.2 0.4 0.6 0.8 1

"result_1.txt"

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08

0 0.2 0.4 0.6 0.8 1

"error_1.txt"

• x=−x,x0= 1.0 に対して,h= 0.5から順にhを1/2 倍したときの数値解と,真の解との誤差.

0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 0.2 0.4 0.6 0.8 1

"result_3.txt"

0 0.5 1 1.5 2 2.5

0 0.2 0.4 0.6 0.8 1

"error_3.txt"

ex08.tex,v 1.1 2007-05-27 16:19:24+09 naito Exp

★ ルンゲクッタ型公式

• x =x,前進オイラー法, 改良オイラー法,ホインの３次公式,ルンゲクッタ法で t= 1.0 での誤差を 示したもの.

1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1

10 100 1000 10000 100000 1e+06 1e+07

Euler Improved Euler Heun 3 Runge-Kutta

• x = −x, 前進オイラー法, 改良オイラー法, ホインの３次公式, ルンゲクッタ法で t = 4π での (x(t))²+ (x(t))² と1.0との誤差を示したもの.

1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1 100

10 100 1000 10000 100000 1e+06

Euler Improved Euler Heun 3 Runge-Kutta

• x=−sin(x), 前進オイラー法, 改良オイラー法, ホインの３次公式,ルンゲクッタ法でt= 4πでの

12(x(t))²−cos(x(t))と ¹₂(x(0))²−cos(x(0))との誤差を示したもの.

1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 0.0001 0.01 1

10 100 1000 10000 100000

Euler Improved Euler Heun 3 Runge-Kutta

• x=−x, (x(t))²+ (x(t))² の値の変化の様子. 前進オイラー法,改良オイラー法,ホインの３次公式, ルンゲクッタ法. (t= 0.01)

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14

0 2 4 6 8 10 12 14

Euler

0 5e-07 1e-06 1.5e-06 2e-06 2.5e-06 3e-06 3.5e-06

0 2 4 6 8 10 12 14

i-Euler

前進オイラー法 改良オイラー法

-1.2e-06 -1e-06 -8e-07 -6e-07 -4e-07 -2e-07 0

0 2 4 6 8 10 12 14

Heun 3

-1.8e-11 -1.6e-11 -1.4e-11 -1.2e-11 -1e-11 -8e-12 -6e-12 -4e-12 -2e-12 0

0 2 4 6 8 10 12 14

Rk

ホインの３次公式 ルンゲクッタ

ex08.tex,v 1.1 2007-05-27 16:19:24+09 naito Exp

• x=−sin(x), ¹₂(x(t))²−cos(x(t))の値の変化の様子. 前進オイラー法,改良オイラー法,ホインの ３次公式,ルンゲクッタ法. (t= 0.01)

0 0.01 0.02 0.03 0.04 0.05 0.06

0 2 4 6 8 10 12 14

Euler

0 2e-07 4e-07 6e-07 8e-07 1e-06 1.2e-06 1.4e-06 1.6e-06 1.8e-06 2e-06

0 2 4 6 8 10 12 14

i-Euler

前進オイラー法 改良オイラー法

-3.5e-07 -3e-07 -2.5e-07 -2e-07 -1.5e-07 -1e-07 -5e-08 0

0 2 4 6 8 10 12 14

Heun 3

-1.5e-11 -1e-11 -5e-12 0 5e-12 1e-11

0 2 4 6 8 10 12 14

Rk

ホインの３次公式 ルンゲクッタ

● 実習内容

• (5)から(8)を改良オイラー法,ホインの３次公式, ルンゲクッタ法を用いて解きなさい.