10. 常微分方程式 (4) 実用ルーチンのベンチマーク

実用ルーチンの選択

XLPackに収録されている常微分方程式用のプログラムはすべてステップ幅自動調節機能を持っている. 密出力の必要性や特記事項を参考に選択するとよい. スティフではない問題にスティフな問題用のルーチンを使用しても解を得ることができるが効率が非常に悪くなる.

XLPackに収録されているルーチン (スティフではない問題用)

XLPackに収録されているルーチン (スティフな問題用)

ここでは関数の呼び出し回数でしか評価していないが, 総じて, 高精度公式を使っているDop853とアダムス・バシュフォース・モルトン法のDeabmが速い. また, 補外法のOdexも悪くないようである. ステップ幅の制御機能を持たない通常のルンゲ・クッタ法も比較してみたが関数の変化が激しいような実践問題では苦しいようである.

スティフではない問題

剛体の回転に対するオイラー方程式

剛体の回転に対するオイラー方程式は次のとおり.

  
  I1y1' = (I2 - I3)y2y3
  I2y2' = (I3 - I1)y3y1
  I3y3' = (I1 - I2)y1y2 + f(x)

 
y₁, y₂, y₃ は回転ベクトルの座標, I₁, I₂, I₃ は主慣性モーメントである. 第3座標には次の外部強制力が加わっている.

  
  f(x) = 0.25 sin2x  (3π ≦ x ≦ 4π)
         0  (その他のx)

 
初期値は I₁ = 0.5, I₂ = 2, I₃ = 3, y₁(0) = 1, y₂(0) = 0, y₃(0) = 0.9 である.

これを, Derkf, Dopri5, Dverk, Dop853, Deabm および Odex を用いて Tol(要求精度) = 10^-6, 10^-7, …, 10^-15 と変化させて計算したときに得られた t = 20 における y₁, y₂, y₃ の相対誤差の平均に対する関数評価回数を両対数プロットしたものである. 誤差は右にいくほと少ない(精度がよい). 比較のためオイラー法(Euler)と4次のルンゲ・クッタ法(RK4)も表示した. こちらは h = 0.01, 0.005, …, 10^-5 と変化させて計算した.

他の問題でもそうであるが, Tol = 10^-13, 10^-14, 10^-15 などは過剰な要求だと思われ, ある程度のところからは結果の精度は上がらない. ルーチンによって途中で計算を終了するものと構わず計算を続けるものがあるので注意が必要だろう.

この問題は4次のルンゲ・クッタ法は少し変な動きを示しているが精度10^-11付近まで問題なく計算できている. オイラー法でもステップ幅を小さくしていけば時間はかかるが正しく計算できそうである.

ローレンツ方程式

ローレンツ方程式はカオス的ふるまいをする方程式である.

  
  y1' = -σy1 + σy2
  y2' = -y1y3 + ry1 - y2
  y3' = y1y2 - by3

 
ザルツマンの値 σ = 10, r = 28, b = 8/3 を使うと非周期的になる. 初期値は y₁(0) = -8, y₂(0) = 8, y₃(0) = 27 である. 次のようなローレンツ・アトラクタとよばれる図形になる.

上と同様に, Derkf, Dopri5, Dverk, Dop853, Deabm および Odex の結果を示す. 比較のためオイラー法(Euler)と4次のルンゲ・クッタ法(RK4)も表示した. t = 16 における相対誤差の平均をプロットした.

この問題は精度が出にくく, 10^-7～10^-8付近で限界のようである.

4次のルンゲ・クッタ法では問題なく計算できたが, オイラー法ではステップ幅を小さくしても意味のある値は得られなかった.

アレンストーフ軌道

