10.3 スティフな方程式, 微分代数方程式

10.3.1 解法の安定性

初期値問題を解いている間に数値解が異常な値になったり発散したりすることがある. その原因は, 微分方程式自体の問題である場合と, 数値解法の問題である場合がある.

10.3.1.1 微分方程式の安定性

「10.2.1.1 1 段法の誤差」の説明図に示されるように, 微分方程式の解曲線が時間が進むにつれて離れていく場合, 初期値のわずかな差が大きく拡大する現象が生じることがある. これは不安定な微分方程式の例である.

逆に, 微分方程式の解曲線が時間が進むにつれて近づいていく場合, むしろ誤差が縮小していくので安定な微分方程式の例である.

関数がどちらの傾向を示すかは微分方程式の \(\partial f/\partial y\) の符号により決まる.

10.3.1.2 数値解法の安定性

安定な微分方程式であっても, 使用する数値解法によって不安定現象を示すことがある. これは数値解法の安定性の問題である.

10.3.1.2.1 ステップ幅の制約

いま, 微分方程式の導関数 \(\partial f/\partial y\) が一定値 \(\lambda\) で, \(\lambda < 0\) とする (これは安定な微分方程式の初期値問題である). また, 各ステップで加えられる誤差はステップによらず一定であるとする. その場合, 数値解法の安定性は次の線形テスト問題を調べることでわかる. \[ dy/dt = \lambda y \space (t > 0, \lambda (定数) < 0), \space y(0) = y_0 \] この問題の解は, \[ y(t) = y_0 e^{\lambda t} \] となるから, \(t \to \infty\) のとき \(y(t) \to 0\) となる.

(1) オイラー法

上のテスト問題をオイラー法で解くと,
\[
y_n = y_0(1 + h\lambda)^n
\] となるから, \(n \to \infty\) のとき \(y_n \to 0\) とならなければならない. そのためには,
\[
|1 + h\lambda| < 1 \] でなければならない. これを安定条件という. \(\lambda\) を複素数として \(z = h\lambda\) とおき, 安定条件 \(|R(z)| < 1\) (ただし, \(R(z) = 1 + z\)) を複素平面上に図示すると次のようになる.

色付き部分 \(\{ z \in C; |R(z)| < 1 \}\) を絶対安定領域という. \(\lambda\) が実数であれば \(-2 < z = h\lambda < 0\) である. 区間 \([-2, 0]\) を絶対安定区間という. このとき, ステップ幅は \(0 < h < -2/\lambda\) としなければならない.

(2) 後退オイラー法

後退オイラー法では \(R(z)\) は次のようになる.
\[
R(z) = \frac{1}{1 – z}
\] 絶対安定領域を図示すると次のようになる.

このように, 左半平面が絶対安定領域に含まれるときに A 安定であるといい, ステップ幅 h に対する制限がなくなる.

(3) ルンゲ・クッタ法

s 段 p 次のルンゲ・クッタ法では \(R(z)\) は次のようになる.
\[
R(z) = 1 + z + \frac{z^2}{2!} + \dots + \frac{z^p}{p!} + \gamma_{p+1}z^{p+1} + \dots + \gamma_sz^s
\] ここで, \(\lambda_{p+1}, \dots\) はルンゲ・クッタ法のパラメータにより決まる.

p = s = 2, 3, 4 について絶対安定領域を図示すると次のようになる.

次数が上がると安定領域も広くなることがわかる. 絶対安定区間はそれぞれ, [-2, 0], [-2.513, 0], [-2.785, 0] である. ルンゲ・クッタ法の陽解法は A 安定ではない.

(4) 多段法 (予測子・修正子法)

多段法も陽解法は A 安定ではない.

多段法の陰解法では A 安定な解法は 2 次以下である. すなわち, A 安定な多段法で最も高次なものは台形則 (アダムス・モルトン法の k = 1) である.

台形則では \(R(z)\) は次のようになる.
\[
R(z) = \frac{1+\frac{z}{2}}{1-\frac{z}{2}}
\] この絶対安定領域は左半平面全部となり, 台形則は A 安定である.

数値実験 (21)

次の初期値問題をオイラー法により解く.
\[
y’ = -10y, \space y(0) = 1
\] この方程式の解は, \(y(t) = e^{-10t}\) である. ステップ幅を 0.09, 0.19, 0.2, 0.21 と変えて計算した結果, 次のようになった.

