「ルンゲ=クッタ法」カテゴリーアーカイブ

ケプラー問題に対する陽的解法と陰的解法

2体のケプラー問題を数値的に解きます。
ここでは、適切な変換をして求めるのではなく、刻み幅制御されたプログラムで無理やり計算します。

そして、
陽的解法であるルンゲ=クッタ=フェールベルグ法

陰的解法であるルンゲ=クッタ=ガウス=ルジャンドル法
の離心率に対する計算量の違いを調べてみます。

Kepler問題

二次元で二体の問題を考えます。運動方程式は
\(
\begin{align}
\frac{d^2 x}{dt^2}&=-\frac{x}{(x^2+y^2)^{3/2}} \\
\frac{d^2 y}{dt^2}&=-\frac{y}{(x^2+y^2)^{3/2}}
\end{align}
\)
\(t=[0,20],~~x(0)=1-e,~~x'(0)=0,~~y(0)=0,~~y'(0)=\sqrt{\frac{1+e}{1-e}}\)
です。ここで、\(e\)は離心率を表します。

この問題に対する解は良く知られていて、
\(
\displaystyle (x+e)^2+\frac{y^2}{1-e^2}=1
\)

で表され、
\(
0\le e\lt 1
\)

の時、楕円となります。

関数の評価回数の離心率の依存性


楕円の軌道を持つ範囲において、
計算は離心率\(e\)が1に近づくほど難しくなります。
なぜなら、原点付近を通過する際に、質点の導関数の変化が大きくなるからです。

使用したプログラムの説明は
陽的解法はhttps://slpr.sakura.ne.jp/qp/runge-kutta-ex/
陰的解法はhttps://slpr.sakura.ne.jp/qp/irkgl-program/
です。実際のプログラムも置いておきます。

離心率を変化させたときの軌道はこんな感じです。

さて、計算速度の評価ですが、関数が何回評価されたかで比較します。
念の為、陰的解法はLU分解もあるので単純な比較は難しいことを注記しておきます。

図の一番上は関数の評価回数の離心率依存性を表しています。縦軸は評価回数、横軸は\(1-e\)を表します。真ん中の図は一番上と同じですが、軸を対数にとっています。一番下は、\(t=20\)まで計算した時に\(t=0\)の初期エネルギーとの相対誤差\(|(E(t=20)-E(t=0))/E(t=0)|\)を表しています。

特徴的な振る舞いは、離心率に対して陰的解法の評価回数は線形の振る舞い、陽的解法は指数で振る舞っている点です。
これは、質点が原点の近くを通るような問題の際に違いが際立つ事を示しています。
また、エネルギーの保存に関しても陰的解法の方が良いことが分かるでしょう。

一応注意しておきますが、ここでいう陰的解法はルンゲ=クッタ=ガウス=ルジャンドル法の振る舞いです。一般的な陰的解法については話していないことに注意してください。

プログラム


きっかけ

きっかけとしてはtwitterで流れてきまして、やってみよう、と思いました。

陰的ルンゲ=クッタ法の高速化

陰的ルンゲ=クッタ法の基本的な考えは
陰的ルンゲ=クッタ法
にて説明しました。
基本的には
\(
y(x+h)=y(x)+\Delta y
\)

の右辺\(y(x)+\Delta y\)を計算します。
しかし、陰的ルンゲ=クッタ法の方法を数値計算で行おうとすると望ましくない部分が現れます。
それは、

  1. \(y(x)=O(h^0),~~\Delta y=O(h^1)\)を同時に扱わなければならず、桁落ちが激しい
  2. 関数の評価回数が多い
  3. ヤコビアン、LU分解の計算コストが非常に高い

という点です。
簡易ニュートン法を用いる事を前提にしておくと、上の問題は若干解決することが出来ます。

本稿では陰的ルンゲ=クッタ法を発展させ、数値計算的に陰的ルンゲ=クッタ法のアルゴリズムを工夫し、どのように計算量を減らすか?に焦点を当てます。

目次

  1. 計算方法
  2. プログラム
  3. プログラムの評価
  4. 結論
  5. RADAU5の使い方
  6. 参考文献

注意
[2,5]の同著者の新しい論文では、本稿の計算方法([1]に従う方法)ではなく、
もう一段階変換してからニュートン法を利用しています。
正直な所、私自身が追いきれなかっとのと、複素数が入ってくるプログラムでしたので、[1]の方法で止めておきます。
2つの計算方法の収束の早さや精度を比較を[2]で行っていますので、気になる方はそちらをご覧ください。

計算方法


解きたい問題は\(N_{\text{eq}}\)本の連立一次微分方程式

です。これを\(s\)段のルンゲ=クッタ法で解くことを考えます。
座標や添え字は以下のように決めました。

すると次のステップの値は

と求められます。\(z_{np}^{[i]}\)を求める事が問題になります。ここで、\(d_j,~~(j=1,2,\cdots,s)\)は既知で

のように、Butcherテーブルから求められます。
具体的に、3段6次の陰的公式であるGauss-Legendreの場合、

と求められます。
\(z_{np}^{[i]}\)の具体的な形は

で計算されます。\(z_{np}^{[i]}\)をベクトルとして表したいので、

と変換します。
この非線形方程式はニュートン法によって求める事が出来て、以下の通り、\(k\)回の繰り返しで解が収束するまで計算されます。

ここで、解が収束した結果が求めたい\(z_{np}^{[i]}\)となります。すなわち、

です。行列\(J’\)は

と求められます\((n,m=1,2,\cdots,N_{\text{eq}}, p,q=1,2,\cdots,s)\)。ここで、\(\delta_{i,j}\)はクロネッカーのデルタ、\(a_{i,j}\)はButcherテーブルの値、\(J\)はヤコビアンで、

と計算されます。また、ベクトル\(\mathbf{w}_{np}^{(k)}\)は、

と定義します。

初期値の推定


初期値は\(i=0\)の初めのステップでは

として求め、それ以降では過去の結果を多項式補間して求めるとかなりよい精度です[1]。


ニュートン法の停止


ニュートン法は以下の条件が満たされた時、終了します。

まず、\(i=0\)の初めのステップでは必ず\(k\ge 2\)まで行い、
その後、以下の条件が満たされた時、終了します。

ここで、\(\eta_k^{[i]}\)は

