月別アーカイブ: 2019年1月

プログラミングと数値計算, ルンゲ=クッタ法

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

コメントする

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
４倍精度 q

プログラムの評価


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

  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