(2)常微分方程式のルンゲ・クッタ法による解法

電気 303/ 電情 303 数値解析 (12) 微分方程式の数値解法 (2)

常微分方程式の

ルンゲ・クッタ法による解法

• 解くべき微分方程式はdx

dt =f(x, t)

• 取り扱うのは離散変数法(時刻(t₀, t₁, . . . , t_N) で微分方程式の数値解を計算)

• h_k = t_k+1−t_kをステップ幅という. h_kがk によらず一定のときにはhと書く.

• x(t_k)のことをxkと書く.

• 常微分方程式の数値解法:1段法 と 多段法.

段は英単語stepの訳.

• 1段法の一般形:

explicit(陽的)な公式: ξ_k+1 =ξ_k+h_kΨ(t_k,ξ_k, h_k) implicit(陰的)な公式:

ξ_k+1 =ξ_k+hkΨ(tk, t_k+1,ξ_k,ξ_k+1, hk)

• まず, Euler法より精密なΨの構成法である Runge-Kutta型公式について述べ(おもにex-

plicitな公式),続いて多段法について述べる.

• 参考文献は前回と同じなので再掲しない.

• 以下では, x(t)は dx

dt =f(x, t)の解とする.

• Runge-Kutta型公式は代表的な1段法.

• 平均値の定理から, ∃Φ(tk,xk, h_k), xk+1 = xk+h_kΦ(tk,xk, h_k). Φ(tk,xk, h_k)は未知な ので, 既知のΨ(t,x(t), hk)で近似する.

• Ψ(t,x(t), hk)−Φ(t,x(t), hk)を局所離散化誤 差といい, τ(t, r)と書く.

• ある公式が∃C > 0, ∃ρ > 0, ∀r : 0 < r < ρ,

∀t ∈ [a, b], kτ(t, r)k∞ ≤ Cr^pを満たすとき, その公式をp次精度の公式, あるいはp次の 公式という.

区間の右端では,たとえばx(t+r)−x(t)をx(t)−x(t−r)で読み換えるな どといった,若干の読み換えが必要.

• ステップ幅が一定であるとき(hとする),p次 精度の公式に対し, あるL >0が存在して,

∀k,kxk−ξ_kk ≤ e^L(b−a) L Ch^p

となることが示せる([齊藤], pp. 225–226).

• Runge-Kutta型の公式は, いくつかの点にお けるf(x, t)の値を組み合わせてΨを構成す ることでp次精度(p≥1)の公式を得る方法 の総称.

• 色々なバリエーションがある.

• Heun法はRunge-Kutta法の一種.

• Ψを構成するために関数fの値をs個の点に おいて評価するとき, これをs段数(s-stage) のRunge-Kutta法と呼ぶ([齊藤])

• s段数のRunge-Kutta法を構成する目的は, その公式がp次精度となるようにすること.

• s ≥pである必要があることが示せるが([齊 藤]), s=pとは限らない.

• 一般論の前に, 代表的な公式を紹介する(見 方がわかればよく, 覚える必要はない).

• 差分方程式ξ_k+1=ξ_k+h_kΨ(tk,ξ_k, h_k)にお ける関数Ψをp次精度となるようにしたい.

⊲ 2段数のRunge-Kutta法

• 2段数のRunge-Kutta法は, k1 = f(ξ_k, t_k), k₂ =f(ξ_k+h_kβk₁, t_k+αh_k)とし, Ψ(t_k,ξ_k) = c₁k₁+c₂k₂とする公式である. α,β,c₁,c₂を うまく選ぶと2次精度となる.

• 2段数のexplicitなRunge-Kutta法において, α > 0とし, β = α, c₁ = 1− _2α¹ , c₂ = _2α¹ とすると, 2次精度が達成できることが示せ る. この意味で, 2段数2次精度のexplicitな

Runge-Kutta法には,無限のバリエーション

がある([齊藤], pp. 227–229).

• よく知られた選び方を次ページに示す.

2段数Runge-Kutta法 α β c₁ c₂ Heun法 1 1 ¹₂ ¹₂ 改良Euler法 ¹₂ ¹₂ 0 1

開し,その低次の項が零になるように係数合わせをすること によって得られる. 計算は繁雑ではあるが,機械的にできる.

すべてのRunge-Kutta法を特徴付けることにはあまり意味