h = 0.2 を境に, これ以上ステップ幅を大きくすると発散した. 0.2 は上の (1) で説明した絶対安定区間の上限にあたり, これを超えると不安定になることが確かめられた. なお, 0.1 < h < 0.2 では振動しながら収束した.

A 安定な解法である後退オイラー法を使って同じ例題を同じように計算すると次のようになった.

後退オイラー法を使用するとステップ幅に関係なく計算できることが確認された.

注 – 「10.1.5.2 ミルン法」で説明したミルン法の不安定性はここで説明した安定性とは別の問題で, ステップ幅を小さくしても不安定性はおさまらない.

10.3.1.2.2 ステップ幅の制約 (連立方程式の場合)

連立方程式の場合, 線形テスト問題で \(y\) をベクトルに拡張して扱う.
\[
\frac{d\boldsymbol{y}}{dt} = \boldsymbol{J}\boldsymbol{y}, \space \boldsymbol{y}(0) = \boldsymbol{y_0}
\] \(\boldsymbol{J}\) は \(\boldsymbol{y}\) の元数と同じ次元数の正方行列となる. 安定な微分方程式の初期値問題であるための条件は, \(\boldsymbol{J}\) のすべての固有値 \(\lambda_i\) の実部が負であることである.
\[
Re(\lambda_i) < 0 \space (すべての i について) \] また, \(z = h\lambda_i\) として, 複素平面上で \(z\) が絶対安定領域に含まれていればその \(h\) について安定である.

10.3.2 スティフな方程式

スティフな方程式 (硬い方程式ともいう) は化学反応や電気回路などの問題で現れる. 非常に急な変化をする解を持ち, 通常の解法では解きにくいものをいう.

10.3.2.1 スティフな方程式の例

はじめに計算例を示す.

数値実験 (22)

次の 2 つの初期値問題を考える.

例 1

\[
\frac{d\boldsymbol{y}}{dt} =
\begin{pmatrix}
-2 & 1 \\
2 & -3 \\
\end{pmatrix}
\boldsymbol{y} +
\begin{pmatrix}
-cos(t) \\
3 cos(t) – sin(t) \\
\end{pmatrix}
, \space \boldsymbol{y}(0) =
\begin{pmatrix}
1 \\
2 \\
\end{pmatrix}
\]

例 2

\[
\frac{d\boldsymbol{y}}{dt} =
\begin{pmatrix}
-2 & 1 \\
1998 & -1999 \\
\end{pmatrix}
\boldsymbol{y} +
\begin{pmatrix}
-cos(t) \\
1999 cos(t) – sin(t) \\
\end{pmatrix}
, \space \boldsymbol{y}(0) =
\begin{pmatrix}
1 \\
2 \\
\end{pmatrix}
\] これらは解析解を求めることができて, それぞれの一般解は次のように表される.

例 1 の一般解

\[
\boldsymbol{y}(t) = C_1 e^{-t}
\begin{pmatrix}
1 \\
1 \\
\end{pmatrix}
+ C_2 e^{-4t}
\begin{pmatrix}
1 \\
-2 \\
\end{pmatrix}
+
\begin{pmatrix}
0 \\
cos(t) \\
\end{pmatrix}
\]

例 2 の一般解

\[
\boldsymbol{y}(t) = C_1 e^{-t}
\begin{pmatrix}
1 \\
1 \\
\end{pmatrix}
+ C_2 e^{-2000t}
\begin{pmatrix}
1 \\
-1998 \\
\end{pmatrix}
+
\begin{pmatrix}
0 \\
cos(t) \\
\end{pmatrix}
\] 定数 \(C_1\) と \(C_2\) は初期値により決まり, 例 1 と 2 の初期値では両方共 \(C_1 = 1, C_2 = 0\) となるため, 解析解は一致し次のようになる.
\[
\boldsymbol{y}(t) = e^{-t}
\begin{pmatrix}
1 \\
1 \\
\end{pmatrix}
+
\begin{pmatrix}
0 \\
cos(t) \\
\end{pmatrix}
\] この関数をグラフに描くと次のようになる.

これら (例 1, 例 2) をルンゲ・クッタ・フェールベルグ法 (XLPack に収録されているプログラム Derkfa) を使って目標精度 \(10^{-7}\) で解を計算し, 得られた解の相対精度 (左目盛) と関数評価回数 (右目盛) をプロットすると次のようになった.

