10.4 常微分方程式の実用プログラムのベンチマーク

10.4.1 実用プログラムの選択

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

XLPack に収録されているプログラム (スティフではない問題用)

※1 それぞれのリバースコミュニケーション版も使用可能(ルーチン名に”_r”のサフィックスが付く).
※2 Deabm では指定しない限り補間値が使用できればそれを使用する.

XLPack に収録されているプログラム (スティフではない 2 階常微分方程式 (y” = f(t, y)) 用)

※1 それぞれのリバースコミュニケーション版も使用可能(ルーチン名に”_r”のサフィックスが付く)

XLPack に収録されているプログラム (スティフな問題用)

※1 それぞれのリバースコミュニケーション版も使用可能(ルーチン名に”_r”のサフィックスが付く)
※2 Debdf では指定しない限り補間値が使用できればそれを使用する.

10.4.2 スティフでない問題

本節の例題の詳細は参考文献[5]を参照のこと.

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

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

剛体の回転に対するオイラー方程式は次のとおり.
\[
\begin{align}
& I_1y_1′ = (I_2 – I_3)y_2y_3 \\
& I_2y_2′ = (I_3 – I_1)y_3y_1 \\
& I_3y_3′ = (I_1 – I_2)y_1y_2 + f(x) \\
\end{align}
\] \(y_1, y_2, y_3\) は回転ベクトルの座標, \(I_1, I_2, I_3\) は主慣性モーメントである. 第 3 座標には次の外部強制力が加わっている.
\[
\begin{equation}
f(x) =
\begin{cases}
0.25 sin^2x & (3\pi \le x \le 4\pi) \\
0 & (その他の x) \\
\end{cases}
\end{equation}
\] 初期値は \(I_1 = 0.5, \space I_2 = 2, \space I_3 = 3, \space y_1(0) = 1, \space y_2(0) = 0, \space y_3(0) = 0.9\) である.

これを, Derkfa, Dopri5a, Dverka, Dop853a, Deabm および Odexa を用いて \(Tol \space (要求精度) = 10^{-6}, 10^{-7}, \dots, 10^{-15}\) と変化させて計算した. 比較のためオイラー法 (Euler) と 4 次ルンゲ・クッタ法 (RK4) も表示した. こちらは \(h = 0.01, \space 0.005, \space \dots, \space 10^{-5}\) と変化させて計算した.

x 軸は t = 20 における \(y_1, \space y_2, \space y_3\) の相対誤差の平均で, 右にいくほど小さい (精度がよい).　y 軸はそのときの関数評価回数である. 両軸とも対数目盛でプロットした.

他の問題でもそうであるが, \(Tol = 10^{-13}, 10^{-14}, 10^{-15}\) などは過剰な要求だと思われ, ある程度のところからは結果の精度は上がらない. プログラムによって途中で計算を終了するものと構わず計算を続けるものがある.

この問題は 4 次ルンゲ・クッタ法は少し変な動きを示しているが精度 \(10^{-11}\) 付近までは問題なく計算できている. オイラー法でもステップ幅を小さくしていけば正しく計算できるようであるが, 時間がかかりすぎる (関数評価回数が多すぎる).

(2) ローレンツ方程式

ローレンツ方程式はカオス的ふるまいをする方程式である.
\[
\begin{align}
& y_1′ = -\sigma y_1 + \sigma y_2 \\
& y_2′ = -y_1y_3 + ry_1 – y_2 \\
& y_3′ = y_1y_2 – by_3 \\
\end{align}
\] ザルツマンの値 \(\sigma = 10, \space r = 28, \space b = 8/3\) を使うと非周期的になる. 初期値を \(y_1(0) = -8, \space y_2(0) = 8, \space y_3(0) = 27\) として計算すると次のようなローレンツ・アトラクタとよばれる図形が得られる.

上と同様に, Derkfa, Dopri5a, Dverka, Dop853a, Deabm および Odexa を用いて \(Tol \space (要求精度) = 10^{-6}, 10^{-7}, \dots, 10^{-15}\) と変化させて計算した. 比較のためオイラー法 (Euler) と 4 次ルンゲ・クッタ法 (RK4) も表示した. こちらは \(h = 0.01, \space 0.005, \space \dots, \space 10^{-5}\) と変化させて計算した.