がないので,以下ではよく知られたパラメータの選び方のみ をを紹介する. 以下で述べる,より高い段数のRunge-Kutta 法についても,「なぜそんなパラメータを選ぶか」が一見し て明らかでないと思われるが,それらは「辻褄が合うように 係数合わせをして,うまくいったものの例」に過ぎない.

⊲ 4段数のRunge-Kutta法

• 4段数のRunge-Kutta法は, kj =f(ξ_k+h_kP

l<jβ_jlkl, t_k+α_jh_k) (ただし j = 1, . . . ,4,P

l<jβ_jl =α_j)とし, Ψ(tk,ξ_k) = P4

j=1c_jkj(ただしP4

j=1c_j = 1)とする公式で ある. (j = 1のときはP

l<j· · · の項は· · · の 部分にかかわらず零と定義する).

Rungeの1/6公式 α₁ α₂ α₃ α₄ c₁ c₂ c₃ c₄

0 ¹₂ ¹₂ 1 ¹₆ ¹₃ ¹₃ ¹₆ β₂₁= ¹₂, β₃₂= ¹₂, β₄₃= 1

記載のないパラメータは零

この公式は, RK4, RK41,古典的Runge-Kutta法など,色々の名前で呼ばれる.

Kuttaの1/8公式 α₁ α₂ α₃ α₄ c₁ c₂ c₃ c₄

0 ¹₃ ²₃ 1 ¹₈ ³₈ ³₈ ¹₈ β₂₁= ¹₃, β₃₁=−¹₃, β₃₂= 1 β₄₁= 1, β₄₂=−1, β₄₃ = 1

記載のないパラメータは零

• 続いて, Runge-Kutta法の一般形とその表記 法について述べる.

• Ψ(tk,ξ_k) =

s

X

j=1

c_jkj, Ps

j=1c_j = 1 という構 成を考える. kjは後述. j = 1, . . . , sに注意.

Explicit(陽的):P

l<jβ_jl=α_jとしたとき, kj =f(ξ_k+h_kX

l<j

β_jlkl, t_k+α_jh_k), Implicit(陰的):Ps

l=1β_jl=α_jとしたとき, kj =f(ξ_k+h_k

s

X

l=1

β_jlkl, t_k+α_jh_k)

• Euler法は1段数のRunge-Kutta法であると 解釈できる.

• s段数のRunge-Kutta型公式は, s個のパラ メータα₁, . . . , α_s,s個のパラメータc₁, . . . , c_s およびs²個のパラメータ(βil)によって特徴 付けられる. ただし,α_j =Ps

l=1β_jlである.

• Runge-Kutta型公式のパラメータを簡潔に表 示するために, Butcher配列という形式が用 いられる. これは, 次のような表である.

α₁ β₁₁ β₁₂ · · · β_1s α₂ β₂₁ β₂₂ · · · β_2s ... ... ... ... αs β_s1 β_s2 · · · βss

c₁ c₂ · · · cs

• explicitな公式のButcher配列では,右下がり の対角線とその上の要素は零である.

• s = p = 2, s = p = 3, s = p = 4の場合 は,パラメータを含んだRunge-Kutta型公式 の一般形(のButcher配列)を書き下すことが できる([Butcher]).

• s = p = 2とした場合のexplicitな公式の一 般形についてはすでに述べた.

• s=p= 3,s=p= 4については,パラメータ 値に関する場合分けが必要で繁雑なので,こ の講義では述べない([Butcher]参照).

• 以下, Butcher配列の例を示す([齊藤]).

前進Euler法

1 後退Euler法

1 2段数 (s =p= 2) [齊藤]

Heun法 改良Euler法 Crank–Nicolson法 0 0 0

1 1 0

1 2

1 2

0 0 0

1 2

1 2 0 0 1

0 0 0 1 ¹₂ ¹₂

1 2

1 2

Heun法 Kutta法

0 0 0 0

1 4

1

4 0 0

2

3 −²₉ ⁸₉ 0

1

4 0 ³₄

0 0 0 0

1 2

1

2 0 0

1 −1 2 0

1 6

2 3

1 6

0 0 0 0 0

1 2

1

2 0 0 0

1

2 0 ¹₂ 0 0 1 0 0 1 0

1 6

1 3

1 3

1 6

0 0 0 0 0

1 3

1

3 0 0 0

2

3 −¹₃ 1 0 0 1 1 −1 1 0

1 8

3 8

3 8

1 8

