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

● 前回の講義のまとめ

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

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

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

★ ルンゲ・クッタ型公式

• 常微分方程式の初期値問題(1) の数値解を,オイラー法よりも高次の公式を得る方法として, 以下の 方法を考える.

ki=f(tn+hci, xn+h

s

X

j=1

aijkj), i= 1, . . . , s,

xn+1=xn+h

s

X

i=1

biki.

(2)

ここで,{aij}^s_i,j=1,{bi}^s_i=1,{ci}^s_i=1は,指定の次数の公式となるような条件の下で適切にきめる定数 である. この形の公式をルンゲ・クッタ型公式と呼び,sをその段数と呼ぶ.

• いま,i≤j の時に aij= 0 を仮定すると,{ki}^s_i=1 は, k1 から順に計算が可能である. この条件を満 たすとき,陽的ルンゲ・クッタ型公式と呼ぶ.

なお,i < j の時に aij = 0となる公式を,半陰的ルンゲ・クッタ型公式,このような条件がない公式 を陰的ルンゲ・クッタ型公式と呼ぶ.

• 前進オイラー法は「１段陽的ルンゲ・クッタ型公式」と考えられ,後退オイラー法は「１段陰的ルン ゲ・クッタ型公式」と考えられる.

【前進オイラー法】 s= 0,a11= 0,b1= 1,c1= 0とおくと, k1=f(tn, xn),

xn+1=xn+hk1=xn+hf(tn, xn) となり, 前進オイラー法を得ることができる.

【後退オイラー法】 s= 0,a11= 1,b1= 1,c1= 1とおくと, k1=f(tn+h, xn+hk1),

xn+1=xn+hk1=xn+hf(tn+1, xn+1) となり, 後退オイラー法を得ることができる.

• ルンゲ・クッタ型公式で得た数値解{xn}がtn+1≤T に対して,

|xn+1−x(tn+1)| ≤Ch,T|xn−x(tn)|+Ch^p+1+O(h^p+2)

を満たすとき,その公式をp次公式であるという. ただし, Ch,T は hについてp次以下であるとす る. この意味で, 前進オイラー法および後退オイラー法は「1段 1次」の公式である.

★ ルンゲ・クッタ型公式の次数条件

• s段のルンゲ・クッタ型公式がp次あるための条件を求める. 簡単のため, x^′(t) =f(x(t)),

x(0) =x0

(3)

の形の場合を考える.

• (1)の代わりに,補助的にτ(t) =tを導入し,

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

x(0) =x0, τ(0) = 0

と考え,X(t) = (τ(t), x(t))についての微分方程式と考えれば, (1)は(3)の形に帰着される. この時, ルンゲ・クッタ型スキームは,

ℓi = 1, i= 1, . . . , s, ki =f(tn+h

s

X

j=1

aijℓi, xn+h

s

X

j=1

aijkj), i= 1, . . . , s,

xn+1 =xn+h

s

X

i=1

biki,

tn+1 =tn+h

s

X

i=1

biℓi, と書ける. したがって,

s

X

j=1

aij=cj, (4)

s

X

j=1

bi= 1 (5)

を仮定することにより, (3) についてのルンゲ・クッタ型公式から, (1)についてのルンゲ・クッタ型 公式を導出することができる. よって,以下では(4), (5)を仮定して,次数条件を計算する.

• はじめに,

Xi=xn+hX

j

aijkj

とおく. この時,

d

dhXi^α=X

j

aijkj^α+hX

j

aij

d dhk^αj, d

dhk^αj =X

β

fβ^α

d dhX_j^β,

が成り立つ. したがって,

Xj^α

h=0=x(tn)^α, k^αj

h=0=f^α, d

dhX_i^α _h=0

=X

j

aijf^α=cif^α, d

dhk^αj

_h=0=X

β

fβ^αcif^β =ci

X

β

fβ^αf^β, がなりたつ. 一方,ki をhについてテイラー展開すれば,

k^α_i = k_i^α|_h=0+h d dhk_i^α

_h=0

+O(h²) となるので,