で定義され、特に\(\theta_k\)は収束割合(convergion rate)と呼ばれます。
また、\(\kappa=0.1\sim 0.01,~~Ntol=\text{min}(0.03,\sqrt{Rtol})\)です[5]。ここで、\(Rtol\)は相対誤差を意味します。
実際に組んでみますと、理論が違うせいなのか、\(Ntol=\text{min}(0.03,\sqrt{Rtol})\)では十分なほど収束はしませんでした。なので、実際に組む時には\(Ntol=10^{-9}\)にしてしまっています。もしかしたら、\(10^{-9}\)では足らず、もっと必要かもしれません。

\(i\ge 1\)のステップでは\(k=1\)の時は前ステップの\(\eta\)の値を使い、


で判定します。ここで、倍精度演算ならば\(Uround=5\times 10^{-16}\)です。

\(i\ge 1,~~k\ge 2\)では

に従い、計算します。

ニュートン法の繰り返しの\(k\ge 2\)において、どこかで

が満たされる様であれば、刻み幅が大きすぎて収束しないことを意味します。なので、刻み幅を小さくする必要があります。

誤差判定


区間\(i\)の\(n\)番目の方程式の解の誤差は、以下の連立方程式を解いて得られます[1]。

ここで、\(\gamma_0\)はButcherテーブルの行列Aの実固有値で、ガウス=ルジャンドル陰的ルンゲクッタの3段6次であれば、
\(
\gamma_0=0.215314423116112178244733530380696
\)

を得ます。そして、上の方程式を解いた後に、計算結果を棄却するか判定するために、量

を計算します。
もしも\(||\text{err}^{[i]}||\lt 1\)であれば、その刻み幅で計算した結果を採用します。
しかし、\(||\text{err}^{[i]}||\ge 1\)であれば、以下の連立方程式を解きます([1]の”Hump”問題を参照)。

上を計算し、\(||\text{err}^{[i]}||\)を再計算した結果、もしも\(||\text{err}^{[i]}||\lt 1\)であれば、その刻み幅で計算した結果を採用し、そうでなければ刻み幅を次の節に従って変更します。

刻み幅制御


刻み幅制御をするためには, 2つ新しい刻み幅の推定値である\(h_{i+1}^{(1)},h_{i+1}^{(2)}\)を計算します。

その結果、刻み幅が小さくなるか大きくなるかに従って、どちらの刻み幅を採用するか決定します[2]。

プログラム


Fortran90で書いたプログラムはこちら。LU分解と連立方程式を解くため、LAPACKを使います。
それにリンクしてコンパイル、実行をしてください。
モジュールを使用していますが、これは関数の呼び出し回数を計測するためだけにグローバル変数として使っているため、消してもプログラムに何の影響もありません。

追記)
色々計算してみました。その結果、\(tol=10^{-8}\)より小さい値は使わない方が良さそうです。どうもこれ以上の精度にしてしまうと誤差の溜まり具合が増えてしまう感じがします。

ある配列で定義したグリッド上の値

計算結果を、ある配列で定義したグリッド上で欲しい場合、以下のプログラムで行うことが出来ます。
下に載せたプログラムと、上のコードの中にあるサブルーチン(grk, irkgl, dirk6, discrete_h, Jacobian)を一緒にコンパイルしてください。

終点だけの結果が欲しい場合

終点だけの結果が欲しい場合、不要なwork配列などは省略できます。
下に載せたプログラムと、上のコードの中にあるサブルーチン(grk, irkgl, dirk6, discrete_h, Jacobian)を一緒にコンパイルしてください。

収束判定の余地


上のプログラムは収束判定を[1]と同じにしているため、過剰に評価しているパラメータになっているかもしれません。

その余地としては、計算回数を減らすために重要な順に

  1. ニュートン法の収束判定 … Ntol
  2. 刻み幅の安全係数 … fac
  3. 刻み幅の離散化 … discrete_h内のimax
  4. ヤコビアンの更新条件 … if(Newt.le.2.or.theta.lt.1d-3)Jup=0の箇所

です。現状のプログラムでは安全のために過剰評価気味にしています。

4倍精度とGECP


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

また、倍精度、4倍精度に対応しているので、それらのプログラムを使って
陰的ルンゲクッタ法を書いたものが以下のものです。
irkgl_dge.f90
irkgl_qge.f90

倍精度 d
4倍精度 q

プログラムの評価



5つの問題について、評価した結果を上に載せます。評価の良し悪しは、
連立方程式の右辺の関数が何回評価されたかによって決めました。
比較対象は

  1. 陰的解法である本稿のプログラム(自作)
  2. 陽的解法であるルンゲ=クッタ=フェールベルグの4,5次(自作)
  3. 陰的解法であるRADAU IIAに基づくプログラム[3]

です。横軸に要求した精度、縦軸に実際に評価された回数を載せました。
まず、自作の陽的解法、陰的解法を比べますと、硬くない方程式であるEq.1,2,3では
陽的解法の方が10倍近く早いです。
しかし、硬い方程式であるEq.4,5は10倍から1000倍ほど陰的公式の方が早いという結果が得られました。望み通り、陰的解法が動いていることが確認できます。

…さて、専門家が書いたRADAU5ですが、めちゃくちゃ早いです。硬い方程式であるEq.4,5でも、自作したやつの1/10の回数で大体終わっています。
しかも、硬くない方程式であるEq.1,2,3ですら、自作の陽的解法よりも少ない回数で計算を終えています。本当にどういう事なんでしょうね…。
上には上がいるものです。

結論


ちゃんとした陰的解法が欲しいのであれば、自作せず、専門家のプログラムを使いましょう。

RADAU5の使い方


RADAU5のFORTRANプログラムはhttps://www.unige.ch/~hairer/software.htmlにあります。
使い方に手間取ったので、どうやって使うのかメモしておきます。

  1. http://www.unige.ch/~hairer/testset/testset.htmlに移動し、Van der Pol方程式のメインプログラムをダウンロードする。場所は
    ・VDPOL … driver for RADAU5,
    と書かれている所をクリックするとhttp://www.unige.ch/~hairer/testset/stiff/vdpol/driver_radau5.fに飛ぶ。これを保存し、driver_radau5.fという名前で保存する。
  2. https://www.unige.ch/~hairer/software.htmlに移動
  3. Stiff Differential Equations and Differential-Algebraic Problems
    という項目の
    ・RADAU5
    ・DC_DECSOL
    ・DECSOL
    のリンク先のプログラムをダウンロード。それぞれradau5.f, dc_decsol.f, decsol.fという名前で保存する。
  4. 合計4つのプログラムをダウンロードしたら、コンパイルを
    gfortran decsol.f dc_decsol.f radau5.f driver_radau5.f

    でコンパイルし、実行する。

