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

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

微分方程式
\(
\displaystyle \frac{d^2}{dt^2}x=t
\)

を倍精度実数、刻み幅制御で解く場合のコードはこれです。複素数の場合は本文参照。

もくじ

  1. ルンゲクッタ法の系統
  2. 陽的ルンゲクッタ法の段数と次数について
  3. 誤差について
  4. Butcher tableによるルンゲクッタ法の記述
  5. 埋め込まれた陽的ルンゲクッタ法
  6. ルンゲ-クッタ-フェールベルグ法による刻み幅の自動制御(適応刻み幅制御)
  7. ルンゲ-クッタ-フェールベルグ法による刻み幅の自動制御のプログラム(実数/複素数 で変数を 自動/指定 で解く 1/2 階微分方程式のプログラム)
  8. 4倍精度ルーチン
  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によるプログラムです。ほぼ上の説明をそのままプログラミングしたものです。

実例を8つ載せます。

  • 実数で、時間の刻みを自動で決める1階微分方程式のプログラム
  • 実数で、時間の刻みを指定する1階微分方程式のプログラム
  • 複素数で、時間の刻みを自動で決める1階微分方程式のプログラム
  • 複素数で、時間の刻みを指定する1階微分方程式のプログラム
  • 実数で、時間の刻みを自動で決める2階微分方程式のプログラム
  • 実数で、時間の刻みを指定する2階微分方程式のプログラム
  • 複素数で、時間の刻みを自動で決める2階微分方程式のプログラム
  • 複素数で、時間の刻みを指定する2階微分方程式のプログラム

です。
時間の指定する、というのは、指定した時刻の間は自動で刻み幅を制御しながらやりますが、出力するのは指定した時間のみ、という意味です。

