● 前回の講義のまとめ常微分方程式の初期値問題

(1)

● 前回の講義のまとめ

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

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

(1)

の数値解を構成する計算スキームを考える. 以下では,

f (t, p)

は滑らかであると仮定する.

★ 後退オイラー法

•

差分化を

dx

dt (t k ) ←→ x k − x k+1

− h

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

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

によって

{ x k }

を与える計算スキームを得ることができる. これを後退オイラー法と呼ぶ.

•

後退オイラー法に対しても,

p

に関するリプシッツ条件と

^∂f _∂t , ^∂f _∂p

の有界性があれば,局所離散化誤差

e k

と, 大域的離散化誤差

E k = | x(t k ) − x k |

は

e k = O(h ² ), E k = O(h)

をみたす.

•

一般に数値計算スキームにおいて,

x k

（または

{ x j } ^k j=0

）から

x k+1

を与える計算が,何らかの方程式を解くことなく行えるとき,陽的方法とよび,そうでないとき陰的解法と呼ぶ.

•

後退オイラー法は陰的解法の代表例であり,

x k

から

x k+1

を計算する際には,何らかの代数方程式を解かなければならない.

【例１】

x ^′ (t) = λx(t)

の時,後退オイラー法の計算式は

x k+1 = x k + λhx k+1

となる. したがって, これを

x k+1

について解いて,

x k+1 = 1

(1 − λh) x k

とすればよい.

【例２】

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

を考える. これは,

y ^′ (t) = − sin(x(t)), x ^′ (t) = y(t)

と書き換えて考えると,

y k+1 = y k − h sin(x k+1 ), x k+1 = x k + hy k+1

を解くことになる. この式を

X k = (x k , y k )

が与えられたとき,

X k+1 = (x k+1 , y k+1 )

を求める問題と考え,

Φ(x, y) = (x k + hy k+1 , y k − h sin(x))

と考え,

X = Φ(X)

を逐次近似で求めると考えればよい.

または, この場合には,上の式から

x k+1

を消去して,

y k+1 = y k − h sin(x k + hy k+1 )

を解くと考えて,

Φ(y) = y k − h sin(x k + hy)

とおき,

y = Φ(y)

を逐次近似で求めると考えてもよい.

(2)

•

このように,後退オイラー法は,前進オイラー法と比較して,実際の計算が複雑となっているが,多くの後退解法は,それに対応する前進解法と比較して,求めた数値解の性質がよいことが多い.

★ テイラー展開法

•

オイラー法よりも高次の近似を実現する計算スキームとして,

x k+1 = x k + hf (t k , x k ) + h ²

2! F 2 (t k , x k ) + · · · + h ^p

p! F p (t k , x k )

によって

{ x k }

を与えるテイラー法を考えることができる. ここで,

F j (t, x)

は,

^d _dt

^j−1j

f (t, x(t))

の

x ^′ (t)

などに,微分方程式の条件

x ^′ = f (t, x)

を代入したものである.

•

テイラー法は,

x(t)

のテイラー展開

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

2! x ^′′ (t k ) + · · · h ^p

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

を利用した近似スキームである. （オイラー法は,テイラー法の

p = 1

に相当するスキームである.）

したがって,

f (t, p)

に必要なだけの滑らかさがあれば,

e k = O(h ^p+1 ), E k = O(h ^p )

が成り立つ.

•

しかし,テイラー法を実際に適用するには

f (t, p)

の高次導関数を計算する必要があるため, 実用上は用いられない.

★ 修正オイラー法

•

オイラー法で用いた

x ^′ (t)

の近似である前進差分

x ^′ (t k ) ∼ x(t k+1 ) − x(t k ) h

だけでなく,後退差分

x ^′ (t k ) ∼ x(t k− 1 ) − x(t k )

− h

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

x ^′ (t k ) ∼ x(t k+1 ) − x(t k− 1 ) 2h

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

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

2 x ^′′ (t k ) + h ³

3! x ⁽³⁾ (t k ) + O(h ⁴ ), x(t k− 1 ) = x(t k ) − hx ^′ (t k ) + h ²

2 x ^′′ (t k ) − h ³

3! x ⁽³⁾ (t k ) + O(h ⁴ ),

の上の式から下を引けば,

x(t k+1 ) − x(t k− 1 ) = 2hx ^′ (t k ) + h ³

3 x ⁽³⁾ (t k ) + O(h ⁵ )

(3)

となる. この

h ³

以上の項を無視すれば,