参考文献


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

[2]E.~Hairer and G.~Wanner. Stiff diferential equations solved by Radau methods, J. Comput. Appl. Math., 111:93-111, 1999.

[3]E. Hairer, Fortran and Matlab Codes https://www.unige.ch/~hairer/software.html

[4]10. 常微分方程式 (2)https://www.ktech.biz/jp/archives/1003, K Technologies Sites (2014)

[5]Nicola Guglielmi and E.~Hairer, User’s Guide for code RADAU5 – Version 2.1 (packed in “radar5-v2.1.tar”) http://www.unige.ch/~hairer/software.html, 2005

陰的ルンゲ=クッタ法

もし利用されるのであれば、実際のプログラムは
陰的ルンゲ=クッタ法の高速化
にして下さい。このページのプログラムは上のページのプログラムに比べて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

ルンゲ=クッタ法の説明と刻み幅制御

ルンゲ=クッタ法(Runge-Kutta method、RK法)とは?
僕の知る限りの知識で紹介します。

特に良く使われる陽的ルンゲ=クッタ法は、
・実装が簡単
・良いアルゴリズムではない
という手法です。

良いアルゴリズムである陰的ルンゲ=クッタ法は、
陰的ルンゲ=クッタ法
をご覧ください。

もくじ

  1. ルンゲ=クッタ法の系統
  2. 陽的ルンゲ=クッタ法の段数と次数について
  3. 誤差について
  4. Butcher tableによるルンゲ=クッタ法の記述
  5. 埋め込まれた陽的ルンゲ=クッタ法
  6. ルンゲ=クッタ=フェールベルグ法による刻み幅の自動制御(適応刻み幅制御)
  7. ルンゲ=クッタ=フェールベルグ法による刻み幅の自動制御のプログラム(実数/複素数 で 1/2 階微分方程式を解くプログラム)
  8. 不連続な点を含む場合
  9. 刻み幅制御のベンチマーク(振り子)
  10. (追記)ルンゲ=クッタ=ドルマンド=プリンス法
  11. 陽的ルンゲ=クッタ法の導出
  12. 参考文献

理論はいいから4次ルンゲ=クッタ法の数値計算ではどうやるかだけ教えてくれ!という方は4次ルンゲ=クッタ法へどうぞ。

ルンゲ=クッタ法の系統


ルンゲ=クッタ法は微分方程式の数値計算解を得るための手法のことを指します。

通常の参考書で微分方程式を解くために良く紹介されているのは、オイラー法と中点法、4次ルンゲ=クッタ法でしょうか。
オイラー法も中点法も4次ルンゲ=クッタ法も、”陽的ルンゲクッタ法” と呼ばれる枠組みの1つです。

オイラー法は正確には “陽的1段1次ルンゲ=クッタ法” と呼ばれ、
中点法は “陽的2段2次ルンゲ=クッタ法”、
4次ルンゲクッタ法(RK4)は、”陽的4段4次ルンゲ=クッタ法” と呼ばれています。

“段”と”次”とはなんなんでしょう?それは、
計算の大変さ(段)と、計算の正確さ(次)
です。”“の値が小さければ小さいほど計算時間が少なくて済みますし、”“の値が高ければ高いほど計算が正確です。

オイラー法は1という計算コストで正確さ1が得られますし、RK4は4という計算コストで正確さ4が得られます。

4次ルンゲ=クッタ法が使われる理由

理由は実装が簡単でそれなりの精度を持つから。です。

陽的ルンゲ=クッタ法において、pという計算コスト(p段)で、pより大きな正確さq(q次)を得ることはできません
Derivation of Runge–Kutta methodsによれば、
\(q\)次の正確さ(q次のオーダー)を得たい場合、最低限必要な段数\(p_{\mbox{min}}(q)\)は

という関係にあります。

ここで注目するべきは4次の時までは計算コストに比例して計算精度が上がっていきます。
しかし、4次以上では、計算コストの増加と計算の正確さが見合わなくなっていきます。故に
計算効率が良いのは4次だろう、
と予想できます。
また、重要な理由として、4次ルンゲ=クッタ法に現れる係数が5次以降と比べて圧倒的にシンプルであることが挙げられます。4次では\(0,1/2,1/6\)程度の係数だけが使われ、プログラム作成時の入力ミスがほぼ生まれません。しかし、5次では\(28561/56430, -7200/2197\)といった係数が数多く出てきます。
これらの理由から,4次ルンゲ=クッタ法(RK4)が数値計算科学の世界でよく使われるのです。

陽的ルンゲ=クッタ法に限って言えばプログラムの実装が非常に簡単であることが挙げられます。陰的ルンゲ=クッタ法と呼ばれるアルゴリズムもあり、これは陽的ルンゲ=クッタよりも優れていますが、計算量が多くなり、若干複雑なアルゴリズムになります。陰的ルンゲ=クッタ法を詳しく知りたい方は陰的ルンゲ=クッタ法をご覧ください。

Q. オイラー法もものすごく細かい分点を取れば、その計算効率はRK4と同じなんじゃないの?
A. 刻み幅の乗数で効いてくるのでそうではありません。高次を使っても計算が信頼できるのであれば、大きなステップで進めるほうが早いです。例え、目標精度への計算時間が同じだとしても、計算機の有効桁数によって否定されてしまいます。
RK4で典型的にとられる時間ステップの間隔は、おおよそ\(10^{-2}\sim 10^{-4}\)程度であり、RK4のエラーのオーダーは\(O(h^5)\)です。
そして、科学計算で使う際の有効桁数は倍精度型で16桁です。
1ステップ当りの誤差は\(h\)の減少に伴い、解が\(h^4\)の早さで収束していく、と言えます。
だから16桁の計算では\(h=10^{-1}\to h=10^{-4}\)に変化させる時、誤差は\(O(h^5)=c 10^{-5}\to c 10^{-20}\)と変化します。
おおよそ\(c\approx 1\)と見積もれば、(有効桁数16桁を多少超えてしまいますが、)有効桁数いっぱいまで正しい値が出るであろうと期待できます。

