陰的ルンゲ-クッタ法 #

1. 陽的/陰的ルンゲ-クッタ法の概要 #

下記の微分方程式を、適当な初期条件の元で数値的に解くことを考えます。 \begin{align}\label{rkeq} \frac{dy}{dx}=f(x,y),~~~y(x_\text{in})=y_\text{in} \end{align}

$p$段のルンゲ-クッタ法を用いると、微小量$h$だけ離れた関数の値$y(x_\text{in}+h)$は、 \begin{align} \label{rkeq_y0h} y(x_\text{in}+h)=y_\text{in}+h\sum_{i=1}^p b_i k^{(i)} \end{align} で与えられます。ここで、$k^{(i)}$は \begin{gather} k^{(i)}\equiv f\left(x_\text{in}+c_i h, y_\text{in} + h\sum_{j=1}^p a_{i,j}k^{(j)}\right) \end{gather} を満たし、係数$a_{i,j}, b_{i}, c_{i}$はブッチャー配列によって与えられる既知の係数です。

係数$a_{i,j}$の値によって、陽的解法か陰的解法かが区別されます¹。 ルンゲ-クッタ法の解き方の大きな区別として、下記のように分類されます。

その他、半陰的ルンゲ-クッタ法, 対角陰的ルンゲ-クッタ法 (DIRK) ($j\gt i$に対して、$a_{i,j}=0$)という方法があり、上記の中間くらいの解法も有ります。本稿では示しませんので、各自お調べいただければと思います。

本稿では、完全陰的ルンゲ-クッタ法（以降、単に陰的ルンゲ-クッタ法といいます）の計算はどのように行われているかについて示します。

1.1. ルンゲ-クッタ法の段数と次数 #

同じルンゲ-クッタ法でも陽的/陰的の違いがあることを指摘しましたが、更に段数と次数という違いもあることも指摘しておきます。

つまり、“ルンゲ-クッタ法を使用します” といっても性質を把握できるわけではなく、段数、次数の区別も存在するわけです。

ざっくりと下記のように分けられます。

• 段数: 計算の大変さ（段数が多いほど大変）
• 次数: 計算の正確さ（次数が大きいほど正確）²

たとえ陽的/陰的、段数、次数が同一でも、まだ一意にルンゲ-クッタ法を指定することはできません。 これはルンゲ-クッタ法は条件の数が未知数の数より少ないので、無数の指定方法が存在するためです。

そのため、“ルンゲ-クッタ法を使用する"と言う際は、幾つかの有名なブッチャー配列についている特別な名前を指定するか、ブッチャー配列そのものを全て指定するのが良いでしょう。

2. 陰的ルンゲ-クッタ法の計算方法 #

これから本稿で説明する内容は、p段q次の陰的ルンゲ-クッタ法です。

2.1. 概要 #

陰的ルンゲ-クッタ法とはいえ、式\eqref{rkeq_y0h}を解けばよいことは変わりません。

陽的解法とは異なり、$p$段陰的ルンゲ-クッタ法で大変なのは下記の係数$k^{(i)},~(i=1,2,\cdots,p)$を求める箇所であり、ここをどのように解くかが数値計算上の問題となります。 \begin{gather} k^{(i)}=f\left(x_\text{in}+c_i h, y_\text{in} + h\sum_{j=1}^p a_{i,j}k^{(j)}\right),~~~(i=1,2,\cdots,p) \end{gather}

ここで、新たに関数$g^{(i)},~~(i=1,2,\cdots,p)$を定義して、$g^{(i)}$のゼロ点を探す問題にします。
つまり、下記の方程式を同時に満たす$\mathbf{k}\equiv (k^{(1)}, k^{(2)}, \cdots, k^{(p)})$を探します。