制限三体問題の例である. 質量がμとμ'(= 1 – μ)の2つの天体が平面内回転運動をしており, その2つに比べて質量が無視できる第3の物体が同一平面内を運動する.

  
  y1'' = y1 + 2y2' - μ'(y1 + μ)/D1 - μ(y1 - μ')/D2
  y2'' = y2 - 2y1' - μ'y2/D1 - μy2/D2
  D1 = ((y1 + μ)2 + y22)3/2
  D2 = ((y1 - μ')2 + y22)3/2
  μ = 0.012277471, μ' = 1 - μ

 
次の初期値を用いると, t = 17.0652165601579625588917206249 で周期解となることが知られている.

  
  y1(0) = 0.994, y1'(0) = 0
  y2(0) = 0, y2'(0) = -2.00158510637908252240537862224

 
これはアレンストーフ軌道とよばれている. 1周期分を以下に示す.

2階常微分方程式なので変数変換して2倍の元数の1階の連立方程式にして解く. 上と同様に, Derkf, Dopri5, Dverk, Dop853, Deabm および Odex の結果を示す. 比較のためオイラー法(Euler)と4次のルンゲ・クッタ法(RK4)も表示した. 1周期後の相対誤差の平均をプロットした.

4次のルンゲ・クッタ法の結果を見ると本来Derkf(RKF45)より少し悪い程度のはずだがはるかに多くの(100倍程度の)関数計算を必要としている. これは, 関数の値の変化が急激なところと緩やかなところの差が大きいためと考えられる. このような問題はステップ幅の自動調節機能がないと苦しいようである.

アレンストーフ軌道 (2)

アレンストーフ軌道は固定座標系では次のようになる.

  
  y1'' = μ'(a1(x) - y1)/D1 + μ(b1(x) - y1)/D2
  y2'' = μ'(a2(x) - y2)/D1 + μ(b2(x) - y2)/D2
  D1 = ((y1 - a1(x))2 + (y2 - a2(x))2)3/2
  D2 = ((y1 - b1(x))2 + (y2 - b2(x))2)3/2
  a1(x) = -μcos(x), a2(x) = -μsin(x)
  b1(x) = μ'cos(x), b2(x) = μ'sin(x)
  μ = 0.012277471, μ' = 1 - μ

 
初期値は次のようになる.

  
  y1(0) = 0.994, y1'(0) = 0
  y2(0) = 0, y2'(0) = -2.00158510637908252240537862224 + 0.994

 
この軌道は重なってわかりにくいが, 1周期分を示すと次のようになる.

この方程式は y” = f(t, y) の形の2階常微分方程式なので, この形専用の補外法ルーチンOdex2が使える. 比較対象としては通常のOdexの他にDoprinとRknint(これらはXLPackには含まれていない)を使う. Doprinは7(6)次埋め込み型ニュストレム法[3], Rknintは6(4)次 または 12(10)次埋め込み型ニュストレム法[4]のステップ幅の自動調節機能付きルーチンである.

比較のために4次のルンゲ・クッタ法(RK4), 4次のニュストレム法(NY4)および5次のニュストレム法(NY5)の結果も示した.

Odex2は予定通りOdexの約半分の関数計算で解を求めている. 埋め込み型ニュストレム法と比較しても遜色ない.

上の回転座標系の場合と同じで, ステップ幅の自動調節機能がないRK4, NY4, NY5は多くの関数計算を必要とした.

文献[3] E. Hairer, S.P. Norsett and G. Wanner, “Solving Ordinary Differential Equations I”, Springer Series in Computational Mathematics, Springer-Verlag (1987) [文献[1]の旧版]

文献[4] R.W. Brankin, I. Gladwell, J.R. Dormand, P.J. Prince & W.L. Seward (1989): Algorithm 670. A Runge-Kutta-Nystrom code. ACM Trans. Math. Softw., Vol.15, p.31-40.

天文力学の問題

座標と質量がそれぞれ (x_i, y_i), m_i = i (i = 1, …, 7) の7個の星の平面内の運動である.

  
  xi'' = Σ mj(xj) - xi)/rij (Σは j ≠ i)
  yi'' = Σ mj(yj) - yi)/rij (Σは j ≠ i)
  rij = ((xj) - xi)2 + (yj) - yi)2)3/2 (i, j = 1, ..., 7)

 
初期値は次のとおりである.

  
  x1(0) = 3, x2(0) = 3, x3(0) = -1, x4(0) = -3, 
  x5(0) = 2, x6(0) = -2, x7(0) = 2,
  y1(0) = 3, y2(0) = -3, y3(0) = 2, y4(0) = 0, 
  y5(0) = 0, y6(0) = -4, y7(0) = 4,
  x6'(0) = 1.75, x7'(0) = -1.5, y4'(0) = -1.25, y5'(0) = 1, 
  xi'(0) = yi'(0) = 0 (その他のi)

 
t = 0, 0.1, 0.2, …, 3 をプロットすると次のようになる.

2階常微分方程式なので変数変換して1階の2倍の元数の14元連立方程式にして解く. 上と同様に, Derkf, Dopri5, Dverk, Dop853, Deabm および Odex の結果を示す. 比較のためオイラー法(Euler)と4次のルンゲ・クッタ法(RK4)も表示した. t = 3 における7つの座標の相対誤差の平均をプロットした.