下のモジュールRKmodと共にコンパイルし、実行してください。

  • 実数で、時間の刻みを自動で決める1階微分方程式のプログラム

    下のプログラムは一階微分方程式
    \(
    \displaystyle \frac{dy}{dx}=y\cos(x)
    \)
    を初期条件
    \(
    \left. y(x)\right|_{x=0}=1
    \)
    の条件下で解くものです。
    解析解は\(y(x)=\exp(\sin(x))\)となります。

    上のプログラムを実行すると、
    “fort.10″に
    位置, 関数\(y(x)\), 刻み幅, 積分\(\int_{x_0}^{x}y'(x)dx\)
    の順にファイルが出力されます。
    gnuplotで関数を出力すると、
    expsin1
    が出力されます。紫色の線は解析解\(\exp(\sin(x))\)を表しています。関数の変化が激しいところで刻み幅が自動的に小さくなっているのが分かるでしょう。

  • 実数で、時間の刻みを指定する1階微分方程式のプログラム

    下のプログラムは一階微分方程式
    \(
    \displaystyle \frac{dy}{dx}=y\cos(x)
    \)
    を初期条件
    \(
    \left. y(x)\right|_{x=0}=1
    \)
    の条件下で解くものです。解析解は\(y(x)=\exp(\sin(x))\)となります。
    ただし、出力の間隔を指定します。

    上のプログラムを実行すると、
    “fort.11″に
    位置, 関数\(y(x)\), 積分\(\int_{x_0}^{x}y'(x)dx\)
    の順にファイルが出力されます。
    gnuplotで関数を出力すると、
    expsin2
    が出力されます。紫色の線は解析解\(\exp(\sin(x))\)を表しています。
    ここで、時間の間隔は等間隔に荒く指定しました。ですが、点と点の間は自動で刻み幅制御が行われているので、精度は指定したものとなっています。等間隔でなくても構いません。

  • 複素数で、時間の刻みを自動で決める1階微分方程式のプログラム

    下のプログラムは一階微分方程式
    \(
    \displaystyle \frac{dy}{dx}=\left(\frac{1}{2x}+i\right)y
    \)
    を初期条件
    \(
    \left. y(x)\right|_{x=1}=\exp(i)
    \)
    の条件下で解くものです。解析解は\(y(x)=\exp(ix)\sqrt{x}\)となります。

    上のプログラムを実行すると、
    “fort.10″に
    位置, 関数\(Re\{y(x)\}\), 関数の虚部\(Im\{y(x)\}\), 刻み幅、積分の実部\(Re\{\int_{x_0}^{x}y'(x)dx\}\), 積分の虚部\(Im\{\int_{x_0}^{x}y'(x)dx\}\)
    の順にファイルが出力されます。
    gnuplotで関数を出力すると、
    expsqrt1
    が出力されます。水色、オレンジの点は数値計算で解いた解の実部、虚部を表し、紫、緑の実線は解析解の実部、虚部を表します。

  • 複素数で、時間の刻みを指定する1階微分方程式のプログラム

    下のプログラムは一階微分方程式
    \(
    \displaystyle \frac{dy}{dx}=\left(\frac{1}{2x}+i\right)y
    \)
    を初期条件
    \(
    \left. y(x)\right|_{x=1}=\exp(i)
    \)
    の条件下で解くものです。解析解は\(y(x)=\exp(ix)\sqrt{x}\)となります。

    上のプログラムを実行すると、
    “fort.10″に
    位置, 関数\(Re\{y(x)\}\), 関数の虚部\(Im\{y(x)\}\), 積分の実部\(Re\{\int_{x_0}^{x}y'(x)dx\}\), 積分の虚部\(Im\{\int_{x_0}^{x}y'(x)dx\}\)
    の順にファイルが出力されます。
    gnuplotで関数を出力すると、
    expsqrt2
    が出力されます。水色、オレンジの点は数値計算で解いた解の実部、虚部を表し、紫、緑の実線は解析解の実部、虚部を表します。

  • 実数で、時間の刻みを自動に決める2階微分方程式のプログラム

    下のプログラムは二階微分方程式
    \(
    \displaystyle \frac{d^2y}{dx^2}= \frac{dy}{dx}\cos{x}-y\sin(x)
    \)
    を初期条件
    \(
    \begin{align}
    \left. y(x)\right|_{x=0}=1 \\
    \left. \frac{dy(x)}{dx}\right|_{x=0}=1
    \end{align}
    \)
    の条件下で解くものです。解析解は\(y(x)=\exp(\sin{x}),\frac{dy(x)}{dx}=\cos{x}\exp(\sin{x})\)となります。

    上のプログラムを実行すると、
    “fort.10″に
    位置, 関数\(y(x)\), 導関数\(\frac{dy(x)}{dx}\)
    の順にファイルが出力されます。
    gnuplotで関数を出力すると、
    expsinx_d2
    が出力されます。水色、オレンジの点は数値計算で解いた,解\(y(x)\)、導関数\(\frac{dy(x)}{dx}\)を表し、紫、緑の実線は解析解、とその導関数を表します。

  • 実数で、時間の刻みを指定する2階微分方程式のプログラム

    下のプログラムは二階微分方程式
    \(
    \displaystyle \frac{d^2y}{dx^2}= \frac{dy}{dx}\cos{x}-y\sin(x)
    \)
    を初期条件
    \(
    \begin{align}
    \left. y(x)\right|_{x=0}=1 \\
    \left. \frac{dy(x)}{dx}\right|_{x=0}=1
    \end{align}
    \)
    の条件下で解くものです。解析解は\(y(x)=\exp(\sin{x}),\frac{dy(x)}{dx}=\cos{x}\exp(\sin{x})\)となります。

    上のプログラムを実行すると、
    “fort.10″に
    位置, 関数\(y(x)\), 導関数\(\frac{dy(x)}{dx}\)
    の順にファイルが出力されます。
    gnuplotで関数を出力すると、
    expsinx_d2_d
    が出力されます。水色、オレンジの点は数値計算で解いた,解\(y(x)\)、導関数\(\frac{dy(x)}{dx}\)を表し、紫、緑の実線は解析解、とその導関数を表します。

  • 複素数で、時間の刻みを自動で決める2階微分方程式のプログラム

    下のプログラムは二階微分方程式
    \(
    \displaystyle \frac{d^2y}{dx^2}=\left(-\frac{1}{4x^2}+\frac{i}{x}-1\right)y
    \)
    を初期条件
    \(
    \begin{align}
    \left. y(x)\right|_{x=1}=\exp(i) \\
    \left. \frac{dy(x)}{dx}\right|_{x=1}=\left(\frac{1}{2}+i\right)\exp(i)
    \end{align}
    \)
    の条件下で解くものです。解析解は\(y(x)=\exp(ix)\sqrt{x},\frac{dy(x)}{dx}=\left(\frac{1}{2x}+i\right)\exp(ix)\sqrt{x} \)となります。

    上のプログラムを実行すると、
    “fort.10″に
    位置, 関数\(Re\{y(x)\}\), 関数の虚部\(Im\{y(x)\}\), 微分の実部\(Re\{\frac{dy(x)}{dx}\}\), 微分の虚部\(Im\{\frac{dy(x)}{dx}\}\), 刻み幅
    の順にファイルが出力されます。
    gnuplotで関数を出力すると、
    sqrtxexpix_d2
    が出力されます。水色、オレンジの点は数値計算で解いた,解\(y(x)\)の実部、解\(Re\{\frac{dy(x)}{dx}\}\)の実部を表し、紫、緑の実線は解析解の実部、とその導関数の実部を表します。

  • 複素数で、時間の刻みを指定する2階微分方程式のプログラム

    下のプログラムは2階微分方程式
    \(
    \displaystyle \frac{d^2y}{dx^2}=\left(-\frac{1}{4x^2}+\frac{i}{x}-1\right)y
    \)
    を初期条件
    \(
    \begin{align}
    \left. y(x)\right|_{x=1}=\exp(i) \\
    \left. \frac{dy(x)}{dx}\right|_{x=1}=\left(\frac{1}{2}+i\right)\exp(i)
    \end{align}
    \)
    の条件下で解くものです。解析解は\(y(x)=\exp(ix)\sqrt{x},\frac{dy(x)}{dx}=\left(\frac{1}{2x}+i\right)\exp(ix)\sqrt{x} \)となります。

    上のプログラムを実行すると、
    “fort.10″に
    位置, 関数\(Re\{y(x)\}\), 関数の虚部\(Im\{y(x)\}\), 微分の実部\(Re\{\frac{dy(x)}{dx}\}\), 積分の虚部\(Im\{\frac{dy(x)}{dx}\}\)
    の順にファイルが出力されます。
    gnuplotで関数を出力すると、
    sqrtxexpix_d2e
    が出力されます。水色、オレンジの点は数値計算で解いた,解\(y(x)\)の実部、解\(Re\{\frac{dy(x)}{dx}\}\)の実部を表し、紫、緑の実線は解析解の実部、とその導関数の実部を表します。

4倍精度ルーチン

置いておきます。
メインプログラム

program main
  use qRKmod
  implicit none
  integer::i,N,Nx,info
  real*16::dx,tol,xbound,x0,tx
  real*16,allocatable::x(:)
  complex*32,allocatable::y(:),ss(:)
  complex*32::rkfc
  external::rkfc
 
  N=2
  allocate(y(1:N))
 
  x0=1q0; dx=1000q0; Nx=2
  allocate(x(0:Nx-1))
  do i=0,Nx-1
     x(i)=x0+dble(i)*dx
  enddo

  y(1)=exp(cmplx(0q0,1q0,kind=16))
  y(2)=cmplx(0.5q0,1q0,kind=16)*y(1)

  call qrk_preparation("rkf45")

  tol=1q-18
  write(10,'(3e30.17e4)')x(0),real(y(1)),imag(y(1))
  do i=1,Nx-1
     info=0; tx=x(i-1); xbound=x(i)
     call qrkf45_e(rkfc,tx,y,xbound,info,tol)
     if(info.eq.-1)write(6,'(A,f10.5,A,f10.5)')"strange point between ",x(i-1)," and ",x(i)
 
     write(10,'(3e30.17e4)')x(i),real(y(1)),imag(y(1))
  enddo
  write(6,*)x(0),y(1)
 
  call qrk_deallocation("rkf45")
 
  stop
end program main
!----------------
function rkfc(N,x,y,s)
  implicit none
  integer,intent(in)::N,s
  real*16,intent(in)::x
  complex*32::y(1:N)
  complex*32::rkfc
 
  rkfc=0q0
  if(s.eq.1)then
     rkfc=y(2)
  elseif(s.eq.2)then
     rkfc=(-1q0/(4q0*x*x)+cmplx(0q0,1q0,kind=16)/x-1q0)*y(1)
  else
     write(6,*)"unknown s, program stop"
     stop
  endif
 
  return
end function rkfc
スポンサーリンク

不連続な点を含む場合


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

さて、ここで微分方程式
\(
\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\)位までですかね。これは、機械的な誤差があることによって不安定な平衡点からほんの少しだけ外れてしまったからです。だからカオスとかを考えるときとかは注意しなければなりません。

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


ルンゲクッタ法の導出は煩雑です。単に複雑なだけです。
導出過程について詳しく述べられているページは、早川様が記述された以下の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)
で無料で公開されています。


「ルンゲクッタ法の説明と刻み幅制御」への1件のフィードバック

コメントを残す

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