d dt A B C = A B C d dt x = Ax, A 0 B 0 C 0 = mm 0 mm 0 mm AP = PΛ P AP = Λ P A = ΛP P d dt x = P Ax d dt (P x) = Λ(P x) d dt P x =

第

₃

章 常微分方程式の数値積分法

本章の目標は MATLAB に装備されている数値積分法の仕組みを理解すること にある。Runge-Kutta 法の次数、Butcher 表の見かた、陽的、陰的、埋め込み型、 可変ステップ、硬い方程式といった概念を理解し、適切な数値積分アルゴリズム を選択する判断力を身に付けることを目指す。

3.1

Taylor

展開の式を使わないのはなぜか？

数値積分は Taylor 展開と比べてどこがいいのか？ y(x0+ h) = y(x0) + h· y0(x0) + h2 2!y 00_(x 0) +· · · Taylor 展開では、高階微分の項が複雑になる。数値積分公式は、微分項の計算 なしに積分ができるすぐれものの近似式である。

3.2

分類

1 段階法 多段階法

陽公式 Euler, Runge-Kutta Adams

陰公式 Implicit Euler, Implicit Runge-Kutta Gear

陽公式 左辺 yn+1 、右辺 yn など (n+1 未満) 陰公式 左辺 yn+1 右辺も yn+1。 まず陽公式で yn+1を予測し、それを陰公式に代入して修正する、といった使い方 をする。

3.3

紙と鉛筆の限界を知る

線形なら対角化すれば解析的に解ける。      d dt[A] = −2k[A] d dt[B] = 2k[A]− k[B] d_{[C] = k[B]}

d dt    A B C    =    −2 0 0 2 −1 0 0 1 0       A B C    ,    A0 B0 C0    =    1 mM 0 mM 0 mM    d dtx = Ax AP = PΛ P−1AP = Λ P−1A = ΛP−1 P−1 d dtx = P −1_Ax d dt(P −1_{x) = Λ(P}−1_x) d dtP −1_{x =}    0 0 0 0 −1 0 0 0 −2    P−1x P−1x =    C1 C2exp(−t) C3exp(−2t)    =    1 1 1 2 1 0 1 0 0       A B C    よって A + B + C = C1 2A + B = C2exp(−t) A = C3exp(−2t)

   A0 B0 C0    =    1 mM 0 mM 0 mM    より C3 = 1、C2 = 2、C1 = 1 とわかる。 A + B + C = 1 2A + B = 2 exp(−t) A = exp(−2t) この連立方程式を解いて、    A B C    =    exp(−2t) 2 exp(−t) − 2 exp(−2t) 1− 2 exp(−t) + exp(−2t)    mM

3.4

1

段階法

3.4.1

陽的な

1

段階法

: Euler

法

最も単純な数値積分アルゴリズムは Euler 法です。定義としては、Taylor 展開 の式 y(x0+ h) = y(x0) + h· y0(x0) + h2 2!y 00_(x 0) +· · · から 2 次以降の項を取り除いた次の式が Euler 法です。 y(x0+ h) = y(x0) + h· y0(x0) Euler 法の近似精度は、Taylor 展開式の 1 次までの精度と同等です。 これだと一般的すぎるので、生物の場合にあてはめて考えてみましょう。時刻 t における物質 S の濃度を S(t) とし、そこから微小時間 ∆t だけ進んだ時刻 t + ∆t における物質 S の濃度を S(t + ∆t) で表します。これを上の Euler 法の定義に当て はめると、 S(t + ∆t) = S(t) + ∆tS0(t) となります。 S0_{(t) というのは S の時間微分}dS _{のことですので、}

S(t + ∆t) = S(t) + ∆tdS dt という式が得られます。dS dt は反応速度式で表されます。ここではテストケース として dS dt =−k · S としてみましょう。 S(t + ∆t) = S(t)− ∆t · k · S これを増分の形にすると、 ∆S(t) = S(t + ∆t)− S(t) = −∆t · k · S よって Euler 法を用いる場合、微小時間 ∆t の間に S がどの程度増減するかは −∆t · k · S で計算できることがわかりました。

3.4.2

演習