• B= (βjl)1≤j≤s,1≤l≤sをButcher配列の中央部 にある行列, α= (α₁, . . . , α_s)^T を左側にある ベクトル, c = (c1, . . . , c_s)を下側にあるベク トルとし, 対応するButcher配列を(B,α,c) と表記する.

• explicit な Runge-Kutta法に関し, p ≥ 5な らp段数ではp次精度は実現不可,p≥7なら p+ 1段数ではp次精度は実現不可,p≥8な らp+ 2段数ではp次精度は実現不可という 事実が知られている([Hairer et al], Vol. 1).

• 上記は指定された下限より大きいsに対して 公式の存在を保証するわけではない.

• 5次精度については6段数, 6次精度について は7段数, 8次精度については11段数, 10次 精度については17段数の公式が知られてい る([Hairer et al], Vol. 1).

• 浮動小数点演算の誤差のため, 高次精度の公 式がつねに実用的とは限らない.

• 今までの議論では, ステップ幅の決め方は議 論されていなかった.

• 後に述べるように, Butcher配列(B,α,c)を 持つp+1次精度の公式と, Butcher配列(B,α,c) を持つp次精度の公式が存在すれば, これを 使って誤差を見積ることができる.

• いくつかの(p, p+ 1)の組み合わせについて, 上述のようなButcher配列が存在することが 知られている.

• このようになっているとき, s段数p+ 1次精 度の公式にp次精度の公式が埋め込まれてい ると表現する([齊藤]).

• 上記のようなRunge-Kutta法の組み合わせ を埋め込み型Runge-Kutta法と呼ぶ([齊藤]).

その用途は, 数値解に含まれる誤差の見積り である.

• (B,α,c)に対応するRunge-Kutta法が生成 する数値解を(ξ_k), (B,α,c)に対応するRunge- Kutta法が生成する数値解を(ξ_k)とする.

• 詳細は省くが,kξ_k−ξ_kkの値から,数値解に 含まれる誤差を見積ることができる(ただし, 近似を含むので,精密な評価ではない).

• 誤差の見積りが一定以下になるまでステップ 幅をどんどん小さくしてゆく手法をステップ 幅の自動調整という.

• ステップ幅の自動調整により, 数値解の正し さを, ある程度保証することができる(近似 が含まれるので厳密ではない). 詳細につい ては[齊藤]を参照.

• Scilabに実装されているRKF45は埋め込み 型Runge-Kutta法で, 6段数5次精度の公式 に4次精度の公式が埋め込まれている.

• この講義では今までほとんど述べなかったが, implicitなRunge-Kutta法は硬い問題(後述) やdifferential-algebraic equationを解くのに 用いられる. ただし, 内部反復を必要とする ため, 必ずしも使いやすいわけではない.

• implicitなRunge-Kutta法に興味がある者は [Hairer et el.]や[Butcher]などを参照せよ.

• 線形微分方程式 _dtx =Axにおいて, 行列Aの 固有値がすべて相異なり,それら実部がすべて負 であるものとする.

• Aの固有値のうち実部が最小のものをλm, 最大 のものをλ_M とする.

• |ReλM|/|Reλm|を硬度比という.

• 硬度比が大きい問題を硬い問題という. 硬い問題 は数値的に解きにくいことが知られている.

• 多段法は, ξ_k+1を構成するために, ξ_kだけで なく, 過去の数値解の系列(ξ_k−L, . . . ,ξ_k)に 関連した情報を使う方法の総称(復習).

• Runge–Kutta法では,複数の点での関数値の 線形結合により高精度の公式を得ていたが, 多段法は, (ξk−L, . . . ,ξ_k)を使うことで類似し た効果を得ることを意図している.

• 記法の便宜上, ξ_k+Lを(ξ_k, . . . ,ξ_k+L−1)に関 連した情報を使って構成するという形で問題 を考え直す.

• 簡単のために, この講義では, 多段法ではス テップ幅はhに固定されているものとする

(ステップ幅可変の場合については[Hairer et

al.], Vol. 1を参照).

• Runge–Kutta型公式と比較すると,f(ξ_k, t_k), . . . , f(ξ_k+L, t_k+L) の線形結合で近似関数を 構成するのが簡単:

ξ_k+L =ξ_k+h PL

l=0β_l⁽⁰⁾f(ξ_k+l, t_k+l)

• β_l⁽⁰⁾ (0≤l ≤L) はパラメータとして残す.

• ξ_kを左辺に移項して・・・

