● 講義資料

(1)

● 講義資料

▼ 講義予定

•

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

–

テイラー展開法

–

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

– Runge-Kutta Method

の安定性

● 前回の講義のまとめ

★ 前進オイラー法

•

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

x n+1 = x n + hf(t n , x n )

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

•

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

dx

dt

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

dx

dt (t) = lim

h → 0

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

を

t = t n

で考えれば,

h > 0

が小さいと仮定して,

dx

dt (t n ) = lim

h → 0

x(t n + h) − x(t n )

h ∼ x(t n+1 ) − x(t n ) h

と書き直すことができる. （これを前進差分と呼ぶ.）これを, 微分方程式

x ^′ = f (t, x)

に適用すれば,

x(t n+1 ) − x(t n )

h ∼ f (t n , x n ) x(t n+1 ) ∼ x(t n ) + hf(t n , x(t n ))

となる. したがって,前進オイラー法の計算スキーム

x n+1 = x n + hf(t n , x n )

に相当する式を得ることができる.

•

前進オイラー法のスキームの右辺は,直前の

x n

の値があれば直接計算可能である. その意味で陽的解法

(explicit method)

となっている.

•

前進オイラー法のスキームは,直前の

x n

の値のみで,次の

x n+1

を計算可能である. その意味で１段階法となっている.

(2)

★ 後退オイラー法

•

前進オイラー法での差分化を行った点を

t = t n

ではなく

t = t n+1

での値と考える. すなわち

dx

dt (t n+1 ) = lim

h → 0

x(t n + h) − x(t n )

h ∼ x(t n+1 ) − x(t n ) h

と書き直すことができる. これを,微分方程式

x ^′ = f (t, x)

に適用すれば,

x(t n+1 ) − x(t n )

h ∼ f (t n+1 , x n+1 ) x(t n+1 ) ∼ x(t n ) + hf (t n+1 , x(t n+1 ))

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

x n +1 = x n + hf(t n +1 , x n +1 )

を後退オイラー法

(backward Euler method)

と呼ぶ.

•

後退オイラー法のスキームは,直前の

x n

の値のみで,次の

x n+1

を計算可能である. その意味で１段階法となっている.

•

後退オイラー法のスキームの右辺は,直前の

x n

の値だけからは直接計算できない. その意味で陰的解法

(implicit method)

となっている.

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

–

微分方程式

x ^′ (t) = x(t)

に後退オイラー法を適用すると,

x n+1 = x n + hx n+1

を得る. ここで,

x n

から

x n+1

を得るためには, この式を代数的に

x n+1

について解いて,

x n +1 = 1 1 − h x n

として計算する必要がある.

–

微分方程式

x ^′ (t) = sin(x(t))

に後退オイラー法を適用すると,

x n+1 = x n + h sin(x n+1 )

を得る. ここで,

x n

から

x n+1

を得るために, この式を代数的に

x n+1

について解くことはできない. したがって,逐次近似

y 0 = x n ,

y k+1 = x n + h sin(y k )

を用いて,その極限として

x n+1

を求める必要がある. （なお,

h

が十分に小さければ, この逐次近似が収束することを証明できる.）

(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 n+1 ) = x(t n ) + hx ^′ (t n ) + h ²

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

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

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

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

2 f ^′ (x(t n ))f (x(t n )) + O(h ³ )

と計算できる. ここで,

f

がリプシッツ条件

|f (p) − f (q)| ≤ L|p − q|

をみたし,

f

はなめらかであると仮定すると,

|x(t n +1 ) − x n +1 | ≤ (1 + Lh)|x(t n ) − x n | + M h ²

が成り立つ.

•

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

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

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

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

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

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

2 f ^′ (x(t n+1 ))f (x(t n+1 )) + O(h ³ )

より,

|x(t n+1 ) − x n+1 | ≤ (1 − Lh) ⁻ ¹ |x(t n ) − x n | + M h ²

が成り立つ.

•

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

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

の形の不等式を得ることができた. この評価を局所離散化誤差と呼び,ある時刻

t = t n

において,数値解の値

x n

と厳密解の値

x(t n )

が一致しているとき, １ステップ進んだ際の誤差が

O(h ² )

であることを示している.

(4)

•

この局所離散化誤差の不等式から