1. グルコース濃度の変動に関する常微分方程式 d[Glucose]_dt = −k[Glucose] をオ イラー法で解きなさい。積分ステップ幅 ∆t = 1.0sec とし、t = 0 ∼ 5.0

sec の範囲について解くこと。反応定数 k = 0.2sec−1、t = 0 sec のとき

[Glucose]=1.0mM とする。 2. (1) の常微分方程式の解析解を求め、t = 5.0 sec 時点までの数値解との時系 列平均誤差 ε = 1 n n ∑ i=1 |Xanalytical i − Xinumerical| X_ianalytical を計算しなさい。 3. 積分ステップ幅 ∆t を 1.0 sec より大きくした場合と小さくした場合では、解 析解との時系列平均誤差がどのように変化するか仮説を立てなさい。理由 も述べること。 4. 仮説を検証しなさい。

3.4.3

陰的な

1

段階法

:

陰的

Euler

法

3.4.4

ステップ刻み幅を決めるには「絶対安定」 の概念が必要

3.5

多段階法

多段階法とは以下のような一般式で書ける数値積分法のことである。 α0yn+ α1yn−1+· · · + αkyn−k = h(β0fn+ β1fn−1+· · · + βmfn−m) ynは tn時点における y の値を意味する。また、yn−k = y(tn− kh)、f = y0で ある。

3.5.1

利点

Euler と比べてどこがいいのか？→高次の項を入れて局所離散誤差を小さくで きる。RK と異なり、導関数はステップあたり 1 回計算すれば済む。

3.5.2

陽的な多段階法

: Adams-Bashforth

法

陽的な多段階法の代表例である。あまり使ったことはないが、Gear 法を理解す る目的で足がかりとなるので紹介する。 Adams 法は yn= α1yn−1+ h(β1fn−1+ β2fn−2) という形の多段階法である。係数 α, β を Taylor 展開から定める。まず上の式を 書き換える。

y(tn) = α1y(tn−1) + h(β1f (tn−1, y(tn−1)) + β2f (tn−2, y(tn−2)))

tn−1 = t− h より、Taylor 展開する。 y(tn) ' α1 { y(tn)− h d dty(tn) + h2 2 d2 dt2y(tn) } + hβ1 { f (tn)− h d dtf (tn) + h2 2 d2 dt2f (tn) } + hβ2 { f (tn)− 2h d dtf (tn) + 4h2 2 d2 dt2f (tn) }

h の次数が同じ項ごとに整理する。 h0 _{: y(t} n) = α1y(tn) h1 _{: 0 =} −α 1hy0(tn) + hβ1y0(tn) + hβ2y0(tn) h2 : 0 = α1h 2 2 y00(tn)− h 2_β 1y00(tn)− 2h2β2y00(tn) これで、α, β に関する連立方程式を立てることができる。      1 = α1 0 = −α1+ β1+ β2 0 = 1 2α1− β1− 2β2 これを解くと、2 次の Adams-Bashforth 法の係数は α1 = 1, β1 = 3₂, β2 = 1₂ と求 まる。よって Adams-Bashforth 法は yn= yn−1+ h ( 3 2fn−1+ 1 2fn−2 ) と書ける。化学反応系 dS dt =−k · S にあてはめると、 ∆S = ∆t ( −3 2k· St=t−∆t− 1 2k· St=t−2∆t ) となる。

3.5.3

陰的な多段階法

: Adams-Moulton

法

3.5.4

予測子修正子法

予測子修正子 (predictor-corrector) 法とは、陽公式と陰公式を組み合わせて両者 のいいとこどりを図る方法のことである。Adams-Bashforth 法と Adams-Moulton 法を組み合わせる例が有名である。MATLAB の ode113 がこれにあたる。 長所 短所 陽公式 過去の値を使って、次のステップの値が出せる 安定性で陰公式に劣る 陰公式 安定性が優れている 現時点の値がないと計算できない 陽公式で予測値を計算し (予測子)、予測値をもとに現時点の値を陰公式で計算 (修正子)

3.6

安定性

3.6.1

絶対安定

線形テスト問題 多段階法 安定性多項式の根がすべて 1 未満。 安定性多項式 Runge-Kutta 法