• ξ_k+L−ξ_k =h PL

l=0β_l⁽⁰⁾f(ξ_k+l, t_k+l)

• 左辺を{ξ_k+l}0≤l≤L をすべて使う形に拡張す

ると,線形多段法の一般形が得られる(次式).

L

X

l=0

α_lξ_k+l =h

L

X

l=0

β_lf(ξ_k+l, t_k+l)

!

• 先の式の{αl}0≤l≤L, {βl}0≤l≤Lはパラメータ で, これらの取り方によって色々な公式が得 られる.

• 上述一般形の中で, 等号の右辺にξ_k+Lが含 まれるものをimplicit(陰的)な公式,右辺に ξ_k+Lが含まれないものをexplicit(陽的)な 公式という. これらの特徴は1段法と同様.

• 多段法を非線形にまで(形式的に)拡張するこ とも可能, この形で一般形が述べられている 本もあるので([山本]), 念のため書いておく.

PL

l=0αlξ_k+l=hΨ tk, . . . , t_k+L,ξ_k, . . . ,ξ_k+L;h

• 以下で,代表的な多段法をいくつか紹介する が, その前に注意点を述べる.

• 段数Lの多段法を使うときには, ξ₀, . . . , ξ_L−1を(1段 法等によって)事前に計算する.

• (線形)多段法のメリットは{αl}_0≤l≤Lや{βl}_0≤l≤Lをう まく選ぶと高次精度の公式が得られること.一方で,安 定性などに関して注意すべき点も多い.

• 多段法が安定で収束するための必要十分条件として,以 下のものが知られている([森]): (1)多項式α₀+α₁z+

· · ·+α_Lz^Lの根の絶対値がすべて1以下で,絶対値1の

根は単根, (2)PL

l=0α_l= 0, (3)Pl

l=0lα_l−PL

l=0β_l= 0.

• 多段法において, explicitな公式で仮の解を求

め, それをimplicitな公式に代入してから修

正する方法を, 予測子修正子法という. また, この方法で使われるexplicitな公式を予測子, implicitな公式を修正子という.

• explicitな解法とimplicitな解法の組み合わ せは色々.

• 微分方程式^dx

dt =f(x, t), x(t_k+L−1) =x_k+L−1 は次の積分方程式と等価.

xk+L−xk+L−1 =Rtk+L

tk+L−1 f(ϕ(τ, tk,xk), τ)dτ

• しかし,この積分はϕ(τ, tk,xk)が未知だから 計算できない.

• そこで,ϕ(τ, tk,xk)を時間τに関する1変数 関数と見倣し, Lagrange補間多項式で近似す ることを考える.

• 近似に(t_k,xk), . . . ,(t_k+L−1,x_k+L−1)を使うと, xk+L−xk+L−1 ≃ PL−1

l=0 β_lf(xk+l, t_l)という 形に, (tk,xk), . . . ,(tk+L,xk+L)を使うと,xk+L− x_k+L−1 ≃PL

l=0β_lf(xk+l, t_l)という形になる.

• 前ページのxjをξ_jで置き換えると(k ≤j ≤ k+L), 多段法による常微分方程式の数値解 法の公式が得られる. これを, Adams型公 式という.

• 右辺に(ξ_k+L, t_k+L)を含まない公式をAdams–

Bashforth公式,これを含む公式をAdams–

Moulton公式という.

• Adams–Bashforth 公式 はexplicit, Adams–

Moulton 公式はimplicit である.

• 右辺の係数{βl}0≤l≤Lは, Lagrange補間多項 式を積分することで得られる. この係数はL の値だけで決まり, 問題によって変わるわけ ではないので, 事前に計算するか, 文献等の 値をそのまま用いればよい.

例として, L=2の場合のAdams–Bashforth公式を計算してみる. tk

においてf(ξk, tk),tk+1においてf(ξk+1, tk+ 1)という値を取るLa- grange補間多項式は, t−tk+1

tk−tk+1

f(ξk, tk) + t−tk

tk+1−tk

f(ξk+1, tk+1)で あり, tk+1 −tk = hであることに注意すると,

Z tk+2

tk+1

t−tk+1

tk−tk+1

=

−1

2h(t−tk+1)²

tk+2 tk+1=−h

2, Z tk+2

tk+1

t−tk

tk+1−tk

= 1

2h(t−tk)²

tk+2 tk+1=3h

2 である. よって, L=2の場合のAdams–Bashforth公式は,