t = T = N h

においての誤差

|x(T ) − x N | ≤ Che ^{T L} ,

前進オイラー法,

|x(T ) − x N | ≤ Ch,

後退オイラー法,

を得ることができる. すなわち,オイラー法は誤差が

h

に関して１次の公式であることがわかる.

★ オイラー法と解の存在定理

•

通常,１階常微分方程式の初期値問題の解の存在定理の証明は,

x 0 (t) = x 0 , x n+1 (t) = x 0 +

Z t 0

f (x n (s)) ds

によって定義される関数列

{x n }

が,区間上で一様収束し,その極限

x

（これが連続関数になる）が

x(t) = x 0 + Z t

0

f (x(s)) ds

をみたすことを示す. （一様収束しないと,極限がこの積分方程式を満たすことを示せない.

また,この証明の中で,

f

がリプシッツ条件を満たすことが本質的に使われる.）

なお,

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

でこの近似列を計算すると

x n (t) = X n

k=0

t ⁿ n!

である.

•

一方,前進オイラー法の誤差に関する不等式

|x(T ) − x N | ≤ Che ^{T L}

は,

h → 0

としたとき,

x N

が

x(T )

に各点ごとに収束することを示している. 単純にこの不等式を見ているだけでは, 各点収束しか示していないので, 数値解

(h

に依存している)の構成から解の存在を言うことはできないが,より詳しく計算すれば,

f

が有界の時,

h k → 0

となる列をうまくとれば,数値解

{x ^h _n

^k

}

が一様収束する部分列を取り出すことができる. さらに,

f

がリプシッツ条件をみたせば,収束先は一意的となる. よって,オイラー法によって１階常微分方程式の初期値問題の解の存在定理を示すことが可能であることに注意する必要がある.

なお,

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

の時には,オイラー法で構成される数値解は

x n = (1 + h) ⁿ

であり,

t = nh

と考えれば,

x n = (1 + h) ^t/h

となるので,

h → 0

の時,

x n → x(t) = e ^t

が成り立つ.

(5)

● 講義資料

▼ 改良オイラー法による数値計算

以下の図は,改良オイラー法で種々の常微分方程式の初期値問題の数値解を計算した例である.

（特に書いていない限り

h = 0.01

である.）

0 0.5 1 1.5 2 2.5 3

0 0.2 0.4 0.6 0.8 1

x’ = x

0 0.2 0.4 0.6 0.8 1 1.2 1.4

0 0.2 0.4 0.6 0.8 1

x’ = -x

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

-0.2 0 0.2 0.4 0.6 0.8 1 1.2

0 0.5 1 1.5 2 2.5 3 3.5 4 x’ = (1-x^2)

0 0.5 1 1.5 2 2.5 3 3.5 4

0 2 4 6 8 10

x’ = sin(x)

x ^′ (t) = (1 − x(t) ² ), x(0) = 0 x ^′ (t) = sin(x(t)), x(0) = 1

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

0 2 4 6 8 10 12 14 16 x’’ = -x

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 x’’ = -x

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

(6)

-3 -2 -1 0 1 2 3

0 2 4 6 8 10 12 14 16

x’’ = -sin(x)

-2 -1 0 1 2

x’’ = -sin(x)

x ^′′ (t) = − sin(x(t)), x(0) = 2, x ^′ (0) = 0

-3 -2 -1 0 1 2 3

0 5 10 15 20 25

x’’ = (1-x^2)x’ - x

-3 -2 -1 0 1 2 3

-3 -2 -1 0 1 2 3 x’’ = (1-x^2)x’ - x

x ^′′ (t) = (1 − x(t) ² )x ^′ (t) − x(t), x(0) = 1.5, x ^′ (0) = 1

-10 -5 0 5 10

X’’ = -X/|X|^3

X ^′′ (t) = −X(t)/|X (t)| ³ , X (0) = (1, 1, 0), X ^′ (0) = (1, 0, 0)

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2

0 2 4 6 8 10

Double Pendulum

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 Double Pendulum

二重振り子,

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

(7)

▼ ホインの３次公式による数値計算

以下の図は,ホインの３次公式で種々の常微分方程式の初期値問題の数値解を計算した例である.

h = 0.01

である.）