どちらの例題もほぼ目標精度を満たして解が得られた. 図の Err1 (y1, y2) は例 1, Err2 (y1, y2) は例 2 の相対誤差を左目盛で表す. しかし, 右目盛で表される関数評価回数 NFun1 (例 1) と NFun2 (例 2) では NFun2 が異常に多いことがわかる (t = 5 において NFun1 = 695 に対して NFun2 = 20183).

例 1 と例 2 の方程式は同じ初期値と同じ解を持つが, y の係数行列の固有値は例 1 では -1 および -4 であるのに対して, 例 2 では -1 および -2000 になる. そのため 2 つの例では挙動が異なる.

「10.3.1 解法の安定性」の議論から, これらの微分方程式を解くために必要なステップ幅 h を求めると, 例 1 ではオイラー法ならば h < 0.5, 4 次のルンゲ・クッタ法ならば h < 0.696 とすればよい. しかし, 例 2 ではオイラー法ならば h < 0.001, 4 次のルンゲ・クッタ法ならば h < 0.00139 としなければならない.
実際に試してみると, h をこれらの制約条件よりも大きくすると発散して解くことができないことが確かめられた.

上の Derkfa の計算例でも同様のことが起きていてステップ幅が非常に小さく調整されるため, 異常に多くの関数評価回数を必要としたものと思われる.

ステップ幅に制限のない後退オイラー法で同じ問題を計算すると次のようになった.

ステップ幅は 0.01 とした. 後退オイラー法は 1 次の公式であるためこのステップ幅では精度が悪いが, 大きいステップ幅でも計算できることが確かめられた. 例 1 と例 2 を比べても, 例 2 が特に解きにくいということはない.

上の例のように安定性の問題により通常の解法ではステップ幅に制約があるために, 計算に非常に時間がかかったり, 条件によっては振動や発散をしたりするものをスティフな方程式という. 上の例もそうであるが, 一見急な変化をするようには見えないものであっても, 一般解に指数項を含んでいたりするとスティフな方程式となることがある.

10.3.2.2 スティフな方程式の定義

次のような連立常微分方程式を考える.
\[
\frac{d\boldsymbol{y}}{dt} = \boldsymbol{A}\boldsymbol{y} + \phi(t), \space \boldsymbol{y}(0) = \boldsymbol{y_0}
\] 行列 \(\boldsymbol{A}\) の固有値 \(\lambda_i\) はすべて異なりその実数部は負であるとする.

硬度比 (stiffness ratio) を次のように, 実数部が最大の固有値 \(\lambda_{max}\)と実数部が最小の固有値 \(\lambda_{min}\) の実数部の値の比と定義する.
\[
\frac{|Re(\lambda_{max})|}{| Re(\lambda_{min})|}
\] この硬度比が大きい (\(> 約10^4\)) ときにスティフな方程式ということが多いようである.

しかし, 硬度比が小さくてもスティフなケースや, 硬度比が大きくても初期値によってはスティフではなく解けるケースがあり, 文献 [3] では次のような定義も与えている.