x 軸は t = 16 における \(y_1, \space y_2, \space y_3\) の相対誤差の平均で, 右にいくほど小さい (精度がよい).　y 軸はそのときの関数評価回数である. 両軸とも対数目盛でプロットした.

この問題は精度が出にくく, t が大きい領域 (例えば, t = 50) では高精度のプログラムでも計算が難しいとされる (文献[5]). 本例 (t = 16) でも最大でも \(10^{-7} \sim 10^{-8}\) 程度の精度しか得られない. オイラー法ではステップ幅を小さくしても意味のある値は得られなかった.

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

制限三体問題の例である. 質量が \(\mu\) と \(\mu’ (= 1 – \mu)\) の 2 つの天体 (月と地球を想定) が平面内回転運動をしており, その 2 つに比べて質量が無視できる第 3 の物体が同一平面内を運動するものとする. 月と地球の座標がそれぞれ (1, 0) と (0, 0) に固定されているような座標系で表すものとする. 第 3 の物体の x 座標を \(y_1\), y 座標を \(y_2\) で表すと方程式は次のようになる.
\[
\begin{align}
& y_1^{\prime\prime} = y_1 + 2y_2′ – \mu'(y_1 + \mu)/D_1 – \mu(y_1 – \mu’)/D_2 \\
& y_2^{\prime\prime} = y_2 – 2y_1′ – \mu’y_2/D_1 – \mu y_2/D_2 \\
& D_1 = ((y_1 + \mu)^2 + y_2^2)^{3/2} \\
& D_2 = ((y_1 – \mu’)^2 + y_2^2)^{3/2} \\
& \mu = 0.012277471, \mu’ = 1 – \mu \\
\end{align}
\] 次の初期値を用いると, t = 17.0652165601579625588917206249 で周期解になることが知られている.
\[
\begin{align}
& y_1(0) = 0.994, \space y_1′(0) = 0 \\
& y_2(0) = 0, \space y_2′(0) = -2.00158510637908252240537862224 \\
\end{align}
\] これはアレンストーフ軌道とよばれており, 次のような軌道になる.

緑の丸は地球, 赤の丸は月, 青の丸は第 3 の物体を表す.

この方程式は 2 階常微分方程式なので 1 階の 4 元連立常微分方程式に変数変換して解く. 上と同様に, Derkfa, Dopri5a, Dverka, Dop853a, Deabm および Odexa を用いて \(Tol \space (要求精度) = 10^{-6}, 10^{-7}, \dots, 10^{-15}\) と変化させて計算した. 比較のためオイラー法 (Euler) と 4 次ルンゲ・クッタ法 (RK4) も表示した. こちらは \(h = 0.01, \space 0.005, \space \dots, \space 10^{-5}\) と変化させて計算した.

x 軸は 1 周期後 (t = 17.065216560158) における \(y_1, \space y_2\) の相対誤差の平均で, 右にいくほど小さい (精度がよい).　y 軸はそのときの関数評価回数である. 両軸とも対数目盛でプロットした.

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

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

同じ制限三体問題の方程式を絶対座標系で表すと次のようになる.
\[
\begin{align}
& y_1^{\prime\prime} = \mu'(a_1(x) – y_1)/D_1 + \mu(b_1(x) – y_1)/D_2 \\
& y_2^{\prime\prime} = \mu'(a_2(x) – y_2)/D_1 + \mu(b_2(x) – y_2)/D_2 \\
& D_1 = ((y_1 – a_1(x))^2 + (y_2 – a_2(x))^2)^{3/2} \\
& D_2 = ((y_1 – b_1(x))^2 + (y_2 – b_2(x))^2)^{3/2} \\
& a_1(x) = -\mu cos(x), \space a_2(x) = -\mu sin(x) \\
& b_1(x) = \mu’cos(x), \space b_2(x) = \mu’sin(x) \\
& \mu = 0.012277471, \space \mu’ = 1 – \mu \\
\end{align}
\] 初期値は次のようになる.
\[
\begin{align}
& y_1(0) = 0.994, \space y_1′(0) = 0 \\
& y_2(0) = 0, \space y_2′(0) = -2.00158510637908252240537862224 + 0.994 \\
\end{align}
\] この座標系では 3 つ全部の天体が動き, 軌道は次のようになる.