\begin{equation} \left\{ \begin{aligned} &g^{(1)}(\mathbf{k})\equiv k^{(1)}-f\Bigl(x_\text{in}+c_1 h, y_\text{in}+h\sum_{j=1}^{p}a_{1,j}k^{(j)}\Bigr)=0 \\ &g^{(2)}(\mathbf{k})\equiv k^{(2)}-f\Bigl(x_\text{in}+c_2 h, y_\text{in}+h\sum_{j=1}^{p}a_{2,j}k^{(j)}\Bigr)=0 \\ &\hspace{6em}\vdots \\ &g^{(p)}(\mathbf{k})\equiv k^{(p)}-f\Bigl(x_\text{in}+c_p h, y_\text{in}+h\sum_{j=1}^{p}a_{p,j}k^{(j)}\Bigr)=0 \end{aligned}\label{gk_irk} \right. \end{equation}

また、略記のために$\mathbf{g}=(g^{(1)}, g^{(2)}, \cdots, g^{(p)})$を定義すると、 \eqref{gk_irk}は

\begin{align}\label{gk0_irk} \mathbf{g}(\mathbf{k})=\mathbf{0} \end{align}

と書くことができます。ここで、$\mathbf{k}=(k^{(1)}, k^{(2)}, \cdots, k^{(p)})$と置きました。

今、これを満たす$\mathbf{k}$を求めたい訳です。

さて、非線形の連立方程式\eqref{gk0_irk}を解けば良いですが、非線形な方程式なのでそのまま解くのは難しいです。 なので、ある点周りでは線形方程式として扱えると仮定して解いていきましょう。

本稿では、メジャーな方法であるニュートン-ラフソン法を想定して非線形方程式を解きますが、参考文献³のような他の手法でも構いません。

まず、方程式\eqref{gk0_irk}には自明な解$\mathbf{k}=\mathbf{k}_0=(k_0^{(1)},k_0^{(2)},\cdots,k_0^{(p)})$が存在します。それは、$h=0$の場合で、 \begin{align} k_0^{(i)}=f(x_\text{in}, y_\text{in}),~~(\text{if }h=0) \end{align} です。

もちろん、$h=0$では元の方程式と初期条件をそのまま代入しただけで、微分方程式をまるで解いたことになりません。 ですが、$h\to 0$の場合を考えるうえで非常に役に立ちます。

$h\to 0$の解は、$h=0$の解から大きく離れないことが想像できます。 つまり、$h=0$の解$\mathbf{k}_0$と、$h$が有限の解$\mathbf{k}$は、ランダウの記号を用いて

\begin{align} \label{gk_irk2} |\mathbf{k}-\mathbf{k}_0| \sim O(h^{1}) \end{align}

の範囲にあると仮定します。別の言い方をすれば、点$\mathbf{k}=\mathbf{k}_0$の周りで線形であると考えて、解もこの線形の範囲内に存在すると考えるのです。 逆に言えば、もし$h$が大きくて線形の範囲から逸脱する場合、おおよそ線形の範囲になるまで$h$を小さくします。

仮定\eqref{gk_irk2}の下では、式\eqref{gk0_irk}は例えばニュートン-ラフソン法を用いることで解を求めることができます。 つまり、ニュートン-ラフソン法の$l$回目の反復において更新される解を$\mathbf{k}_l$と書くと、

\begin{align} \label{nr_irk} \mathbf{k}_{l+1}=\mathbf{k}_{l}-J_{l}^{-1}\mathbf{g}_{l},~~(l=0,1,\cdots) \end{align}

であり、解が\eqref{gk_irk2}の範囲にある限り、反復$l$を十分大きくすれば解\eqref{gk0_irk}が得られます。 \eqref{nr_irk}の右辺は、通常は$J_{l}$の逆行列をあらわに求めることはせず、$J_{l}$のLU分解を利用して解きます⁴。

その後、$\mathbf{k}$を式\eqref{rkeq_y0h}に代入して$y(x_\text{in}+h)$を得ることができます。

また、式\eqref{nr_irk}において$g_{l}\equiv g(\mathbf{k}_l)$と置き、 また$J_{l}\equiv J(\mathbf{k}_{l})$は点$\mathbf{k}_{l}$における$\mathbf{g}_l$のヤコビ行列

\begin{align} J_l= \left( \begin{array}{*4{>{\displaystyle}c}} \frac{\partial g^{(1)}}{\partial k^{(1)}}&\frac{\partial g^{(1)}}{\partial k^{(2)}}&\cdots &\frac{\partial g^{(1)}}{\partial k^{(p)}} \\ \frac{\partial g^{(2)}}{\partial k^{(1)}}&\frac{\partial g^{(2)}}{\partial k^{(2)}}&\cdots &\frac{\partial g^{(2)}}{\partial k^{(p)}} \\ \vdots&\vdots&\ddots&\vdots \\ \frac{\partial g^{(p)}}{\partial k^{(1)}}&\frac{\partial g^{(p)}}{\partial k^{(2)}}&\cdots &\frac{\partial g^{(p)}}{\partial k^{(p)}} \\ \end{array} \right) \end{align}

と定義しました。$J$の行列要素を$J_{i,j}$と書けば、

\begin{eqnarray} \label{jac} J_{i,j} &=& \frac{\partial}{\partial k^{(j)}} \left[k^{(i)}-f\Bigl(x_\text{in}+c_i h, y_\text{in} + h\sum_{s=1}^{p}a_{i,s}k^{(s)}\Bigr)\right] \\ &=& \delta_{i,j}-ha_{i,j}\frac{\partial f(x,y)}{\partial y}_{\substack{\begin{aligned}x&=x_\text{in}+c_ih,\\ y&=y_\text{in}+h\sum\limits_{r=1}^{p}a_{i,r}k^{(r)}\end{aligned}}} \end{eqnarray}

と求められます。ここで$\delta_{i,j}$はディラックのデルタ関数で、下記の関数です。

\begin{align} \delta_{i,j}= \left\{ \begin{aligned} &1,\hspace{2em} (i=j) \\ &0,\hspace{2em} (i\ne j) \end{aligned} \right. \end{align}

2.1.1. ブッチャー配列 #

ブッチャー配列（Butcher tableau）は、ルンゲ-クッタ法\eqref{rkeq_y0h}で現れる多くの係数をひとまとめにした配列です。
例えば下記のような配列（表）で記載されます。

ここで示した係数($a_{i,j}, b_{i}, c_{i}$)は\eqref{rkeq_y0h}で説明した係数、 $\tilde{b}_i$は同じブッチャー配列で異なる次数を与える$b_i$であり、 微小量$h$を制御する場合に利用されます。

ルンゲ-クッタ法では、自由度の数が制限の数より多いため、同じルンゲ-クッタ法でも複数のブッチャー配列を取り得ます。 どのブッチャー配列を利用するかによって、同じルンゲ-クッタ法でも、その性質が大きく変わります。

例えば、時間発展する物理学の問題を解く場合、解きたい問題が満たすべき性質を考慮してブッチャー配列を選び、利用します⁵。 その性質は、例えば下記のような性質です。

• 段数、次数（計算の大変さ、正確さ）
• 陰的解法、陽的解法（硬い方程式に適する、適さない）
• 安定性（数値的に安定に解けるか？）⁵
• シンプレクティック性（保存されるべき保存量が保存される解き方か？）
• 時間反転対称性（時間を遡るように解いたとき、元の値に戻ることが保証される解き方か？）
• 計算コスト（計算が大変か？）

様々なブッチャー配列の中でも、特に優秀な性質を持つ配列を利用したルンゲ-クッタ法が、本稿で紹介するガウス-ルジャンドル配列を用いた方法です。

2.1.2. ガウス-ルジャンドル配列 #

下記のブッチャー配列で与えられる3段6次のガウス-ルジャンドル配列は、数値計算的に非常に喜ばしい性質を持ったアルゴリズムです。

このガウス-ルジャンドル配列の持つ性質は、下記の通りです⁵。

• p段で2p次の精度を持つ
• 陰的解法なので、硬い方程式に適する
• A-安定である
• シンプレクティック性を保つ（たとえ途中で刻み幅を変更しても保つ）
• 時間反転対称性を持つ

ただし先に示した通り$\mathbf{k}$を求めるために非線形方程式を解かねばならず、計算コストが非常に掛かるため、必要に応じて採用しましょう。

3. 連立一次微分方程式の場合 #

今まで、一次の微分方程式だけを考えていましたが、$N$本の微分方程式の場合を考えておきましょう。 解きたい問題は下記のものです。

\begin{align} \frac{d\mathbf{y}}{dx}=\mathbf{f}(x,\mathbf{y}),\hspace{2em}\mathbf{y}(x_\text{in})=\mathbf{y}_\text{in} \\ \end{align} この時、$x_\text{in}$から微小区間$h$進んだ位置の関数値$\mathbf{y}(x_\text{in}+h)$は、 \begin{align} \mathbf{y}(x_\text{in}+h) = \mathbf{y}_\text{in} + h\sum_{i=1}^{p}b_i \mathbf{k}^{(i)} \end{align} となります。ここで、$\mathbf{y}, \mathbf{f}$は$N$次元のべクトルで、下記のように定義しました。

\begin{align} \mathbf{y}= \left( \begin{array}{*4{>{\displaystyle}c}} y_1\\ y_2\\ \vdots \\ y_N \end{array} \right), \hspace{3em} \mathbf{f}= \left( \begin{array}{*4{>{\displaystyle}c}} f_1(x,\mathbf{y})\\ f_2(x,\mathbf{y})\\ \vdots \\ f_N(x,\mathbf{y}) \end{array} \right) \end{align} ルンゲ-クッタ法の係数、関数はそれぞれの段で$N$次元のベクトルとなります。$i$段目の係数、関数をそれぞれ$\mathbf{k}^{(i)}, \mathbf{g}^{(i)}$と書いて \begin{align} \mathbf{k}^{(i)}&\equiv\mathbf{f}\left(x_\text{in}+c_i h, \mathbf{y}_\text{in} + h\sum_{j=1}^{p} a_{i,j}\mathbf{k}^{(j)}\right) \\ \mathbf{g}^{(i)}(\mathbf{k}; h) &\equiv \mathbf{k}^{(i)} - \mathbf{f}\left(x_\text{in}+c_i h, \mathbf{y}_\text{in} + h\sum_{j=1}^{p} a_{i,j}\mathbf{k}^{(j)}\right)=\mathbf{0} \label{Neqgi_irk} \end{align} と定義します。具体的に係数を書きますと、下記のようになります。 \begin{align} \mathbf{k}^{(i)}\equiv \left( \begin{array}{*4{>{\displaystyle}c}} k_{1}^{(i)}\\ k_{2}^{(i)}\\ \vdots \\ k_{N}^{(i)} \\ \end{array} \right), \hspace{3em} \mathbf{k}\equiv \left( \begin{array}{*4{>{\displaystyle}c}} \mathbf{k}^{(1)}\\ \mathbf{k}^{(2)}\\ \vdots \\ \mathbf{k}^{(p)} \\ \end{array} \right) =\left( \begin{array}{*4{>{\displaystyle}c}} k_{1}^{(1)}\\ k_{2}^{(1)}\\ \vdots \\ k_{N}^{(1)} \\ k_{1}^{(2)} \\ k_{2}^{(2)} \\ \vdots \\ k_{N}^{(p)} \end{array} \right) \end{align}

$\mathbf{k}^{(i)}$は$N$次元、$\mathbf{k}$は$N\times p$次元のベクトルとなります。

式\eqref{Neqgi_irk}は$i$段目だけです。そのため、全ての係数$\mathbf{k}$についてゼロ点を探したいので \begin{align} \label{Neqg_irk} \mathbf{g}(\mathbf{k}; h) \equiv \mathbf{k} - \mathbf{F}(\mathbf{k}; h)=\mathbf{0} \end{align} と書いて$N\times p$次元のベクトル$\mathbf{g}$を定義し、この関数のゼロ点を探すことが目的となります。 ここで、式\eqref{Neqg_irk}において、 \begin{align} \mathbf{F}(\mathbf{k}; h)\equiv \left( \begin{array}{*4{>{\displaystyle}c}} \mathbf{f}\Bigl(x_\text{in}+c_1 h, \mathbf{y}_\text{in} + h\sum\limits_{j=1}^{p} a_{1,j}\mathbf{k}^{(j)}\Bigr) \\ \mathbf{f}\Bigl(x_\text{in}+c_2 h, \mathbf{y}_\text{in} + h\sum\limits_{j=1}^{p} a_{2,j}\mathbf{k}^{(j)}\Bigr) \\ \vdots \\ \mathbf{f}\Bigl(x_\text{in}+c_p h, \mathbf{y}_\text{in} + h\sum\limits_{j=1}^{p} a_{p,j}\mathbf{k}^{(j)}\Bigr) \\ \end{array} \right) =\left( \begin{array}{*4{>{\displaystyle}c}} f_1\Bigl(x_\text{in}+c_1 h, \mathbf{y}_\text{in} + h\sum\limits_{j=1}^{p} a_{1,j}\mathbf{k}^{(j)}\Bigr) \\ f_2\Bigl(x_\text{in}+c_1 h, \mathbf{y}_\text{in} + h\sum\limits_{j=1}^{p} a_{1,j}\mathbf{k}^{(j)}\Bigr) \\ \vdots \\ f_N\Bigl(x_\text{in}+c_1 h, \mathbf{y}_\text{in} + h\sum\limits_{j=1}^{p} a_{1,j}\mathbf{k}^{(j)}\Bigr) \\ f_1\Bigl(x_\text{in}+c_2 h, \mathbf{y}_\text{in} + h\sum\limits_{j=1}^{p} a_{2,j}\mathbf{k}^{(j)}\Bigr) \\ f_2\Bigl(x_\text{in}+c_2 h, \mathbf{y}_\text{in} + h\sum\limits_{j=1}^{p} a_{2,j}\mathbf{k}^{(j)}\Bigr) \\ \vdots \\ f_N\Bigl(x_\text{in}+c_p h, \mathbf{y}_\text{in} + h\sum\limits_{j=1}^{p} a_{p,j}\mathbf{k}^{(j)}\Bigr) \\ \end{array} \right) \end{align} を定義しました。

ニュートン-ラフソン法を使用するための初期値は、$h\to 0$が解近傍の$O(h)$の範囲にある解と考えて \begin{align}\label{Nkin_irk} \mathbf{k}^{(i)}_0 = \mathbf{f}(x_\text{in}, \mathbf{y}_\text{in}),~~~(i=1,2,\cdots,p) \end{align} と書くことができます。

$l$回目の繰り返しにおけるニュートン-ラフソン法は、下記のように書くことができます。 \begin{align} \mathbf{k}_{l+1} = \mathbf{k}_{l} - {J_{l}}^{-1} \mathbf{g}_{l} \end{align} ここで$J$はヤコビ行列であり、多次元の場合は下記のように書けます。 \begin{align} J = \frac{\partial \mathbf{g}}{\partial \mathbf{k}} &= \begin{pmatrix} \frac{\partial \mathbf{g}^{(1)}}{\partial \mathbf{k}^{(1)}} & \frac{\partial \mathbf{g}^{(1)}}{\partial \mathbf{k}^{(2)}} & \ldots & \frac{\partial \mathbf{g}^{(1)}}{\partial \mathbf{k}^{(p)}} \\ \frac{\partial \mathbf{g}^{(2)}}{\partial \mathbf{k}^{(1)}} & \frac{\partial \mathbf{g}^{(2)}}{\partial \mathbf{k}^{(2)}} & \ldots & \frac{\partial \mathbf{g}^{(2)}}{\partial \mathbf{k}^{(p)}} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial \mathbf{g}^{(p)}}{\partial \mathbf{k}^{(1)}} & \frac{\partial \mathbf{g}^{(p)}}{\partial \mathbf{k}^{(2)}} & \ldots & \frac{\partial \mathbf{g}^{(p)}}{\partial \mathbf{k}^{(p)}} \end{pmatrix} =\begin{pmatrix} \frac{\partial g_1^{(1)}}{\partial k_1^{(1)}} & \frac{\partial g_1^{(1)}}{\partial k_2^{(1)}} & \ldots & \frac{\partial g_1^{(1)}}{\partial k_N^{(1)}} & \frac{\partial g_1^{(1)}}{\partial k_1^{(2)}} & \ldots && \frac{\partial g_1^{(1)}}{\partial k_N^{(p)}} \\ \frac{\partial g_2^{(1)}}{\partial k_1^{(1)}} & \frac{\partial g_2^{(1)}}{\partial k_2^{(1)}} & \ldots & \frac{\partial g_2^{(1)}}{\partial k_N^{(1)}} & \frac{\partial g_2^{(1)}}{\partial k_1^{(2)}} & \ldots && \frac{\partial g_2^{(1)}}{\partial k_N^{(p)}} \\ \vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ \frac{\partial g_N^{(1)}}{\partial k_1^{(1)}} & \frac{\partial g_N^{(1)}}{\partial k_2^{(1)}} & \ldots & \frac{\partial g_N^{(1)}}{\partial k_N^{(1)}} & \frac{\partial g_N^{(1)}}{\partial k_1^{(2)}} & \ldots && \frac{\partial g_N^{(1)}}{\partial k_N^{(p)}} \\ \frac{\partial g_1^{(2)}}{\partial k_1^{(1)}} & \frac{\partial g_1^{(2)}}{\partial k_2^{(1)}} & \ldots & \frac{\partial g_1^{(2)}}{\partial k_N^{(1)}} & \frac{\partial g_1^{(2)}}{\partial k_1^{(2)}} & \ldots && \frac{\partial g_1^{(2)}}{\partial k_N^{(p)}} \\ \vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ \frac{\partial g_N^{(p)}}{\partial k_1^{(1)}} & \frac{\partial g_N^{(p)}}{\partial k_2^{(1)}} & \ldots & \frac{\partial g_N^{(p)}}{\partial k_N^{(1)}} & \frac{\partial g_N^{(p)}}{\partial k_1^{(2)}} & \ldots && \frac{\partial g_N^{(p)}}{\partial k_N^{(p)}} \\ \end{pmatrix} \end{align}

ヤコビ行列の$u$行$v$列の要素$J_{u,v}$は、下記の通りに定義し計算できます。 \begin{align} J_{u,v}\equiv J^{(i,j)}_{n,m} &\equiv \frac{\partial g_n^{(i)}}{\partial k_m^{(j)}} \\ &=\delta_{i,j}\delta_{n,m}-ha_{i,j} \frac{\partial f_n}{\partial y_m}\Biggl|_{\substack{\begin{aligned}x&=x_\text{in}+c_ih,\\ \mathbf{y}&=\mathbf{y}_\text{in}+h\sum\limits_{r=1}^{p}a_{i,r}\mathbf{k}^{(r)}\end{aligned}}} \end{align} ここで、添え字の範囲は \begin{align} n&=1,2,\cdots,N,~~~m=1,2,\cdots,N\\ i&=1,2,\cdots,p,~~~j=1,2,\cdots,p\\ u&=1,2,\cdots,Np,~~~v=1,2,\cdots,Np \end{align} です。

4. 計算結果とプログラム #

実際に計算した結果を示したいわけですが、陰的ルンゲ-クッタ法を使用して微分方程式を解く際に、ブッチャー配列をガウスルジャンドルに設定するためのライブラリがどうやらなかなか無いようです。 そのため作りました。Pythonによるプログラムは以下からダウンロードできます。

(MIT license)
IRK 3-stage 6-order Gauss-Legendre Butcher tableau
https://slpr.sakura.ne.jp/sikinote/supplement/python/numeric/IRK/IRK_GL36.zip

中には、下記のプログラムが入っています。

• main.py: 可読性を重視
• main_opt.py: 高速化を重視（main.pyより10%程度早いです）

実際に6次の精度を持っているかだけを確認し、陰的ルンゲ-クッタ法の確認としましょう。

微分方程式 \begin{align} \frac{d^2y}{dx^2}=-y, \end{align} を初期条件$y(0)=0,~y’(0)=1$の下で$x=100$まで解いていきます。この時、刻み幅を変化させたときに解析解$y=\sin(x)$との相対誤差を確かめてみます。 実際に計算すると、下記のようになります。

この時、今回の問題に対し計算が正常にできていると思われる刻み幅$h$は、おおよそ$[1\sim 0.02]$位だと分かります。 これより刻み幅$h$が小さいと精度が16桁（倍精度計算の10進数のおおよその桁数）を超えてしまうため、これ以上小さくしても意味がないです。つまり、倍精度演算にとって性能が良すぎるため、ただただ計算量が増えるという領域になります。

を$10^{-1}$倍すると、誤差が$10^{-6}$倍となるため、$h$について6次の精度を持っていることが分かります。

5. 付録 #

5.1. 数値計算におけるガウス-ルジャンドル #

数値計算の分野でガウス-ルジャンドルと言うといくつか存在しますので、調査する際は気を付けましょう。一例を挙げます。

• 陰的ルンゲ-クッタ法のブッチャー配列の一つとしてガウス-ルジャンドル配列
• 数値積分のアルゴリズムとしてガウス-ルジャンドル求積法
• 円周率を求めるアルゴリズムとしてガウス-ルジャンドル法

5.2. 陰的解法の効率化 #

第2, 3章で陰的ルンゲ-クッタ法の計算方法について説明しました。ここからは数値計算の効率化に関する話です。

陰的ルンゲ-クッタ法の理論は、4章までで終わっています。
単にルンゲ-クッタ法を利用したいだけの方は、以降の話は全く必要ありません。

以降は、数値計算の高速化の話です。

陰的解法は面倒です。 特に陰的解法では次元が少し増えるだけで、ニュートン法とヤコビ行列の要素計算で特に重い計算となります。

これらの計算効率を向上させるために、過去の研究者たちがさまざまな工夫を凝らしてきました。 ここでは、その中からニュートン法の実行に関係する下記の効率的な数値計算法を紹介します。

• 節5.1) ニュートン法の逆行列を含む項を、LU分解を利用して計算する
• 節5.2) ヤコビ行列の要素$J_{u,v}$を近似（簡易ニュートン法）