非斉次項の変動量 M(t) を次のように定義する.
\[
M(t) = \frac{\|\phi'(t)\|}{\|\phi(t)\|}, \space ただし \space \phi(t) = 0 \space のときは \space M(t) = 0 \space とする
\] 次に減衰比 (damping ratio) を次のように定義する.
\[
\frac{|Re(\lambda_{max})|}{M(t)}
\] このとき次のように定義する.

– 減衰比が大きな問題を, その t の近傍でスティフな系という.
– 硬度比が大きいが, 減衰比がそれよりずっと小さい問題は疑似的にスティフな系という.

上の例 2 は硬度比は 2000 で, 減衰比は t により変動し \(1 \sim 約1999^2\) となっている. これはスティフな系といえる.

10.3.2.3 スティフな方程式の解法

スティフな方程式を解くためにはステップ幅に制約がない解法, すなわち, A 安定な解法を使う必要がある.

A 安定な解法としては後退オイラー法や陰的台形則があるが, 実用的にはもっと高次の公式が必要であり種々の解法が研究されてきた. 以下, 代表的な解法として, 後退微分公式 (BDF) および陰的ルンゲ・クッタ (IRK) 法について説明する. また, ローゼンブロック法と補外法についても少し触れる.

10.3.2.3.1 後退微分公式 (BDF)

アダムス・モルトン法では \((t_i, f_i) \space (i = n – k + 1, \dots, n + 1)\) を補間データとして \(f\) を補間したが, \(y\) そのものを補間する方法がある. すなわち, \((t_i, y_i) \space (i = n – k + 1, \dots, n + 1)\) を補間データとして補間多項式 \(q(t)\) を求め, それが少なくても格子点の 1 か所 (\(i = n + 1 – r\) の点) で微分方程式を満たすようにする.
\[
q'(t_{n+1-r}) = f(y_{n+1-r}, t_{n+1-r})
\] ここで, \(r = 1\) とすると陽解法になり, \(k = 1\) と \(k = 2\) の場合はそれぞれオイラー法と中点則になる. しかし, \(k \ge 3\) の場合には不安定になり微分方程式の解法として使うことはできない.

\(r = 0\) とすると陰解法になるが \(k \le 6\) で安定な解法が得られる. これは後退微分公式 (Backward differentiation formula (BDF) または ギヤ法 (Gear method)) とよばれ, k 次の公式になる.
\[
\sum_{j=1}^k (1/j) \nabla^j y_{n+1} = hf_{n+1}
\] ただし,
\[
\nabla^0 f_n = f_n, \space \nabla^{j+1} f_n = \nabla^j f_n – \nabla^j f_{n-1} \space (後退差分)
\] である.

\(k = 1 \sim 6\) について整理すると次のようになる.
\[
\begin{align}
& k = 1 : y_{n+1} – y_n = h f_{n+1} \space (後退オイラー法に一致) \\
& k = 2 : (3/2)y_{n+1} – 2y_n + (1/2)y_{n-1} = h f_{n+1} \\
& k = 3 : (11/6)y_{n+1} – 3y_n + (3/2)y_{n-1} – (1/3)y_{n-2} = h f_{n+1} \\
& k = 4 : (25/12)y_{n+1} – 4y_n + 3y_{n-1} – (4/3)y_{n-2} + (1/4)y_{n-3} = hf_{n+1} \\
& k = 5 : (137/60)y_{n+1} – 5y_n + 5y_{n-1} – (10/3)y_{n-2} + (5/4)y_{n-3} – (1/5)y_{n-4} = h f_{n+1} \\
& k = 6 : (147/60)y_{n+1} – 6y_n + (15/2)y_{n-1} – (20/3)y_{n-2} + (15/4)y_{n-3} – (6/5)y_{n-4} + (1/6)y_{n-5} = h f_{n+1} \\
\end{align}
\] BDF の絶対安定領域を \(k = 1 \sim 4\) について下図に示す.

図の曲線の外側が安定領域である. \(k = 1\) (後退オイラー法) と \(k = 2\) は A 安定であるが, それ以外では少し左半平面にはみだしている部分があり A 安定ではない. \(k\) が大きくなるほど安定領域は狭くなり, \(k \ge 7\) では原点でも不安定になる.

スティフな方程式を解く際に, 実際には多くの場合に完全に A 安定である必要はない. そこで, 条件を緩めた代わりに高次な解法として BDF がよく使われている. \(4 \sim 6\) 次の BDF は A 安定の条件を緩めた A(α) 安定である.

原点を頂点として, 左側に非有界であり, 実軸に対して上に α° 下に α° の扇型の領域が絶対安定領域に含まれる場合, A(α) 安定であるという. A(90°) 安定は A 安定に等しい. 方程式の係数行列の固有値がすべてこの中に入っているならば, 実質的に A 安定と同じである. 3 次の BDF は A (86°) 安定, 4 次の BDF は A(73°) 安定, 5 次の BDF は A(52°) 安定, 6 次の BDF は A(18°) 安定である.

左半平面がほとんど絶対安定領域に含まれるが, 原点付近を含まない虚軸近傍だけが絶対安定領域に含まれない場合を硬安定という. 虚軸に近いところに係数行列の固有値ない限り, これも実質的に A 安定と同じである. \(4 \sim 6\) 次の BDF は硬安定である.

数値実験 (23)

数値実験 (22) の例 2 のスティフな方程式を BDF を使って解く.

2 次 ～ 6 次の BDF を用いて h = 1 ～ 0.001 と変化させて計算したときの \(t = 5\) における \(y_2\) の相対誤差 (右に行くほど精度がよい) に対する関数評価回数を対数プロットしたものである. なお, 出発値は 3 次の IRK 法で計算した.

2 次 ～ 6 次まで問題なく解を求めることができた.

10.3.2.3.2 陰的ルンゲ・クッタ (IRK) 法

s 段ルンゲ・クッタ法の一般形は次のようなものであった.
\[
\begin{align}
& k_i = f(t_n + c_i h, y_n + h \sum_{j=1}^s a_{ij}k_j) \space (i = 1, 2, \dots, s) \\
& y_{n+1} = y_n + h \sum_{i=1}^s b_ik_i \\
\end{align}
\] \(a_{ij}, b_i\) および \(c_i\) はパラメータで, \(b_i\) および \(c_i\) は次式を満たす.
\[
\begin{align}
& c_i = \sum_{i=1}^s a_{ij} \space (i = 1, 2, \dots, s) \\
& \sum_{i=1}^s b_i = 1 \\
\end{align}
\] \(a_{ij} = 0 \space (i \le j)\) の場合, 陽解法になる. それ以外の場合, 陰的 (implicit) ルンゲ・クッタ (IRK) 法という. IRK 法は A 安定なものを作れるが, 1 ステップごとに連立非線形方程式を解かなければならないため計算の手間がかかるという難点がある.

(1) ガウス法

ガウス積分則に基づく最大次数の IRK 法をガウス法という. s 段ガウス法の次数は 2s 次である. \(c_1, c_2, \dots, c_s\) は区間 [0, 1] に変数変換されたルジャンドル多項式
\[
\frac{d^s}{dx^s} (x^s(x – 1)^s)
\] のゼロ点で与えられる.

クンツマン-ブッチャー法とよばれる4次と6次の公式の係数を以下に示す. この方法は A 安定である.

2 段 4 次クンツマン-ブッチャー法
\[
\begin{array}{c|cc}
1/2-√3/6 & 1/4 & 1/4-√3/6 \\
1/2+√3/6 & 1/4+√3/6 & 1/4 \\
\hline
& 1/2 & 1/2 \\
\end{array}
\]

3 段 6 次クンツマン-ブッチャー法
\[
\begin{array}{c|ccc}
1/2-√15/10 & 5/36 & 2/9-√15/15 & 5/36-√15/30 \\
1/2 & 5/36+√15/24 & 2/9 & 5/36-√15/24 \\
1/2+√15/10 & 5/36+√15/30 & 2/9+√15/15 & 5/36 \\
\hline
& 5/18 & 4/9 & 5/18 \\
\end{array}
\]

(2) ラダウ法とロバット法

ラダウとロバットの積分則に基づく IRK 法である. \(c_1, c_2, \dots, c_s\) は次の多項式のゼロ点で与えられる.
\[
\begin{align}
& \frac{d^{s-1}}{dx^{s-1}} (x^s(x – 1)^{s-1}) \space (Ⅰ型: ラダウ右) \\
& \frac{d^{s-1}}{dx^{s-1}} (x^{s-1}(x – 1)^s) \space (Ⅱ型: ラダウ左) \\
& \frac{d^{s-2}}{dx^{s-2}} (x^{s-1}(x – 1)^{s-1}) \space (Ⅲ型: ロバット) \\
\end{align}
\] Ⅰ型のラダウⅠA 法 (\(c_1 = 0\) になる), Ⅱ型のラダウⅡA 法 (\(c_s = 1\) になる) の次数は s 段のとき 2s – 1 次である. 両方とも A 安定である. 1 段ラダウⅡA 法は陰的オイラー法に一致する.

2 段 3 次ラダウⅡA 法
\[
\begin{array}{c|cc}
1/3 & 5/12 & -1/12 \\
1 & 3/4 & 1/4 \\
\hline
& 3/4 & 1/4 \\
\end{array}
\]

3 段 5 次ラダウⅡA 法
\[
\begin{array}{c|ccc}
(4-√6)/10 & (88-7√6)/360 & (296-169√6)/1800 & (-2+3√6)/225 \\
(4+√6)/10 & (296+169√6)/1800 & (88+7√6)/360 & (-2-3√6)/225 \\
1 & (16-√6)/36 & (16+√6)/36 & 1/9 \\
\hline
& (16-√6)/36 & (16+√6)/36 & 1/9 \\
\end{array}
\]

Ⅲ型の方法はロバットⅢA, ⅢB, ⅢC の３つがある. s 段のとき 2s – 2 次となる. いずれも A 安定である.

以下にロバットⅢA 法の例を示す。

3 段 4 次ロバットⅢA 法
\[
\begin{array}{c|ccc}
0 & 0 & 0 & 0 \\
1/2 & 5/24 & 1/3 & -1/24 \\
1 & 1/6 & 2/3 & 1/6 \\
\hline
& 1/6 & 2/3 & 1/6 \\
\end{array}
\]

4 段 6 次ロバットⅢA 法
\[
\begin{array}{c|cccc}
0 & 0 & 0 & 0 & 0 \\
(5-√5)/10 & (11+√5)/120 & (25-√5)/120 & (25-13√5)/120 & (-1+√5)/120 \\
(5+√5)/10 & (11-√5)/120 & (25+13√5)/120 & (25+√5)/120 & (-1-√5)/120 \\
1 & 1/12 & 5/12 & 5/12 & 1/12 \\
\hline
& 1/12 & 5/12 & 5/12 & 1/12 \\
\end{array}
\]

数値実験 (24)

数値実験 (22) の例 2 のスティフな方程式を IRK 法を使って解く.

2 段 4 次クンツマン-ブッチャー法 (ガウス法), 3 段 6 次クンツマン-ブッチャー法 (ガウス法), 2 段 3 次ラダウⅡA 法, 3 段 5 次ラダウⅡA 法, 3 段 4 次ロバットⅢA 法, 4 段 6 次ロバットⅢA 法を用いて \(h = 1 \sim 0.001\) と変化させて計算したときの \(t = 5\) における \(y_2\) の相対誤差 (右に行くほど精度がよい) に対する関数評価回数を対数プロットしたものである.

IRK 法ではステップごとに n x s 元非線形連立方程式 (この例では 2 x s 元) を解かなければならず, BDF に比べると s 倍大きい方程式になる. 関数 f() の計算が比較的容易な場合には非線形連立方程式を解く計算量は無視できなくなる. そのため, 図のように関数評価回数だけを計算量の指標にはできず, 他の指標 (例えば CPU 時間など) を使用すべきかもしれない.

10.3.2.3.3 対角陰的ルンゲ・クッタ法とローゼンブロック法

IRK 法ではステップごとに n x s 元非線形連立方程式を解かなければならないが, 場合によっては大きな方程式になってしまうのでこの負担を軽減する工夫が検討されている.

\(a_{ij} = 0 \space (i < j)\) で \(a_{ii} \ne 0\) (少なくとも 1 つの i について) としたものを対角陰的 (diagonally implicit) ルンゲ・クッタ (DIRK) 法という. 半陰的 (semi-implicit) ルンゲ・クッタ法ともいう. さらに, すべての対角成分が同一の値である場合, 単純対角陰的 (singly diagonally implicit) ルンゲ・クッタ (SDIRK) 法という.

DIRK 法の場合, s 個の n 元非線形方程式を逐次解くことにより解を求めることができるようになる. n 元非線形方程式をニュートン法で解くとき \(a_{ii}\) を含む係数の線形連立方程式を解く必要があるが, SDIRK 法であれば LU 分解を 1 回で済ませて結果を繰り返し使うことができるようになる.

DIRK 法は次のように表される.
\[
\begin{align}
& k_i = f(y_n + h \sum_{j=1}^{i-1} a_{ij}k_j + a_{ii}k_i) \space (i = 1, 2, \dots, s) \\
& y_{n+1} = y_n + h \sum_{i=1}^s b_ik_i \\
\end{align}
\] 工夫を進めて次のように変えたのがローゼンブロック法 (Rosenbrock method) である.
\[
\begin{align}
& (I – h \gamma_{ii}J)k_i = f(y_n + h \sum_{j=1}^{i-1} \alpha_{ij}k_j) + hJ \sum_{j=1}^{i-1} \gamma_{ij}k_j \space (i = 1, 2, \dots, s) \\
& y_{n+1} = y_n + h \sum_{i=1}^s b_ik_i \\
\end{align}
\] \(\alpha_{ij}, \gamma_{ij}, b_i\) は係数, \(J = f'(y_n)\) である. DIRK 法を線形化したものとみることができるので線形陰的 (linearly implicit) ルンゲ・クッタ法ともよばれる. 各ステップでは \(k_i\) を未知数とし \((I – h\gamma_{ii}J)\) を係数とする線形連立方程式を解くことで計算を進められるようになる.

これに属する種々の公式が提案されており, 公式により A 安定あるいは A(α) 安定である.

10.3.2.3.3 補外法

補外法をスティフな方程式に適用するために GBS アルゴリズムの線形陰的な (ローゼンブロックタイプの) 拡張が研究されている.
\[
\begin{align}
& (I – h_jJ)(y_1 – y_0) = h_jf(t_0, y_0) \space [線形陰的オイラー法] \\
& (I – h_jJ)(y_{i+1} – y_i) = -(I + h_jJ)(y_i – y_{i-1}) + 2h_jf(t_i, y_i) \space [線形陰的中点則] \\
& T_{j1} = (1/2)(y_{n_{j-1}} + y_{n_{j+1}}) \\
\end{align}
\] \(J\) はヤコビ行列の近似で, \(J = 0\) とおくと GBS アルゴリズムに一致する.

線形陰的中点則や線形陰的オイラー法を使った補外法もいくつか実装されており, α が 90° に近い A(α) 安定である.

10.3.3 微分代数方程式

次のように常微分方程式と代数方程式が連立した方程式を微分代数方程式 (Differential algebraic equations (DAE)) という.
\[
\begin{align}
& x’ = f(t, x, z) \\
& 0 = g(t, x, z) \\
\end{align}
\] 代数方程式は制約条件を表し, この制約のもとで常微分方程式を解くとみることができる. 代数制約式 \(g()\) を解析的に微分と消去を繰り返し \(z’ = h(t, x, z)\) のような形に変換すれば, すべての未知数について陽的な連立常微分方程式が得られる (方程式が特異でなければ可能である). この変換に必要な微分回数を微分代数方程式の指数という.

10.3.3.1 単振り子の例

代数方程式の例として単振り子の計算を行う.

数値実験 (25)

下図のように, 長さ \(l[m]\) の糸に重さ \(m[kg]\) のおもりがついた振り子を考える. 糸と鉛直線の角度が \(\theta[rad]\) で重力加速度を \(g[m/s^2]\) とする.

この振り子の運動は \(\theta\) を使って次の 2 階常微分方程式で表される.
\[
\theta^{\prime\prime} = -(g/l)sin(\theta), \space \theta(0) = \theta_0
\] \(\theta_0\) で手を離すと振り子は周期運動を開始し, その周期 \(T\) は次のようになることがわかっている.
\[
T = 4\sqrt{\frac{l}{g}}K sin(\frac{\theta_0}{2})
\] ただし, \(K\) は第 1 種完全楕円積分である.

この問題は特別な形をした 2 階常微分方程式なので (10.2.5.2 参照) ニュストレム法 (XLPack の Dopn43) を使って解くことができる. もちろん, 1 階連立形に変換して Derkfa を使って解くこともできて, 同じ結果が得られる.　初期値を \(\theta(0) = \pi/4\) とし, t = 0 ～ 5 における解を求め 0.1 ごとにプロットすると次のようになる.

振り子の動きをアニメーションで表示すると次のようになる.

さて, 同じ問題を \(q_1 = l sin(\theta), q_2 = -l cos(\theta)\) として直交 (x – y) 座標系で表す.

この場合の方程式は次のようになる. \(z\) は張力に関する変数である.
\[
\begin{align}
& mq_1^{\prime\prime} = -zq_1 \\
& mq_2^{\prime\prime} = -zq_2 – mg \\
& q_1^2 + q_2^2 = l^2 \\
\end{align}
\] この方程式は \(\theta\) を使って変数変換し z を消去すれば上の \(\theta\) を使った常微分方程式に戻すことができるが, ここではこのまま直交座標系で計算を行う方法を考える.

\(v_1 = q_1′, \space v_2 = q_2′, \space \lambda = z/m\) と変数変換して 1 階連立形に書き直すと次のようになる.
\[
\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 \\
\end{align}
\] これは微分代数方程式の一般形になっている. 最後の式が代数制約式であり, 位置の制約条件を表す (物理的にはおもりの軌道が半径 l の円周上にあることを示す). これを 1 回微分すると次のようになる.
\[
q_1q_1′ + q_2q_2′ = q_1v_1 + q_2v_2 = 0
\] これは速度の制約条件を表す. もう 1 回微分すると次のようになる.
\[
q_1v_1′ + q_2v_2′ + v_1^2 + v_2^2 = -\lambda – gq_2 + v_1^2 + v_2^2 = 0
\] これは加速度の制約条件を表す. これをもう 1 回微分すると \(\lambda’\) の式が得られるのでこの微分代数方程式の指数は 3 である.

10.3.3.2 微分代数方程式の解法 (1)

微分代数方程式は指数が高いほど解きにくいので, 制約式の解析的微分を繰り返し指数を下げて解くのが 1 つの方法である.

単振り子の制約式を 2 回微分して得られた \(\lambda\) を元の連立方程式に代入すると次の 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}
\] この連立常微分方程式は普通に解くことができて, 初期値を \(q_1(0) = l sin(\pi/4), \space q_2(0) = -l cos(\pi/4)\) とし, Derkfa を使って計算すると \(\theta\) を使った上図と等価の x-y 座標で表した解が得られる.

