● 講義資料

● 講義資料

▼ 講義予定

• 常微分方程式の初期値問題の数値解法 – Runge-Kutta Methodの安定性 – Runge-Kutta Methodと数値積分

● 前回の講義のまとめ

★ Runge-Kutta法

• 次の計算スキームをRunge-Kutta 法と呼び,s をその段数と呼ぶ.

ki =f(x0+h

s

X

j=1

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

x1=x0+h

s

X

i=1

biki,

• 特に,i≤j に対しては, aij= 0が成り立つものを陽的 Runge-Kutta法と呼ぶ. すなわち, ki =f(x0+h

i−1

X

j=1

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

x1=x0+h

s

X

i=1

biki, である. 具体的に書けば,











k1=f(x0),

k2=f(x0+ha21k1), x1=x0+h(b1k1+b2k2)

s= 2,















k1=f(x0),

k2=f(x0+ha21k1),

k3=f(x0+ha31k1+ha32k2), x1=x0+h(b1k1+b2k2+b3k3)

s= 3,

などとなる.

• Runge-Kutta法の代表的なものは,以下の通りである.

改良オイラー法 （２段２次）

b1= 0, b2= 1, a21= 1 2,











k1=f(x0), k2=f(x0+h

2k1), x1=x0+hk2

ホインの２次公式 （２段２次）

b1=b2=1

2, a21= 1,











k1=f(x0), k2=f(x0+hk1), x1=x0+h

2(k1+k2) ホインの３次公式 （３段３次）

b1= 1/4, b2= 0, b3= 3/4, a21= 1/3, a31= 0, a32= 2/3























k1=f(x0), k2=f(x0+h

3k1), k3=f(x0+2h

3 k2), x1=x0+h

4(k1+ 3k3) クッタの３次公式 （３段３次）

b1= 1/6, b2= 2/3, b3= 1/6, a21= 1/2, a31=−1, a32= 2























k1=f(x0), k2=f(x0+h

2k1), k3=f(x0−hk1+ 2hk2), x1=x0+h

6(k1+ 4k2+k3) Runge-Kutta法 （４段４次公式）

a21= 1/2, a31= 0, a32= 1/2, a41= 0, a42= 0, a43= 1, b1= 1/6, b2= 1/3, b3= 1/3, b4= 1/6































k1=f(x0), k2=f(x0+h

2k1), k3=f(x0+h

2k2), k4=f(x0+h

2k3), x1=x0+h

6(k1+ 2k2+ 2k3+k4)

通常「Runge-Kutta法」と言うときには,上記の４段４次のRunge-Kutta法を指す.

• s段のRunge-Kutta型公式がp次となる条件を求めるためには, 以下の方法が考えられる.

（多くの書籍では,この方法を組み合わせ論的に整理した証明を用いている）

▼ テイラー展開を書く

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

2 x^′′(t) +h³

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

4!x⁽⁴⁾(t) +O(h⁵), ここで,x[0, T]−→Rⁿ とすれば,t∈[0, T]を固定したとき,

f(x(t)) =f ∈Rⁿ,

f^′(x(t)) =f^′∈Hom(Rⁿ,Rⁿ), であることに注意して,

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

dtf(x(t)) =f^′(x(t))x^′(t) =f^′f, となる. よって,

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

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

=x(t) +hf+h²

2 f^′f+O(h³) が成り立つ.

▼ ki の展開式を計算する ki をhの関数として展開すると,以下が成り立つ.

ki(h) =f(x0+hX

aijkj(h)), ki^′(h) =f^′(x0+hX

aijkj(h)) hX

aijkj^′(h) +X

aijkj(h) , これらのh= 0での値は,

ki(0) =f, k_i^′(0) =f^′X

aijkj

=X aij

f^′f となる. これを書き直せば,

ki(0) =f, k^′_i(0) =X

aij

f^′f

▼ テイラー展開に代入する これらのki(0),k^′_i(0)を使って, x1=x0+hX

biki(h)

=x0+hX

biki(0)

+h²X

bik^′_i(0)

+O(h³)

=x0+hX bi

f+h²X biaij

f^′f +O(h³) が成り立つ.

▼ テイラー展開に代入する この計算を比較すれば,

hの項 X

bi= 1, h² の項 X

biaij =1 2, をみたせば,局所離散化誤差をh³ とすることができる.

• より高次の条件式は以下の通りである.

３次の条件式







Xbiaijaik = 1 3, Xbiaijajℓ =1

6

４次の条件式























Xbiaijaikaiℓ=1 4, Xbiaijajkaiℓ= 1

8, Xbiaijajkajℓ= 1 12, Xbiaijajkakℓ = 1

24

• これらの条件式を各sについて具体的に書けば以下のようになる.

▼ ２段２次公式

b1+b2= 1, b2a21= 1

2

▼ ３段３次公式

b1+b2+b3= 1,

b2a21+b3(a31+a32) = 1 2, b2a²₂₁+b3(a31+a32)²= 1 3, b3a32a21=1

6

▼ ４段４次公式

b1+b2+b3+b4= 1,

a21b2+a31b3+a32b3+a41b4+a42b4+a43b4=1 2, a²₂₁b2+ (a31+a32)²b3+ (a41+a42+a43)²b4=1

3, a21a32b3+ (a21a42+a31a43+a32a43)b4= 1

6, a³₂₁b2+ (a31+a32)³b3+ (a41+a42+a43)³b4=1

4,

a32a21(a31+a32)b3+b4(a42a21+a43(a31+a32))(a41+a42+a43) = 1 8, a²₂₁a32b3+ (a²₂₁a42+ (a31+a32)²a43)b4= 1

12, a21a32a43b4= 1

24

これらの条件を満たせば,s段p次公式,すなわち,

|x(tn)−xn| ≤Ch^p となる公式となる.

• なお,１次の条件式P

bi= 1 が成り立つことは,計算スキームが意味を持つための必要条件 であるので,適切性の条件と呼ばれる. 逆に,Pbi = 1さえ成り立っていれば, aij の数値の 入力ミスがあっても, より低い次数の数値解をつくることができてしまうことに注意が必要 である.

• 陽的Runge-Kutta型公式については,以下の事実が知られている.

– s段でs+ 1次以上の公式をつくることはできない.

– s≥5の時,s段s次公式をつくることはできない.

• なお, Runge-Kutta法の係数{aij},{bi},{ci},ただし,ci=P^s

j=1aij を c1 a11 · · · a1s

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

b1 · · · bs

と表記することが多く,これをButcher 行列と呼ぶ.

• 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

● 講義資料

▼

Runge-Kutta

以下の図は,s段陽的ルンゲ・クッタ型公式の絶対安定領域を示している.

-3 -2 -1 0 1 2 3

-3 -2 -1 0 1 2 3

Stable Region for Explicit Runge-Kutta Methods s=1 s=2 s=3 s=4

以下の図は, 前進 Euler 法 (λh = −2.0), Runge-Kutta 法 (λh = −3.0) で, y^′(t) = −λy(t), y(0) = 1.0を解いた例.

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

0 20 40 60 80 100 120 140

Forward Euler. z = 2.5

0 1e+06 2e+06 3e+06 4e+06 5e+06 6e+06

0 20 40 60 80 100 120 140 160

Runge-Kutta. z = 3.0

▼ 補間法

以下の図は,f(x) = 1

1 + 25x² を多項式補間した計算例 等間隔補間

n= 2 n= 4 n= 6

-1 -0.5 0.5 1

0.2 0.4 0.6 0.8 1

-1 -0.5 0.5 1

-0.4 -0.2 0.2 0.4 0.6 0.8 1

-1 -0.5 0.5 1

0.2 0.4 0.6 0.8 1

n= 10 n= 14 n= 18

-1 -0.5 0.5 1

0.5 1 1.5 2

-1 -0.5 0.5 1

-0.5 0.5 1 1.5 2 2.5

-1 -0.5 0.5 1

-2 -1 1 2

Chebyshev補間

n= 2 n= 4 n= 6

-1 -0.5 0.5 1

-0.2 0.2 0.4 0.6 0.8 1

-1 -0.5 0.5 1

0.2 0.4 0.6 0.8 1

-1 -0.5 0.5 1

0.2 0.4 0.6 0.8 1

n= 8 n= 10 n= 20

-1 -0.5 0.5 1

0.2 0.4 0.6 0.8 1

-1 -0.5 0.5 1

0.2 0.4 0.6 0.8 1

-1 -0.5 0.5 1

0.2 0.4 0.6 0.8 1

● 実習内容

1. （★★★）Runge-Kutta型公式から得られる種々の数値積分公式を利用して,以下の積分の 近似値を計算しなさい. また,刻み幅hを変化させたとき,近似値と真の値との誤差のグラフ をつくりなさい.

(a) Z 2

1

dx x (b)

Z 1 0

dx 1 +x²

2. （★★）以下の関数を, 等間隔のサンプル点をとって, 2n次多項式で近似しなさい. （近似 多項式のグラフと,関数のグラフの両方を書いて比較しなさい）

(a) f(x) = cos(x) (区間は [−π/2, π/2]) (b) f(x) = 1

1 + 25x² (区間は[−1,1])

3. （★）10次以下のLegrandre多項式の零点（の近似値）をすべて計算しなさい.

4. （★）10次以下のChebyshev 多項式の零点（の近似値）をすべて計算しなさい.

5. （★）等間隔のサンプル点を用いて, f(x) = 1

1 + 25x² を区間[−1,1上で2n次多項式で近 似しなさい. （近似多項式のグラフと,関数のグラフの両方を書いて比較しなさい）

6. （★）チェビシェフ多項式の零点をサンプル点として,f(x) = 1

1 + 25x² を区間[−1,1上で 2n次多項式で近似しなさい. （近似多項式のグラフと,関数のグラフの両方を書いて比較し なさい）