ξ_k+2−ξ_k+1=h

−1

2f(ξk, tk) +3

2f(ξk+1, tk+1)

.

• BDF法はBackward Difference Formula法の 略. この方法はimplicitな解法の一種.

• BDF法は硬い系向けの解法として広く使わ れている([三井]).

• ξ_k+Lを(ξ_k, . . . ,ξ_k+L−1)に関連した情報から 構成するという多段法の問題設定に戻る.

• ξ_k+Lを未定のままにしたまま,

(tk,ξ_k), (tk+1,ξ_k+1), . . . , (tk+L,ξ_k+L)を通る

L次のLagrange補間多項式を構成する. こ

れをp(t)とする.

• ステップ幅が等間隔だったから,t_k+L−1 =t_k+l− h, t_k+L−2 =t_k+l−2h, . . .等に注意する.

• p(t)が微分方程式の解の近似関数となるよう にξ_k+Lを決めたい.

• 決め方はひと通りというわけではないが,p(t) の導関数のt=t_k+Lにおける値が微分方程式

dx

dt =f(x, t)の右辺と一致するように決める というのはひとつの考え方である.

• d dtp(t)

t=tk+L

=f(ξ_k+L, t)となるようにξ_k+L を定める解法をBDF法という.

• Lagrange補間多項式を微分し, t=t_k+Lにお ける値を評価することは,代数的にできる(ス テップ幅が等間隔であることに注意).

• いくつかのLに対するBDF法は・・・

L= 1 : −ξ_k+ξ_k+1=hf(ξ_k+1, tk+1) L= 2 : ¹₂ξ_k−2ξk+1+³₂ξ_k+2=hf(ξk+2, tk+2)

L= 3 : −¹₃ξ_k+³₂ξ_k+1−3ξk+2+¹¹₆ξ_k+3=hf(ξk+3, tk+3)

• BDF法はL >6に対して不安定となること

が知られている.

微分方程式 d

dt x₁ x₂

!

= 0 −1 1 0

! x₁ x₂

!

, x₁(0) x₂(0)

!

= 1

0

!

を解く. 解はx₁(t) = cost, x₂(t) = sint. この微分 方程式は硬い系ではないので, BDF法などを使う 利点はないが,すべての解法を試してみる.

dx=[0 -1;1 0]*x;

endfunction t0=0;

//時間の刻みは問題に合わせて要調整. t=0:.1:100;

x0=[1;0];

y=ode(x0,t0,t,f); //デフォルト

y=ode("adams",x0,t0,t,f); //Adams法 y=ode("stiff",x0,t0,t,f); //BDF法

y=ode("rk",x0,t0,t,f); //Runge-Kutta法 y=ode("rkf",x0,t0,t,f); //RKF45

デフォルト解法はAdams型予測子・修正子法とBDF 法の自動切り換え

通り.

解法 誤差 計算時間

デフォルト 5.33×10⁻⁵ 0.031 Adams 5.78×10⁻⁵ 0.111 BDF法 7.61×10⁻⁴ 0.062 Runge-Kutta法 9.10×10⁻⁹ 0.811 RKF45 1.66×10⁻⁵ 0.062

以下は文献[三井]に掲載されている硬い系である.

d dt

x₁ x₂

!

= −2 1

1998 −1999

! x₁ x₂

!

+ −cost

1999 cost−sint

! ,

初期値をx₁(0) = 1, x₂(0) = 2とする. 解はx₁(t) = e^−t, x₂(t) =e^−t+ costである(真の解はMathematicaによって求 めた).

dx=[-2 1;1998 -1999]*x+..

[-cos(t);1999*cos(t)-sin(t)];

endfunction

//注意: 記号.. は「次の行に続く」という意味 t0=0;

t=0:.5:100;

x0=[1;2];

y=ode(x0,t0,t,f); //デフォルト

y=ode("adams",x0,t0,t,f); //Adams法 y=ode("stiff",x0,t0,t,f); //BDF法

y=ode("rk",x0,t0,t,f); //Runge-Kutta法 y=ode("rkf",x0,t0,t,f); //RKF45

通り.

解法 誤差 計算時間

デフォルト 3.3×10⁻⁷ 0.031 Adams 2.1×10⁻⁶ 0.484 BDF法 4.9×10⁻⁷ 0.109 Runge-Kutta法 7.5×10⁻¹ 29.547 RKF45 8.9×10⁻⁴ 2.028