微分代数方程式の制約式を微分して得られた常微分方程式は「不変性を持つ常微分方程式」といわれ, 制約条件が陽には示されていないが連立常微分方程式の中に含まれている. このような方程式には安定性の問題があり, 計算を進めていったときに誤差の累積により制約条件が満たされなくなっていくことがある.

途中で誤差が蓄積してきたと仮定して \(q_2\) の初期値を \(-l cos(\pi/4) + 0.05\) として (本来このような初期値は許されないが) 安定性の問題をシミュレートしてみる. Derkfa で計算し x-y 平面での振り子の軌跡をプロットすると, 下図のように制約条件から外れた軌跡を描いた. 赤線は制約条件を満たす本来の軌跡である.

これを安定化する方法の一例として, \(x’ = f(t, x)\) の代わりに次の方程式を解く方法がある. \(h(x)\) は不変方程式で, \(h(x) = 0\) のときに制約条件を満たす. 従って, 制約条件を満たしていればこの方程式は元の方程式と同じ解を持つ.
\[
x’ = f(t, x) – \gamma F(x)h(x)
\] 詳細は省略するが, 今考えている単振り子の問題の場合, 具体的には次のようにするとよい (文献[4]を参照のこと).
\[
\begin{align}
& h(x) =
\begin{pmatrix}
q_1^2 + q_2^2 – l^2 \\
q_1v_1 + q_2v_2 \\
\end{pmatrix} \\
& F(x) = D(HD)^-1 \\
& D^T = H =
\begin{pmatrix}
2q_1 & 2q_2 & 0 & 0 \\
v_1 & v_2 & q_1 & q_2 \\
\end{pmatrix} \\
\end{align}
\] \(\gamma = 10\) として上図と同様に計算すると次のようになった.