その他にも下記の効率化が有りますが、ここでは詳細は延べません。更なる効率化については参考書⁶をご覧ください。

• 解く問題を、差分だけに対する問題に変換してから陰的ルンゲ-クッタ法を適用する
• ニュートン法の初期値を、過去の結果を補外した点とする（$h\to 0$の解よりも初期値が解に近いことを期待）
• ニュートン法の終了条件を、過去の計算結果を利用して決定する
• 埋め込まれた次数$\tilde{b}$を利用して、要求精度を満たす最大の刻み幅$h$で計算する

5.2.1. LU分解 #

まずはニュートン法の反復計算\eqref{nr_irk}で生じる逆行列の計算を考えます。つまり、下記の右辺第二項です。

\begin{align} \mathbf{k}_{l+1}=\mathbf{k}_{l}-J_{l}^{-1}\mathbf{g}_{l},~~(l=0,1,\cdots) \end{align}

右辺第二項は下記の連立一次方程式の解$\mathbf{x}$を求める事と同じです。

\begin{align} A\mathbf{x}=\mathbf{b},~~~\mathbf{x}=A^{-1}\mathbf{b} \end{align}

ここで、$A$は$N\times N$の行列、$\mathbf{x}, \mathbf{b}$は$N\times 1$の列ベクトルです。 数値計算の基本ですが、逆行列$A^{-1}$をあらわに求めるのは以下の理由で避けられます。