x^α_n+1=x^α_n+hX

i

bi

k^α_i|_h=0+h d dhk^α_i

_h=0

+O(h³)

=x(tn)^α+h X

i

bi

!

f^α+h² X

i

bici

! X

β

fβ^αf^β+O(h³)

(6)

が成り立つ. 一方,x(t)をテイラー展開すれば,

x(tn+1)^α=x(tn)^α+hf^α+h² 2

X

β

fβ^αf^β+O(h³) (7)

となる. よって, (6) と(7)から,

x(tn+1)^α−x^α_n+1=x(tn)^α−x^α_n +h f^α(x(tn))−X

i

bif^α(xn)

!

+h² 2



 X

β

(fβ^αf^β)(x(tn))−2X

i

bici(fβ^αf^β)(xn)



+O(h³)

(8)

が成り立つ.

• このことから,次を示すことができる. ただし,以下ではf およびその高階微分について,リプシッツ 条件を仮定する.

【１次の条件式】 (5)が成り立つならば,

|x(tn+1)^α−x^αn+1| ≤(1 +Lh)|x(tn)^α−x^αn|+O(h²)

が成り立つ. したがって, (5)は,ルンゲ・クッタ型公式が 1 次公式となるための必要十分条件 である.

【２次の条件式】

s

X

i=1

bici= 1

2 (9)

が成り立つならば,

|x(tn+1)^α−x^α_n+1| ≤

1 +Lh+Lh² 2

|x(tn)^α−x^α_n|+O(h³)

が成り立つ. したがって, (5)と (9)は, ルンゲ・クッタ型公式が2 次公式となるための必要十 分条件である.

★ 具体的なルンゲ・クッタ型公式

• 以下では,ルンゲ・クッタ型公式の係数を

c1 a11 · · · a1s

... ... ... cs as1 · · · ass

b1 · · · bs

と表記する. これをButcher 行列と呼ぶ.

• 代表的なルンゲ・クッタ型公式は以下の通りである.

1段 1次公式

前進オイラー法 後退オイラー法 0

1

1 1

1

2段 2次公式

改良オイラー法 ホインの２次公式 0

1/2 1/2

0 1

0

1 1

1/2 1/2

3段 3次公式

ホインの３次公式 クッタの３次公式 0

1/3 1/3

2/3 0 2/3

1/4 0 3/4

0 1/2 1/2

1 −1 2

1/6 2/3 1/6 4段 4次公式

ルンゲ・クッタの公式 ルンゲ・クッタ・ジルの公式 0

1/2 1/2

1/2 0 1/2

1 0 0 1

1/6 1/3 1/3 1/6

0

1/2 1/2

1/2 √¹

2−¹₂ 1−√¹ 2

1 0 −^√1₂ 1 +^√1₂

1/6 ¹₃(1−√¹

2) ¹₃(1 +√¹

2) 1/6

• これらの中でも特に有名な公式がルンゲ・クッタ法（ルンゲ・クッタの公式）であり,それを標準的 な記法で書けば

k1=f(tn, xn),

k2=f(tn+h/2, xn+hk1/2), k3=f(tn+h/2, xn+hk2/2), k4=f(tn+h, xn+hk3), xn+1=xn+h

6(k1+ 2k2+ 2k3+k4) となる.

• なお, 5 段以上のルンゲ・クッタ型公式については, 段数と等しい次数の公式を得ることができない ことが知られている.

★ ルンゲ・クッタ型公式と数値積分

• ルンゲ・クッタ法を

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

x(0) = 0 (10)

に適用することを考える. この時, (10)の解は x(t) =

Z t 0

f(s)ds (11)

であるので, (10)の数値解を構成することは,積分(11)の近似値を求めることに他ならない.

• 前進オイラー法を(10)に適用すると,

xn+1=xn+hf(tn) となり,これは,

Z tn+1

tn

f(s)ds∼hf(tn)

となる近似である. すなわち,積分を「分割区間の端点の高さ×分割幅」の長方形で近似したことに なる.

• ホインの２次公式を(10)に適用すると,