制約条件を満たす方向に安定化しているのがわかる.

このように, 微分代数方程式は, 問題の定式化にいくつかのやり方がある場合には指数の小さなものを選び, 指数が大きければ解析的に微分を繰り返すことにより指数を下げることにより常微分方程式に変換して解くのがよい. なお, 問題によっては制約条件を満たすように安定化を図った方がよいことがある.

10.3.3.3 微分代数方程式の解法 (2)

微分代数方程式を直接離散化すると非常にスティフな方程式になる. しかし, 微分代数方程式を直接離散化して解くことができるプログラムもある. 指数を下げるなどの再定式化には手間がかかるので, これを使って解けるのであればその方がよい. しかし, 指数 1 にしか使えないプログラムも多いので, この場合でも指数が小さくなるように定式化しておくのがよいことは変わらない.

単振り子の方程式を再掲する.
\[
\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}
\] 最後の式については, 微分を繰り返すと指数 2 と 1 の式が得られる.
\[
\begin{align}
& 0 = q_1v_1 + q_2v_2 \space [指数 2 の式]\\
\end{align}
\] \[
\begin{align}
& 0 = -\lambda – gq_2 + v_1^2 + v_2^2 \space [指数 1 の式]\\
\end{align}
\] XLPack の Radaua は指数 1 ～ 3 の微分代数方程式の問題を直接離散化して解くことができる. また, Rodasa と Seulexa は指数 1 の問題であれば直接離散化して解くことができる.

上の指数 1 ～ 3 の 3 つの式を使用した数値実験結果が 10.4.4 (1) に掲載されている.