緑の丸は地球, 赤の丸は月, 青の丸は第 3 の物体を表す. わずかではあるが地球も動いている.

この方程式は \(y{\prime\prime} = f(t, y)\) の形の 2 階常微分方程式なので, 1 階常微分方程式に変換せずにこのままの形で解くことができるこの形専用のプログラムを使うことができる.

Dopn43, Dopn64, Dopn86, Dopn1210 および Odex2a を用いて \(Tol \space (要求精度) = 10^{-6}, 10^{-7}, \dots, 10^{-15}\) と変化させて計算した. 比較のために 4 次のニュストレム法 (NY4) および 5 次のニュストレム法 (NY5) も表示した. こちらは \(h = 0.01, \space 0.005, \space \dots, \space 10^{-5}\) と変化させて計算した.

(3) の Derkfa と RK4 の関係と同様に, Dopn43 は NY4 および NY5 と同程度の性能が予想されたが実際にはかなり少ない関数評価回数で解を求めている. やはり, ステップ幅の自動調節機能の効果であると思われる.

(5) 天文力学の問題

座標と質量がそれぞれ \((x_i, y_i), \space m_i = i \space (i = 1, \dots, 7)\) の 7 個の星の平面内の運動である.
\[
\begin{align}
& x_i^{\prime\prime} = \sum_{j \ne i} m_j(x_j) – x_i)/r_{ij} \\
& y_i^{\prime\prime} = \sum_{j \ne i} m_j(y_j) – y_i)/r_{ij} \\
& r_{ij} = ((x_j) – x_i)^2 + (y_j) – y_i)^2)^{3/2} \space (i, j = 1, \dots, 7) \\
\end{align}
\] 初期値は次のとおりである.
\[
\begin{align}
& x_1(0) = 3, x_2(0) = 3, x_3(0) = -1, x_4(0) = -3, x_5(0) = 2, x_6(0) = -2, x_7(0) = 2, \\
& y_1(0) = 3, y_2(0) = -3, y_3(0) = 2, y_4(0) = 0, y_5(0) = 0, y_6(0) = -4, y_7(0) = 4, \\
& x_6′(0) = 1.75, x_7′(0) = -1.5, x_i'(0) = = 0 \space (その他の i) \\
& y_4′(0) = -1.25, y_5′(0) = 1, y_i'(0) = 0 \space (その他の i) \\
\end{align}
\] t = 0 ～ 3 を表示すると次のようになる.

この区間では 6 回の衝突が起こっており, そのときに軌道と加速度が急激に変化する.

この方程式は 7 元の 2 階常微分方程式なので 1 階の 14 元連立常微分方程式に変数変換して解く. 上と同様に, Derkfa, Dopri5a, Dverka, Dop853a, Deabm および Odexa を用いて \(Tol \space (要求精度) = 10^{-6}, 10^{-7}, \dots, 10^{-15}\) と変化させて計算した. 比較のためオイラー法 (Euler) と 4 次ルンゲ・クッタ法 (RK4) も表示した. こちらは \(h = 0.01, \space 0.005, \space \dots, \space 10^{-5}\) と変化させて計算した.

x 軸は t = 3 における 7 つの座標の相対誤差の平均で, 右にいくほど小さい (精度がよい).　y 軸はそのときの関数評価回数である. 両軸とも対数目盛でプロットした.

RK4 の結果を見ると (3) と同様にステップ幅の自動調節機能の効果が大きいことがわかる.

よく見るとこの方程式は \(y^{\prime\prime} = f(t, y)\) の形の 2 階常微分方程式なので, この形専用のプログラムを使うことができて, 7 元の 2 階常微分方程式のまま変数変換なしで計算することができる.