この方程式は y” = f(t, y) の形の2階常微分方程式なので, この形専用の補外法ルーチンOdex2が使える. 比較対象としては通常のOdexの他にDoprinとRknintを表示した. Doprinは7(6)次埋め込み型ニュストレム法, Rknintは6(4)次 または 12(10)次埋め込み型ニュストレム法のステップ幅の自動調節機能付きルーチンである(アレンストーフ軌道 (2) を参照).

比較のために4次のルンゲ・クッタ法(RK4), 4次のニュストレム法(NY4)および5次のニュストレム法(NY5)の結果も示した.

この問題でもステップ幅の自動調節機能がないRK4, NY4, NY5は多くの関数計算を必要としており, 自動調節機能の効果がわかる.

感染モデル(遅延微分方程式)

遅延微分方程式(DDE)の例として伝染性疾患の古典的モデルを示す.

人口の中で y₁(t)は感染可能な人(susceptible), y₂(t)は感染者(infected), y₃(t)は治癒者など免疫があり感染不可能な人(removed)の数を示す値とする. 単位時間あたりに感染する人の数は積y₁(t)y₂(t)に比例し, 新しく治癒した人の数は感染者の数に比例するとすれば次のモデルが得られる.

  
  y1'(t) = -y1(t)y2(t)
  y2'(t) = y1(t)y2(t) -  y2(t)
  y3'(t) = y2(t)

 
初期値は次のとおりとする.

  
  y1(0) = 5, y2(0) = 0.1, y3(0) = 1

 
このモデルでは次のように全員が感染してそして治癒して免疫を獲得するのでそれ以上何も起こらなくなる.

もし免疫を持つ人が一定期間τ後に免疫の効果がなくなり再び感染可能になるとすると周期的に感染症が発生すると予想される. また, 潜伏期間τ₂も導入することにすれば次のモデルが得られ, 過去の値が必要な遅延微分方程式となり通常のプログラムでは解くことができない.

  
  y1'(t) = -y1(t)y2(t - τ2) + y2(t - τ)
  y2'(t) = y1(t)y2(t - τ2) -  y2(t)
  y3'(t) = y2(t) - y2(t - τ)

 
τ = 10, τ₂ = 1 としてRetardを使って計算すると次のように周期的発生の図が得られた. なお, 初期値は上と同じであるが, 遅延微分方程式ではt = -τ ～ 0 の値も必要になるのでここでは初期値に同じとした.

ここでは関数の呼び出し回数でしか評価していないが, Debdf (BDF), Radau5 (5次IRK), Seulex (補外法) は同程度であった. Radau(IRK 5次～13次可変)は高次の公式が選択されるのか他より速いようである. Rodas (ローゼンブロック法)は遅いようにみえるが計算量は多くない可能性がある(今回は計算量による評価は行っていない).

スティフな問題

ファン・デル・ポル振動子

常微分方程式

  
  y'' + αy' + y = 0

 
の解は, α > 0 ならば減衰し, α < 0 ならば不安定になる.

yが小さいときは α > 0 に, yが大きいときは α < 0 となるようにαを変化させるために α = ε(y² – 1) (ε > 0) とし, さらに y₁‘ = y₂ とおくと次の連立常微分方程式が得られる.

  
  y1' = y2
  y2' = ε(1 - y12)y2 - y1  (ε > 0)

 
この方程式では, 小さい振動は増幅され,　大きい振動は減衰される. 従って, すべての解がそこへ収束するような安定な周期解(極限周期軌道)の存在が予想される. ε = 1, 初期値 y₁(0) = 2, y₂(0) = -2 を始点として計算すると次のように周期軌道に入った.

各ルーチンの比較のために, ε = 10^-6, 初期値 y₁(0) = 2, y₂(0) = 0 として計算する. Debdf, Radau5, Radau, Rodas および Seulex を用いて Tol(要求精度) = 10^-6, 10^-7, …, 10^-15 と変化させて計算したときに得られた t = 11 における y₁ および y₂ の相対誤差の平均に対する関数評価回数を両対数プロットする. 誤差は右にいくほと少ない(精度がよい).

オレゴネータ

3次元の極限周期軌道を持つ化学反応モデルで, 次の式で表される.

  
  y1' = 77.27(y2 + y1(1 - 8.375×10-6y1 - y2))
  y2' = (1/77.27)(y3 - (1 + y1)y2)
  y3' = 0.161(y1 - y3)

 
初期値は y₁(0) = 1, y₂(0) = 2, y₃(0) = 3 である. t = 360 のときの値を求める. 次のように解の絶対値のオーダーが急激に変化するスティフな方程式である.