これに対し、オイラー法で同じことをするには\(h\)を\(10^{-16}\)にしなくてはいけません。
\(t\)の値が\(10^{-16}\)変わった時に、桁落ちの問題を回避できるほど関数\(x\)の値に変化が生まれるか?
が問題になり、まぁそんな急激な変化は生まれないでしょう、と予想できます。これでは桁落ちの問題を回避するほどの変化は到底望めません。

よって計算の効率と有効桁数の限界から、RK4なのです。

また、あまりに高次の方法を使うとルンゲ現象に代表される不安定性といったことが起こるかもしれません。
高次は高精度という意味とイコールではないことに注意しましょう。この高次≠高精度については等間隔の分点における数値積分の時に書いたので気になる方はどうぞ。。

誤差について


4次ルンゲ=クッタ法の、1ステップ当りの誤差は\(h^5\)に比例,すなわち\(O(h^5)\)です。
しかし、通常は細かいルンゲ=クッタ法を何度も繰り返して計算します。
区間\([a,b]\)を刻み幅\(h\)の4次ルンゲ=クッタ法で\(N\)回のステップを繰り返し計算する場合、誤差は
\(
\displaystyle N\cdot O(h^5)=\frac{b-a}{h}\cdot O(h^5)=O(h^4)\)

となります。よって、\(N\)回繰り返すような計算では、オーダーが1つ落ちる事に注意しましょう。


スポンサーリンク


陽的ルンゲ=クッタ法の段数と次数について


さて、ここまで“段”は計算コスト、で“次”は計算の正確さ、という曖昧な表現でしたが、その表現をちゃんと知りましょう。
段と次を知るためにはルンゲ=クッタ法の計算方法を知る必要があります。
具体例を載せます。
\(
\displaystyle \frac{dx}{dt}=f(t,x)
\)

の、\(t_{n}\rightarrow t_{n}+h\ \ (=t_{n+1})\)における初期値問題に対する、
陽的1段1次ルンゲ=クッタ法(オイラー法)の計算スキームは、
\(
x_{n+1}=x_{n}+h\cdot f(t_{n},x_{n})
\)
です。

陽的4段4次ルンゲ=クッタ法(RK4)の計算スキームは、
\(
\begin{align}
k_1&=f(t_n, x_n) \\
k_2&=f(t_n+h/2, x_n+h k_1/2) \\
k_3&=f(t_n+h/2, x_n+h k_2/2) \\
k_4&=f(t_n+h, x_n+h k_3) \\
x_{n+1}&=x_{n}+{(k_1+2k_2+2k_3+k_4)}h/6
\end{align}
\)
として与えられます。

一般的に、陽的s段のルンゲ=クッタ法とは
\(
\begin{align}
g_i&=x_n+h\sum_{j=1}a_{i,j}k_j\ \ \ (j\lt i, \ i=1,2,…,s) \\
k_i&=f(t_n+c_ih, g_i) \\
x_{n+1}&=x_n+h\sum_{i=1}^s b_ik_i
\end{align}
\)
として書けます。
ここで行列形式で与えられる係数\(a_{i,j}, b_{i},c_{i}\)によって、そのs段ルンゲ=クッタ法が持つ次数が決められます。段数はここから由来します。

点\((t_n, x(t_n))\)周りで関数をテーラー展開し、その関数が点\((t_n+h\ \ (=t_{n+1}), x(t_{n+1}))\)で作る点を近似解とするのがルンゲ=クッタ法です。
故に、\(x(t_{n+1})\)は、
\(
\begin{align}
x(t_{n+1})=x(t_n)+\left.\frac{h}{1!}\frac{dx}{dt}\right|_{t=t_n}+\left.\frac{h^2}{2!}\frac{d^2x}{dt^2}\right|_{t=t_n}+\left.\frac{h^3}{3!}\frac{d^3x}{dt^3}\right|_{t=t_n}+\left.\frac{h^4}{4!}\frac{d^4x}{dt^4}\right|_{t=t_n}+…
\end{align}
\)
と書けます。
ここで、テイラー展開としてどの程度一致させて\(x(t_n+h)\)を決定するか?を表すのが次数に当たります。

言葉で書くなら、

局所的な精度は,高次項が解のテイラー(Taylor)展開とどの程度マッチしているかによって測定される.オイラー法は一次精度なので,誤差は一次高く\(h^2\)というベキ乗から始まる. ~NDSolveの”ExplicitRungeKutta”メソッドより

ともあります。

Butcher tableによるルンゲ=クッタ法の記述


行列形式で与えられるルンゲ=クッタ法での係数\(a_{i,j}, b_{i},c_{i}\)は何なのか?
具体的に記述してみましょう。
オイラー法(1段1次)はもっとも単純で、係数は
\(
\begin{align}
a_{1,1}&=0 \\
b_{1}&=1 \\
c_{1}&=0
\end{align}
\)
です。これを一般的な表記法の式に当てはめれば、
\(
\begin{align}
g_1&=x_n+h a_{1,1}k_1 \\
k_1&=f(t_n+c_1h, g_1) \\
x_{n+1}&=x_n+h b_1k_1
\end{align}
\)
となります。

中点法は、
\(
\begin{align}
a_{1,1}&=0 \\
a_{1,2}&=0 \\
a_{2,1}&=1/2 \\
a_{2,2}&=0 \\
b_{1}&=0 \\
b_{2}&=1 \\
c_{1}&=0 \\
c_{2}&=1/2
\end{align}
\)
という組で与えられます。

この係数行列の組はまとめてButcher tableと呼ばれる表記をするのが便利です。

これは、\(a_{i,j}, b_{i},c_{i}\)を
Buchertable
としてまとめて書く表記法です。

再び、オイラー法はButcher tableで書くと
euler_butcher
とまとめて書くことができます。
中点法は
midpoint_butcher
RK4は
rk4_butcher
です。

高次のルンゲ=クッタ法(10,12,14次)


4次、5次…とずっとあるわけです。
こんなページがありました[3]。
High-Order Explicit Runge-Kutta Methods
この上のページには

  • 17段10次(8次が埋め込まれてる)
  • 25段12次(10次が埋め込まれてる)
  • 35段14次(12次が埋め込まれてる)

といったButcher tableにおける係数の値が書かれています。埋め込まれてる、の意味は次の節で説明します。
ただし、上のページのbutcher tableは
highorderbutcher2
となっているので注意が必要です。

埋め込まれた陽的ルンゲ=クッタ法


