はじめに
みなさん、Runge-Kutta使ってますか？＞挨拶
Runge-Kuttaは死ぬほど亜種がありますが、いわゆる古典的なRunge-Kutta法は4次のアルゴリズムになっています。これは実装が簡単で精度が高いので、数値計算の教科書には必ず書いてあるようなメジャーな手法ですが、これが4次精度になっていることを確認するのは相当面倒くさいです1。
数値計算屋は、一生に一度は古典Runge-Kuttaの精度確認をするべきですが、私はこの歳になるまでそれをサボっていました。しばらく前から、なぜかTwitterなどでRunge-Kutta法が話題になったので、良い機会だと思って精度確認をしようとしたのですが、実際に「そら」で計算しようとすると、途中で計算がとっちらかって嫌になります。計算に腕力が必要なのもそうですが、見通しよく計算しないとごちゃごちゃになるでしょう。本稿では、古典4次Runge-Kuttaが4次精度であることを、何が面倒なのかを見ながら示そうと思います。
Runge-Kutta法とは
以下のような微分方程式を考えましょう。
$$
\frac{dy}{dt} = f(y)
$$
時刻$t$における物理量$y$の時間微分係数が、$y$自身の関数$f(y)$で与えられています。
一般には$f(y,t)$と時間依存しても良いのですが、ここでは時間微分が時間に非依存である場合のみを考えます2。
ある時刻$t=0$における初期値$y(0)$が与えられた時に、任意の時刻$t$における$y(t)$の値が我々の知りたい量です。
さて、もしこの微分方程式が求積できればそれが答えです。
$$
y(t) = y(0) + \int_0^t f(y) dt
$$
しかし、一般にはこの積分は解析的に実行できないため、なんらかの近似をします。
具体的には、小さい時間刻み$h$をとって、その時間だけ時間発展させます。
$$
y(t+h) = y(t) + \int_t^{t+h} f(y) dt
$$
この積分を近似することで、「次のステップ」の$y$の値を更新できます。
最も単純な近似は、$h$が小さければそのあいだ微分値$f(y)$が変化しない、とするものでしょう。
すると積分の外に出すことができます。
$$
y(t+h) = y(t) + h f(y(t))
$$
本当は時刻$t$から$t+h$まで時々刻々と変化する微分係数$f(y)$を、時刻$t$の値で代表させた、と考えることができます。これが一次のEuler法です。これは非常に精度が悪く、あっという間に計算が破綻します。
現在時刻$t$における微分係数を使って積分を評価するのはさすがに乱暴に過ぎました。そこで、一次のオイラー法で、時刻を半分だけ進めた場所を推定し、その場所での微分値を使うことにしましょう。以後、簡単のため$y(t)$を$y$、$f(y)$を$f$と略記します。
\begin{align}
k_1 &= f \\
k_2 &= f\left(y + \frac{h k_1}{2} \right) \\
y(t+h) &= y(t) + h k_2
\end{align}
これは中点法と呼ばれ、二次精度になります。
ここでは一度「途中の点」を推定し、その場所での微分係数を考えることで精度を上げることができました。このような推定を繰り返すほど、精度が上げられるような気がします。このような「途中の点」を4点とるのが4次のRunge-Kutta法です。
広く使われる4次の古典Runge-Kutta法のアルゴリズムはこんな感じです。
\begin{align}
k_1 &= f \\
k_2 &= f\left(y + \frac{h k_1}{2} \right) \\
k_3 &= f\left(y + \frac{h k_2}{2} \right) \\
k_4 &= f\left(y + h k_3 \right) \\
y(t+h) &= y(t) + \frac{h}{6} \left( k_1 + 2k_2 +2k_3 + k_4\right)
\end{align}
ここで、最後の線形和は、先程の積分を$h$の4次まで厳密に近似しています。
$$
\int_t^{t+h} f(y) dt = \frac{h}{6} \left( k_1 + 2k_2 +2k_3 + k_4\right) + O(h^5)
$$
要するに線形多段法とは、$t$から$t+h$までの微分係数の平均を、最も良く「直線」で近似しましょう、という方法論です。このあたりの議論はそれこそ数値計算の教科書に多数書いてあると思われるので、そちらを参照してください。
4次精度の証明
厳密なテイラー展開
さて、多くの教科書では、中点法くらいまでは精度確認をするかもしれませんが、4次のRunge-Kuttaは非常に計算が煩雑になるために、「これは4次精度になる」とだけ書いてスルーします。これを真面目に計算しましょう、というのが本稿の目的です。
まず、「4次まで厳密」とは、$y(t+h)$をテイラー展開した時の係数が4次まで一致する、という意味です。とりあえず展開してみましょう。
$$
y(t+h) = y + h y' + \frac{h^2}{2} y''+ \frac{h^3}{6} y''' + \frac{h^4}{24} y''' + O(h^5)
$$
普通、時間微分はドットで表すことが多いですが、ここでは趣味でプライムを使っています。この微分係数$y', y'', y'''$を$f$を使って表現していきます。まず$h$の一次の項ですが、$y$が満たす微分方程式から、$y' = f$です。
次に$y''$ですが、これは$f$の時間微分になります。ここで、$f(y(t))$と、中に$y$自身を含んでいるため、微分すると合成関数の微分公式から中身の$y$の微分が出てきます。
\begin{align}
y'' &= \frac{d}{dt} f\\
&= \frac{df}{dy} \frac{dy}{dt} \\
&= f_y f
\end{align}
ここで$f$の$y$微分を$f_y$と書きました。
次に$y'''$を求めるため、$y''$を微分します。すると、関数の積の微分なので
\begin{align}
y''' &= \frac{d}{dt} y''\\
&= f_y' f + f_y f' \\
&= f_{yy} f^2 + f f_y^2
\end{align}
となります。$f$を微分すると$f f_y$が、$f_y$を微分すると$f_{yy} f$が出てくるのに注意しましょう。こうして、微分すればするほど$f_{yyy}$やら$f^2$やらがわらわら出てきます。これが Runge-Kuttaの精度確認のめんどくさポイントその1 です。まぁここは
- $f$を微分すると$f_y f$になる
- $f_y$ を微分すると$f_{yy} f$になる
- $f_{yy}$ を微分すると$f_{yyy} f$になる
と、下の添字の$y$が増えて$f$がかかる、と覚えてひたすら計算すれば、面倒ですがまだなんとかなります。こうして、4次までのテイラー展開は以下のように書けます。
\begin{align}
y(t+h) &= y \\
&+ hf \\
&+ \frac{h^2}{2} f_y f\\
&+ \frac{h^3}{6} (f_{yy} f^2 + f_y^2 f)\\
&+ \frac{h^4}{24} (f^3 f_{yyy} + 4 f^2 f_y f_{yy} + f f_{y}^3) \\
&+ O(h^5)
\end{align}
これを、後の便利のために以下のように整理しておきましょう。
\begin{align}
\frac{y(t+h) - y}{h} &= f \\
&+ \frac{h}{2} f_y f\\
&+ \frac{h^2}{6} (f_{yy} f^2 + f_y^2 f)\\
&+ \frac{h^3}{24} (f^3 f_{yyy} + 4 f^2 f_y f_{yy} + f f_{y}^3) \\
&+ O(h^4)
\end{align}
Runge-Kuttaの公式を見るとわかるように、左辺は$(k_1+2k_2+2k_3+k_4)/6$です。この対応を見るのが本稿の目的です。
Runge-Kutta法の精度確認
ではいよいよ4次の古典Runge-Kutta法の精度を確認しましょう。もういちど公式を確認します。
\begin{align}
k_1 &= f \\
k_2 &= f\left(y + \frac{h k_1}{2} \right) \\
k_3 &= f\left(y + \frac{h k_2}{2} \right) \\
k_4 &= f\left(y + h k_3 \right) \\
\frac{y(t+h)-y(t)}{h} &= \frac{\left( k_1 + 2k_2 +2k_3 + k_4\right)}{6}
\end{align}
つまり、我々が確認すべき式はこんな感じです。
\begin{align}
\frac{\left( k_1 + 2k_2 +2k_3 + k_4\right) }{6} &= f \\
&+ \frac{h}{2} f_y f\\
&+ \frac{h^2}{6} (f_{yy} f^2 + f_y^2 f)\\
&+ \frac{h^3}{24} (f^3 f_{yyy} + 4 f^2 f_y f_{yy} f + f_{yy}^3) \\
&+ O(h^4)
\end{align}
さて、Runge-Kuttaは4次精度ですが、微分係数は3次まで合えば良いので、3次まで展開して係数を比較すれば良いことになります。とりあえず$k_1$から$k_4$までを素直にテイラー展開してみましょう。
\begin{align}
k_1 &= f \\
k_2 &= f + \left(\frac{hk_1}{2}\right) f_y +
\frac{1}{2} \left(\frac{hk_1}{2}\right)^2 f_{yy} +
\frac{1}{6} \left(\frac{hk_1}{2}\right)^3 f_{yyy} + O(h^4)\\
k_3 &= f + \left(\frac{hk_2}{2}\right) f_y +
\frac{1}{2} \left(\frac{hk_2}{2}\right)^2 f_{yy} +
\frac{1}{6} \left(\frac{hk_2}{2}\right)^3 f_{yyy} + O(h^4) \\
k_4 &= f + (h k_3) f_y +
\frac{1}{2} (h k_3)^2 f_{yy} +
\frac{1}{6} (h k_3)^3 f_{yyy} +O(h^4)
\end{align}
ここで、$k_2$は$k_1$に、$k_3$は$k_2$に依存するため、順番に代入していく必要があります。これが Runge-Kuttaの精度確認のめんどくさポイントその2 です。そして、最後に$h$の0次、$h$の1次・・・と集めていって、厳密な展開と係数を比較するのですが、特に難しい計算はなく、ただただ展開していくだけなので苦行です。これはちょっと手でやるのは嫌です。でも我々にはMathematicaがあるじゃないか！
というわけで、ここはMathematicaにやらせましょう。テイラー展開を済ませてしまえば、あとは微分演算はありません。これは$f_y$は永遠に$f_y$のままというこです。つまり、$f$、$f_y$、$f_{yy}$といった変数を、互いに独立な変数として扱って良い、ということです。たとえばMathematicaで$f_y$をfy、$f_{yy}$をfyyといった感じで表現しましょう。
すると、3次までのテイラー展開のテンプレートは以下のように表現できます。
t3 := f + x fy + x^2/2 fyy + x^3/6 fyyy
このxを、例えば$hk_1/2$で置換してやれば、$k_2$を得ます。Mathematicaでの操作はこんな感じになります。
t3 := f + x fy + x^2/2 fyy + x^3/6 fyyy
k1 := f
k2 := t3 /. x -> h k1/2
k3 := t3 /. x -> h k2/2
k4 := t3 /. x -> h k3
最後に、$(k_1+2k_2+2k_3+k_4)/6$の$h$に関する4次までの展開を見てやります。
Series[(k1 + 2 k2 + 2 k3 + k4)/6, {h, 0, 3}]
実行結果はこんな感じです。
厳密なテイラー展開と係数が一致しているのがわかりますね。
Jameson-Baker Runge-Kutta法の精度確認
Runge-Kuttaには様々な亜種がありますが、そのうちの一つにLow-storage Runge-Kuttaというものがあります。これは中間変数を保存せずに実装できる線形多段法です。その中で、こちらの記事で紹介されていたJameson-Baker法の精度を計算してみましょう。
Jameson-Baker法は、先の記法を用いるとこんな感じになります。
\begin{align}
k_1 &= f \\
k_2 &= f\left(y+ \frac{h k_1}{4}\right) \\
k_3 &= f\left(y+ \frac{h k_2}{3}\right) \\
k_4 &= f\left(y+ \frac{h k_3}{2}\right) \\
\frac{y(t+h) - y}{h} &= k_4\\
\end{align}
最後に$k_1$から$k_4$の線形和を計算する必要があった古典Runge-Kuttaに比べ、$k_1$は$k_2$の計算だけに、$k_2$は$k_3$の計算だけに使われるため、どんどん上書きできます。したがって中間変数の値を保存する必要がなく、省メモリとなります。
こちらの精度確認も、Mathematicaを使えば一発です。
t3 := f + x fy + x^2/2 fyy + x^3/6 fyyy
k1 := f
k2 := t3 /. x -> h k1/4
k3 := t3 /. x -> h k2/3
k4 := t3 /. x -> h k3/2
Series[k4, {h, 0, 3}]
実行結果はこんな感じになります。
厳密な値と微妙に係数が違うのがわかるかと思います。なお、Jameson-Bakerは、問題が線形の場合、つまり$f_y = f_{yy} = 0$の場合には4次精度になります。
Out[6] /. {fyy -> 0, fyyy -> 0}
これは、$f_y = f_{yy} = 0$の時の厳密な展開と一致します。
まとめ
「数値計算屋なら一生に一度は4次のRunge-Kuttaの精度確認をしろ」という格言3に従い、これまでサボってきた精度確認をしました。中点法くらいならすぐですが、4次は以下の二点で面倒です。
- そもそも厳密な展開係数がちょっと面倒
- アルゴリズムの係数が入れ子構造になっているため、順番にテイラー展開していかないといけないのが面倒
前者はまだ異なる次数の項目が混ざらないので楽ですが、後者は例えば$h^2$の項が$k_2$、$k_3$、$k_4$、$k_2^2$、$k_3^2$、$k_4^2$のそれぞれから出てくるのですこぶる面倒です。正直これを手でやる気はしなかったので4Mathematicaを使いました。テイラー展開さえ手でやってしまえば、あとは置換して$h$のオーダーごとに整理するだけなので、たった6行です。Mathematica素晴らしい！
っていうか、もしMathematicaがなくても、手でまともに計算するより、簡単な数式処理のスクリプトを書いたほうが早いかもしれません。それくらい面倒でした。まぁ、Mathematicaを使ったにせよ、ちゃんと精度を確認できたことには違いないので、これでよしとします。
参考文献
- Runge-Kutta法についてのノート(早川尚男) 京大の早川先生によるノート。一般のRunge-Kutta法の満たすべき条件を丁寧に導出してあります。
- とても簡単な4段4次のRunge-Kutta法について この記事をきっかけにRunge-Kuttaの精度確認の重い腰を上げました。また、Low-storage Runge-Kutta法の存在もこの記事で知りました。
- 岩波講座 応用数学〈第13巻〉〔方法3〕 微分方程式の数値解法 I 三井斌友著 早川先生によると格調高い本とのことです。
- 数値解析の基礎 (篠原 能材 著) 今読み返してみたら、Runge-Kuttaの面倒な局所打ち切り誤差の計算を全段真面目に書いてありました。
謝辞
早川先生にはノートを教えていただきました。aokomoriutaさんに計算ミスを指摘していただきました。ありがとうございます。