Dopn43, Dopn64, Dopn86, Dopn1210 および Odex2a を用いて \(Tol \space (要求精度) = 10^{-6}, 10^{-7}, \dots, 10^{-15}\) と変化させて計算する. 比較のために 4 次のニュストレム法 (NY4) および 5 次のニュストレム法 (NY5) も表示した. こちらは \(h = 0.01, \space 0.005, \space \dots, \space 10^{-5}\) と変化させて計算した.

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

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

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

人口の中で \(y_1(t)\) は感染可能な人 (susceptible), \(y_2(t)\) は感染者 (infected), \(y_3(t)\) は治癒者など免疫があり感染不可能な人 (removed) の数を示す値とする. 単位時間あたりに感染する人の数は積 \(y_1(t)y_2(t)\) に比例し, 新しく治癒した人の数は感染者の数に比例するとすれば次のモデルが得られる.
\[
\begin{align}
& y_1′(t) = -y_1(t)y_2(t) \\
& y_2′(t) = y_1(t)y_2(t) – y_2(t) \\
& y_3′(t) = y_2(t) \\
\end{align}
\] 初期値は次のとおりとする.
\[
y_1(0) = 5, \space y_2(0) = 0.1, \space y_3(0) = 1
\] このモデルでは次のように全員が感染してそして治癒して免疫を獲得するのでそれ以上何も起こらなくなる.

もし免疫を持つ人が一定期間 \(\tau\) 後に免疫の効果がなくなり再び感染可能になるとすると周期的に感染症が発生すると予想される. また, 潜伏期間 \(\tau_2\) も導入することにすれば次のモデルが得られ, 過去の値が必要な遅延微分方程式となり通常のプログラムでは解くことができない.
\[
\begin{align}
& y_1′(t) = -y_1(t)y_2(t – \tau_2) + y_2(t – \tau) \\
& y_2′(t) = y_1(t)y_2(t – \tau_2) – y_2(t) \\
& y_3′(t) = y_2(t) – y_2(t – \tau) \\
\end{align}
\] XLPack の Retarda はドルマン・プリンス法プログラム Dopri5a を遅延微分方程式用に修正したものである. Dopri5a では密出力のために 1 ステップ前までの過去の値を保持するが, Retarda はずっと前までの値を保持するようになっている.

\(\tau = 10, \space \tau_2 = 1\) として Retarda を使って計算すると, 次のように周期的発生の図が得られた. 初期値は上と同じであるが, 遅延微分方程式の定義では \(t = -\tau \sim 0\) の値も必要になる. ここでは初期値と同じ値をとるものとした.

10.4.3 スティフな問題

本節の例題の詳細は参考文献[6]を参照のこと.

Debdf (BDF), Seulexa (補外法) および Radaua (5 次 ～ 13 次可変 IRK) は問題によるが, 同程度の結果であった. Rodasa (ローゼンブロック法) はほかに比べて遅いようにみえるが, ここでは関数評価回数でしか評価していないためであり総計算量は多くない可能性がある.

(1) ファン・デル・ポル振動子

常微分方程式
\[
y^{\prime\prime} + \alpha y’ + y = 0
\] の解は, \(\alpha > 0\) ならば減衰し, \(\alpha < 0\) ならば不安定になる. y が小さいときは \(\alpha > 0\) に, y が大きいときは \(\alpha < 0\) となるように \(\alpha\) を変化させるために \(\alpha = \epsilon(y^2 - 1) \space (\epsilon > 0)\) とし, さらに \(y_1′ = y_2\) とおくと次の連立常微分方程式が得られる.
\[
\begin{align}
& y_1′ = y_2 \\
& y_2′ = \epsilon(1 – y_1^2)y_2 – y_1 \space (\epsilon > 0) \\
\end{align}
\] この方程式では, 小さい振動は増幅され, 大きい振動は減衰される. そのため, すべての解がそこへ収束するような安定な周期解 (極限周期軌道) が存在する. \(\epsilon = 1\), 初期値 \(y_1(0) = 2, y_2(0) = -2\) を始点とした場合, 次のように周期軌道に入る.

