● 講義資料

● 講義資料

▼ 講義予定

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

– Runge-Kutta Methodの安定性

● 前回の講義のまとめ

★ 前進オイラー法

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

xn+1=xn+hf(tn, xn)

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

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

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

dx

dt(t) = lim

h→0

x(t+h)−x(t) h をt=tn で考えれば,h >0 が小さいと仮定して,

dx

dt(tn) = lim

h→0

x(tn+h)−x(tn)

h ∼x(tn+1)−x(tn) h

と書き直すことができる. （これを前進差分と呼ぶ.）これを, 微分方程式x^′ =f(t, x)に適 用すれば,

x(tn+1)−x(tn)

h ∼f(tn, xn) x(tn+1)∼x(tn) +hf(tn, x(tn)) となる. したがって,前進オイラー法の計算スキーム

xn+1=xn+hf(tn, xn) に相当する式を得ることができる.

• 前進オイラー法のスキームの右辺は,直前のxnの値があれば直接計算可能である. その意味 で陽的解法(explicit method)となっている.

• 前進オイラー法のスキームは,直前のxn の値のみで,次のxn+1 を計算可能である. その意 味で１段階法となっている.

★ 後退オイラー法

• 前進オイラー法での差分化を行った点をt=tnではなくt=tn+1 での値と考える. すなわち dx

dt(tn+1) = lim

h→0

x(tn+h)−x(tn)

h ∼ x(tn+1)−x(tn) h と書き直すことができる. これを,微分方程式x^′ =f(t, x)に適用すれば,

x(tn+1)−x(tn)

h ∼f(tn+1, xn+1) x(tn+1)∼x(tn) +hf(tn+1, x(tn+1))

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

xn+1=xn+hf(tn+1, xn+1) を後退オイラー法(backward Euler method)と呼ぶ.

• 後退オイラー法のスキームは,直前のxn の値のみで,次のxn+1 を計算可能である. その意 味で１段階法となっている.

• 後退オイラー法のスキームの右辺は,直前のxnの値だけからは直接計算できない. その意味 で陰的解法(implicit method)となっている.

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

– 微分方程式x^′(t) =x(t)に後退オイラー法を適用すると, xn+1=xn+hxn+1

を得る. ここで, xn から xn+1 を得るためには, この式を代数的に xn+1 について解 いて,

xn+1= 1 1−hxn

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

– 微分方程式x^′(t) = sin(x(t))に後退オイラー法を適用すると, xn+1=xn+hsin(xn+1)

を得る. ここで, xn からxn+1 を得るために, この式を代数的にxn+1 について解くこ とはできない. したがって,逐次近似

y0=xn,

yk+1=xn+hsin(yk)

を用いて,その極限としてxn+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(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³), xn+1=xn+hf(xn)

x(tn+1)−xn+1=x(tn)−xn+h(f(x(tn))−f(xn)) +h²

2 f^′(x(tn))f(x(tn)) +O(h³) と計算できる. ここで,f がリプシッツ条件

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

をみたし,f はなめらかであると仮定すると,

|x(tn+1)−xn+1| ≤(1 +Lh)|x(tn)−xn|+M h² が成り立つ.

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

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

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

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

2 f(x^′(tn+1))f(x(tn+1)) +O(h³), xn+1=xn+hf(xn+1)

x(tn+1)−xn+1=x(tn)−xn+h(f(x(tn+1))−f(xn+1)) +h²

2 f^′(x(tn+1))f(x(tn+1)) +O(h³) より,

|x(tn+1)−xn+1| ≤(1−Lh)⁻¹|x(tn)−xn|+M h² が成り立つ.

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

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

の形の不等式を得ることができた. この評価を局所離散化誤差と呼び,ある時刻t=tn にお いて,数値解の値 xn と厳密解の値 x(tn)が一致しているとき, １ステップ進んだ際の誤差が O(h²)であることを示している.

• この局所離散化誤差の不等式からt=T =N h においての誤差

|x(T)−xN| ≤Che^{T L}, 前進オイラー法,

|x(T)−xN| ≤Ch, 後退オイラー法,

を得ることができる. すなわち,オイラー法は誤差がh に関して１次の公式であることがわ かる.

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

• オイラー法は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)

は,局所離散化誤差がO(h³)となるように思える. この方法を修正オイラー法と呼ぶ.

• 修正オイラー法を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)

-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

二重振り子,x1(0) = 2,x2(0) = 0,x^′₁(0) =x^′₂(0) = 0

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

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

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

-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

二重振り子,x1(0) = 2,x2(0) = 0,x^′₁(0) =x^′₂(0) = 0

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

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

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

-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

二重振り子,x1(0) = 2,x2(0) = 0,x^′₁(0) =x^′₂(0) = 0

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

以下の図は,種々のルンゲ・クッタ型公式で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の常微分方程式の初期値問題の数値解を,改良オイラー法,ホインの ２次公式,ホインの３次公式,ルンゲ・クッタ法で計算しなさい.