0 0.5 1 1.5 2 2.5 3

0 0.2 0.4 0.6 0.8 1

x’ = x

0 0.2 0.4 0.6 0.8 1 1.2 1.4

0 0.2 0.4 0.6 0.8 1

x’ = -x

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

-0.2 0 0.2 0.4 0.6 0.8 1 1.2

0 0.5 1 1.5 2 2.5 3 3.5 4 x’ = (1-x^2)

0 0.5 1 1.5 2 2.5 3 3.5 4

0 2 4 6 8 10

x’ = sin(x)

x ^′ (t) = (1 − x(t) ² ), x(0) = 0 x ^′ (t) = sin(x(t)), x(0) = 1

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

0 2 4 6 8 10 12 14 16 x’’ = -x

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 x’’ = -x

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

(8)

-3 -2 -1 0 1 2 3

0 2 4 6 8 10 12 14 16

x’’ = -sin(x)

-2 -1 0 1 2

x’’ = -sin(x)

x ^′′ (t) = − sin(x(t)), x(0) = 2, x ^′ (0) = 0

-3 -2 -1 0 1 2 3

0 5 10 15 20 25

x’’ = (1-x^2)x’ - x

-3 -2 -1 0 1 2 3

-3 -2 -1 0 1 2 3 x’’ = (1-x^2)x’ - x

x ^′′ (t) = (1 − x(t) ² )x ^′ (t) − x(t), x(0) = 1.5, x ^′ (0) = 1

-10 -5 0 5 10

X’’ = -X/|X|^3

X ^′′ (t) = −X(t)/|X (t)| ³ , X (0) = (1, 1, 0), X ^′ (0) = (1, 0, 0)

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2

0 2 4 6 8 10

Double Pendulum

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 Double Pendulum

二重振り子,

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

(9)

▼ ルンゲ・クッタ法による数値計算

以下の図は,ルンゲ・クッタ法で種々の常微分方程式の初期値問題の数値解を計算した例である.

h = 0.01

である.）

0 0.5 1 1.5 2 2.5 3

0 0.2 0.4 0.6 0.8 1

x’ = x

0 0.2 0.4 0.6 0.8 1 1.2 1.4

0 0.2 0.4 0.6 0.8 1

x’ = -x

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

-0.2 0 0.2 0.4 0.6 0.8 1 1.2

0 0.5 1 1.5 2 2.5 3 3.5 4 x’ = (1-x^2)

0 0.5 1 1.5 2 2.5 3 3.5 4

0 2 4 6 8 10

x’ = sin(x)

x ^′ (t) = (1 − x(t) ² ), x(0) = 0 x ^′ (t) = sin(x(t)), x(0) = 1

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

0 2 4 6 8 10 12 14 16 x’’ = -x

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 x’’ = -x

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

(10)

-3 -2 -1 0 1 2 3

0 2 4 6 8 10 12 14 16

x’’ = -sin(x)

-2 -1 0 1 2

x’’ = -sin(x)

x ^′′ (t) = − sin(x(t)), x(0) = 2, x ^′ (0) = 0

-3 -2 -1 0 1 2 3

0 5 10 15 20 25

x’’ = (1-x^2)x’ - x

-3 -2 -1 0 1 2 3

-3 -2 -1 0 1 2 3 x’’ = (1-x^2)x’ - x

x ^′′ (t) = (1 − x(t) ² )x ^′ (t) − x(t), x(0) = 1.5, x ^′ (0) = 1

-10 -5 0 5 10

X’’ = -X/|X|^3

X ^′′ (t) = −X(t)/|X (t)| ³ , X (0) = (1, 1, 0), X ^′ (0) = (1, 0, 0)

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2

0 2 4 6 8 10

Double Pendulum

-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2 Double Pendulum

二重振り子,

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

(11)

▼ 種々の方法による数値計算

以下の図は,種々のルンゲ・クッタ型公式で

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

を計算し,

t = 1

に相当する値を真の解の

t = 1

での値とを比較したものである. （相対誤差を表示している.）

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

100 1000 10000 100000 1e+06 1e+07 1e+08 1e+09

x’ = x, relative error at t = 1 Forward Euler Backward Euler Improved Euler Heun 2 Heun 3 Kutta 3 Runge-Kutta Runge-Kutta-Gill

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