xn+1 =xn+h

2(f(tn) +f(tn+1)) となり,これは,

Z tn+1

tn

f(s)ds∼ h

2(f(tn) +f(tn+1))

となる近似である. すなわち,積分を「分割区間の両方の端点の高さの平均×分割幅」で近似したこ とになり,これは, 4点(tn,0), (tn, f(tn)), (tn+1, f(tn+1)), (tn+1,0)でつくられる台形の面積に等し い. 積分をこのように近似して計算する方法を台形公式と呼ぶ.

• 台形公式は,区間[0, h]上の関数f(t)を, (0, f(0)), (h, f(h))を通る直線（１次関数）f1(t)で近似し, Z h

0

f(s)ds∼ Z h

0

f1(s)ds と近似したことに他ならない.

● 講義資料

▼ 講義予定

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

– ルンゲ・クッタ型公式の安定性 – シンプレクティック解法

● 講義資料

★ ルンゲクッタ型公式の安定性

• x^′=−x,x(0) = 1.0に対して,hを動かしたときの例.

-0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

0 5 10 15 20 25 30 35 40 45 50

Euler Exp(-t)

0 0.2 0.4 0.6 0.8 1

0 5 10 15 20 25 30 35 40 45 50

I Euler Exp(-t)

オイラー法(t= 1.5) 改良オイラー法(t= 1.5)

0 0.2 0.4 0.6 0.8 1

0 5 10 15 20 25 30 35 40 45 50

Heun 3 Exp(-t)

0 0.2 0.4 0.6 0.8 1

0 5 10 15 20 25 30 35 40 45 50

RK Exp(-t)

ホインの３次公式(t= 1.5) ルンゲクッタ法(t= 1.5)

-3000 -2000 -1000 0 1000 2000 3000 4000

0 5 10 15 20 25 30 35 40 45 50

Euler Exp(-t)

0 2000 4000 6000 8000 10000 12000 14000 16000 18000

0 5 10 15 20 25 30 35 40 45 50

I Euler Exp(-t)

オイラー法(t= 2.5) 改良オイラー法(t= 2.5)

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

0 5 10 15 20 25 30 35 40 45 50

Heun 3 Exp(-t)

0 0.2 0.4 0.6 0.8 1

0 5 10 15 20 25 30 35 40 45 50

RK Exp(-t)

ホインの３次公式(t= 2.5) ルンゲクッタ法(t= 2.5)

-200000 -100000 0 100000 200000 300000 400000

0 5 10 15 20 25 30 35 40 45 50

Euler Exp(-t)

0 1e+07 2e+07 3e+07 4e+07 5e+07 6e+07 7e+07

0 5 10 15 20 25 30 35 40 45 50

I Euler Exp(-t)

オイラー法(t= 3.5) 改良オイラー法(t= 3.5)

-2e+07 -1e+07 0 1e+07 2e+07 3e+07 4e+07 5e+07

0 5 10 15 20 25 30 35 40 45 50

Heun 3 Exp(-t)

0 200000 400000 600000 800000 1e+06 1.2e+06 1.4e+06

0 5 10 15 20 25 30 35 40 45 50

RK Exp(-t)

ホインの３次公式(t= 3.5) ルンゲクッタ法(t= 3.5)

-4 -3 -2 -1 0 1 2 3 4

-4 -3 -2 -1 0 1 2 3

p=4 p=3 p=2 p=1

ルンゲクッタ型公式の安定領域 （kが大きいほど,安定領域は大きくなる）

★ シンプレクティック解法

• x^′′=−x, x(0) = 1.0,x^′(0) = 0.0に種々の方法を適用した例. h= 0.10,相空間ではt= 200πまで を表示.

-1.5 -1 -0.5 0 0.5 1 1.5

-1.5 -1 -0.5 0 0.5 1 1.5

x’(t), p(t)

x(t), q(t) Harmonic Ossirator

Euler

0 5 10 15 20 25

0 10 20 30 40 50 60

Hamiltonian Error

t Harmonic Ossirator

Euler

オイラー法

