= f ( t;x ( t )) dxdt 2. 常微分方程式 数値解析

数値解析

2. _{常微分方程式} dx

dt = f (t, x(t))

石渡哲哉

(

)

この章の目標

常微分方程式の初期値問題 (IVP):

dx

dt = f(t, x(t)), t > a, (1)

x(a) = x₀ (2)

の近似計算法を学び、誤差解析等について理解する。

対象:

• ^{上記のような} 1 _{階の正規形} ODE に対する方法を紹介するが、1 _階

の連立系にしてもまったく同じように適用可能である。スカラー 関数 x, f をベクトル値の関数に読み替えるだけで良い。また、高 階 ODE _は、1 階連立系に書きなおすことができることにも注意。

• f は十分滑らかとし、解の存在と一意性は保証されているとする。

また、このとき解 x _も f の滑らかさに応じて十分に滑らかであ ることにも注意。

代表的な方法

• ^{オイラー法} (Euler method)

• ^ホイン法 (Heun method), _中点法 (_{改良オイラー法})

• ^{ルンゲ・クッタ法} (Runge-Kutta method)

2.1 _{オイラー法} (Euler method)

dx

dt = f(t, x(t))

の左辺を次で近似 (h > 0 _とする)_： x(t + h) − x(t)

h ⇐= dx

dt = lim

h→0

x(t + h) − x(t) h

式を整理すると,

x(t + h) = x(t) + h · f(t, x(t))

となり, x(t) _の値から x(t + h) の値を逐次的に計算することができる.

t_n := a + nh _とし, x(t_n) _{の近似値を} X_{n と書くことにする}. (IVP) を次にように離散化する方法をオイラー法という.

オイラー法のスキーム

{ X_n+1 = X_n + h · f(t_n, X_n) (n = 0, 1, 2, · · · ) X₀ = x₀

h > 0 はステップサイズ、刻み幅、離散化パラメータなどと呼ばれる.

近似の妥当性

Q h → 0 _{としたとき}, X_{n は元の解} x(t) _{に収束するだろうか？}^*1

*1 より詳しく, (t_n, X_n) をつないでできた折れ線が x(t) _{の解曲線に収束するか}, _と いう観点もあり.

■ 解の収束性について (_誤差解析)

閉区間 [a, b] (T := b − a > 0) _上で解 x(t) _{が存在しているとする}. h = T /N (N ∈ R) _とし, _{この区間で離散化誤差} e_n := x(t_n)−X_n (n = 0, 1, 2, · · · , N) _{を調べる。}

定理 (_誤差解析)

離散化誤差に対して以下が成り立つ:

h _や N _{に依存しない定数} C > 0 _{が存在して} max

0≤n≤N |e_n| ≤ Ch.

Remark:

• 上の結果をランダウの記号を使って |e_n| = O(h) _{と書くことも}

ある.

• h ↘ 0 _のとき, _誤差が h _の 1 _{次のオーダーで} 0 _{に収束する}.

• 今回の誤差解析は離散化誤差のみ考慮している. _{丸め誤差の影響}

についても考察する必要がある場合は, _{初期値の設定}, _{各ステップ}

の計算時に丸め誤差が混入するとして解析を行う.

オイラー法はシンプルであり、プログラミングも容易であるが、収束が 遅く、実用上はあまり使われない。

次節以降において、オイラー法の改良や、ODE _{の数値計算でよく使わ}

れているルンゲ・クッタ法について紹介する.

2.2 _{次数と収束}

以下の 1 段階法のスキームを考える:

X_n+1 = X_n + h · Φ(t_n, X_n) (n = 0, 1, 2, · · · , N − 1) (3) x(t) _{を真の解として、}

x(t_n+1) = x(t_n) + h · Φ(t_n, x(t_n)) + O(h^k+1)

を満たすとき, _スキーム (3) _は k _次の / _次数 k _の (_{精度をもつ}) _離散化

である, _という.

このとき, 以下の収束定理が成り立つ:

収束定理 Φ はリプシッツ連続とする. 1 _{段階法のスキーム} (3) _が k _次精度を

もつとすると, _{以下が成立する}: max

0≤n≤N |e_n| ≤ O(h^k).

上の定理より, 近似の精度を上げれば収束も速くなることが分かる.

よって, 比較的シンプルなスキームで^*2, なるべく高精度のものを作り たい, _{ということになる}.

*2 計算の手間が少なければ, それだけ丸め誤差の混入の機会が減る.

2.3 _{オイラー法の改良}

ODE x^′(t) = f(t, x(t)) _を t_{n から} t_{n+1 まで積分すると、}

x(t_n+1) − x(t_n) =

∫ _tn+1

t_n

f(s, x(s))ds i.e.

x(t_n+1) = x(t_n) +

∫ t_n+1 t_n

f(s, x(s))ds

オイラー法などの 1 _段階法:X_n+1 = X_n + h · Φ(t_n, X_n) _では、積

分

∫ t_n+1 t_n

f(s, x(s))ds _を h · Φ(t_n, X_n) でどう近似するかということ に相当する。

オイラー法では、この積分を 1 _点 t = t_{n での値を利用して}, _高さ f(t_n, X_n), _幅 h の長方形で近似している.(_区分求積)

この近似を、よりよいものに変える.

代表的なものを２つ紹介： ホイン法、中点法（修正オイラー法)