3.7

Runge-Kutta

法

1 段階法に分類される数値積分法の中では、Runge-Kutta 法とその亜種がシス テム生物学の数理モデルに用いられる。

3.7.1

利点

LM と同じく、高次の項を入れて局所離散誤差を小さくしたかったことが発明 の動機。LM は段数分だけ前のステップの値を用意しなければならないが、RK は 1 ステップ前の値から計算できる。しかし、次数分だけ導関数を計算する必要が ある点は劣る。安定領域は Adams 系より広い？

3.7.2

s

段の陽的

Runge-Kutta

法

                       k1 = f (x0, y0) k2 = f (x0+ c2h, y0+ ha21k1) k3 = f (x0+ c3h, y0+ h(a31k1+ a32k2)) .. . ks = f (x0+ csh, y0+ h(as1k1+ as,s−1ks−1)) y1 = y0+ h(b1k1 +· · · + bsks) ただし

c2 = a21 c3 = a31+ a32 .. . cs = as1+· · · + as,s−1 これを s 段の陽的 Runge-Kutta 法と呼ぶ。

3.7.3

Butcher

表

s 段の陽的 Runge-Kutta 法の係数を以下の形式で表記したものを Butcher 表と 呼ぶ。今後、さまざまな数値積分法のアルゴリズムを説明する際に標準的に用い られる形式なので、慣れて欲しい。 k1の係数 0 k2 c2 a21 k3 c3 a31 a32 .. . ... ... . .. ks cs as1 as2 · · · as,s−1 y1 b1 b2 · · · bs−1 bs

3.7.4

Runge-Kutta

法の係数を決める

(RK2)

2 段 2 次の Runge-Kutta 法は以下のようになる。          k1 = f (x0, y0) k2 = f (x0+h₂, y0+h₂k1) y1 = y0+ hk2

3.7.5

Runge-Kutta

法の次数

一般に、 ||y(x0+ h)− y1|| ≤ Khp+1

となる Runge-Kutta 公式のことを「p 次の Runge-Kutta 法」と呼ぶ。Taylor 級 数とは、p 次の項まで一致する。

3.7.6

Taylor

展開との精度比較