x k+1 = x k− 1 + 2hf (t k , x k )

というスキームを得る. これを修正オイラー法と呼ぶ.

•

修正オイラー法は,初期値

x 0

だけではなく,

x 1

（出発値と呼ぶ）を必要とする. 通常は

x 1

を求めるためにオイラー法を用いる.

•

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

e k = O(h ³ ), E k = O(h ² )

が成り立つ.

•

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

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

(2)

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

x k+1 = x k − 1 − 2hx k (3)

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

t ² + 2ht − 1 = 0

の２つの解

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

h ² + 1

を用いて,

x k = c 1 η ₁ ^k + c 2 η ₂ ^k (4)

と書くことができる. いま,

η 2 < − 1

より,

c 2 6 = 0

ならば,

| x k | → ∞

となり, (1)の解の挙動

x(t) → 0

とは大きく異なることになる.

ここで,初期条件

x 0

と出発値

x 1

を用いて,

c 1 , c 2

は

c 2 = x 1 − x 0

η 2 − η 1

= x 0 − x 1

2 √ h ² + 1 , c 1 = x 0 − c 1

と書け,

x 0 6 = x 1

であるかぎりは,任意の

h > 0

に対して

c 2 6 = 0

である. したがって

(4)

より,

c 2 6 = 0

である限りは,

k

が大きいときに

x k

の挙動は

η ₂ ^k

に支配されることになり,

| x k | → ∞

かつ

x k

は振動する.

•

また,

c 2 = 0

となる初期値を与えたときにも, 実際に計算機内部で計算される初期値, その計算誤差から,

c 2 = ǫ

に対応するものとなり,

c 2 6 = 0

に対応する初期値を与えたときと同様の振る舞いを示す.

•

修正オイラー法を用いて

(3)

の数値解を構成すると,真の解とは全く異なる振動しながら「発散」する解があらわれる. このような現象を数値的不安定性と呼び,本来「安定」な数学的対象を「離散化」または「数値化」したことにより発生する不安定性である. 最も安定な対象であると思われる

x ^′ (t) = − x(t)

でさえも不安定となってしまうため,修正オイラー法も実用上は用いられない.

(4)

•

このような数値的不安定な現象は, ３項間（または,それ以上の）線形漸化式にしたがって計算を行う際にもあらわれる. その代表例がフィボナッチ数列

a n+1 = a n + a n − 1

である. これは, その特性方程式

t ² − t − 1 = 0

の２つの解

α = 1 + √

5 2 , β = 1 − √ 5

2

が

α > 1, − 1 < β < 0

となることに起因する. この場合,一般項は

a n = C 1 α ⁿ + C 2 β ⁿ

と書けるが,

C 1 = 0

となる初期条件は

a 1 /a 0 = β

となり,実際の計算機内部では

C 1 = ǫ

と与えたことに他ならない. したがって,

a n

の漸近的な挙動は

a n ∼ α ⁿ

となってしまう.

● 講義資料

▼ 講義予定

•

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

–

改良オイラー法とルンゲ・クッタ法

–

一般的なルンゲ・クッタ型公式

–

ルンゲ・クッタ型公式と数値積分の関係

● 講義資料

★ 改良オイラー法

• x ^′ = x, x 0 = 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, x 0 = 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"

(5)

★ ルンゲ・クッタ型公式

• 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

(6)

• x ^′′ = − sin(x),

前進オイラー法, 改良オイラー法, ホインの３次公式,ルンゲ・クッタ法で

t = 4π

での