Debdf, Radaua (次数自動可変), Rodasa および Seulexa の比較を行う. ここでは, \(\epsilon = 10^{-6}\), 初期値 \(y_1(0) = 2, y_2(0) = 0\) とする. \(Tol (要求精度) = 10^{-6}, 10^{-7}, \dots, 10^{-15}\) と変化させて計算したときに得られた \(t = 11\) における \(y_1\) および \(y_2\) の相対誤差の平均に対する関数評価回数を両対数プロットする. 相対誤差の平均は右にいくほと少ない (精度がよい).

(2) オレゴネータ

3 次元の極限周期軌道を持つ化学反応モデルで, 次の式で表される.
\[
\begin{align}
& y_1′ = 77.27(y_2 + y_1(1 – 8.375×10^{-6}y_1 – y_2)) \\
& y_2′ = (1/77.27)(y_3 – (1 + y_1)y_2) \\
& y_3′ = 0.161(y_1 – y_3) \\
\end{align}
\] 初期値は \(y_1(0) = 1, \space y_2(0) = 2, \space y_3(0) = 3\) である.

次のように解の絶対値のオーダーが急激に変化するスティフな方程式である.

極限周期軌道が存在する.

上と同様に, Debdf, Radaua (次数自動可変), Rodasa および Seulexa の比較を行う. \(Tol (要求精度) = 10^{-6}, 10^{-7}, \dots, 10^{-15}\) と変化させて計算したときに得られた \(t = 360\) における \(y_1, \space y_2\) および \(y_3\) の相対誤差の平均に対する関数評価回数を両対数プロットする. 相対誤差の平均は右にいくほと少ない (精度がよい).

(3) ブラセレータ

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

  
            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)\) とおけば次式が得られる.
\[
\begin{align}
& u’ = A + u^2v – (B + 1)u \\
& v’ = Bu – u^2v \\
\end{align}
\] 拡散を伴う場合には拡散項を追加した偏微分方程式になる. 1 次元の場合は次のようになる.
\[
\begin{align}
& \frac{\partial u}{\partial t} = A + u^2v – (B + 1)u + \alpha \frac{\partial^2u}{\partial x^2} \\
& \frac{\partial v}{\partial t} = Bu – u^2v + \alpha \frac{\partial^2v}{\partial x^2} \\
\end{align}
\] ここで, \(0 \le x \le 1, \space A = 1, \space B = 3, \space \alpha = 1/50\) とし, 境界条件, 初期条件を次のようにする.
\[
\begin{align}
& u(0, t) = u(1, t) = 1, \space v(0, t) = v(1, t) = 3 \\
& u(x, 0) = 1 + \frac{1}{2}sin(2\pi x), \space v(x, 0) = 3 \\
\end{align}
\] x に関する 2 階微分を \(x_1 = i/(N + 1) \space (1 \le i \le N), \space \Delta x = 1/(N + 1)\) の N 点格子上の差分で置き換えると次の常微分方程式系が得られる.
\[
\begin{align}
& u’ = 1 + u_i^2v_i – 4u_i + \frac{\alpha}{(\Delta x)^2}(u_{i-1} – 2u_i + u_{i+1}) \\
& v’ = 3u_i – u_i^2v_i + \frac{\alpha}{(\Delta x)^2}(v_{i-1} – 2v_i + v_{i+1}) \\
& u_0(t) = u_{N+1}(t) = 1, \space v_0(t) = v_{N+1}(t) = 3 \\
& u_i(0) = 1 + \frac{1}{2}sin(2\pi x_i), \space v_i(0) = 3 \space (i = 1, \dots, N) \\
\end{align}
\] N = 500 とすると, 最大固有値約 -20000 を持つスティフな 1000 元の微分方程式系になる. 端付近 (x が 0 または 1 に近いところ) と中央付近 (x が 0.5 付近)の値をプロットすると次のようになった.

上と同様に, Debdf, Radaua (次数自動可変), Rodasa および Seulexa の比較を行う. \(Tol (要求精度) = 10^{-6}, 10^{-7}, \dots, 10^{-15}\) と変化させて計算したときに得られた \(t = 10\) における上図の 6 点の値の相対誤差の平均に対する関数評価回数を両対数プロットする. 相対誤差の平均は右にいくほと少ない (精度がよい).