ホイン法: _{積分を台形で近似}:

∫ _tn+1

t_n

f(s, x(s))ds ≈ h

2 (f(t_n, x(t_n)) + f(t_n+1, x(t_n+1))

≈ h

2 (f(t_n, X_n) + f(t_n+1, X_n+1))

ただし, _今 X_{n+1 の値は未知なので}, これをオイラー法を利用して近似 値 X˜_n+1 = X_n + hf(t_n, X_n) _{で代用する}.

ホイン法のスキーム













X_n+1 = X_n + h · Φ(t_n, X_n) (n = 0, 1, 2, · · · ) where Φ(t_n, X_n) = ¹₂(k₁ + k₂)

k₁ = f(t_n, X_n)

k₂ = f(t_n+1, X_n + hk₁) X₀ = x₀

中点法: 積分を区間の中点での値を用いた長方形で近似:

∫ t_n+1 t_n

f(s, x(s))ds ≈ hf(t_n + h

2 , x(t_n + h

2 ))≈ hf(t_n + h

2 , X_n+1/2)

ただし, 今 X_n+1/2 の値は未知なので, これを刻み幅を h/2 としたオイ ラー法を利用して近似値 X˜_n+1/2 = X_n + ^h₂ f(t_n, X_n) _{で代用する}.

中点法のスキーム













X_n+1 = X_n + h · Φ(t_n, X_n) (n = 0, 1, 2, · · · ) where Φ(t_n, X_n) = k₂

k₁ = f(t_n, X_n)

k₂ = f(t_n + ^h₂, X_n + ^h₂ k₁) X₀ = x₀

以上の２つの方法は 2 _{次精度をもつ} 1 _{段階法となる}.(_各自)

=⇒ ^{先の定理より} 2 次収束することも分かる.

これらの方法のように、[t_n, t_n+1] _{の複数の点を用い}, _{それらの点での} f _の値を (_{何らかの意味での}) 平均をとることで積分の近似精度を上げる ことができることがある. この方向での高次化の方法にルンゲ・クッタ 法が知られている. (これまで紹介した方法は、すべてこの方法の一種と なっている.)

2.4 _{ルンゲ・クッタ法} (Runge-Kutta method)

s _{段のルンゲ・クッタ法}































X_n+1 = X_n + h · Φ(t_n, X_n) (n = 0, 1, 2, · · · )

where Φ(t_n, X_n) = b₁k₁ + b₂k₂ + · · · b_sk_s k₁ = f(t_n, X_n) (c₁ = 0, a₁₁ = 0)

k₂ = f(t_n + c₂h, X_n + ha₂₁k₁)

k₃ = f(t_n + c₃h, X_n + h(a₃₁k₁ + a₃₂k₂)) ...

k_s = f(t_n + c_sh, X_n + h(a_s1k₁ + a_s2k₂+

· · · + a_ss−1ks−1)) X₀ = x₀

※ a_ij, b_i, c_i : ルンゲ・クッタ法の係数パラメータ

Remarks:

• ^{オイラー法は} s = 1, _{ホイン法は} s = 2, a₂₁ = 1, b₁ = b₂ = 1/2, c₂ = 1 のルンゲクッタ法とみることができる.

• s = 1, 2, 3, 4 に対しては、係数パラメータをうまく選ぶことで s

段 s 次のルンゲクッタ法のスキームを作ることができる.

• s ≥ 5 _では、s _段 s 次のルンゲクッタ法のスキームは構成できな いことが示されている. ^*3

*3 陰的ルンゲ・クッタ法であれば可能.

実用上は、次の 4 段のルンゲ・クッタ法のスキームがよく使われている. s = 4, a₂₁ = a₃₂ = 1/2, a₃₁ = a₄₁ = a₄₂ = 0, a₄₃ = 1

b₁ = b₄ = 1/6, b₂ = b₃ = 1/3 c₂ = c₃ = 1/2, c₄ = 1

この設定にすると, 4 次精度を持つことを示すことができる. 4 _段 4 _次 RK _{法のスキーム}





















X_n+1 = X_n + h · Φ(t_n, X_n) (n = 0, 1, 2, · · · )

where Φ(t_n, X_n) = ¹₆(k₁ + 2k₂ + 2k₃ + k₄) k₁ = f(t_n, X_n)

k₂ = f(t_n + ¹₂h, X_n + ¹₂hk₁) k₃ = f(t_n + ¹₂h, X_n + ¹₂hk₂) k₄ = f(t_n + h, X_n + hk₃)

X₀ = x₀

【参考】3 _段 3 _次 RK _{法のスキーム}

















X_n+1 = X_n + h · Φ(t_n, X_n) (n = 0, 1, 2, · · · ) where Φ(t_n, X_n) = ¹₆(k₁ + 4k₂ + k₃)

k₁ = f(t_n, X_n)

k₂ = f(t_n + ¹₂h, X_n + ¹₂hk₁)

k₃ = f(t_n + h, X_n − hk₁ + 2hk₂) X₀ = x₀