まず、ステップの長さを半分にした Euler 法を書く y1 = y0+ hf ( x0 + h 2, y0+ h 2f0 ) 2 変数関数の Taylor 展開公式を用いてこれを展開する。 f(x0+ h₂, y0+ h₂f0 ) = f (x0, y0) + ( h 2 ∂ ∂x + h 2f0 ∂ ∂y ) f (x0, y0) + _2!1 ( h 2 ∂ ∂x+ h 2f0 ∂ ∂y )2 f (x0, y0) +· · · = f (x0, y0) + h₂ (fx+ f fy) (x0, y0) + 1₂ ( h2 4 ∂2 ∂x2 + h2 2 f0 ∂ ∂x ∂ ∂y+ h2 4 f 2 0 ∂2 ∂y2 ) f (x0, y0) +· · · = f (x0, y0) + h₂ (fx+ f fy) (x0, y0) + h₈2 (fxx+ 2f fxy + f2fyy) (x0, y0) +· · · ∴ y1 = y0+ hf (x0, y0) + h 2 2 (fx+ fyf ) (x0, y0) + h₈3 (fxx+ 2fxyf + fyyf2) (x0, y0) +· · · これと厳密解の精度を比較するため、厳密解の Taylor 展開を行う。 y(x0+ h) = y0 + hy0(x0, y0) + h 2 2 y00(x0, y0) + h3 6 y000(x0, y0) +· · · = y0 + hf (x0, y0) + h 2 2(fx+ fyf )(x0, y0) + h₆3(fxx+ 2fxyf + fyyf2+ fxfy+ fy2f )(x0, y0) +· · · 数値解と厳密解の差をとる。1 階微分の項までは両者一致するので、2 階の項の 差がものを言う。 y(x0+ h)− y1 = h3 24{(fxx+ 2fxyf + fyyf 2_{+ 4(f} xfy + fy2f )}(x0, y0) +· · · よって、 ||y(x0+ h)− y1|| ≤ Kh3 となるので、前述の 2 段の Runge-Kutta 法の式は 2 次の精度である。

3.7.7

素材？

dy dx = f (x, y) df = ∂f ∂xdx + ∂f ∂ydy df dx = ∂f ∂x + ∂f ∂y dy dx = y00 y00 = fx+ fyf dy00 = ∂y 00 ∂xdx + ∂y00 ∂y dy dy00 dx = ∂y00 ∂x + ∂y00 ∂y dy dx = ∂ ∂x(fx+ fyf ) + ∂ ∂y(fx+ fyf ) dy dx = {fxx+ (fxyf + fyfx)} + {fxy + (fyyf + fy2)}f y000 = fxx+ 2fxyf + fyyf2+ fxfy+ fy2f

3.8

ステップ幅の変更

3.8.1

埋め込み型公式

「公式名 p(ˆp)」は次数 p、誤差評価に用いる式の次数 ˆp であることを示す (例:

Dormand-Prince 5(4), Bogacki-Shampine 3(2) がそれぞれ ode45, ode23 に対応)。

また、「Dormand-Prince 法は陽的な Runge-Kutta(5,4) の式である」のようにも書 かれる。 埋め込み型＝可変ステップ

3.8.2

RK3(2)

の係数を求める

RK3 の係数 まず、3 段 3 次公式の係数を求める。

         k1 = f (xn, yn) (1) k2 = f (xn+ c2h, yn+ ha21k1) (2) k3 = f (xn+ c2h, yn+ h(a31k1+ a32k2)) (3) yn+1 = yn+ h(b1k1+ b2k2 + b3k3) (4) k2をテイラー展開する。今ここで求めたいのは 3 次公式なので、h3が掛かる項 までテイラー展開して、厳密解と比較しなければならない。k2には上の (4) 式で h が掛かるので、h2の項まで展開すれば事足りる。 k2 = f (xn, yn) + _1!h ( c2_∂x∂ + a21k1_∂y∂ ) f + h_2!2 ( c2_∂x∂ + a21k1_∂y∂ )2 f + O(h3) = f (xn, yn) + h(c2fx+ a21f fy) + h 2 2 (c 2 2fxx + 2c2a21f fxy + a221f2fyy) + O(h3) k3も同様にテイラー展開する。 k3 = f (xn, yn) + _1!h ( c3_∂x∂ + (a31k1+ a32k2)_∂y∂ ) f + h_2!2 ( c3_∂x∂ + (a31k1+ a32k2)_∂y∂ )2 f + O(h3) = f (xn, yn) + h(c3fx+ a31f fy+ a32k2fy) + h₂2{c2 3fxx+ 2c3(a31f + a32k2)fxy+ (a31f + a32k2)2fyy+ O(h3) k2が掛かっている項を f で書き直しておきたい。そこで a32k2を a32k2 = a32 { f + h(c2fx+ a21f fy) + O(h2) } あるいは a32k2 = a32(f + O(h)) と表す。これらを項ごとの h の次数に応じて k3のテイラー展開式に代入する。 k3 = f (xn, yn) + h[c3fx+ a31f fy+ a32{f + h(c2fx+ a21f fy) + O(h2)} fy] + h2 2 [c 2 3fxx+ 2c3{a31f + a32(f + O(h))}fxy +{a31f + a32(f + O(h))}2fyy] さてこれで RK 公式のテイラー展開が完了した。次にやることは厳密解もテイ ラー展開し、h の次数ごとに各項を RK 公式と厳密解とで比較することである。厳 密解のテイラー展開は、次の通りである。

y(x + h) = y(x) + hdy dx + h2 2! d2_y dx2 + h3 3! d3_y dx3 + O(h 4 ) y(x + h)− y(x) = hf +h 2 2 (fx+ fyf ) + h3 6 (fxx+ 2fxyf + fyyf 2_{+ f} yfx+ fy2f ) + O(h 4₎ d2_y dx = fx+ fyf となるのは、全微分の公式 df dx = ∂f ∂x + ∂f ∂y df dx による。 さて、RK 公式と厳密解のテイラー級数を h の次数ごとに比較しよう。 厳密解 RK 公式 1 次 : hf = h(b1+ b2+ b3)f 2 次 : h 2 2 (fx+ fyf ) = b2h 2 (c2fx+ a21f fy) + b3h2(c3fx+ a31f fy+ a32f fy) 3 次 : h 3 6 (fxx+ 2fxyf + fyyf 2 + fyfx+ fy2f ) = 1 2b2h 3 (c2₂fxx+ 2c2a21f fxy+ a221f 2 fyy) + b3h3a32(c2fx+ a21f fy)fy + 1 2b3h 3_{c2 3fxx +2c3a31f fxy + 2c3a32f fxy +(a31+ a32)2f2fyy} h が 1 次の項の比較より、 b1+ b2+ b3 = 1 次いで、h2_{の項を比較すると、} fxについて : h2 2 fx = b2h 2_c 2fx+ b3h2c3fx 1 2 = b2c2+ b3c3 f fyについて : h2 2 fyf = b2h 2_a 21f fy+ b3h2(a31+ a32)f fy 1 2 = b2a21+ b3(a31+ a32) 最後に h3_{の項を比較する。}

fxxについて : h3 6 fxx = 1 2b2h 3_c2 2fxx+ 1 2b3h 3_c2 3fxx 1 6 = 1 2b2c 2 2+ 1 2b3c 2 3 f fxyについて : h3 6 (2fxyf ) = 1 2b2h 3_(2c 2a21f fxy) + 1 2b3h 3_(2c 3a31f fxy + 2c3a32f fxy) 1 3 = b2c2a21+ b3c3(a31+ a32) fyyf2について : h3 6 fyyf 2 ₌ 1 2b2h 3_a2 21f 2_f yy+ 1 2b3h 3_(a 31+ a32)2f2fyy 1 6 = 1 2b2a 2 21+ 1 2b3(a31+ a32) 2 fyfxについて : h3 6 fyfx = b3h 3_a 32c2fxfy 1 6 = b3a32c2 f_y2f について : h 3 6 f 2 yf = b3h3a32a21f fy2 1 6 = b3a32a21 これで係数間の条件式が次のように揃った。                          c2 = a21 c3 = a31+ a32 b1+ b2 + b3 = 1 b2c2+ b3c3 = 1 2 b2c22+ b3c23 = 1 3 b3a32c22 = 1 6 求めたい係数 8 個 (a21, a31, a32, b1, b2, b3, c2, c3) に対して条件式 6 本なので、不定 パラメータを c2 = u, c3 = v のようにおく。 { b2u + b3v = 1₂ b2u2 + b3v2 = 1₃ より { b2u2+ b3uv = 1₂u b2u2+ b3v2 = 1₃

(uv− v2)b3 = 1 2u− 1 3 b3 = 1 2u− 1 3 uv− v2 = 3u− 2 6v(u− v) また同様に { b2uv + b3v2 = 1₂v b2u2+ b3v2 = 1₃ より (uv− u2)b2 = 1 2v− 1 3 b2 = 3v− 2 6u(v− u) また、b3a32c2 = 1 6より、 3u− 2 6v(u− v)a32u = 1 6 a32 = v(u− v) u(3u− 2) 残りの変数は次のように書ける。 b1 = 1− b2− b3 a21 = c2 = u a31 = c3− a32= v− a32 埋め込み型公式 (RK2) の係数決定 埋め込み型公式とは、Runge-Kutta 式の係数 b1, b2, b3,· · · の値だけを変えた式 ˆ yn+1 = yn+ h(ˆb1k1+· · · + ˆbsks) のことである。通常、埋め込み型公式は Runge-Kutta 式より 1 次多いか少ない次数にしておき、両者の差がユーザー定義の誤差

限界値 εtolを越えないようにステップ幅 h を調節する。これを数式で書くと以下 のようになる。 |y1i− ˆy1i| ≤ εtol (中略) RK3(2) 公式では ˆb3 = 0 なので、 ˆ_b 1 = 1− ˆb2 b2u + b3v = 1₂ に ˆb3 = 0 を代入して、 ˆ_b 2u = 1₂ ˆ_{b2 =} 1 2u = 1 2c2 ˆ_{b1 = 1−} 1 2c2 ルンゲの 2 次公式の係数をこの埋め込み型公式 (RK2) にあてはめる。ルンゲの 2 次公式のブッチャー表は 0 1 2 1 2 0 1 であるから、c2 = 1₂, a21= 1₂, ˆb1 = 0, ˆb2 = 1 が得られ、RK3 のブッチャー表は 0 1 2 1 2 c3 a31 a32 b1 b2 b3 となる。この形の RK3 には、たとえば Ralston の公式がある (REF.FIXME)。 0 1 2 1 2 3 4 0 3 4 2 9 1 3 4 9 よって、Ralston の公式にしたい場合、a31 = 0, a32 = 3₄, b1 = 2₉, b2 = 1₃, b3 = 4 9, c3 = 3 4 とすればよい。Ralston の式以外にも RK3 公式は複数提案されている。 これらは誤差係数の議論から導出されているが、本書の目的である「プライマー」 としての役割を越えるため取り扱わない。(REF.FIXME) を参照されたい。

3.8.3

ステップ幅変更アルゴリズム

|y1i− ˆy1i| ≤ εtol まず、 err = v u u t 1 n n ∑ i=1 ( y1i− ˆy1i εtol )2 err が 1 を上回らないようにしたい。 err ' Chq+1 1 ' Chq+1_opt より err ' (_hh opt) q+1 hopt ' h(_err1 ) 1 q+1 許容誤算限界値を下回るステップ幅 hoptを得る。ステップ幅を変更するアルゴ リズムは以下の通りである。 1. err が 1 以下であれば、ステップ幅 h は変更せずに用いる。 2. err が 1 を上回ったら、hoptを計算し、新しい h として採用する。

3.9

硬い方程式

3.9.1

Gear

法または後退微分法

(BDF: Backward Differential

Formula)

陰的な多段階法の代表例であり、後述する「硬い方程式」に強い。

α0yn+ α1yn−1+· · · + αkyn−k = hfn

という形の多段階法である。陽的 Adams 法同様、Taylor 展開で係数を求める。 Gear の”Numerical initial value problems” pp.214 に”stiffly stable method”とし て登場する。

L.H.S. ' α0y(tn) + α1 { y(tn)− h d dty(tn) + h2 2 d2 dt2y(tn) } + α2 { y(tn)− 2h d dty(tn) + 4h2 2 d2 dt2y(tn) } R.H.S. = hd dty(tn) h の次数が同じ項ごとに整理する。 h0 _{: α}

0y(tn) + α1y(tn) + α2y(tn) = 0

h1 _{: −α} 1hy0(tn)− 2hα2y0(tn) = hy0(tn) h2 : α1h 2 2 y00(tn) + α22h 2_y00_(t n) = 0      α0 + α1+ α2 = 0 −α1− 2α2 = 1 1 2α1+ 2α2 = 0 これを解くと、2 次の Gear 法の係数は α0 = 3₂, α1 =−2, α2 = 1₂ と求まる。 よって 2 次の Gear 法は 3 2yn− 2yn−1+ 1 2yn−2 = hfn 化学反応系 dS dt =−k · S にあてはめると、 3 2St− 2St−∆t+ 1 2St−2∆t =−∆t · k · St となる。

3.9.2

硬い方程式

(stiff equation)

の安定性

硬度比 (stiffness ratio) 連立常微分方程式のヤコビ行列の固有値から、その系の「硬さ」を推定できる。 硬度比 = | 実部最大の固有値 | | 実部最小の固有値 |

A 安定 離散変数法の絶対安定領域が、複素平面の左半面を含むこと。1 段階、2 段階の BDF は A 安定である。 硬安定 (stiffly stable) • 左半面における原点近傍 • 左半面を走る虚軸の平行線のうち、その左側全体が絶対安定領域になっって いるものが存在する 1∼6 段階の BDF は硬安定である。

3.10

Further reading

三井斌友「常微分方程式の数値解法」岩波書店 Gear, Springer’s yellow book http://www.scholarpedia.org/article/Backward_differentiation_formulas by Gear