なお, この問題の場合, 連立方程式の元数が大きくヤコビ行列が上下 2 の帯幅の帯行列であることがわかっているので, Mljac = Mujac = 2 (Debdf の場合は Ml = Mu = 2) と設定しておくことができる. そうしないと計算量が膨大になってしまう.

10.4.4 微分代数方程式

(1) 単振り子

「10.3.3.1 微分代数方程式」において取り上げた単振り子を直交座標系で表した例について再び考える.

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

この振り子の運動は次の常微分代数方程式で表される. \(v_1, v_2\) は \(q_1, q_2\) 方向の速度である.
\[
\begin{align}
& q_1′ = v_1 \\
& q_2′ = v_2 \\
& v_1′ = -\lambda q_1 \\
& v_2′ = -\lambda q_2 – g \\
& 0 = q_1^2 + q_2^2 – l^2 \space [指数 3] \\
\end{align}
\] この式は指数 3 である. 最後の式をを微分して得られる指数 2 および指数 1 の式は次のようであった.
\[
\begin{align}
& 0 = q_1v_1 + q_2v_2 \space [指数 2] \\
& 0 = -\lambda – g q_2 + v_1^2 + v_2^2 \space [指数 1] \\
\end{align}
\] さらに, \(\lambda\) [指数 1] を元の連立方程式に代入すると次の 4 元連立常微分方程式 (指数 0) が得られた.
\[
\begin{align}
& q_1′ = v_1 \\
& q_2′ = v_2 \\
& v_1′ = (gq_2 – v_1^2 – v_2^2)q_1 \\
& v_2′ = (gq_2 – v_1^2 – v_2^2)q_2 – g \\
\end{align}
\] 以上の式を使って \(t = 0 (q_1 = sin(\pi/4), q_2 = -cos(\pi/4), v_1 = 0, v_2 = 0)\) から始めて 2 周期後の位置を求めその相対誤差と関数計算回数 (上の 5 本の式をセットで計算して 1 回とカウントする) をプロットする. なお, \(Tol = 10^{-4}, 10^{-5}, \dots, 10^{-12}\) と変えて計算した.

プログラムは Radaua, Rodasa および Seulexa を使用した. これらは Mu’ = φ(u) 型の問題を解くことができるように作られていて, 質量マトリックス M が単位行列の場合 (M = I) には常微分方程式を陰解法で解き, そうでない場合には微分代数方程式 (DAE) を解く. この場合 M を次のように設定した.
\[
M =
\begin{pmatrix}
1 \\
& 1 \\
& & 1 \\
& & & 1 \\
& & & & 0 \\
\end{pmatrix}
\] Radaua は指数 2 と 3 に適用するときはパラメータ Nind1, Nind2, Nind3 を設定する必要があり, ここではそれぞれ 4, 0, 1 と設定した.

さらに, 比較のために 4 元連立常微分方程式 (指数 0) を Derkfa で解いた場合 (安定化なし) の結果もプロットした.

(i) 指数 1 の場合

Radaua は指数 3 まで, Rodasa および Seulexa は指数 1 に適用できるが, それぞれ解を求めることができた.

(ii) 指数 2 の場合

Seulexa は指数 2 は適用範囲外であるが解を求めることができた. ただし, 計算量は多かった. Rodasa は低精度の領域では解を求めることができたが計算量は非常に多かった.

(iii) 指数 3 の場合

指数 3 で使えるのは Radaua であるが, Derkfa (指数 0) に比べると計算量は多い. また, 次数 (5, 9, 13 次) を自動選択するようにすると不安定になったため, ここでは 5 次に固定した.

やはり指数 3 では計算が難しくなるようである. できるだけ指数が低くなるように定式化した方がよさそうである.

(4) トランジスター回路

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