上と同様に, Debdf, Radau5, Radau, Rodas および Seulex の結果を示す.

ブラセレータ

ブラセレータは次の反応式で表される仮想的な化学反応系である.

  
            k1
       A  ----->  X
            k2
   B + X  ----->  Y + D (2分子反応)
            k3
  2X + Y  ----->  3X    (自己触媒3分子反応)
            k4
       X  ----->  E

 
AとBは一定値であるとし, また, 反応速度k_iをすべて1にする. u(t) = X(t), v(t) = Y(t) とおけば次式が得られる.

  
  u' = A + u2v - (B + 1)u
  v' = Bu - u2v

 
拡散を伴う場合には拡散項を追加した偏微分方程式になる. 1次元の場合は次のようになる.

  
  ∂u/∂t = A + u2v - (B + 1)u + α∂2u/∂x2
  ∂v/∂t = Bu - u2v + α∂2v/∂x2

 
ここで, 0 ≦ x ≦ 1, A = 1, B = 3, α = 1/50, 境界条件, 初期条件は次のとおりとする.

  
  u(0, t) = u(1, t) = 1, v(0, t) = v(1, t) = 3
  u(x, 0) = 1 + (1/2)sin(2πx), v(x, 0) = 3

 
xに関する2階微分を x₁ = i/(N + 1) (1 ≦ i ≦ N), Δx = 1/(N + 1) のN点格子上の差分で置き換えると次の常微分方程式系が得られる.

  
  u' = 1 + ui2vi - 4ui + α/(Δx)2(ui-1 - 2ui + ui+1)
  v' = 3ui - ui2vi + α/(Δx)2(vi-1 - 2vi + vi+1)

  u0(t) = uN+1(t) = 1, v0(t) = vN+1(t) = 3
  ui(0) = 1 + (1/2)sin(2πxi), vi(0) = 3  (i = 1, ..., N)

 
N = 500 とすると, 最大固有値約-20000を持つスティフな1000元の微分方程式系になる. 端付近(xが0または1に近いところ)と中央付近(xが0.5付近)の値をプロットすると次のようになった.

上と同様に, Debdf, Radau5, Radau, Rodas および Seulex の結果を示す. 誤差は t = 10 における上図の6点の値の相対誤差の平均値とした.なお, この問題の場合, 連立方程式の元数が大きくヤコビ行列が上下2の帯幅の帯行列であることがわかっているので, Mljac = Mujac = 2 (Debdfの場合は Ml = Mu = 2) と設定しておく必要がある. そうしないと計算量が膨大になってしまう.

微分代数方程式

単振り子

微分代数方程式の項で例として使用した単振り子を直交座標で再び取り上げる.

下図のように, 長さl[m]の糸に重さm[kg]のおもりがついた振り子を考える. おもりの位置は(q₁, q₂)とする. g[m/s²]は重力加速度である.

この振り子の運動は次の常微分代数方程式で表された. v₁, v₂はq₁, q₂方向の速度である.

  
  q1' = v1
  q2' = v2
  v1' = -λq1
  v2' = -λq2 - g
  0 = q12 + q22 - l2  [指数3]

 
この式は指数3である. 指数2および指数1のとき最後の式は次のようであった.

  
  0 = q1v1 + q2v2  [指数2]
  0 = -λ - gq2 + v12 + v22  [指数1]

 
指数1の式のλを最初の4つの式に代入して得られる常微分方程式(= 指数0)は次のとおりである.

  
  q1' = v1
  q2' = v2
  v1' = (gq2 - v12 - v22)q1
  v2' = (gq2 - v12 - v22)q2 - g

 
以上の式を使って t = 0 (q₁ = sin(π/4), q₂ = -cos(π/4), v₁ = 0, v₂ = 0) から始めて2周期後の位置を求めその相対誤差と関数計算回数(上の5本(または4本)の式をセットで計算して1回とカウントする)をプロットする. なお, Tol = 1.0e-4, 1.0e-5, …, 1.0e-12 と変えて計算した.

Radau5, Radau, Rodas, SeulexおよびDasslを使用した(Radaupは対象としなかったがRadau5とRadauの結果に準じると考えてよい). 比較のために, 変換した常微分方程式をDerkfで解いた場合および同じく安定化して解いた場合も載せた.

Radau5, Radau, Rodas, Seulexは Mu’ = φ(u) 型の問題を解くように作られていて, この場合Mを次のように設定した.

  
      ( 1             )
      (    1          )
  M = (       1       )
      (          1    )
      (             0 )

 
