陰的ルンゲ=クッタ法

もし利用されるのであれば、実際のプログラムは
陰的ルンゲ=クッタ法の高速化
にして下さい。このページのプログラムは上のページのプログラムに比べて100倍くらい遅いです。


陰的ルンゲ=クッタ法とはどのように計算するのか? を主眼において説明します。

陽的ルンゲ=クッタ法は
ルンゲ=クッタ法の説明と刻み幅制御
4次ルンゲ=クッタ法
をご覧ください。

  1. 陽的解法と陰的解法
  2. 一階の微分方程式
  3. 複数の微分方程式
  4. Runge-Kutta-Gauss-Legendre法のButcherテーブル
  5. プログラム
  6. 刻み幅制御
  7. 参考文献

陽的解法と陰的解法


陽的ルンゲ=クッタ法のButcherテーブルは以下の通り

テーブルの右上と対角要素の値がゼロです。

しかし、陰的ルンゲ=クッタ法のButcherテーブルは

であり、すべての要素に値が入っています。

陽的解法は\(k_i\)を求めるためには\(k_{j\lt i}\)の要素が分かっていれば順繰りに計算できますが、
陰的解法は\(k_i\)を求めるために\(k_j\)の全ての要素が分かっていなければなりません。すなわち、陰的解法では1ステップ進めるために自己無撞着の連立方程式を解かなくてはいけないことを意味しています。

陰的解法の利点は、硬い方程式をスムーズに解くことが出来る点や、シンプレクティック性を持つ公式が存在する点です。

では、陰的公式の数値解法について述べて行きます。
急に任意の数の連立微分方程式を考えるのは複雑なので、方程式の数が一つだけの問題をまず考えて、その後一般化して考えます。

一階の微分方程式


問題

を初期条件\(y(x_0)=y_0\)の下で数値的に解き、\(y(x_0+h)\)の値を数値的に得ることを考えます。
\(s\)段の陰的ルンゲ=クッタ法を用いることを考えると、解は(1.2b)の形で書くことが出来ます。

ここで、係数\(k_i,~~(i=1,2,\cdots,s)\)は

を満たします。式(1.3)を式変形して、関数のゼロ点を探す問題に焼きなおします。すなわち、

の通り、関数\(g_i\)のゼロ点を探すのです。

関数\(g_i\)のゼロ点を与える式(1.4)の解を\(\mathbf{k}^{(\infty)}\)と書き、
解の近傍に存在する\(\mathbf{k}^{(0)}\)が刻み幅\(h\)のオーダーに収まる範囲を考えます。すなわち、

の範囲を考えます。\(\mathbf{k}^{(0)}\)を初期値とすれば、\(\mathbf{k}^{(0)}\)は\(h\)のリーディングオーダーの範囲で

と与えられます。ここで、関数\(y(x)\)もオーダー\(y(x_0+h)-y(x_0)\in O(h^1)\)の範囲にあることを仮定しています。

この想定の下ではNewton-Raphson法を用いることが出来て、

を繰り返して解が収束するまで計算を繰り返せば良く、初期値は\(\mathbf{k}^{(0)}\)です。
ここで略記した表記は

を意味します。また、行列Jは

で表される行列です。

\(J\)の行列要素を求める事はさておいて、Newton-Raphson法を用いるためには

解いて\(\mathbf{z}\)を求めなければなりません。\(J\)の逆行列を求める事は数値的にメリットは無いので、行列のLU分解からこの問題を連立方程式を解く問題に変えます。
式(1.12)のように変形します。

行列が下三角行列\(L\)と上三角行列\(L\)の積で書けるとします。すなわち、行列\(J\)を

の通りに分解するのです。連立方程式

を解くために、LU分解より、

です。\(U\mathbf{z}=\mathbf{w}\)と置けば、

をまず解いて、\(\mathbf{w}\)を得ます。その後、

を解くことで\(\mathbf{z}\)を得ます。

行列要素の計算

続いて、行列要素を計算することを考えます。

もう少し計算進めると、行列要素は

となります。ここで、\(\delta_{n,m}\)はKroneckerのデルタです。

簡易Newton法

行列\(J\)の要素をNewton法の繰り返し度に毎回計算するのは大変です。
なので、\(O(h^1)\)の範囲で関数の符号などが大きく変わらないと想定します。

となります。
すると、行列要素\(J\)は

と表されます。すなわち、Newton法の繰り返し毎に行列\(J\)を計算する必要はなく、

ということです。

複数の微分方程式


続いて以下のN本の微分方程式を解くことを考えます。

