error_g1.eps

(1)

● 講義資料

▼ 講義予定

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

● 前回の講義のまとめ

★ 前進オイラー法 • 最も単純な離散変数法は, 前進オイラー法 xn+1= xn+ hf (tn, xn) である. （単に「オイラー法」と言ったときには, この前進オイラー法を指す.） • 微分方程式を離散変数法で数値的に解く場合には, 微分 dx dt を何らかの形で離散化する必要がある. そのため, 以下のように差分化を行う. dx dt(t) = limh→0 x(t + h) − x(t) h を t = tn で考えれば, h > 0 が小さいと仮定して, dx dt(tn) = limh→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 を計算可能である. その意味で１段階法となっている.

(2)

★ 後退オイラー法 • 前進オイラー法での差分化を行った点を t = tnではなく t = tn+1 での値と考える. すなわち dx dt(tn+1) = limh→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+ h sin(xn+1) を得る. ここで, xn から xn+1 を得るために, この式を代数的に xn+1 について解くことはできない. したがって, 逐次近似 y0= xn, yk+1= xn+ h sin(yk) を用いて, その極限として xn+1 を求める必要がある. （なお, h が十分に小さければ, この逐次近似が収束することを証明できる.）

(3)

★ オイラー法の誤差 • 今後, 考える微分方程式の初期値問題は, 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) + h2 2 x ′′ (tn) + O(h3), = x(tn) + hf (x(tn)) + h2 2 f (x ′ (tn))f (x(tn)) + O(h3), xn+1= xn+ hf (xn) x(tn+1) − xn+1= x(tn) − xn+ h(f (x(tn)) − f (xn)) + h2 2 f ′ (x(tn))f (x(tn)) + O(h3) と計算できる. ここで, f がリプシッツ条件 |f (p) − f (q)| ≤ L|p − q| をみたし, f はなめらかであると仮定すると, |x(tn+1) − xn+1| ≤ (1 + Lh)|x(tn) − xn| + M h2 が成り立つ. • 後退オイラー法の場合には, x(tn) = x(tn+1) − hx′(tn+1) +h 2 2 x ′′ (tn+1) + O(h3), = x(tn+1) − hf (x(tn+1)) + h2 2 f (x ′ (tn+1))f (x(tn+1)) + O(h3), xn+1= xn+ hf (xn+1) x(tn+1) − xn+1= x(tn) − xn+ h(f (x(tn+1)) − f (xn+1)) + h2 2 f ′ (x(tn+1))f (x(tn+1)) + O(h3) より, |x(tn+1) − xn+1| ≤ (1 − Lh)−1|x(tn) − xn| + M h2 が成り立つ. • 前進オイラー法, 後退オイラー法のいずれの場合にも, ǫn+1= Aǫn+ M h2, ǫn= |x(tn) − xn| の形の不等式を得ることができた. この評価を局所離散化誤差と呼び, ある時刻 t = tn において, 数値解の値 xn と厳密解の値 x(tn)が一致しているとき, １ステップ進んだ際の誤差が O(h2₎_{であることを示している.}

(4)

• この局所離散化誤差の不等式から t = T = Nh においての誤差 |x(T ) − xN| ≤ CheT L, 前進オイラー法, |x(T ) − xN| ≤ Ch, 後退オイラー法, を得ることができる. すなわち, オイラー法は誤差が h に関して１次の公式であることがわかる. ★ より高次の計算スキーム • オイラー法は h に関して 1 次の公式であり, 実用上は誤差が大きくて使い物にならない. したがって, h に関して, より高次の公式, すなわち, |x(T ) − xN| = O(hp), p > 1 となる公式を探す必要がある. • 最初に思いつくのは, テイラー展開 x(tn+1) = x(tn) + hx ′ (tn) + h2 2 x ′′ (tn) + O(h3), = x(tn) + hf (x(tn)) + h2 2 f (x ′ (tn))f (x(tn)) + O(h3), の hk _{の項までを取り入れたスキームである. たとえば, k = 2 と考え,} xn+1= xn+ hf (xn) + h2 2 f ′ (xn)f (xn) と計算すると, 局所離散化誤差が O(h3₎ _{となることが期待できる. 実際, この計算スキーム} は, (hk _{の項までのテイラー展開に相当するものを使えば, ）局所離散化誤差は O(h}k+1₎_となり, k 次の公式となることが証明できる. （この方法をテイラー展開法と呼ぶ.）しかし, 計算スキームの中に f の高階微分が入り, プログラム中に書き込む情報が多くなってしまうため, 実用上は使われることはない. • 別のアイデアとして, 前進オイラー法・後退オイラー法の時のテイラー展開の計算をみると, x(tn+1) = x(tn) + hx ′ (tn) + h2 2 x ′′ (tn) + h3 6 x (3)_{(t) + O(h}4_), x(tn−1) = x(tn) − hx ′ (tn) + h2 2 x ′′ (tn) − h3 6 x (3)_{(t) + O(h}4_), となっているので, これを加えれば, x(tn+1) − x(tn−1) = 2hx′(tn) + h3 3 x (3)_{(t) + O(h}5₎ となり, xn+1= xn−1+ 2hf (xn) は, 局所離散化誤差が O(h3₎_{となるように思える. この方法を修正オイラー法と呼ぶ.}

(5)

• 修正オイラー法を 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 の特性方程式 t2_{+ 2ht − 1 = 0} の解の配置にある. この方程式の解 α1, α2 は, α1< −1, 0 < α2< 1 であり, 一般解は xn= A1αn1 + A2αn2 と書ける. ここで, xn→ 0であるための必要十分条件は A1= 0であるが, x0= C, x1= Cα1 以外の場合には A16= 0 となる. そのため, |xn| → ∞となることがわかる. なお, ほぼ同じ状況がフィボナッチ数列 an+1= an+ an−1 に関しても生じる.

(6)

● 講義資料

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

以下の図は, 改良オイラー法で種々の常微分方程式の初期値問題の数値解を計算した例である. （特に書いていない限り 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)2_{), 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

(7)

-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)2_)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)|3_{, 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′1(0) = x ′ 2(0) = 0

(8)

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

以下の図は, ホインの３次公式で種々の常微分方程式の初期値問題の数値解を計算した例である. （特に書いていない限り 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)2_{), 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

(9)

(10)

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

以下の図は, ルンゲ・クッタ法で種々の常微分方程式の初期値問題の数値解を計算した例である. （特に書いていない限り 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)2_{), 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

(11)

(12)

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

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