抵抗 R を流れる電流は \(I = U/R\), コンデンサ C を流れる電流は \(I = C dU/dt\) を満たす. トランジスタでは 2 から 3 に流れる電流が 99 倍増幅されて 4 から 3 に流れる電流になるものとし, これらの電流は電位差 \(U2 – U3\) に依存し, 関数 \(f(U2 – U3)\) で表されものとする. キルヒホッフの法則によりある点に出入りする電流の和はゼロであるから, 点 \(1, 2, \dots, 5\) において方程式を立てると次のように指数 1 の微分代数方程式が得られる.
\[
\begin{align}
& \frac{U_e(t)}{R_0} – \frac{U_1}{R_0} + C_1(U_2′ – U_1′) = 0 \\
& \frac{U_b}{R_2} – U_2(\frac{1}{R_1} + \frac{1}{R_2}) + C_1(U_1′ – U_2′) -0.01f(U_2 – U_3) = 0 \\
& f(U_2 – U_3) – \frac{U_3}{R_3} – C_2U_3′ = 0 \\
& \frac{U_b}{R_4} – \frac{U_4}{R_4} + C_3(U_5′ – U_4′) -0.99f(U_2 – U_3) = 0 \\
& -\frac{U_5}{R_5} + C_3(U_4′ – U_5′) = 0 \\
\end{align}
\] ここで \(f(U)\) と \(U_e(t)\) は次の値を使用するものとする.
\[
\begin{align}
& f(U) = 10^{-6}(e^{\frac{U}{0.026}} – 1) \\
& U_e(t) = 0.4 sin(200 \pi t) \\
\end{align}
\] これを \(Mu’ = \phi(u)\) 型の問題として表すと, M, u’ および \(\phi(u)\) は次のようになる.
\[
M =
\begin{pmatrix}
-C_1 & C_1 \\
C_1 & -C_1 \\
& & -C_2 \\
& & & -C_3 & C_3 \\
& & & C_3 & -C_3 \\
\end{pmatrix}
\] \[
u’ =
\begin{pmatrix}
U_1′ \\
U_2′ \\
U_3′ \\
U_4′ \\
U_5′ \\
\end{pmatrix}
\] \[
\phi(u) =
\begin{pmatrix}
-\frac{U_e(t)}{R_0} + \frac{U_1}{R_0} \\
-\frac{U_b}{R_2} + U_2(\frac{1}{R_1} + \frac{1}{R_2}) + 0.01f(U_2 – U_3) \\
-f(U_2 + U_3) – \frac{U_3}{R_3} \\
-\frac{U_b}{R_4} + \frac{U_4}{R_4} + 0.99f(U_2 – U_3) \\
\frac{U_5}{R_5} \\
\end{pmatrix}
\] この式による \(t = 0 \sim 0.2\) における入力信号 \(U_e(t)\) に対する出力信号 \(U_5(t)\) は次のようになる. 初期値は, \(U_1 = 0, \space U_2 = 3, \space U_3 = 3, \space U_4 = 6, \space U_5 = 0 \) とした.

トランジスタ回路による増幅作用が確認できる.

\(Mu’ = \phi(u)\) 型の DAE に適用できる Radaua, Rodasa および Seulexa を用いて計算し比較を行う. \(Tol (要求精度) = 10^{-6}, 10^{-7}, \dots, 10^{-15}\) と変化させて計算したときに得られた \(t = 10\) における上図の 6 点の値の相対誤差の平均に対する関数評価回数を両対数プロットする. 相対誤差の平均は右にいくほと少ない (精度がよい).

なお, \(y_1 = U_1 – U_2, \space y_2 = U_3, \space y_3 = U_4 – U_5, \space z_1 = U_2, \space z_2 = U_5\) とおくと, 次のように微分代数方程式の一般形 (\(y’ = f(t, y, z), \space 0 = g(t, y, z)\)) に変形することもできる.
\[
\begin{align}
& y_1′ = \frac{U_e(t) – y_1 + z_1}{R_0C_1} \\
& y_2′ = \frac{f(z_1 – y_2) -\frac{y_2}{R_3}}{C_2} \\
& y_3′ = \frac{z_2}{R_5C_3} \\
& 0 = \frac{U_e(t) – y_1 – z_1}{R_0} + \frac{U_b}{R_2} – z_1(\frac{1}{R_1} + \frac{1}{R_2}) – 0.01f(z_1 – y_2) \\
& 0 = \frac{U_b – y_3 – z_2}{R_4} – \frac{z_2}{R_5} – 0.99f(z_1 – y_2) \\
\end{align}
\]