• 逆行列を求める場合の計算量は、通常$O(N^3)$となり非常に重い。
• 行列に特異性（逆行列が存在しない or 困難）がある場合は、計算誤差が溜まりやすい

逆行列を求めない代わりに、下記の方法で連立一次方程式が解かれます⁷。

※1. 疎行列の場合、疎行列LU分解と呼ばれる方法があります。
※2. コレスキー分解（正定値対称行列に対する解法）もあります。
※3. 反復解法として、GMRES法や共役勾配法（正定値対称行列に対する解法）などが有ります。

本稿ではLU分解を利用することを考えます。LU分解とは、行列$A$を下三角行列$L$と上三角行列$U$の積に分解する事を言います。 もし$A$が$L$, $U$に分解することができれば、連立方程式は前進代入と後退代入を用いて容易に解くことができます。

つまり、$A$のLU分解が実行出来た場合、 \begin{align} A\mathbf{x}=\mathbf{b} \to LU\mathbf{x}=\mathbf{b} \end{align} と書くことができます。ここで、新たに$\mathbf{y}\equiv U\mathbf{x}$を定義して、$\mathbf{y}$について連立方程式を解けば \begin{align} L\mathbf{y}=\mathbf{b} \to \text{get } \mathbf{y} \end{align} となります。$L$は下三角行列なので前進代入（計算量$O(N^2)$）で容易に解けます。更に、$\mathbf{x}$について解けば、 \begin{align} U\mathbf{x}=\mathbf{y} \to \text{get } \mathbf{x} \end{align} となります。$U$は上三角行列なので後退代入（計算量$O(N^2)$）で容易に解けます。