関数\(g\)のゼロ点を探す問題に焼きなおすと、

となります。
Newton法で繰り返すことを考えると、

です。それぞれ

を満たします。行列\(J\)は、

であり、\(\mathbf{k}\)の初期値は

です。行列要素\(J\)は

です。簡易Newton法では

として得られます。

ルンゲ=クッタ=ガウス=ルジャンドル法のButcherテーブル


ガウス=ルジャンドル求積法に基づくルンゲ=クッタ法のButcherテーブルは以下のものです。

2段4次

3段6次

このガウス=ルジャンドル求積法に基づく陰的ルンゲ=クッタ法は凄まじく素晴らしい特性を持ちます。

  • A安定である[4]
  • シンプレクティックである[3,4]
  • 対称である[4]
  • 硬い方程式に向く[1,3,4]
  • 刻み幅を途中で変更してもシンプレクティック性が保たれる[3]

至れり尽くせりの性質を持ちますね。とりあえずこれを使っておけば問題ないということです。

プログラム


Fortran90による、N本連立方程式を陰的ルンゲ=クッタ法で解くプログラムです。
簡易Newton法を用いています。

lapackを用いるので、mklかlapackにリンクしてコンパイルしてください。

3段6次 陰的Runge-Kutta

2段4次 陰的Runge-Kutta

さて、実際に解析解と刻み幅を変えて精度の比較をしましょう。
正しく実装できていれば刻み幅を1桁減らすと解析解との一致桁数が4桁、6桁合うようになるはずです。
横軸を刻み幅、縦軸を解析解との相対誤差にとったものが以下のグラフです。
問題はばねの問題で、\(x=100\)で比較しました。

確かに刻み幅を1桁減らすと解析解との一致桁数が4桁、6桁一致するようになっています。
上記計算は倍精度で行っています。なので、3段6次公式は倍精度演算にとってはオーバースペック気味です。
なぜなら、刻み幅が0.05辺りで桁あふれのような振る舞いが見えてしまっているからです。4倍精度演算になればちょうどいいくらいでしょう。

ここ以降で私が書いている部分はあまり良くない説明とプログラムです。
もし刻み幅制御のプログラムが欲しければ、
陰的ルンゲ=クッタ法の高速化
を参照してください。

刻み幅制御


3段6次の公式には別の次数が埋め込まれていますが、それは2次のオーダーです。ちょっと6次のオーダーに対して2次のオーダーを用いるのはもったいない気がします。

ここでは、良くない方法ですが、別の次数を得るために別の陰的公式を用いることを考えます。
別の陰的公式として、LobattoⅢCと呼ばれる3段4次の公式を採用しました。
端点の値が評価されるという性質があるからです。この性質が無いと不連続点をうまく検知できませんでした。
刻み幅制御のやり方は陽的Runge-Kutta-Fehlbergの方法と同じにしましたので、ここでは書きません。
詳しくはルンゲクッタ法の説明と刻み幅制御 -シキノートを参照してください。
陰的公式の正しい刻み幅制御の方法は文献[1]に載っていますので、本格的に知りたい方はそちらを参照してください。

刻み幅を制御した後、実際に返される値は3段6次の結果です。
懸念事項があります。3段6次のGauss-Legendre求積法に基づく公式はシンプレクティック積分ですが、
刻み幅を変えてしまうとシンプレクティック性が失われているのではないか?という事項です。
私が見聞きしている限り、シンプレクティック性を保つためには等間隔刻み幅をもちいなければならないはずです。

@adhara_mathphysさんによるコメントを戴きまして、Gauss-Legendre求積法に基づいた陰的Runge-Kutta公式は刻み幅の変化に依らずシンプレクティック性が保たれるようです。


またシンプレクティック性を保つ積分法は、全力学的エネルギーが正準変数の二次形式で書けるならどんなに時間発展してもエネルギーが保存されるという性質を持ちます。ばねの問題であればこの性質を満たし全力学的エネルギーは保存しますが、ケプラー問題は満たさないため、シンプレクティック性は保ちますが、全力学的エネルギーは保存しません。

3段6次と3段4次を利用した刻み幅制御陰的Runge-Kutta法

大規模計算にはメモリ等を考えていないので向いていません。ご注意ください。

3段6次:Gauss–Legendre methods(実際に返される値)
3段4次:LobattoⅢC (誤差を見積もるためだけに利用)

下のプログラムは
\(
\frac{dy}{dx}=y \cos(x),~~y(0)=1
\)

を解くものです。解析解は
\(
y(x)=e^{sin(x)}
\)

となります。


用途に合わせたその他プログラム