注 – これらのルーチンは M = I(単位行列) のときは常微分方程式を解く陰解法として動作し, M <> I であれば微分代数方程式を解く.

Radau5とRadauは指数2と3に適用するときはパラメータNind1, Nind2, Nind3を設定する必要があり, それぞれ 4, 1, 0 および 4, 0, 1 とした.

指数1の場合

Radau5とRadauは Tol = 1.0e-6以下ではアンダーフローで停止した. Dasslは Tol = 1.0e-7以下では収束せず停止した. Derkf(常微分方程式)の場合, 安定化を行うと高精度の領域ではかえって計算量が増えたが, 精度が低いときには効果があった.

指数2の場合

Seulexは指数2は適用範囲外であるが解を求めることができたものの計算量が非常に多くなった. RodasとDasslは低精度の領域では解を求めることができたが計算量は非常に多かった.

指数3の場合

指数3で使えるのはRadau5とRadauであるが常微分方程式に比べると計算量は非常に多い.

トランジスター回路

次のトランジスター増幅回路を考える. 図の中で 1, 2, …, 5 の点における電圧を U₁, U₂, …, U₅とする. U_e(t)は入力信号電圧, U_bは電源電圧である(U_b = 6V). 出力信号は5における電圧U₅として取り出すことができる.

f(U)とU_e(t)は次のとおりである.

  
  f(U) = 10-6(exp(U/0.026) - 1)
  Ue(t) = 0.4・sin(200πt)

 
抵抗Rを流れる電流は I = U/R, コンデンサCを流れる電流は I = C・dU/dt を満たす. トランジスタでは2から3に流れる電流が99倍増幅されて4から3に流れる電流になるものとし, これらの電流は電位差 U2 – U3 に依存し, 関数f(U2 – U3)で表されものとする. キルヒホッフの法則によりある点に出入りする電流の和はゼロであるから, 1, 2, …, 5 の点において方程式を立てると次のように指数1の微分代数方程式が得られる.

  
  Ue(t)/R0 - U1/R0 + C1(U2' - U1') = 0
  Ub/R2 - U2(1/R1 + 1/R2) + C1(U1' - U2') -0.01f(U2 - U3) = 0
  f(U2 - U3) - U3/R3 - C2U3' = 0
  Ub/R4 - U4/R4 + C3(U5' - U4') -0.99f(U2 - U3) = 0
  -U5/R5 + C3(U4' - U5') = 0

 
これを Mu’ = φ(u) 型の問題として表すと次のようになる.

  
      ( -C1  C1             )
      (  C1 -C1             )
  M = (         -C2         )
      (             -C3  C3 )
      (              C3 -C3 )

       ( U1' )
       ( U2' )
  u' = ( U3' )
       ( U4' )
       ( U5' )

          ( -Ue(t)/R0 + U1/R0                         )
          ( -Ub/R2 + U2(1/R1 + 1/R2) + 0.01f(U2 - U3) )
  φ(u) = ( -f(U2 + U3) - U3/R3                       )
          ( -Ub/R4 + U4/R4 + 0.99f(U2 - U3)           )
          ( U5/R5                                     )

 
Radau5を使って t = 0 ～ 0.1 を計算すると, 入力信号と出力信号は次のようになった.

なお, y₁ = U₁ – U₂, y₂ = U₃, y₃ = U₄ – U₅, z₁ = U₂, z₂ = U₅ とおくと, 次のように微分代数方程式の一般形 (y’ = f(t, y, z), 0 = g(t, y, z)) に変形できる.

  
  y1' = (Ue(t) - y1 + z1)/(R0C1)
  y2' = (f(z1 - y2) -y2/R3)/C2
  y3' = z2/(R5C3)
  0 = (Ue(t) - y1 - z1)/R0 + Ub/R2 - z1(1/R1 + 1/R2) - 0.01f(z1 - y2)
  0 = (Ub - y3 - z2)/R4 - z2/R5 - 0.99f(z1 - y2)

 
Radau5, Radau, Rodas, SeulexおよびDasslを比較してみる. t = 0 から始めて t = 0.1 におけるU₅の値の相対誤差と関数計算回数をプロットする. なお, Tol = 1.0e-4, 1.0e-5, …, 1.0e-10 と変えて計算した.

Mu’ = φ(u) 型の式を使った時の結果は次のとおりである.

一般形 (y’ = f(t, y, z), 0 = g(t, y, z)) に変形した式を使うと次のようになった.

似たような結果であるが, Dasslは不規則な動きを示した.