-1.5 -1 -0.5 0 0.5 1 1.5

-1.5 -1 -0.5 0 0.5 1 1.5

x’(t), p(t)

x(t), q(t) Harmonic Ossirator

Improved_Euler

0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008

0 10 20 30 40 50 60

Hamiltonian Error

t Harmonic Ossirator

Improved_Euler

改良オイラー法

-1.5 -1 -0.5 0 0.5 1 1.5

-1.5 -1 -0.5 0 0.5 1 1.5

x’(t), p(t)

x(t), q(t) Harmonic Ossirator

Heun_3

-0.003 -0.0025 -0.002 -0.0015 -0.001 -0.0005 0

0 10 20 30 40 50 60

Hamiltonian Error

t Harmonic Ossirator

Heun_3

ホインの３次公式

-1.5 -1 -0.5 0 0.5 1 1.5

-1.5 -1 -0.5 0 0.5 1 1.5

x’(t), p(t)

x(t), q(t) Harmonic Ossirator

Runge-Kutta

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

0 10 20 30 40 50 60

Hamiltonian Error

t Harmonic Ossirator

Runge-Kutta

ルンゲクッタ法

-1.5 -1 -0.5 0 0.5 1 1.5

-1.5 -1 -0.5 0 0.5 1 1.5

x’(t), p(t)

x(t), q(t) Harmonic Ossirator

Symplectic_Euler

-0.06 -0.04 -0.02 0 0.02 0.04 0.06

0 10 20 30 40 50 60

Hamiltonian Error

t Harmonic Ossirator

Symplectic_Euler

シンプレクティックオイラー法

-1.5 -1 -0.5 0 0.5 1 1.5

-1.5 -1 -0.5 0 0.5 1 1.5

x’(t), p(t)

x(t), q(t) Harmonic Ossirator

SV

-0.0025 -0.002 -0.0015 -0.001 -0.0005 0

0 10 20 30 40 50 60

Hamiltonian Error

t Harmonic Ossirator

SV

シュテルマ・ベルレ法

-1.5 -1 -0.5 0 0.5 1 1.5

-1.5 -1 -0.5 0 0.5 1 1.5

x’(t), p(t)

x(t), q(t) Harmonic Ossirator

Ruth

-3e-05 -2e-05 -1e-05 0 1e-05 2e-05 3e-05

0 10 20 30 40 50 60

Hamiltonian Error

t Harmonic Ossirator

Ruth

ルース法

-1.5 -1 -0.5 0 0.5 1 1.5

-1.5 -1 -0.5 0 0.5 1 1.5

x’(t), p(t)

x(t), q(t) Harmonic Ossirator

Implicit_RK_Step_1

-6e-16 -4e-16 -2e-16 0 2e-16 4e-16 6e-16 8e-16

0 10 20 30 40 50 60

Hamiltonian Error

t Harmonic Ossirator

Implicit_RK_Step_1

陰的中点法（１段陰的ルンゲクッタ法）

• x^′′=−x,x(0) = 1.0,x^′(0) = 0.0に種々の方法を適用し,真の解との誤差を表示した例. h= 0.10.

-0.1 -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08

0 50 100 150 200 250 300 350

Symplectic Euler

-0.15 -0.1 -0.05 0 0.05 0.1 0.15

0 50 100 150 200 250 300 350

SV

シンプレクティック・オイラー法 シュテルマ・ベルレ法

-5e-05 -4e-05 -3e-05 -2e-05 -1e-05 0 1e-05 2e-05 3e-05 4e-05 5e-05

0 50 100 150 200 250 300 350

Ruth

-5e-05 -4e-05 -3e-05 -2e-05 -1e-05 0 1e-05 2e-05 3e-05 4e-05 5e-05

0 50 100 150 200 250 300 350

IRK-2

ルース法 陰的中点法

● 実習内容

1. 単振動とは,ハミルトニアン

H(q(t), p(t)) = 1

2p(t)²+1

2q(t)², q(t), p(t)∈R で表される運動である.