¹ ₂ (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

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

(7)

• 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

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

● 実習内容

x ^′ (t) = λx(t), x(0) = x 0 , λ 6 = 0. (5)

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

x ^′′ (t) + ax ^′ (t) + bx(t) = f (x), x(0) = x 0 , x ^′ (0) = v 0 . (7) x ^′′ (t) + sin(x(t)) = 0, x(0) = x 0 , x ^′ (0) = v 0 . (8) x ^′′ (t) + (x(t) ² − 1)x ^′ (t) + x(t) = 0, x(0) = x 0 , x ^′ (0) = v 0 . (9) 1.

上記の常微分方程式の数値解を改良オイラー法, ホインの３次公式,ルンゲ・クッタ法を用いて構成

しなさい. また, これらのうち２階常微分方程式については,その解を相平面上に図示しなさい.

2.

適切な常微分方程式の初期値問題を考えることにより

log 2

の近似値を求めなさい.

● 前回の講義のまとめ 常微分方程式の初期値問題

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

(1)

f (t, p)

•

dx

dt (t k ) ←→ x k − x k+1

− h

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

{ x k }

•

p

∂f ∂t , ∂f ∂p

e k

E k = | x(t k ) − x k |

e k = O(h 2 ), E k = O(h)

•

x k

{ x j } k j=0

x k+1

•

x k

x k+1

x ′ (t) = λx(t)

x k+1 = x k + λhx k+1

x k+1

x k+1 = 1

(1 − λh) x k

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

y ′ (t) = − sin(x(t)), x ′ (t) = y(t)

y k+1 = y k − h sin(x k+1 ), x k+1 = x k + hy k+1

X k = (x k , y k )

X k+1 = (x k+1 , y k+1 )

Φ(x, y) = (x k + hy k+1 , y k − h sin(x))

X = Φ(X)

x k+1

y k+1 = y k − h sin(x k + hy k+1 )

Φ(y) = y k − h sin(x k + hy)

y = Φ(y)

•

•

x k+1 = x k + hf (t k , x k ) + h 2

2! F 2 (t k , x k ) + · · · + h p

p! F p (t k , x k )

{ x k }

F j (t, x)

d dt

f (t, x(t))

x ′ (t)

x ′ = f (t, x)

•

x(t)

x(t k+1 ) = x(t k ) + hx ′ (t k ) + h 2

2! x ′′ (t k ) + · · · h p

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

p = 1

f (t, p)

e k = O(h p+1 ), E k = O(h p )

•

f (t, p)

•

x ′ (t)

x ′ (t k ) ∼ x(t k+1 ) − x(t k ) h

x ′ (t k ) ∼ x(t k− 1 ) − x(t k )

− h

x ′ (t k ) ∼ x(t k+1 ) − x(t k− 1 ) 2h

x(t k+1 ) = x(t k ) + hx ′ (t k ) + h 2

2 x ′′ (t k ) + h 3

3! x (3) (t k ) + O(h 4 ), x(t k− 1 ) = x(t k ) − hx ′ (t k ) + h 2

2 x ′′ (t k ) − h 3

3! x (3) (t k ) + O(h 4 ),

x(t k+1 ) − x(t k− 1 ) = 2hx ′ (t k ) + h 3

3 x (3) (t k ) + O(h 5 )

h 3

x k+1 = x k− 1 + 2hf (t k , x k )

•

x 0

x 1

x 1

•

● 前回の講義のまとめ常微分方程式の初期値問題

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

^∂f _∂t , ^∂f _∂p

e k = O(h ² ), E k = O(h)

{ x j } ^k j=0

x ^′ (t) = λx(t)

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

y ^′ (t) = − sin(x(t)), x ^′ (t) = y(t)

x k+1 = x k + hf (t k , x k ) + h ²

2! F 2 (t k , x k ) + · · · + h ^p

^d _dt

x ^′ (t)

x ^′ = f (t, x)

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

2! x ^′′ (t k ) + · · · h ^p

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

e k = O(h ^p+1 ), E k = O(h ^p )

x ^′ (t)

x ^′ (t k ) ∼ x(t k+1 ) − x(t k ) h

x ^′ (t k ) ∼ x(t k− 1 ) − x(t k )

x ^′ (t k ) ∼ x(t k+1 ) − x(t k− 1 ) 2h

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

2 x ^′′ (t k ) + h ³

3! x ⁽³⁾ (t k ) + O(h ⁴ ), x(t k− 1 ) = x(t k ) − hx ^′ (t k ) + h ²

2 x ^′′ (t k ) − h ³

3! x ⁽³⁾ (t k ) + O(h ⁴ ),

x(t k+1 ) − x(t k− 1 ) = 2hx ^′ (t k ) + h ³

3 x ⁽³⁾ (t k ) + O(h ⁵ )

h ³

e k = O(h ³ ), E k = O(h ² )

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

t ² + 2ht − 1 = 0

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

h ² + 1

x k = c 1 η ₁ ^k + c 2 η ₂ ^k (4)

2 √ h ² + 1 , c 1 = x 0 − c 1

η ₂ ^k

x ^′ (t) = − x(t)

t ² − t − 1 = 0

a n = C 1 α ⁿ + C 2 β ⁿ

a n ∼ α ⁿ

• x ^′ = x, x 0 = 1.0

• x ^′ = − x, x 0 = 1.0

• x ^′ = x,

• x ^′′ = − x,