埋め込まれた“という表現が出てきたのでその説明を行いましょう。
日本語では『埋め込まれた陽的ルンゲ=クッタ法』、英語では『embedded explicit runge-kutta method』と呼ばれるものがあります。
これは、p段q次陽的ルンゲ=クッタ法を作ったら、別の次数の陽的ルンゲ=クッタ法も、係数行列\(a_{i,j}, c_{i}\)を使って作れるじゃありませんか!
というものです。

Butcher tableは、この場合extended Butcher tableと呼ばれ、こういう形式で書かれます。
embedded_rk
この埋め込まれたルンゲ=クッタ法のいいところは、

  1. 計算誤差の評価ができる
  2. 刻み幅を自動的に制御できる、適応刻み幅制御。(応用として。)

という点です。ルンゲ=クッタ法によって得られた解が真の解とどのくらい違っているのか?が評価できるんです。

例えば、4次のルンゲ=クッタ法を使って得られた解\(x^{(4)}(t)\)と5次のルンゲ=クッタ法を使って得られた解\(x^{(5)}(t)\)があったとします。
もしも、\(x^{(4)}(t)\)と解\(x^{(5)}(t)\)の解の差を調べ、その差が無かったらその数値計算解は真の解に限りなく近い、と判断することができ、差が大きかったらその解は真の解から離れていて、数値計算の精度が足らない、と判断することができます。どちらも1つだけの解では出来なかったことです。

精度が足らない場合、刻み幅を小さくすれば精度が上がります。また逆に、精度が十分に足りている場合、刻み幅を大きくし、計算時間を減らすことができます。
これが適応刻み幅制御なのです。

違った次数のルンゲ=クッタ法を、まるで別々に計算してもいいのですが、そうすると計算時間が単純に考えておおよそ2倍になります。
解を評価するために2倍の計算時間が必要というのは良くない計算効率です。
そこで考えられたのが埋め込まれたルンゲ=クッタ法なのです。

具体例を挙げましょう。
一番簡単な埋め込まれたルンゲ=クッタ法は、ホイン法と呼ばれています。
butcher_heun
1行目は2次のオーダーを持ち、2行目は1次のオーダーを持ちます。

また、4次と5次を持つ埋め込まれたルンゲ=クッタ法は、ルンゲ=クッタ=フェールベルグ(Runge-Kutta-Fehlberg)法と呼ばれています。
その埋め込まれたルンゲ=クッタ法は、
butcher_rkf45
と書かれます。1行目は5次のオーダー、2行目は4次のオーダーを持ちます。

ルンゲ=クッタ=フェールベルグ法による刻み幅の自動制御(適応刻み幅制御)


さて、次数の違う2つのルンゲ=クッタ法を用いて、適応刻み幅制御を行いたいと考えます。
刻み幅を制御するにあたって、適当に精度良かったから2倍にしてもまだ大丈夫だろ、とか差が大きすぎるから刻み幅半分にしよう、ということをやってはいけません。
適当にやったら計算時間が余計にかかり、精度が良くない変な結果が得られます。

[5~9]によれば、ルンゲ=クッタ=フェールベルグ法において区間\(i\)での最適な刻み幅\(h’\)は区間\(i\)の誤差評価の結果を使って、
\(
\displaystyle h’=\left(\frac{\varepsilon h}{2|x^{(5)}_{i}-x^{(4)}_{i}|}\right)^{1/4}
\)

と予想できます。もちろん、この\(h’\)は区間\(i\)の最適な幅ですが、関数に劇的な変化は無いだろうとして、次の区間の計算の一番初めに用いる刻み幅を推定するのです。
なので、\(i+1\)番目の計算区間では、計算するときはこの\(h’\)の値を使えばいいんです。
(ちなみに、m次ルンゲ=クッタ法の場合では
\(
\displaystyle h’=\left(\frac{\varepsilon h}{2|x^{(m+1)}_{i}-x^{(m)}_{i}|}\right)^{1/m}
\)

と推測されます。)

詳しい理由は分かりませんが、5次オーダーではなく、4次です。5次のオーダーは誤差評価のためだけに用いられているようです。
ルンゲ=クッタ=フェールベルグ法の計算スキームは[7]に詳しく書かれています。
日本語訳して、その計算スキームを書けば下のようになります。
刻み幅制御rkf45_4order

ルンゲ=クッタ=フェールベルグ法による刻み幅の自動制御(プログラム)


講義のレポート等の宿題で使うのは僕の意向と異なるので使用はご控えください。
研究目的、趣味、確かめの場合はミスがあるかもしれないことを念頭に置いたうえならば使用と改変をご自由にしてください。
このプログラムを使用して生じた責任は取りません。

fortran90によるプログラムです。ほぼ上の説明をそのままプログラミングしたものです。

  • 実数、1階微分方程式の場合
  • 微分方程式
    \(
    \displaystyle \frac{d}{dx}y=y\cos x,~~y(x=0)=1
    \)

    を倍精度実数、刻み幅制御で\(x=10\)まで解く事を考えます。
    解析解は
    \(
    \displaystyle y(x)=exp(sin(x))
    \)

    です。コードは以下の通りです。


  • 実数、2階微分方程式の場合
  • 微分方程式
    \(
    \displaystyle \frac{d^2}{dx^2}y=-\frac{1}{4}y,~~y(x=0)=1,~ y'(x=0)=0
    \)

    を倍精度実数、刻み幅制御で\(x=20\)まで解く事を考えます。コードは以下の通りです。

    ※ここで使われているサブルーチンdrkf45は実数一階微分のプログラム内にあるルーチンと一字一句同一です。


  • 複素数、1階微分方程式の場合
  • 微分方程式
    \(
    \displaystyle \frac{d}{dx}y=y\cos x,~~ y(x=0)=1+i\frac{1}{2}
    \)

    を倍精度実数、刻み幅制御で\(x=10\)まで解く事を考えます。コードは以下の通りです。


  • 複素数、2階微分方程式
  • 微分方程式
    \(
    \displaystyle \frac{d^2}{dx^2}y=-\frac{1}{4}y,~~ y(x=0)=1+i\frac{1}{2},~y'(x=0)=0+i0
    \)

    を倍精度実数、刻み幅制御で\(x=20\)まで解く事を考えます。コードは以下の通りです。

    ※ここで使われているサブルーチンcrkf45は複素数一階微分のプログラム内にあるルーチンと一字一句同一です。

——————————–

等間隔の出力の場合は、以下の通りで実行できます。
サブルーチンはdrkf45は変わっていません。