この運動を,適当な初期条件のもとに「シンプレクティックオイラー法」で計算し, (q, p)平面内で運 動を図示しなさい.

2. オモリの質量 m,長さℓの単振り子の運動の位置エネルギー,運動エネルギーは, つりあいの位置か らの角度をθ(t)としたとき,

T = 1

2mℓ²(θ^′)², U =−mℓgcosθ となる. したがって,ハミルトニアンは

H(q, p) = p²

2mℓ² −mℓgcosq, 運動方程式は,

dq dt = p

mℓ², dp

dt =−mℓgsinq となる.

この運動を,m,ℓ,g を全て1 に取り,適当な初期条件のもとに「シンプレクティックオイラー法」で 計算し, (q, p)平面内で運動を図示しなさい.

● レポート問題

★ レポート問題

以下の問題は評価対象のレポート問題です.

【レポート問題03】 ルンゲ・クッタ法を用いて,常微分方程式の初期値問題

x^′′(t) + sin(x(t)) = 0, x(0) = 2.5, x^′(0) = 0.0

の数値解を計算しなさい. ただし,時間刻み幅をh= 0.01とし,T = 12.0までの解を求め,この数値 解を相平面上に図示しなさい.

【レポート問題04】 ルンゲ・クッタ法を常微分方程式の初期値問題

x^′(t) =f(t), x(0) = 0 に適用して,積分

Z t 0

f(s)ds

の近似値を求めることを考える. この時,ルンゲ・クッタ法から得られる積分の近似式が

Z h 0

f(s)ds∼ h

6(f(0) + 4f(h/2) +f(h))

であることを証明し,右辺の式は,区間[0, h]でf をどのような関数で近似したことなるかをのべな さい.

【レポート問題05】 s段のルンゲ・クッタ型公式が3次公式となるための条件を求めなさい. また,クッ タの３次公式とホインの３次公式が,それらの条件を満たすことを示しなさい.

これらの問題については,A: 10 点, C: 0 点,と採点します. （中間的な点をつける可能性あり）

★ 注意事項

【締め切り】 締め切りのおおよその目安は,２００９年７月１５日（水）とします.

その他の注意事項については,前回のレポートの注意事項と同じです.

★ レポート問題

(Extra)

以下の問題は評価対象のレポート問題ですが,この問題を「評価の点の満点」にはカウントしません.（つ まり,以下の問題での得点は「ボーナスポイント」となります.）

【レポート問題E03】 後退オイラー法を用いて,常微分方程式の初期値問題

x^′′(t) + sin(x(t)) = 0, x(0) = 2.5, x^′(0) = 0.0

の数値解を計算しなさい.ただし,時間刻み幅をh= 0.01とし,T = 12.0までの解を求めなさい.ま た,この数値解を相平面上に図示しなさい.

【レポート問題E04】 シンプレクティック・オイラー法を用いて,常微分方程式の初期値問題

x^′′(t) + sin(x(t)) = 0, x(0) = 2.5, x^′(0) = 0.0

の数値解を計算しなさい. ただし,時間刻み幅をh= 0.01とし,T = 12.0までの解を求め,この数値 解を相平面上に図示しなさい. また, この微分方程式にシンプレクティック・オイラー法を適用可能 であることを示しなさい.

【レポート問題E05】 陰的中点法を用いて,常微分方程式の初期値問題

x^′′(t) + sin(x(t)) = 0, x(0) = 2.5, x^′(0) = 0.0

の数値解を計算しなさい. ただし,時間刻み幅をh= 0.01とし,T = 12.0までの解を求め,この数値 解を相平面上に図示しなさい.

【レポート問題E06】 s段のルンゲ・クッタ型公式が 4次公式となるための条件を求めなさい. また,ル ンゲ・クッタ法がそれらの条件を満たすことを示しなさい.

問題E03, E04については, 10 点満点,問題E05, E06 については,15 点満点で採点します. （中間的な

点をつける可能性あり）Extraな問題の締め切りは２００９年８月半ばとし,提出方法は他の問題と同じと します. （詳細な締め切りは,最終回の講義でお知らせします）