● 講義資料

● 講義資料

▼ 講義予定

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

– Runge-Kutta Methodの安定性

● 前回の講義のまとめ

★ より高次の計算スキーム

• オイラー法はhに関して1次の公式であり,実用上は誤差が大きくて使い物にならない. し たがって,hに関して,より高次の公式,すなわち,

|x(T)−xN|=O(h^p), p >1 となる公式を探す必要がある.

• 最初に思いつくのは,テイラー展開

x(tn+1) =x(tn) +hx^′(tn) +h²

2 x^′′(tn) +O(h³),

=x(tn) +hf(x(tn)) +h²

2 f(x^′(tn))f(x(tn)) +O(h³), のh^k の項までを取り入れたスキームである. たとえば,k= 2 と考え,

xn+1=xn+hf(xn) +h²

2 f^′(xn)f(xn)

と計算すると,局所離散化誤差がO(h³) となることが期待できる. 実際,この計算スキーム は, (h^k の項までのテイラー展開に相当するものを使えば, ）局所離散化誤差はO(h^k+1)と なり,k次の公式となることが証明できる. （この方法をテイラー展開法と呼ぶ.）しかし,計 算スキームの中にf の高階微分が入り,プログラム中に書き込む情報が多くなってしまうた め,実用上は使われることはない.

• 別のアイデアとして,前進オイラー法・後退オイラー法の時のテイラー展開の計算をみると, x(tn+1) =x(tn) +hx^′(tn) +h²

2 x^′′(tn) +h³

6 x⁽³⁾(t) +O(h⁴), x(tn−1) =x(tn)−hx^′(tn) +h²

2 x^′′(tn)−h³

6 x⁽³⁾(t) +O(h⁴), となっているので,これを加えれば,

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

3 x⁽³⁾(t) +O(h⁵) となり,

xn+1=xn−1+ 2hf(xn)

• 修正オイラー法をx^′(t) =−x(t),x(0) = 1に適用してみると, xn+1=xn−1−2hxn, x0= 1, x1= 1 +h となる. この漸化式で定義される数列{xn}は,

|xn| → ∞, n→ ∞

となってしまう. 厳密解はx(t) =e^−tであり,x(t)→0 (t→ ∞)のはずである. （この状況 を,x(t)は安定解であると言う.）

厳密解は安定解であり, xn →0 となるべきところが, 離散化したことにより,数値解が不安 定になってしまう. （この現象を数値的不安定性と呼ぶ.）

したがって,この方法は数値解の構成方法として不適当であることがわかる.

• 修正オイラー法が数値的不安定性を生む原因は,３項間漸化式

xn+1=xn−1−2hxn

の特性方程式

t²+ 2ht−1 = 0 の解の配置にある. この方程式の解α1, α2 は,

α1<−1, 0< α2<1 であり,一般解は

xn=A1αⁿ₁ +A2αⁿ₂

と書ける. ここで,xn→0であるための必要十分条件はA1= 0であるが,x0=C,x1=Cα1

以外の場合にはA16= 0 となる. そのため,|xn| → ∞となることがわかる.

なお,ほぼ同じ状況がフィボナッチ数列

an+1=an+an−1

に関しても生じる.

● 講義資料

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

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

（特に書いていない限り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

-3 -2 -1 0 1 2 3

0 2 4 6 8 10 12 14 16

x’’ = -sin(x)

-2 -1 0 1 2

-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

-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)

-0.5 0 0.5 1 1.5 2

Double Pendulum

-0.5 0 0.5 1 1.5

Double Pendulum

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

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

（特に書いていない限り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

-3 -2 -1 0 1 2 3

0 2 4 6 8 10 12 14 16

x’’ = -sin(x)

-2 -1 0 1 2

-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

-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)

-0.5 0 0.5 1 1.5 2

Double Pendulum

-0.5 0 0.5 1 1.5

Double Pendulum

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

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

（特に書いていない限り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

-3 -2 -1 0 1 2 3

0 2 4 6 8 10 12 14 16

x’’ = -sin(x)

-2 -1 0 1 2

-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

-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)

-0.5 0 0.5 1 1.5 2

Double Pendulum

-0.5 0 0.5 1 1.5

Double Pendulum

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

以下の図は,種々のルンゲ・クッタ型公式で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を小さくとってしまうと,誤差が必ずしも小さくならないことに注意が必要 である.

F ψ α γ d

τ1 f i bi 1 1 1

τ2 f^′(f)

j

i biaij 1 2 2

τ3,1 f^′′(f, f)

k j

i biaijaik 1 3 3

τ3,2 f^′(f^′(f))

k j

i biaijajk 1 6 3

τ4,1 f^′′′(f, f, f)

k ℓ j

i biaijaikaiℓ 1 4 4

τ4,2 f^′(f^′′(f, f))

j k

ℓ

i biaijaikakℓ 3 8 4

τ4,3 f^′′(f^′(f), f)

ℓ k

j

i biaijajkajℓ 1 12 4

τ4,4 f^′(f^′(f^′(f)))

ℓ k j

i biaijajkakℓ 1 24 4

● 実習内容

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