LU分解自体の計算量は、先に示した通り$O(\frac{2}{3}N^3)$ですが、一度LU分解をしてしまえば異なる$\mathbf{b}$について$O(N^2)$で解くことができるため、問題によっては非常に良い方法となります。

5.2.1.1. $PA=LU$の形によるLU分解 #

さて、数値計算プログラムによっては行列$A$のLU分解は、 \begin{align} A=LU \end{align} ではなく、 \begin{align} PA=LU \end{align} と、置換行列$P$を含んで求められる場合があります（例えばpythonのscipy.linalg.lu()をデフォルトで使用する場合）。 この$P$は、“数値的不安定さを解消するためのピポット選択に伴う行の入れ替え” を意味しています。

$PA=LU$によってLU分解がなされた場合、元の方程式を変形して$P$を追加する形にする必要があります。

例えば連立方程式 \begin{align} A\mathbf{x}=\mathbf{b} \end{align} を解きたい時、$PA=LU$によってLU分解されるプログラムを利用するならば、両辺に$P$を掛けて \begin{align} PA\mathbf{x}=P\mathbf{b} \to LU\mathbf{x}=\mathbf{b}’ \end{align} として、$\mathbf{b}’=P\mathbf{b}$を新たに計算して解く必要があります。

上手く変形できない場合も出てくるかもしれません。通常、プログラムにはそういったオプションを指定することができると思いますので、通常のLU分解を行いましょう（scipy.linalg.lu()であれば、オプションpermute_l=Trueに指定して使用すれば良いです）。