その他、必要と思われるプログラムをまとめて公開します。
4つに分けました。

プログラムの名づけ方

倍精度 d
4倍精度 q

LAPACK使用 l
GECP使用 g

openMP使用 m
不使用 n

自動間隔出力 a
等間隔出力 e

例えば、倍精度でLAPACKを用いて、openMPを使用して等間隔出力する場合のプログラム名は
main_dlme.f90という名前にしてあります。

デフォルトではばねの問題
\(
\begin{align}
\frac{d^2}{dt^2}y(t) = -y(t)
\end{align}
\)

を解く問題になっています。
main_dlma.f90
main_dlme.f90
main_dgma.f90
main_dgme.f90
main_qgma.f90
main_qgme.f90

GECPとは

LU分解と連立方程式の解法を、LAPACKを使わないで、
九州大学の渡部 善隆様が公開なさっているGECP(Gaussian Elimination with Complete Pivoting, 一般実行列に対する連立1次方程式の数値解を完全ピボット選択付き Gauss の消去法によって求める Fortran サブルーチン)
に置き換えたものです
http://yebisu.cc.kyushu-u.ac.jp/~watanabe/RESERCH/GECP/index.html
再配布可能なので、上記プログラムで必要な物には組み込んであります。

参考文献


[1]E. Hairer and G. Wanner, ‘Solving Ordinary Differential Equations II’ Springer, 1996

[2]ルンゲ=クッタ法のリスト

[3]J. M. Sanz-Serna, ‘Runge-Kutta schemes for Hamiltonian systems’, https://www.researchgate.net/publication/200175547_Runge-Kutta_schemes_for_Hamiltonian_systems

[4]遠藤 理平, 数値計算法概論 常微分方程式の解法, http://www.natural-science.or.jp/article/20120323225814.php

その他参照したページ


遠藤 理平, 数値計算法概論 常微分方程式の解法 http://www.natural-science.or.jp/article/20120323225814.php

田中正次, 三村和正, 山下茂, 陰的Runge-Kutta法の特性について(常微分方程式の数値解法) https://repository.kulib.kyoto-u.ac.jp/dspace/bitstream/2433/100230/1/0643-4.pdf

渡部 善隆, 一般実行列に対する連立1次方程式の数値解を完全ピボット選択付き Gauss の消去法http://yebisu.cc.kyushu-u.ac.jp/~watanabe/RESERCH/GECP/index.html


「陰的ルンゲ=クッタ法」への5件のフィードバック

  1. こんにちは。いつもありがたく拝見させていただいております。にょいと申します。
    ひとつだけ分からない点があります。初期値k^(0)はどのように考えれば良いのでしょうか?
    お時間がありましたら少しヒントだけでもお教え願えればと存じます。

    1. にょい様

      初期値k^(0)はニュートン法の初期値の事ですかね?
      そうでしたら、式(1.6)に従い決定しています。
      すなわち、
      dy/dx=f(x,y), y(x=x0)=y0
      に対し、初期値は全て
      k^(0)=f(x0,y0)
      と置いています。
      お尋ねしていただいた内容と合っていますか?

      追記)
      どのようにk^(0)を決めているのか?
      という事でしたら、式(1.3)より、
      h→0を考えた時の解がk^(0)=f(x0,y0)だから、
      ということです。

  2. sikino様

    ご返信ありがとうございます。
    ご回答いただいた内容ではっきり分かりました。
    式(1.6)を見逃しておりました。
    どのようにk^(0)を決めているのかということも、
    はっきりと理解することができました。

    (ある問題を自作のNewmark-β法+予測子修正子法で
    解こうとした所、相当に細かい刻みにしても、一部の
    パラメータの場合でしか結果が得られず、困っておりました。
    結果を得るにも長時間かかったり、PCがフリーズしたり・・・
    今度は陰的Runge-Kutta法を試してみたいと思います!)

    1. にょい様

      内容があっていて一安心です。
      Newmark-β法とは初めて聞きました。面白いトピックで調べてみようと思います。
      ありがとうございます。
      解こうとしている問題が硬い方程式であれば
      陰的解法が適しているので、試してみてください。

      もしも精度や正しさが求められるのであれば、
      http://slpr.sakura.ne.jp/qp/irkgl-program/#howtoradau5
      に述べてある、私ではなく専門家が書いたRADAU5
      というプログラムがあるのでそちらを参考にすると良いかもしれません。

      硬い方程式とは、
      例えば2つの時間に関する微分方程式A, Bがあり、
      Aは数秒で変化、Bは数時間で時々変化するような方程式の事です。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です