100 1000 10000 100000 1e+06 1e+07 1e+08 1e+09

x’ = -x, relative error at t = 1 Forward Euler Backward Euler Improved Euler Heun 2 Heun 3 Kutta 3 Runge-Kutta Runge-Kutta-Gill

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

100 1000 10000 100000 1e+06 1e+07 1e+08 1e+09

x’ = x, relative error at t = 1 Forward Euler Backward Euler Improved Euler Heun 2 Heun 3 Kutta 3 Runge-Kutta Runge-Kutta-Gill

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

100 1000 10000 100000 1e+06 1e+07 1e+08 1e+09

x’ = -x, relative error at t = 1 Forward Euler Backward Euler Improved Euler Heun 2 Heun 3 Kutta 3 Runge-Kutta Runge-Kutta-Gill

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

•

誤差がそれぞれの収束次数を応じた傾きで推移していることに注意する必要がある.

•

ある程度以上

h

を小さくとってしまうと,誤差が必ずしも小さくならないことに注意が必要である.

● 実習内容

1.

６回目資料の

1

から

9

の常微分方程式の初期値問題の数値解を,改良オイラー法,ホインの２次公式,ホインの３次公式,ルンゲ・クッタ法で計算しなさい.

● 講義資料

•

–

–

– Runge-Kutta Method

– Runge-Kutta Method

•

x n+1 = x n + hf(t n , x n )

•

dx

dt

dx

dt (t) = lim

h → 0

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

t = t n

h > 0

dx

dt (t n ) = lim

h → 0

x(t n + h) − x(t n )

h ∼ x(t n+1 ) − x(t n ) h

x ′ = f (t, x)

x(t n+1 ) − x(t n )

h ∼ f (t n , x n ) x(t n+1 ) ∼ x(t n ) + hf(t n , x(t n ))

x n+1 = x n + hf(t n , x n )

•

x n

(explicit method)

•

x n

x n+1

•

t = t n

t = t n+1

dx

dt (t n+1 ) = lim

h → 0

x(t n + h) − x(t n )

h ∼ x(t n+1 ) − x(t n ) h

x ′ = f (t, x)

x(t n+1 ) − x(t n )

h ∼ f (t n+1 , x n+1 ) x(t n+1 ) ∼ x(t n ) + hf (t n+1 , x(t n+1 ))

x n +1 = x n + hf(t n +1 , x n +1 )

(backward Euler method)

•

x n

x n+1

•

x n

(implicit method)

–

x ′ (t) = x(t)

x n+1 = x n + hx n+1

x n

x n+1

x n+1

x n +1 = 1 1 − h x n

–

x ′ (t) = sin(x(t))

x n+1 = x n + h sin(x n+1 )

x n

x n+1

x n+1

y 0 = x n ,

y k+1 = x n + h sin(y k )

x n+1

h

•

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 n+1 ) = x(t n ) + hx ′ (t n ) + h 2

2 x ′′ (t n ) + O(h 3 ),

= x(t n ) + hf (x(t n )) + h 2

2 f (x ′ (t n ))f (x(t n )) + O(h 3 ), x n+1 = x n + hf (x n )

x ^′ = f (t, x)

x ^′ = f (t, x)

x ^′ (t) = x(t)

x ^′ (t) = sin(x(t))

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

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

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

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

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

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

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

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

2 f ^′ (x(t n ))f (x(t n )) + O(h ³ )

|x(t n +1 ) − x n +1 | ≤ (1 + Lh)|x(t n ) − x n | + M h ²

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

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

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

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

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

2 f ^′ (x(t n+1 ))f (x(t n+1 )) + O(h ³ )

|x(t n+1 ) − x n+1 | ≤ (1 − Lh) ⁻ ¹ |x(t n ) − x n | + M h ²

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

O(h ² )

|x(T ) − x N | ≤ Che ^{T L} ,

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

t ⁿ n!

|x(T ) − x N | ≤ Che ^{T L}

{x ^h _n

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

x n = (1 + h) ⁿ

x n = (1 + h) ^t/h

x n → x(t) = e ^t

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

x ^′ (t) = (1 − x(t) ² ), x(0) = 0 x ^′ (t) = sin(x(t)), x(0) = 1