5.2.2. 簡易ニュートン法 #

ニュートン法による連立方程式の解法を節5.1で説明しましたが、行列$J$の行列要素$J_{u,v}$の計算も大変です。 \begin{align} J_{u,v}&\equiv J^{(i,j)}_{n,m} \\ &=\delta_{i,j}\delta_{n,m}-ha_{i,j} \frac{\partial f_n}{\partial y_m}\Biggl|_{\substack{\begin{aligned}x&=x_\text{in}+c_ih,\\ \mathbf{y}&=\mathbf{y}_\text{in}+h\sum\limits_{r=1}^{p}a_{i,r}\mathbf{k}^{(r)}\end{aligned}}} \end{align}

さて、仮定\eqref{Nkin_irk}より初期値$\mathbf{k}_{l}$と解$\mathbf{k}_{l+1}$は今、$O(h)$の範囲にあると仮定しています。 そのため、ニュートン法で使用する傾きは同じと仮定しても良いのではないか？という想像ができます。

具体的に低次元の場合で比較してみますと、下記のような違いが生じます。

ニュートン法	簡易ニュートン法

簡易ニュートン法でも、反復回数は増えるかもしれませんが収束しそうです。 そのため、行列要素の計算が辛いことが分かっている場合、簡易ニュートン法を繰り返す方が合計の計算時間が早くなります。

つまり簡易ニュートン法を採用すると、行列要素の計算は下記のように計算することができます。 \begin{align} J^{(i,j)}_{n,m} \approx\delta_{i,j}\delta_{n,m}-ha_{i,j} \frac{\partial f_n}{\partial y_m}\Biggl|_{\substack{\begin{aligned}x=x_\text{in},\\ \mathbf{y}=\mathbf{y}_\text{in}\end{aligned}}} \end{align}

まとめると、簡易ニュートン法を採用することで下記のメリットとデメリットが生じます。

• メリット
  • $\frac{\partial f_n}{\partial y_m}$が$i,j$に依存しなくなったため、$N^2$個の偏微分だけ計算すれば良くなった（元々は$N^2p^2$個だった）
  • $\mathbf{k}$に依存しなくなったので、ニュートン法の繰り返しの度に再計算する必要が無くなった
• デメリット
  • ニュートン法の反復回数が増加する

刻み幅を大きくすると簡易ニュートン法がうまく働く条件から外れてしまい、収束までの反復回数が増えてしまいますが、刻み幅を小さくしすぎるとルンゲ-クッタ法の反復回数が増えてしまいます。

まとめると下記のようなメリット・デメリットがあります。

6. 参考文献 #