スポンサーリンク

不連続な点を含む場合


不連続な点を含む場合、境界条件を指定しないと解くことはできません。

さて、ここで微分方程式
\(
\begin{eqnarray}
\frac{dy}{dx}=
\left\{
\begin{aligned}
0\;\;(x\le 0)\\
1\;\;(0\lt x)
\end{aligned}
\right.
\end{eqnarray}
\)

を初期条件\(y(-0.095)=0\)の下で考えます(意図的に境界条件は考えません)。
これを4次ルンゲ=クッタと適応刻み幅ルンゲ=クッタで解いてみましょう。
プログラム上ではそのまま解くことが出来ます。
実際に解かせてみますと、
rk4_rkf451_disc
となり、2つの結果(赤:4次ルンゲ=クッタ、緑:刻み幅制御ルンゲ=クッタ)は異なってしまいます。これは1階微分の不連続性のため発生します。
不連続点\(x=0\)で関数\(y(x)\)に境界条件を指定しない限り、どちらも正しい解なのです。

さて、なぜこんなことが発生するのでしょうか?以下のように問題を表すことにします。

不連続な点を含む1階の微分方程式を考えます。
ここで不連続、という意味は関数\(y(x)\)の一階微分が、点\(x’\)で
\(
\displaystyle \left. \frac{dy}{dx}\right|^{x=x’+0}_{x=x’-0}=a,\;\;\;(a\ne 0)
\)

であるような点を指しているとします。

上記の例題を考えてみましょう。
上記の例題では\(a=1\)です。微分方程式を解析的に解いてみますと、
\(
\begin{eqnarray}
y(x)=
\left\{
\begin{aligned}
C_0\;\;(x\le 0)\\
x+C_1\;\;(0\lt x)
\end{aligned}
\right.
\end{eqnarray}
\)

となります。ここで\(C_0, C_1\)は定数です。
\(C_0, C_1\)は\(y(x)\)が解きたい問題の境界条件によって決まります。

例えば、\(y(x)\)は全領域に対して繋がっている、という条件を課しましょう。この場合、不連続点\(x’\)で
\(
\displaystyle \left. y(x)\right|^{x=x’+0}_{x=x’-0}=0
\)

という境界条件を満たさなければなりません。この条件を課すと、\(C_1=C_0\)となり、初めて関数\(y(x)\)を一意に決めることが出来ます。

1階微分方程式を解く場合、適応刻み幅制御では関数\(y(x)\)は計算領域内で繋がっている事が課されています。しかし、4次ルンゲ=クッタではその条件は課されません。\(C_1\)の値は初期条件に依存し、一意に関数が決まりません。
どちらが悪いという話ではありません。

通常は適応刻み幅でも、4次ルンゲ=クッタでも\(y(x)\)にどこか連続ではない変な点がある場合、その点で区間を別々に分けて解きます。その後、境界条件に従って値を調節して全体の関数を構成します。

ベンチマーク用


微分方程式の解法がどれくらい正しそうかのベンチマーク問題として振り子(角度が大きい時)を考えましょう(振り子の詳しい解説はこちら)。
以下の\(\omega=1\)としたときの運動方程式
\(
\displaystyle \frac{d^2\theta}{dt^2}=-\sin\theta
\)


\(
t=0,\; \theta=0; \frac{d\theta}{dt}=1.9\cdots (k=0.95)
\)
の初期条件の下解いた場合、
1周期\(T\)は
\(
T=4K(0.95)=10.360044923498004876778\cdots
\)
となります。
この値はwolfram alphaから求めました。
4*EllipticK[0.95*0.95] wolfram alpha

刻み幅制御を行い、45000周期目の値を考えます。45000周期目は時刻
\(
T_{45000}=466202.0215574102\cdots
\)
です。刻み幅制御による精度を\(10^{-12}\)に設定し、数値計算を行わせます。

すると実行結果として”fort.10″に

0.4662020113E+006  -0.2103356901E-001   0.1899883579E+001   0.4224109363E-002
0.4662020155E+006  -0.1300808922E-001   0.1899955473E+001   0.4223843994E-002
0.4662020198E+006  -0.4982881533E-002   0.1899993468E+001   0.4223658015E-002
0.4662020240E+006   0.3042061693E-002   0.1899997567E+001   0.4223520921E-002

というデータが出力されます。
1列目が時刻\(t\)、2列目が\(\theta(t)\),3列目が\(\frac{d\theta(t)}{dt}\),4列目が刻み幅\(h\)です。
1回のステップでの要求精度12桁に対し、最終的な結果は8桁まで正しい値を出しています。
この時、計算回数は77,852,488回\(~10^{8}\)回行われているので、最終的な結果として4桁は少なくとも正しいと考えられます。
4桁合っていればいい状況で8桁もあっているのは、被積分関数が変な振る舞いをしないからでしょう。

また、60000周期で確認してみると(\(T_{60000}=621602.695409880292\cdots\))

0.6216026888E+006  -0.1531918479E-001   0.1899938246E+001   0.4223959920E-002
0.6216026930E+006  -0.7293808996E-002   0.1899986003E+001   0.4223717084E-002
0.6216026973E+006   0.7312355417E-003   0.1899999862E+001   0.4223582630E-002
0.6216027015E+006   0.8756011575E-002   0.1899979827E+001   0.4223462029E-002

です。
1回のステップでの要求精度12桁に対し、最終的な結果は8桁まで正しい値を出しています。
この時、計算回数は103,803,513回\(~10^{8}\)回行われているので、最終的な結果として4桁は少なくとも正しいと考えられます。
・・・まだまだ大丈夫そうですね。

少し特殊な初期条件(\(k=1\))でやってみましょう。
この\(k=1\)は、振り子の質点がちょうど真上に来て静止する非常に不安定な状態です。
何秒間静止していられるか試してみましょう。刻み幅の制御等は上記条件と同じです。
横軸に時間\(t\)、縦軸に\(\theta(t)\)を取った時のグラフです。
ellk1
すぐに破綻しました。正しい値は\(t=20\)位までですかね。これは、機械的な誤差があることによって不安定な平衡点からほんの少しだけ外れてしまったからです。だからカオスとかを考えるときとかは注意しなければなりません。

ルンゲ=クッタ=ドルマンド=プリンス法


フェールベルグ法は初期の頃に考えられた埋め込まれた方法です。
研究が進み、実用上では別の係数の組み合わせが良いことが分かってきました。
その一つが1980年に発見されたドルマンド=プリンス(Dormand-Prince)法です。

これは、7段4(5)次の方法です。
フェールベルグ法は6段4(5)次の方法ですので、次数は同じです。

良く調べていませんが、この違いは、4次の結果を基準にして求めたブッチャー係数(フェールベルグ法)か、5次の結果を基準に求めたブッチャー係数(ドルマンド=プリンス法)か?の違いのようです。

単純に考えて、同じ次数なのにドルマンド=プリンス法の方が段数が増えていて効率が悪いです。
しかし、本来は7段なのですが、7段目に呼び出した結果を取って置けば、次のステップの1段目に同じ値が使えるように設計されているので、プログラム上は6段と(ほぼ)同じ関数の呼び出し回数になります。

プログラムはこんな感じになるでしょう。

適当な刻み幅で出力

等間隔(サブルーチンは上のものと同じなので省略)

program main
  implicit none
  integer::i,N,info,Nx
  double precision::h,tol
  double precision::xa,xb,tx
  double precision,allocatable::y(:),x(:),work0(:)
  double precision,external::grk

  N=1
  allocate(y(1:N),work0(1:N))

  Nx=101
  allocate(x(1:Nx))
  xa=0d0
  xb=10d0
  do i=1,Nx
     x(i)=(i-1)*(xb-xa)/dble(Nx-1)+xa
  enddo

  !initial conditions
  y(1)=1d0  ! x (0)

  tol=1d-8
  write(10,'(2e25.10e3)')x(1),y(1)
  do i=2,Nx
     info=-1
     h=x(i)-x(i-1)
     tx=x(i-1)
     do while(info.le.0)
        call dDP45(grk,tx,h,N,y,x(i),info,tol,work0)
     enddo
     write(10,'(2e25.10e3)')x(i),y(1)
  enddo
 
  stop
end program main

function grk(N,x,y,s)
  implicit none
  integer,intent(in)::N,s
  double precision,intent(in)::x
  double precision,intent(in)::y(1:N)
  double precision::grk
 
  grk=0d0
  if(s.eq.1)then
     grk=y(1)*cos(x)
  else
     write(6,*)"***Error grk"; stop
  endif
 
  return
end function grk

4倍精度ルーチン


4倍精度のサブルーチンです。
計算速度は倍精度の30~50倍かかるので、必要なとき以外使わないようにしましょう。

陽的ルンゲ=クッタ法の導出


ルンゲ=クッタ法の導出は煩雑です。単に複雑なだけです。
導出過程について詳しく述べられているページは、早川様が記述された以下のpdfを参照すると良いと思います。
Runge-Kutta法についてのノート(早川尚男)
計算過程を含め記述されているので分かりやすいです。

参考文献


[1]Derivation of Runge–Kutta methods
[2]NDSolveの”ExplicitRungeKutta”メソッド
[3]High-Order Explicit Runge-Kutta Methods
[4]List of Runge–Kutta methods
[5]Runge-Kutta-Fehlberg Method (RKF45)
[6]Runge-Kutta-Fehlberg method
[7]Lecture:13Runge-Kutta-Fehlberg Method
[8]GPU acceleration of Runge Kutta-Fehlberg and its comparison with Dormand-Prince method
[9]William H. Pressら著『ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ– 』(1993)


↑この本が一番有益だと思います。Fortran版もありますが、英語しかありません。ちなみに、英語で,若干古いバージョンでいいのならば
公式のホームページ
Numerical Recipes in C (1992)
Numerical Recipes in Fortran 77 and Fortran 90 (1992, 1996)
で無料で公開されています。

4次ルンゲ・クッタ法

微分方程式なら任せろーバリバリバリー

ルンゲ=クッタ法ってもともとどういうもの?理論は?刻み幅\(h\)を自動的に制御する方法について知りたい方は、ルンゲ=クッタ法の系統的扱いと刻み幅制御へどうぞ。

4次ルンゲ=クッタ法は微分方程式を数値的に解く手段です。

ルンゲ=クッタ法が良く使われる理由は、ひとえにプログラムの実装のしやすさです。ルンゲクッタ法は非常に良い!というアルゴリズムではありませんが、他の方法よりもシンプルで、プログラムに組み込みやすいのです。

”4次”は問題の解の関数をテイラー展開した場合、4次までは一致するように作られた、という意味です。
例えば微分方程式
\(
\displaystyle \frac{dy}{dx}=g(x,y)
\)

を考えます。
数値計算では初期値(例えばx=0の時、y=1など)を決めて、そこからxをhだけ増やし、微分方程式というルールに従って関数\(y(0+h),y(0+h+h),y(0+h+h+h),\cdots\)を作り上げていきます。
この時、4次ルンゲ=クッタ法で求められる答えの関数というのは
\(
\displaystyle y(x+h)=y(x)+hg(x,y)+\frac{h^2}{2!}\frac{dg}{dx}+\frac{h^3}{3!}\frac{d^2g}{dx^2}+\frac{h^4}{4!}\frac{d^3g}{dx^3}+O(h^5)
\)

という関数になります。
\(
g(x,y)=-xy
\)

(ただし、初期条件\(x=0\)で\(y=1\))
である場合、微分方程式の解析解は
\(
\displaystyle y=e^{-x^2/2}
\)

であるため、4次ルンゲ=クッタ法によって導かれる答えは、
\(
\displaystyle y(x+h)=y(x)\left[1-hx+\frac{h^2}{2}(x^2-1)+\frac{h^3}{6}(-x^3+3x)+\frac{h^4}{24}(x^4-6x^2+3)\right]+O(h^5)
\)

となります。

本題に入りましょう。
4次ルンゲ=クッタ法は6つのステップが必要となります。
初期値を\((x_0,y_0)\)と書くと、(\(y_0=y(x_0)\)です。)

  1. \((x_0,y_0)\)より\(k_a\)を求める
  2. \(k_a\)より\(k_b\)を求める
  3. \(k_b\)より\(k_c\)を求める
  4. \(k_c\)より\(k_d\)を求める
  5. \(k_a,k_b,k_c,k_d\)より\(y(x_0+h)\)を求める
  6. \((x_0+h,y(x_0+h))\)を初期値だと思って手順1に戻る。

という感じです。
3章:連立ルンゲ・クッタ法による微分方程式の解を参考にすると、数値計算での4次ルンゲ=クッタ法の計算スキームは以下のようになります。

解きたい微分方程式を連立1次微分方程式の形で書くと一般的にはこう書けます。
———
\(
\begin{eqnarray}
\left\{
\begin{aligned}
\frac{dy_1}{dx} &= f_1(x,y_1,y_2,\cdots,y_N) \\
\frac{dy_2}{dx} &= f_2(x,y_1,y_2,\cdots,y_N) \\
&\vdots \\
\frac{dy_N}{dx} &= f_N(x,y_1,y_2,\cdots,y_N) \\
\end{aligned}
\right.
\end{eqnarray}
\)
———
まず係数\(k_a\)を求めます。
———
\(
\begin{eqnarray}
\left\{
\begin{aligned}
k_{a1} &= hf_1(x,y_1,y_2,\cdots,y_N) \\
k_{a2} &= hf_2(x,y_1,y_2,\cdots,y_N) \\
&\vdots \\
k_{aN} &= hf_N(x,y_1,y_2,\cdots,y_N) \\
\end{aligned}
\right.
\end{eqnarray}
\)
———
次に\(k_b\)。
———
\(
\begin{eqnarray}
\left\{
\begin{aligned}
k_{b1} &= hf_1(x+\frac{h}{2},y_1+\frac{k_{a1}}{2},y_2+\frac{k_{a2}}{2},\cdots,y_N+\frac{k_{aN}}{2}) \\
k_{b2} &= hf_2(x+\frac{h}{2},y_1+\frac{k_{a1}}{2},y_2+\frac{k_{a2}}{2},\cdots,y_N+\frac{k_{aN}}{2}) \\
&\vdots \\
k_{bN} &= hf_N(x+\frac{h}{2},y_1+\frac{k_{a1}}{2},y_2+\frac{k_{a2}}{2},\cdots,y_N+\frac{k_{aN}}{2}) \\
\end{aligned}
\right.
\end{eqnarray}
\)
———
そして\(k_c\)。
———
\(
\begin{eqnarray}
\left\{
\begin{aligned}
k_{c1} &= hf_1(x+\frac{h}{2},y_1+\frac{k_{b1}}{2},y_2+\frac{k_{b2}}{2},\cdots,y_N+\frac{k_{bN}}{2}) \\
k_{c2} &= hf_2(x+\frac{h}{2},y_1+\frac{k_{b1}}{2},y_2+\frac{k_{b2}}{2},\cdots,y_N+\frac{k_{bN}}{2}) \\
&\vdots \\
k_{cN} &= hf_N(x+\frac{h}{2},y_1+\frac{k_{b1}}{2},y_2+\frac{k_{b2}}{2},\cdots,y_N+\frac{k_{bN}}{2}) \\
\end{aligned}
\right.
\end{eqnarray}
\)
———
最後に\(k_d\)。
———
\(
\begin{eqnarray}
\left\{
\begin{aligned}
k_{d1} &= hf_1(x+h,y_1+k_{c1},y_2+k_{c2},\cdots,y_N+k_{cN}) \\
k_{d2} &= hf_2(x+h,y_1+k_{c1},y_2+k_{c2},\cdots,y_N+k_{cN}) \\
&\vdots \\
k_{dN} &= hf_N(x+h,y_1+k_{c1},y_2+k_{c2},\cdots,y_N+k_{cN}) \\
\end{aligned}
\right.
\end{eqnarray}
\)
———
最後に求めた\(k_a,k_b,k_c,k_d\)を使って\(x+h\)でのそれぞれの関数の値を導きます。、
———
\(
\begin{eqnarray}
\left\{
\begin{aligned}
x &= x+h \\
y_{1} &= y_{1}+(k_{a1}+2k_{b1}+2k_{c1}+k_{d1})/6 \\
y_{2} &= y_{2}+(k_{a2}+2k_{b2}+2k_{c2}+k_{d2})/6 \\
&\vdots \\
y_{N} &= y_{N}+(k_{aN}+2k_{bN}+2k_{cN}+k_{dN})/6 \\
\end{aligned}
\right.
\end{eqnarray}
\)
———
となります。

4次ルンゲ=クッタ法のプログラム


実際に例題を解きましょう。
2つの例題を解きます。


スポンサーリンク


1階微分方程式


1つ目は
\(
\displaystyle \frac{dy}{dx}=-y\sin{x}
\)

です。
一般解は
\(
y=Ae^{\cos{x}}
\)

であり、初期条件\(x=0,y=2\)として解けば解析解は
\(
y=2e^{\cos{x}-1}
\)

です。これを4次ルンゲ=クッタ法を用いて解くには、以下のプログラムで実現できます。

解いて、gnuplotで
plot “fort.10”
数値計算解(赤点)解析解(緑線)と共に出力すると、こんなグラフが得られます。
4次ルンゲクッタ法による計算1

連立1階微分方程式


もう一つの例題は微分方程式
\(
\displaystyle \frac{d^2y}{dx^2}+2\gamma \frac{dy}{dx}+y=0
\)

を解きます。これは物理の世界ではバネの減衰振動を表す運動方程式で(詳しくは減衰振動へ。)、\(0\lt\gamma\lt1\)の場合、解は
\(
y(x)=Ae^{-\gamma x}\cos(x\sqrt{1-\gamma^2}-\alpha)
\)

となります。
\(\gamma=0.15\), 初期条件を今、
\(
\begin{eqnarray}
\left\{
\begin{aligned}
\left.y\right|_{x=0} &= 1 \\
\left.\frac{dy}{dx}\right|_{x=0} &= -0.15 \\
\end{aligned}
\right.
\end{eqnarray}
\)
とすると、解析解は
\(
y(x)=e^{-\gamma x}\cos{(x\sqrt{1-0.15^2})}
\)

となります。
プログラムで計算する際は、まず連立1次微分方程式に焼きなおす必要があります。すなわち、
\(
\begin{eqnarray}
\left\{
\begin{aligned}
\frac{dy}{dx}&=v \\
\frac{dv}{dx}&=-2\gamma v -y \\
\end{aligned}
\right.
\end{eqnarray}
\)
として解けばいいわけです。

その時の変更するべき場所は少しで、mainとrkfdの一部を変えればおしまいです。
変えた場所をコメントで
!+—+
を入れたので確認してください。

gnuplotで数値計算解(赤点)解析解(緑線)のグラフを書くとこうでなります。
4次ルンゲクッタ法による計算2

スポンサーリンク