時間発展演算子の導出

2015年4月24日 sikino 5件のコメント

量子力学における時間発展演算子$\displaystyle \hat{U}_{(t)}=e^{-i\frac{\hat{H}}{\hbar}t}$の導出です。
これは、ハミルトニアンは時間依存しない場合に限ります。時間依存するときは時間順序積という概念が出てきます。

シュレディンガー方程式
$
\begin{align}
i\hbar\frac{\partial}{\partial t}\psi_{(x,t)}=\hat{H}\psi_{(x,t)}
\end{align}
$
は既知のものとして進めます。

形式的な導出でよく説明されるのは以下の手順です。

シュレディンガー方程式の両辺を$i\hbar$で割って、更に波動関数$\psi_{(x,t)}$で割り両辺を時刻$0\sim t$まで積分します。
すなわち、
$
\begin{align}
\int_{\psi_{(0)}}^{\psi_{(t)}}\frac{d\psi}{\psi}&=\int_0^t\frac{\hat{H}}{i\hbar}dt^{\prime} \\
\ln\frac{\psi_{(t)}}{\psi_{(0)}}&=\frac{-i}{\hbar}\hat{H}t \\
\end{align}
$

なので
$
\displaystyle \psi_{(x,t)}=e^{-i\frac{\hat{H}}{\hbar}t}\psi_{(x,0)}
$
証明終了,,

(;’;ﾟ;ж;ﾟ;｀;)ﾌﾞﾎｫｯ
(; ^ω^)は…

何で波動関数で割ってくれちゃってるんですか？ハミルトニアンは演算子ですよ。形式的にしてもひどくない？

ということで、ちゃんと求めましょう。結果は↑のものと同じになります。不思議ですね。

ハミルトニアンが時間依存しないという条件の下求めます。
まず、時間依存するシュレディンガー方程式より、
$
\displaystyle i\hbar\frac{\partial}{\partial t}\psi_{(x,t)}=\hat{H}\psi_{(x,t)}
$
となります。無限小時間$\Delta t$に対して、シュレディンガー方程式は、
$
\displaystyle i\hbar\lim_{\Delta t\rightarrow 0}\frac{\psi_{(x,t+\Delta t)}-\psi_{(x,t)}}{\Delta t}=\hat{H}\psi_{(x,t)}
$
です。以降$\lim$は省略します。

式を変形して、
$
\displaystyle \psi_{(x,t+\Delta t)}=\left(1-i\frac{\hat{H}}{\hbar}\Delta t\right)\psi_{(x,t)}
$
また、もう$\Delta t$だけ時間を進めると、
$
\displaystyle \psi_{(x,t+2\Delta t)}=\left(1-i\frac{\hat{H}}{\hbar}\Delta t\right)^2\psi_{(x,t)}
$
という結果が得られます。

ここで、$\Delta t$をN回作用させることを考えます。
すなわち、$t^{\prime}=N\Delta t$とおき、$\Delta t=\frac{t^{\prime}}{N}$として考えれば、任意の時間$t^{\prime}$に対して、
$
\displaystyle \psi_{(x,t+t^{\prime})}=\left(1-i\frac{\hat{H}}{\hbar}t^{\prime}\frac{1}{N}\right)^N\psi_{(x,t)}
$
という式が成立します。

無限回、微小区間$\Delta t$を動かす操作を考えればいいので、$N\rightarrow \infty$として考えれば、
$
\displaystyle \psi_{(x,t+t^{\prime})}=\lim_{N\rightarrow \infty}\left(1-i\frac{\hat{H}}{\hbar}t^{\prime}\frac{1}{N}\right)^N\psi_{(x,t)}
$
公式
$
\displaystyle \lim_{x\rightarrow \pm\infty}\left(1+\frac{a}{x}\right)^x=e^a
$
より、
$
\displaystyle \psi_{(x,t+t^{\prime})}=e^{-i\frac{\hat{H}}{\hbar}t^{\prime}}\psi_{(x,t)}
$

$t=0$と置き、$t^{\prime}=t$と文字を置き換えれば、
$
\displaystyle \psi_{(x,t)}=e^{-i\frac{\hat{H}}{\hbar}t}\psi_{(x,0)}
$

となり、ちゃんと時間発展演算子が導けました。めでたしめでたし。

ハミルトニアンが時間依存する場合、
$
\displaystyle i\hbar\frac{\partial}{\partial t}\psi_{(x,t)}=\hat{H}_{(t)}\psi_{(x,t)}
$
となります。無限小時間$\Delta t$に対して、シュレディンガー方程式は、
$
\displaystyle i\hbar\lim_{\Delta t\rightarrow 0}\frac{\psi_{(x,t+\Delta t)}-\psi_{(x,t)}}{\Delta t}=\hat{H}_{(t)}\psi_{(x,t)}
$
次の時刻では
$
\displaystyle \psi_{(x,t+2\Delta t)}=
\left(1-i\frac{\hat{H}_{(t+\Delta t)}}{\hbar}\Delta t\right)\left(1-i\frac{\hat{H}_{(t)}}{\hbar}\Delta t\right)\psi_{(x,t)}
$
となります。２乗にはなってくれません。
と、どんどん作用させていくわけです。時間依存しない時と明らかに違ったものになりますが、詳しくは次の機会のお話で。

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

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

2015年4月19日 sikino 5件のコメント

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

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

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

埋め込まれた陽的ルンゲ=クッタ法
ルンゲ=クッタ=フェールベルグ法による刻み幅の自動制御(適応刻み幅制御)
ルンゲ=クッタ=フェールベルグ法による刻み幅の自動制御のプログラム(実数/複素数で１/２階微分方程式を解くプログラム)
不連続な点を含む場合
刻み幅制御のベンチマーク(振り子)
（追記）ルンゲ=クッタ=ドルマンド=プリンス法
陽的ルンゲ=クッタ法の導出
参考文献

理論はいいから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$回繰り返すような計算では、オーダーが１つ落ちる事に注意しましょう。

[adsense1]

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

さて、ここまで“段”は計算コスト、で“次”は計算の正確さ、という曖昧な表現でしたが、その表現をちゃんと知りましょう。
段と次を知るためにはルンゲ=クッタ法の計算方法を知る必要があります。
具体例を載せます。
$
\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}$を

としてまとめて書く表記法です。

再び、オイラー法はButcher tableで書くと

とまとめて書くことができます。
中点法は

RK4は

です。

高次のルンゲ=クッタ法(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は

となっているので注意が必要です。

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

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

Butcher tableは、この場合extended Butcher tableと呼ばれ、こういう形式で書かれます。

この埋め込まれたルンゲ=クッタ法のいいところは、

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

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

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

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

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

具体例を挙げましょう。
一番簡単な埋め込まれたルンゲ=クッタ法は、ホイン法と呼ばれています。

1行目は2次のオーダーを持ち、2行目は1次のオーダーを持ちます。

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

と書かれます。1行目は5次のオーダー、2行目は4次のオーダーを持ちます。

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

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

[5~9]によれば、ルンゲ=クッタ=フェールベルグ法において区間$i$での最適な刻み幅$h’$は区間$i$の誤差評価の結果を使って、
$
\displaystyle h’=\delta h=\left(\frac{\varepsilon h}{2|x^{(5)}_{i}-x^{(4)}_{i}|}\right)^{1/4} h
$
と予想できます。ここで$\varepsilon$はエラーを制御する精度の目安で、おおよそ計算時に要求する相対誤差です。もちろん、この$h’$は区間$i$の最適な幅ですが、関数に劇的な変化は無いだろうとして、次の区間の計算の一番初めに用いる刻み幅を推定するのです。
なので、$i+1$番目の計算区間では、計算するときはこの$h’$の値を使えばいいんです。
(ちなみに、m次ルンゲ=クッタ法の場合では
$
\displaystyle h’=\delta h=\left(\frac{\varepsilon h}{2|x^{(m+1)}_{i}-x^{(m)}_{i}|}\right)^{1/m} h
$
と推測されます。)

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

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

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

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

実数、1階微分方程式の場合

微分方程式
$
\displaystyle \frac{d}{dx}y=y\cos x,~~y(x=0)=1
$
を倍精度実数、刻み幅制御で$x=10$まで解く事を考えます。
解析解は
$
\displaystyle y(x)=exp(sin(x))
$
です。コードは以下の通りです。

▼ここクリックでこの場に展開

program main
implicit none
integer::i,N,info
double precision::x,h,tol,xbound
double precision,allocatable::y(:)
external::grk

N=1 ! Number of differential equations
allocate(y(1:N))

x=0d0; xbound=10d0
y(1)=1d0 !initial condition

h=xbound-x
i=0; info=0; tol=1d-8;
do while(info.le.0)
call drkf45(grk,x,h,N,y,xbound,info,tol)
write(10,'(3e25.10e3)')x,y(1),h
i=i+1
enddo
write(6,*)"Number of referenced times -->",i

stop
end program main

subroutine grk(N,x,y,f)
implicit none
integer,intent(in)::N
double precision,intent(in)::x,y(1:N)
double precision,intent(out)::f(1:N)

! Solve
! d y(1) / dt = y(1) * cos(t)

f(1)=y(1)*cos(x)

return
end subroutine grk

!===============================

subroutine drkf45(grk,x,h,N,y,xbound,info,tol)
! if h < hmin, propagate forcibly with warning.
!
!-----------------
!info = -9 (maybe path the discontinue points)
! = 0 (Running now)
! = 1 (x reach xbound)
!-----------------
!
implicit none
integer,intent(in)::N
double precision,intent(in)::xbound,tol
double precision,intent(inout)::x,h,y(1:N)
integer,intent(inout)::info

integer::i,j,FLAG,key
double precision::R,delta,tx,Sy,err
double precision,allocatable::ty(:),K(:,:),tf(:)
double precision,parameter::hmin=1d-14,hmax=0.5d0
integer,parameter::s=6
double precision::a(1:s,1:s),b1(1:s),b2(1:s),c(1:s),Rc(1:s)
external::grk

c(1:6)=(/0d0, 0.25d0, 0.375d0,&
0.9230769230769230769230769230769230769231d0, 1d0, 0.5d0/)
a(1:6,1:6)=0d0
a(1,1:6)=(/0d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(2,1:6)=(/0.25d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(3,1:6)=(/0.09375d0, 0.28125d0, 0d0, 0d0, 0d0, 0d0/)
a(4,1:6)=(/0.8793809740555302685480200273099681383705d0, &
-3.277196176604460628129267182521620391443d0, &
3.320892125625853436504324078288575329995d0, 0d0, 0d0, 0d0/)
a(5,1:6)=(/2.032407407407407407407407407407407407407d0,-8d0, &
7.173489278752436647173489278752436647173d0, &
-0.2058966861598440545808966861598440545809d0, 0d0, 0d0/)
a(6,1:6)=(/-0.2962962962962962962962962962962962962963d0,2d0, &
-1.381676413255360623781676413255360623782d0, &
0.4529727095516569200779727095516569200780d0,-0.275d0,0d0/)
b2(1:6)=(/0.1185185185185185185185185185185185185185d0, 0.d0,&
0.5189863547758284600389863547758284600390d0, &
0.5061314903420166578061314903420166578061d0, &
-0.18d0, 0.03636363636363636363636363636363636363636d0/)
b1(1:6)=(/0.1157407407407407407407407407407407407407d0, 0d0,&
0.5489278752436647173489278752436647173489d0, &
0.5353313840155945419103313840155945419103d0, -0.2d0, 0d0/)
Rc(1:6)=(/0.002777777777777777777777777777777777777778d0,0d0, &
-0.02994152046783625730994152046783625730994d0, &
-0.02919989367357788410419989367357788410420d0, 0.02d0, &
0.03636363636363636363636363636363636363636d0/)

key=0
allocate(ty(1:N),tf(1:N),K(1:s,1:N))
ty=0d0; tf=0d0; K=0d0

if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
endif

if(h.ge.abs(xbound-x))h=xbound-x

FLAG=1
if(abs(x-xbound).le.hmin)then
info=1
FLAG=0
endif

do while(FLAG.eq.1)
tx=x
do j=1,s
tx=x+c(j)*h
ty(1:N)=y(1:N)
do i=1,j-1
ty(1:N)=ty(1:N)+K(i,1:N)*a(j,i)
enddo
call grk(N,tx,ty,tf)
K(j,1:N)=h*tf(1:N)
enddo

!step 4
R=0d0
do i=1,N
R=R+(Rc(1)*K(1,i)+Rc(3)*K(3,i)+Rc(4)*K(4,i) &
+Rc(5)*K(5,i)+Rc(6)*K(6,i))**2d0
enddo
R=abs(dsqrt(R)/h/dble(N))

Sy=0d0
do i=1,N
Sy=Sy+(y(i)*y(i))
enddo
Sy=dsqrt(Sy)
if(Sy.ge.1d0)then
err=tol*Sy
else
err=tol
endif

!step 5
if(R.le.err.or.key.eq.1)then
x=x+h
y(1:N)=y(1:N)+b1(1)*K(1,1:N)+b1(3)*K(3,1:N) &
+b1(4)*K(4,1:N)+b1(5)*K(5,1:N)
FLAG=0
endif

!step 6
! Avoid zero deviding.
if(R.ge.1d-20)then
delta=(err/(2d0*R))**0.25d0
else
delta=4d0
endif

!step 7
if(delta.le.0.1d0)then
!function changes dramatically.
h=0.1d0*h
elseif(delta.ge.4d0)then
!function changes loosely.
h=4d0*h
else
!function changes moderately.
h=delta*h
endif

!step 8
if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
elseif(abs(h).lt.hmin)then
h=sign(1d0,h)*hmin
key=1
endif

!step 9
if(abs(xbound-x).le.abs(h))then
h=xbound-x
if(abs(h).le.hmin)then
info=1
FLAG=0
endif
end if

if(h.le.0d0.and.xbound-x.ge.0d0)then
info=1
FLAG=0
elseif(h.ge.0d0.and.xbound-x.le.0d0)then
info=1
FLAG=0
endif
enddo

if(key.eq.1)then
write(6,'(A,f10.5,A,f10.5)')&
"Strange point between ",x-h," and ",x
info=-9
endif

deallocate(ty,tf,K)
return
end subroutine drkf45

実数、2階微分方程式の場合

微分方程式
$
\displaystyle \frac{d^2}{dx^2}y=-\frac{1}{4}y,~~y(x=0)=1,~ y'(x=0)=0
$
を倍精度実数、刻み幅制御で$x=20$まで解く事を考えます。コードは以下の通りです。

▼ここクリックでこの場に展開

program main
implicit none
integer::i,N,info
double precision::x,h,tol,xbound
double precision,allocatable::y(:)
external::grk

N=2
allocate(y(1:N))

x=0d0; xbound=20d0
!initial condition
y(1)=1d0 ! y (0) = 1d0
y(2)=0d0 ! y'(0) = 0d0

h=xbound-x
i=0; info=0; tol=1d-8;
do while(info.le.0)
call drkf45(grk,x,h,N,y,xbound,info,tol)
write(10,'(4e25.10e3)')x,y(1),y(2),h
i=i+1
enddo
write(6,*)"Number of referenced times -->",i

stop
end program main

subroutine grk(N,x,y,f)
implicit none
integer,intent(in)::N
double precision,intent(in)::x,y(1:N)
double precision,intent(out)::f(1:N)

! Solve
! d^2 y(1) / dt^2 = - 0.5 * y(1)

f(1)=y(2)
f(2)=-0.25d0*y(1)

return
end subroutine grk

!===============================

subroutine drkf45(grk,x,h,N,y,xbound,info,tol)
! if h < hmin, propagate forcibly with warning.
!
!-----------------
!info = -9 (maybe path the discontinue points)
! = 0 (Running now)
! = 1 (x reach xbound)
!-----------------
!
implicit none
integer,intent(in)::N
double precision,intent(in)::xbound,tol
double precision,intent(inout)::x,h,y(1:N)
integer,intent(inout)::info

integer::i,j,FLAG,key
double precision::R,delta,tx,Sy,err
double precision,allocatable::ty(:),K(:,:),tf(:)
double precision,parameter::hmin=1d-14,hmax=0.5d0
integer,parameter::s=6
double precision::a(1:s,1:s),b1(1:s),b2(1:s),c(1:s),Rc(1:s)
external::grk

c(1:6)=(/0d0, 0.25d0, 0.375d0,&
0.9230769230769230769230769230769230769231d0, 1d0, 0.5d0/)
a(1:6,1:6)=0d0
a(1,1:6)=(/0d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(2,1:6)=(/0.25d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(3,1:6)=(/0.09375d0, 0.28125d0, 0d0, 0d0, 0d0, 0d0/)
a(4,1:6)=(/0.8793809740555302685480200273099681383705d0, &
-3.277196176604460628129267182521620391443d0, &
3.320892125625853436504324078288575329995d0, 0d0, 0d0, 0d0/)
a(5,1:6)=(/2.032407407407407407407407407407407407407d0,-8d0, &
7.173489278752436647173489278752436647173d0, &
-0.2058966861598440545808966861598440545809d0, 0d0, 0d0/)
a(6,1:6)=(/-0.2962962962962962962962962962962962962963d0,2d0, &
-1.381676413255360623781676413255360623782d0, &
0.4529727095516569200779727095516569200780d0,-0.275d0,0d0/)
b2(1:6)=(/0.1185185185185185185185185185185185185185d0, 0.d0,&
0.5189863547758284600389863547758284600390d0, &
0.5061314903420166578061314903420166578061d0, &
-0.18d0, 0.03636363636363636363636363636363636363636d0/)
b1(1:6)=(/0.1157407407407407407407407407407407407407d0, 0d0,&
0.5489278752436647173489278752436647173489d0, &
0.5353313840155945419103313840155945419103d0, -0.2d0, 0d0/)
Rc(1:6)=(/0.002777777777777777777777777777777777777778d0,0d0, &
-0.02994152046783625730994152046783625730994d0, &
-0.02919989367357788410419989367357788410420d0, 0.02d0, &
0.03636363636363636363636363636363636363636d0/)

key=0
allocate(ty(1:N),tf(1:N),K(1:s,1:N))
ty=0d0; tf=0d0; K=0d0

if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
endif

if(h.ge.abs(xbound-x))h=xbound-x

FLAG=1
if(abs(x-xbound).le.hmin)then
info=1
FLAG=0
endif

do while(FLAG.eq.1)
tx=x
do j=1,s
tx=x+c(j)*h
ty(1:N)=y(1:N)
do i=1,j-1
ty(1:N)=ty(1:N)+K(i,1:N)*a(j,i)
enddo
call grk(N,tx,ty,tf)
K(j,1:N)=h*tf(1:N)
enddo

!step 4
R=0d0
do i=1,N
R=R+(Rc(1)*K(1,i)+Rc(3)*K(3,i)+Rc(4)*K(4,i) &
+Rc(5)*K(5,i)+Rc(6)*K(6,i))**2d0
enddo
R=abs(dsqrt(R)/h/dble(N))

Sy=0d0
do i=1,N
Sy=Sy+(y(i)*y(i))
enddo
Sy=dsqrt(Sy)
if(Sy.ge.1d0)then
err=tol*Sy
else
err=tol
endif

!step 5
if(R.le.err.or.key.eq.1)then
x=x+h
y(1:N)=y(1:N)+b1(1)*K(1,1:N)+b1(3)*K(3,1:N) &
+b1(4)*K(4,1:N)+b1(5)*K(5,1:N)
FLAG=0
endif

!step 6
! Avoid zero deviding.
if(R.ge.1d-20)then
delta=(err/(2d0*R))**0.25d0
else
delta=4d0
endif

!step 7
if(delta.le.0.1d0)then
!function changes dramatically.
h=0.1d0*h
elseif(delta.ge.4d0)then
!function changes loosely.
h=4d0*h
else
!function changes moderately.
h=delta*h
endif

!step 8
if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
elseif(abs(h).lt.hmin)then
h=sign(1d0,h)*hmin
key=1
endif

!step 9
if(abs(xbound-x).le.abs(h))then
h=xbound-x
if(abs(h).le.hmin)then
info=1
FLAG=0
endif
end if

if(h.le.0d0.and.xbound-x.ge.0d0)then
info=1
FLAG=0
elseif(h.ge.0d0.and.xbound-x.le.0d0)then
info=1
FLAG=0
endif
enddo

if(key.eq.1)then
write(6,'(A,f10.5,A,f10.5)')&
"Strange point between ",x-h," and ",x
info=-9
endif

deallocate(ty,tf,K)
return
end subroutine drkf45

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

複素数、1階微分方程式の場合

微分方程式
$
\displaystyle \frac{d}{dx}y=y\cos x,~~ y(x=0)=1+i\frac{1}{2}
$
を倍精度実数、刻み幅制御で$x=10$まで解く事を考えます。コードは以下の通りです。

▼ここクリックでこの場に展開

program main
implicit none
integer::i,N,info
double precision::x,h,tol,xbound
complex(kind(0d0)),allocatable::y(:)
external::grk

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

x=0d0; xbound=10d0
y(1)=dcmplx(1d0,0.5d0) !initial condition

h=xbound-x
i=0; info=0; tol=1d-8;
do while(info.le.0)
call crkf45(grk,x,h,N,y,xbound,info,tol)
write(10,'(4e25.10e3)')x,y(1),h
i=i+1
enddo
write(6,*)"Number of referenced times -->",i

stop
end program main

subroutine grk(N,x,y,f)
implicit none
integer,intent(in)::N
double precision,intent(in)::x
complex(kind(0d0)),intent(in)::y(1:N)
complex(kind(0d0)),intent(out)::f(1:N)

f(1)=y(1)*cos(x)

return
end subroutine grk

!===============================

subroutine crkf45(grk,x,h,N,y,xbound,info,tol)
! if h < hmin, propagate forcibly with warning.
!
!-----------------
!info = -9 (maybe path the discontinue points)
! = 0 (Running now)
! = 1 (x reach xbound)
!-----------------
!
implicit none
integer,intent(in)::N
double precision,intent(in)::xbound,tol
double precision,intent(inout)::x,h
complex(kind(0d0)),intent(inout)::y(1:N)
integer,intent(inout)::info
external::grk

integer::i,j,FLAG,key
double precision::R,delta,tx,Sy,err
complex(kind(0d0)),allocatable::ty(:),tf(:),K(:,:)
double precision,parameter::hmin=1d-14,hmax=0.5d0
integer,parameter::s=6
double precision::a(1:s,1:s),b1(1:s),b2(1:s),c(1:s),Rc(1:s)

c(1:6)=(/0d0, 0.25d0, 0.375d0,&
0.9230769230769230769230769230769230769231d0, 1d0, 0.5d0/)
a(1:6,1:6)=0d0
a(1,1:6)=(/0d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(2,1:6)=(/0.25d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(3,1:6)=(/0.09375d0, 0.28125d0, 0d0, 0d0, 0d0, 0d0/)
a(4,1:6)=(/0.8793809740555302685480200273099681383705d0, &
-3.277196176604460628129267182521620391443d0, &
3.320892125625853436504324078288575329995d0, 0d0, 0d0, 0d0/)
a(5,1:6)=(/2.032407407407407407407407407407407407407d0,-8d0, &
7.173489278752436647173489278752436647173d0, &
-0.2058966861598440545808966861598440545809d0, 0d0, 0d0/)
a(6,1:6)=(/-0.2962962962962962962962962962962962962963d0,2d0, &
-1.381676413255360623781676413255360623782d0, &
0.4529727095516569200779727095516569200780d0,-0.275d0,0d0/)
b2(1:6)=(/0.1185185185185185185185185185185185185185d0, 0.d0,&
0.5189863547758284600389863547758284600390d0, &
0.5061314903420166578061314903420166578061d0, &
-0.18d0, 0.03636363636363636363636363636363636363636d0/)
b1(1:6)=(/0.1157407407407407407407407407407407407407d0, 0d0,&
0.5489278752436647173489278752436647173489d0, &
0.5353313840155945419103313840155945419103d0, -0.2d0, 0d0/)
Rc(1:6)=(/0.002777777777777777777777777777777777777778d0,0d0, &
-0.02994152046783625730994152046783625730994d0, &
-0.02919989367357788410419989367357788410420d0, 0.02d0, &
0.03636363636363636363636363636363636363636d0/)

key=0
allocate(ty(1:N),tf(1:N),K(1:s,1:N))
ty=dcmplx(0d0,0d0); K=dcmplx(0d0,0d0)
tf=dcmplx(0d0,0d0)

if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
endif

if(h.ge.abs(xbound-x))h=xbound-x

FLAG=1
if(abs(x-xbound).le.hmin)then
info=1
FLAG=0
endif

do while(FLAG.eq.1)
tx=x
do j=1,s
tx=x+c(j)*h
ty(1:N)=y(1:N)
do i=1,j-1
ty(1:N)=ty(1:N)+K(i,1:N)*a(j,i)
enddo
!do i=1,N
! K(j,i)=h*grk(N,tx,ty,i)
!enddo
call grk(N,tx,ty,tf)
K(j,1:N)=h*tf(1:N)
enddo

!step 4
R=0d0
do i=1,N
R=R+abs(Rc(1)*K(1,i)+Rc(3)*K(3,i)+Rc(4)*K(4,i) &
+Rc(5)*K(5,i)+Rc(6)*K(6,i))**2d0
enddo
R=abs(dsqrt(R)/h/dble(N))

Sy=0d0
do i=1,N
Sy=Sy+abs(y(i)*y(i))
enddo
Sy=dsqrt(Sy)
if(Sy.ge.1d0)then
err=tol*Sy
else
err=tol
endif

!step 5
if(R.le.err.or.key.eq.1)then
x=x+h
y(1:N)=y(1:N)+b1(1)*K(1,1:N)+b1(3)*K(3,1:N) &
+b1(4)*K(4,1:N)+b1(5)*K(5,1:N)
FLAG=0
endif

!step 6
! Avoid zero deviding.
if(R.ge.1d-20)then
delta=(err/(2d0*R))**0.25d0
else
delta=4d0
endif

!step 7
if(delta.le.0.1d0)then
!function changes dramatically.
h=0.1d0*h
elseif(delta.ge.4d0)then
!function changes loosely.
h=4d0*h
else
!function changes moderately.
h=delta*h
endif

!step 8
if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
elseif(abs(h).lt.hmin)then
h=sign(1d0,h)*hmin
key=1
endif

!step 9
if(abs(xbound-x).le.abs(h))then
h=xbound-x
if(abs(h).le.hmin)then
info=1
FLAG=0
endif
end if

if(h.le.0d0.and.xbound-x.ge.0d0)then
info=1
FLAG=0
elseif(h.ge.0d0.and.xbound-x.le.0d0)then
info=1
FLAG=0
endif
enddo

if(key.eq.1)then
write(6,'(A,f10.5,A,f10.5)')&
"Strange point between ",x-h," and ",x
info=-9
endif

deallocate(ty,tf,K)
return
end subroutine crkf45

複素数、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$まで解く事を考えます。コードは以下の通りです。

▼ここクリックでこの場に展開

program main
implicit none
integer::i,N,info
double precision::x,h,tol,xbound
complex(kind(0d0)),allocatable::y(:)
external::grk

N=2
allocate(y(1:N))

x=0d0; xbound=20d0
y(1)=dcmplx(1d0,0.5d0) !initial condition
y(2)=dcmplx(0d0,0d0) !initial condition

h=xbound-x
i=0; info=0; tol=1d-8;
do while(info.le.0)
call crkf45(grk,x,h,N,y,xbound,info,tol)
write(10,'(6e25.10e3)')x,dble(y(1)),dimag(y(1)) &
,dble(y(2)),dimag(y(2)),h
i=i+1
enddo
write(6,*)"Number of referenced times -->",i

stop
end program main

subroutine grk(N,x,y,f)
implicit none
integer,intent(in)::N
double precision,intent(in)::x
complex(kind(0d0)),intent(in)::y(1:N)
complex(kind(0d0)),intent(out)::f(1:N)

f(1)=y(2)
f(2)=-0.25d0*y(1)

return
end subroutine grk

!===============================

subroutine crkf45(grk,x,h,N,y,xbound,info,tol)
! if h < hmin, propagate forcibly with warning.
!
!-----------------
!info = -9 (maybe path the discontinue points)
! = 0 (Running now)
! = 1 (x reach xbound)
!-----------------
!
implicit none
integer,intent(in)::N
double precision,intent(in)::xbound,tol
double precision,intent(inout)::x,h
complex(kind(0d0)),intent(inout)::y(1:N)
integer,intent(inout)::info
external::grk

integer::i,j,FLAG,key
double precision::R,delta,tx,Sy,err
complex(kind(0d0)),allocatable::ty(:),tf(:),K(:,:)
double precision,parameter::hmin=1d-14,hmax=0.5d0
integer,parameter::s=6
double precision::a(1:s,1:s),b1(1:s),b2(1:s),c(1:s),Rc(1:s)

c(1:6)=(/0d0, 0.25d0, 0.375d0,&
0.9230769230769230769230769230769230769231d0, 1d0, 0.5d0/)
a(1:6,1:6)=0d0
a(1,1:6)=(/0d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(2,1:6)=(/0.25d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(3,1:6)=(/0.09375d0, 0.28125d0, 0d0, 0d0, 0d0, 0d0/)
a(4,1:6)=(/0.8793809740555302685480200273099681383705d0, &
-3.277196176604460628129267182521620391443d0, &
3.320892125625853436504324078288575329995d0, 0d0, 0d0, 0d0/)
a(5,1:6)=(/2.032407407407407407407407407407407407407d0,-8d0, &
7.173489278752436647173489278752436647173d0, &
-0.2058966861598440545808966861598440545809d0, 0d0, 0d0/)
a(6,1:6)=(/-0.2962962962962962962962962962962962962963d0,2d0, &
-1.381676413255360623781676413255360623782d0, &
0.4529727095516569200779727095516569200780d0,-0.275d0,0d0/)
b2(1:6)=(/0.1185185185185185185185185185185185185185d0, 0.d0,&
0.5189863547758284600389863547758284600390d0, &
0.5061314903420166578061314903420166578061d0, &
-0.18d0, 0.03636363636363636363636363636363636363636d0/)
b1(1:6)=(/0.1157407407407407407407407407407407407407d0, 0d0,&
0.5489278752436647173489278752436647173489d0, &
0.5353313840155945419103313840155945419103d0, -0.2d0, 0d0/)
Rc(1:6)=(/0.002777777777777777777777777777777777777778d0,0d0, &
-0.02994152046783625730994152046783625730994d0, &
-0.02919989367357788410419989367357788410420d0, 0.02d0, &
0.03636363636363636363636363636363636363636d0/)

key=0
allocate(ty(1:N),tf(1:N),K(1:s,1:N))
ty=dcmplx(0d0,0d0); K=dcmplx(0d0,0d0)
tf=dcmplx(0d0,0d0)

if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
endif

if(h.ge.abs(xbound-x))h=xbound-x

FLAG=1
if(abs(x-xbound).le.hmin)then
info=1
FLAG=0
endif

do while(FLAG.eq.1)
tx=x
do j=1,s
tx=x+c(j)*h
ty(1:N)=y(1:N)
do i=1,j-1
ty(1:N)=ty(1:N)+K(i,1:N)*a(j,i)
enddo
call grk(N,tx,ty,tf)
K(j,1:N)=h*tf(1:N)
enddo

!step 4
R=0d0
do i=1,N
R=R+abs(Rc(1)*K(1,i)+Rc(3)*K(3,i)+Rc(4)*K(4,i) &
+Rc(5)*K(5,i)+Rc(6)*K(6,i))**2d0
enddo
R=abs(dsqrt(R)/h/dble(N))

Sy=0d0
do i=1,N
Sy=Sy+abs(y(i)*y(i))
enddo
Sy=dsqrt(Sy)
if(Sy.ge.1d0)then
err=tol*Sy
else
err=tol
endif

!step 5
if(R.le.err.or.key.eq.1)then
x=x+h
y(1:N)=y(1:N)+b1(1)*K(1,1:N)+b1(3)*K(3,1:N) &
+b1(4)*K(4,1:N)+b1(5)*K(5,1:N)
FLAG=0
endif

!step 6
! Avoid zero deviding.
if(R.ge.1d-20)then
delta=(err/(2d0*R))**0.25d0
else
delta=4d0
endif

!step 7
if(delta.le.0.1d0)then
!function changes dramatically.
h=0.1d0*h
elseif(delta.ge.4d0)then
!function changes loosely.
h=4d0*h
else
!function changes moderately.
h=delta*h
endif

!step 8
if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
elseif(abs(h).lt.hmin)then
h=sign(1d0,h)*hmin
key=1
endif

!step 9
if(abs(xbound-x).le.abs(h))then
h=xbound-x
if(abs(h).le.hmin)then
info=1
FLAG=0
endif
end if

if(h.le.0d0.and.xbound-x.ge.0d0)then
info=1
FLAG=0
elseif(h.ge.0d0.and.xbound-x.le.0d0)then
info=1
FLAG=0
endif
enddo

if(key.eq.1)then
write(6,'(A,f10.5,A,f10.5)')&
"Strange point between ",x-h," and ",x
info=-9
endif

deallocate(ty,tf,K)
return
end subroutine crkf45

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

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

▼ここクリックでこの場に展開

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

N=1
allocate(y(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 drkf45(grk,tx,h,N,y,x(i),info,tol)
enddo
write(10,'(2e25.10e3)')x(i),y(1)
enddo

stop
end program main

subroutine grk(N,x,y,f)
implicit none
integer,intent(in)::N
double precision,intent(in)::x,y(1:N)
double precision,intent(out)::f(1:N)

! Solve
! d y(1) / dt = y(1) * cos(t)

f(1)=y(1)*cos(x)

return
end subroutine grk

!===============================

subroutine drkf45(grk,x,h,N,y,xbound,info,tol)
! if h < hmin, propagate forcibly with warning.
!
!-----------------
!info = -9 (maybe path the discontinue points)
! = 0 (Running now)
! = 1 (x reach xbound)
!-----------------
!
implicit none
integer,intent(in)::N
double precision,intent(in)::xbound,tol
double precision,intent(inout)::x,h,y(1:N)
integer,intent(inout)::info

integer::i,j,FLAG,key
double precision::R,delta,tx,Sy,err
double precision,allocatable::ty(:),K(:,:),tf(:)
double precision,parameter::hmin=1d-14,hmax=0.5d0
integer,parameter::s=6
double precision::a(1:s,1:s),b1(1:s),b2(1:s),c(1:s),Rc(1:s)
external::grk

c(1:6)=(/0d0, 0.25d0, 0.375d0,&
0.9230769230769230769230769230769230769231d0, 1d0, 0.5d0/)
a(1:6,1:6)=0d0
a(1,1:6)=(/0d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(2,1:6)=(/0.25d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(3,1:6)=(/0.09375d0, 0.28125d0, 0d0, 0d0, 0d0, 0d0/)
a(4,1:6)=(/0.8793809740555302685480200273099681383705d0, &
-3.277196176604460628129267182521620391443d0, &
3.320892125625853436504324078288575329995d0, 0d0, 0d0, 0d0/)
a(5,1:6)=(/2.032407407407407407407407407407407407407d0,-8d0, &
7.173489278752436647173489278752436647173d0, &
-0.2058966861598440545808966861598440545809d0, 0d0, 0d0/)
a(6,1:6)=(/-0.2962962962962962962962962962962962962963d0,2d0, &
-1.381676413255360623781676413255360623782d0, &
0.4529727095516569200779727095516569200780d0,-0.275d0,0d0/)
b2(1:6)=(/0.1185185185185185185185185185185185185185d0, 0.d0,&
0.5189863547758284600389863547758284600390d0, &
0.5061314903420166578061314903420166578061d0, &
-0.18d0, 0.03636363636363636363636363636363636363636d0/)
b1(1:6)=(/0.1157407407407407407407407407407407407407d0, 0d0,&
0.5489278752436647173489278752436647173489d0, &
0.5353313840155945419103313840155945419103d0, -0.2d0, 0d0/)
Rc(1:6)=(/0.002777777777777777777777777777777777777778d0,0d0, &
-0.02994152046783625730994152046783625730994d0, &
-0.02919989367357788410419989367357788410420d0, 0.02d0, &
0.03636363636363636363636363636363636363636d0/)

key=0
allocate(ty(1:N),tf(1:N),K(1:s,1:N))
ty=0d0; tf=0d0; K=0d0

if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
endif

if(h.ge.abs(xbound-x))h=xbound-x

FLAG=1
if(abs(x-xbound).le.hmin)then
info=1
FLAG=0
endif

do while(FLAG.eq.1)
tx=x
do j=1,s
tx=x+c(j)*h
ty(1:N)=y(1:N)
do i=1,j-1
ty(1:N)=ty(1:N)+K(i,1:N)*a(j,i)
enddo
call grk(N,tx,ty,tf)
K(j,1:N)=h*tf(1:N)
enddo

!step 4
R=0d0
do i=1,N
R=R+(Rc(1)*K(1,i)+Rc(3)*K(3,i)+Rc(4)*K(4,i) &
+Rc(5)*K(5,i)+Rc(6)*K(6,i))**2d0
enddo
R=abs(dsqrt(R)/h/dble(N))

Sy=0d0
do i=1,N
Sy=Sy+(y(i)*y(i))
enddo
Sy=dsqrt(Sy)
if(Sy.ge.1d0)then
err=tol*Sy
else
err=tol
endif

!step 5
if(R.le.err.or.key.eq.1)then
x=x+h
y(1:N)=y(1:N)+b1(1)*K(1,1:N)+b1(3)*K(3,1:N) &
+b1(4)*K(4,1:N)+b1(5)*K(5,1:N)
FLAG=0
endif

!step 6
! Avoid zero deviding.
if(R.ge.1d-20)then
delta=(err/(2d0*R))**0.25d0
else
delta=4d0
endif

!step 7
if(delta.le.0.1d0)then
!function changes dramatically.
h=0.1d0*h
elseif(delta.ge.4d0)then
!function changes loosely.
h=4d0*h
else
!function changes moderately.
h=delta*h
endif

!step 8
if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
elseif(abs(h).lt.hmin)then
h=sign(1d0,h)*hmin
key=1
endif

!step 9
if(abs(xbound-x).le.abs(h))then
h=xbound-x
if(abs(h).le.hmin)then
info=1
FLAG=0
endif
end if

if(h.le.0d0.and.xbound-x.ge.0d0)then
info=1
FLAG=0
elseif(h.ge.0d0.and.xbound-x.le.0d0)then
info=1
FLAG=0
endif
enddo

if(key.eq.1)then
write(6,'(A,f10.5,A,f10.5)')&
"Strange point between ",x-h," and ",x
info=-9
endif

deallocate(ty,tf,K)
return
end subroutine drkf45

[adsense2]

不連続な点を含む場合

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

さて、ここで微分方程式
$
\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次ルンゲ=クッタと適応刻み幅ルンゲ=クッタで解いてみましょう。
プログラム上ではそのまま解くことが出来ます。
実際に解かせてみますと、

となり、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から求めました。
4EllipticK[0.950.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)$を取った時のグラフです。

すぐに破綻しました。正しい値は$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
double precision::x,h,tol,xbound
double precision,allocatable::y(:),work0(:)
double precision,external::grk

x=0d0; xbound=10d0

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

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

h=xbound-x
info=-1
tol=1d-8
i=0
write(20,'(3e30.16e3)')x,y(1),h
do while(info.le.0)
call dDP45(grk,x,h,N,y,xbound,info,tol,work0)
write(20,'(3e30.16e3)')x,y(1),h
i=i+1
enddo
write(6,*)"the Number, routine called-->",i

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

!===============================

subroutine dDP45(grk,x,h,N,y,xbound,info,tol,work0)
implicit none
integer,intent(in)::N
double precision,intent(in)::xbound,tol
double precision,intent(inout)::x,h,y(1:N),work0(1:N)
integer,intent(inout)::info
double precision,external::grk

! Runge-Kutta-Dormand-Prince method
!
! if h < hmin, propagate forcibly with warning.
!
!-----------------
!work0(1:N) : Just prepare, used and substituted by
! the coefficients K(7,1:N) to reduce computation costs
!info = -1 (Used for first loop;
! you should set info=0 when you try to use this routine)
! = 0 (Running now)
! = 1 (x reach xbound)
! = -2 (Path the discontinue points)
!-----------------
!

integer::i,j,FLAG,key
double precision::R,delta,tx,Sy,err
double precision,allocatable::tmp(:),K(:,:)
double precision,parameter::hmin=1d-12,hmax=0.5d0
integer,parameter::s=7
double precision::a(1:s,1:s),b1(1:s),b2(1:s),c(1:s),Rc(1:s),work(1:N)

c(1:7)=(/0d0, 0.2d0, 0.3d0,0.8d0,&
0.888888888888888888888888888888888888888889d0,&
1.d0, 1.d0/)
a(1:7,1:7)=0d0
a(1,1:7)=(/0d0, 0d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(2,1:7)=(/0.2d0, 0d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(3,1:7)=(/0.075d0, 0.225d0, 0d0, 0d0, 0d0, 0d0, 0d0/)
a(4,1:7)=(/0.977777777777777777777777777777777777778d0, &
-3.733333333333333333333333333333333333333d0, &
3.555555555555555555555555555555555555556d0, &
0.d0, 0d0, 0d0, 0d0/)
a(5,1:7)=(/2.95259868922420362749580856576741350403902d0, &
-11.59579332418838591678097850937357110196616d0, &
9.8228928516994360615759792714525224813290657d0, &
-0.29080932784636488340192043895747599451303155d0, &
0.d0, 0d0, 0d0/)
a(6,1:7)=(/2.846275252525252525252525252525252525252525d0, &
-10.757575757575757575757575757575757575757576d0, &
8.906422717743472460453592529064227177434725d0, &
! ^ Here, mistake compare with
! http://depa.fquim.unam.mx/amyd/archivero/DormandPrince_19856.pdf
0.278409090909090909090909090909090909090909d0, &
-0.273531303602058319039451114922813036020583d0, &
0.d0,0d0/)
a(7,1:7)=(/0.09114583333333333333333333333333333333333333d0, &
0.d0, &
0.4492362982929020664869721473495058400718778077d0, &
0.6510416666666666666666666666666666666666666667d0, &
-0.3223761792452830188679245283018867924528301887d0, &
0.13095238095238095238095238095238095238095238095d0, &
0.d0/)
b2(1:7)=(/0.0899131944444444444444444444444444444444444444d0, &
0.d0, &
0.45348906858340820604971548367774782869122491764d0, &
0.6140625d0, &
-0.271512382075471698113207547169811320754716981d0, &
0.089047619047619047619047619047619047619047619d0, &
0.025d0/)

b1(1:7)=a(7,1:7)

Rc(1:7)=(/0.001232638888888888888888888888888888888888888889d0,&
0.d0, &
-0.004252770290506139562743336328241988619347109913d0, &
0.036979166666666666666666666666666666666666666d0, &
-0.050863797169811320754716981132075471698113207547d0, &
0.0419047619047619047619047619047619047619047619d0, &
-0.025d0/)

key=0
allocate(tmp(1:N),K(1:s,1:N))
tmp=0d0; K=0d0

if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
endif

if(h.ge.abs(xbound-x))h=xbound-x

FLAG=1
if(abs(x-xbound).le.hmin)then
info=1
FLAG=0
endif

do while(FLAG.eq.1)
tx=x

if(info.ne.-1)then
K(1,1:N)=h*work0(1:N)
else
do i=1,N
work(i)=grk(N,x,y,i)
K(1,i)=h*work(i)
enddo
work0(1:N)=work(1:N)
info=0
endif

do j=2,s
tx=x+c(j)*h
tmp(1:N)=y(1:N)
do i=1,j-1
tmp(1:N)=tmp(1:N)+K(i,1:N)*a(j,i)
enddo
do i=1,N
work(i)=grk(N,tx,tmp,i)
K(j,i)=h*work(i)
enddo
enddo

!step 4
R=0d0
do i=1,N
R=R+(Rc(1)*K(1,i)+Rc(3)*K(3,i)+Rc(4)*K(4,i) &
+Rc(5)*K(5,i)+Rc(6)*K(6,i)+Rc(7)*K(7,i))**2
enddo
R=abs(sqrt(R)/h/dble(N))

Sy=0d0
do i=1,N
Sy=Sy+(y(i)*y(i))
enddo
Sy=dsqrt(Sy)
if(Sy.ge.1d0)then
err=tol*Sy
else
err=tol
endif

!step 5
if(R.le.err.or.key.eq.1)then
x=x+h
y(1:N)=y(1:N)+b1(1)*K(1,1:N)+b1(3)*K(3,1:N) &
+b1(4)*K(4,1:N)+b1(5)*K(5,1:N)+b1(6)*K(6,1:N)
FLAG=0
work0(1:N)=work(1:N)
endif

!step 6
! Avoid zero deviding.
if(R.ge.1d-20)then
delta=(err/(2d0*R))**0.2d0
else
delta=4d0
endif

!step 7
if(delta.le.0.1d0)then
!function changes dramatically.
h=0.1d0*h
elseif(delta.ge.4d0)then
!function changes loosely.
h=4d0*h
else
!function changes moderately.
h=delta*h
endif

!step 8
if(abs(h).ge.hmax)then
h=sign(1d0,h)*hmax
elseif(abs(h).lt.hmin)then
h=sign(1d0,h)*hmin
key=1
endif

!step 9
if(abs(xbound-x).le.abs(h))then
h=xbound-x
if(abs(h).le.hmin)then
info=1
FLAG=0
endif
end if

if(h.le.0d0.and.xbound-x.ge.0d0)then
info=1
FLAG=0
elseif(h.ge.0d0.and.xbound-x.le.0d0)then
info=1
FLAG=0
endif
enddo

if(key.eq.1)then
write(6,'(A,f10.5,A,f10.5)')"Strange point between ",x-h," and ",x
info=-2
endif

deallocate(tmp,K)
return
end subroutine dDP45

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

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倍かかるので、必要なとき以外使わないようにしましょう。

▼ここクリックでこの場に展開

subroutine qDP45(grk,x,h,N,y,xbound,info,tol,work0)
implicit none
integer,intent(in)::N
real*16,intent(in)::xbound,tol
real*16,intent(inout)::x,h,y(1:N),work0(1:N)
integer,intent(inout)::info
real*16,external::grk

! Runge-Kutta-Dormand-Prince method
!
! if h < hmin, propagate forcibly with warning.
!
!-----------------
!work0(1:N) : Just prepare, used and substituted by
! the coefficients K(7,1:7) to reduce computation cost
!info = -1 (Used for first loop;
! you should set info=0 when you try to use this routine)
! = 0 (Running now)
! = 1 (x reach xbound)
! = -9 (Path the discontinue points)
!-----------------
!

integer::i,j,FLAG,key
real*16::R,delta,tx,Sy,err
real*16,allocatable::tmp(:),K(:,:)
real*16,parameter::hmin=1q-20,hmax=0.5q0
integer,parameter::s=7
real*16::a(1:s,1:s),b1(1:s),b2(1:s),c(1:s),Rc(1:s),work(1:N)

c(1:7)=(/0q0, 0.2q0, 0.3q0,0.8q0,&
0.888888888888888888888888888888888888888889q0,&
1.q0, 1.q0/)
a(1:7,1:7)=0q0
a(1,1:7)=(/0q0, 0q0, 0q0, 0q0, 0q0, 0q0, 0q0/)
a(2,1:7)=(/0.2q0, 0q0, 0q0, 0q0, 0q0, 0q0, 0q0/)
a(3,1:7)=(/0.075q0, 0.225q0, 0q0, 0q0, 0q0, 0q0, 0q0/)
a(4,1:7)=(/0.977777777777777777777777777777777777778q0, &
-3.733333333333333333333333333333333333333q0, &
3.555555555555555555555555555555555555556q0, &
0.q0, 0q0, 0q0, 0q0/)
a(5,1:7)=(/2.95259868922420362749580856576741350403902q0, &
-11.59579332418838591678097850937357110196616q0, &
9.8228928516994360615759792714525224813290657q0, &
-0.29080932784636488340192043895747599451303155q0, &
0.q0, 0q0, 0q0/)
a(6,1:7)=(/2.846275252525252525252525252525252525252525q0, &
-10.757575757575757575757575757575757575757576q0, &
8.906422717743472460453592529064227177434725q0, &
! ^ Here, mistake compare with
! http://depa.fquim.unam.mx/amyd/archivero/DormandPrince_19856.pdf
0.278409090909090909090909090909090909090909q0, &
-0.273531303602058319039451114922813036020583q0, &
0.q0,0q0/)
a(7,1:7)=(/0.09114583333333333333333333333333333333333333q0, &
0.q0, &
0.4492362982929020664869721473495058400718778077q0, &
0.6510416666666666666666666666666666666666666667q0, &
-0.3223761792452830188679245283018867924528301887q0, &
0.13095238095238095238095238095238095238095238095q0, &
0.q0/)
b2(1:7)=(/0.0899131944444444444444444444444444444444444444q0, &
0.q0, &
0.45348906858340820604971548367774782869122491764q0, &
0.6140625q0, &
-0.271512382075471698113207547169811320754716981q0, &
0.089047619047619047619047619047619047619047619q0, &
0.025q0/)

b1(1:7)=a(7,1:7)

Rc(1:7)=(/0.001232638888888888888888888888888888888888888889q0,&
0.q0, &
-0.004252770290506139562743336328241988619347109913q0, &
0.036979166666666666666666666666666666666666666q0, &
-0.050863797169811320754716981132075471698113207547q0, &
0.0419047619047619047619047619047619047619047619q0, &
-0.025q0/)

key=0
allocate(tmp(1:N),K(1:s,1:N))
tmp=0q0; K=0q0

if(abs(h).ge.hmax)then
h=sign(1q0,h)*hmax
endif

if(h.ge.abs(xbound-x))h=xbound-x

FLAG=1
if(abs(x-xbound).le.hmin)then
info=1
FLAG=0
endif

do while(FLAG.eq.1)
tx=x

if(info.ne.-1)then
K(1,1:N)=h*work0(1:N)
else
do i=1,N
work(i)=grk(N,x,y,i)
K(1,i)=h*work(i)
enddo
work0(1:N)=work(1:N)
info=0
endif

do j=2,s
tx=x+c(j)*h
tmp(1:N)=y(1:N)
do i=1,j-1
tmp(1:N)=tmp(1:N)+K(i,1:N)*a(j,i)
enddo
do i=1,N
work(i)=grk(N,tx,tmp,i)
K(j,i)=h*work(i)
enddo
enddo

!step 4
R=0q0
do i=1,N
R=R+(Rc(1)*K(1,i)+Rc(3)*K(3,i)+Rc(4)*K(4,i) &
+Rc(5)*K(5,i)+Rc(6)*K(6,i)+Rc(7)*K(7,i))**2
enddo
R=abs(sqrt(R)/h/dble(N))

Sy=0q0
do i=1,N
Sy=Sy+(y(i)*y(i))
enddo
Sy=sqrt(Sy)
if(Sy.ge.1q0)then
err=tol*Sy
else
err=tol
endif

!step 5
if(R.le.err.or.key.eq.1)then
x=x+h
y(1:N)=y(1:N)+b1(1)*K(1,1:N)+b1(3)*K(3,1:N) &
+b1(4)*K(4,1:N)+b1(5)*K(5,1:N)+b1(6)*K(6,1:N)
FLAG=0
work0(1:N)=work(1:N)
endif

!step 6
! Avoid zero deviding.
if(R.ge.1q-50)then
delta=(err/(2q0*R))**0.2q0
else
delta=4q0
endif

!step 7
if(delta.le.0.1q0)then
!function changes dramatically.
h=0.1q0*h
elseif(delta.ge.4q0)then
!function changes loosely.
h=4q0*h
else
!function changes moderately.
h=delta*h
endif

!step 8
if(abs(h).ge.hmax)then
h=sign(1q0,h)*hmax
elseif(abs(h).lt.hmin)then
h=sign(1q0,h)*hmin
key=1
endif

!step 9
if(abs(xbound-x).le.abs(h))then
h=xbound-x
if(abs(h).le.hmin)then
info=1
FLAG=0
endif
end if

if(h.le.0q0.and.xbound-x.ge.0q0)then
info=1
FLAG=0
elseif(h.ge.0q0.and.xbound-x.le.0q0)then
info=1
FLAG=0
endif
enddo

if(key.eq.1)then
write(6,'(A,f10.5,A,f10.5)')"Strange point between ",x-h," and ",x
info=-9
endif

deallocate(tmp,K)
return
end subroutine qDP45

subroutine qrkf45(grk,x,h,N,y,xbound,info,tol)
! if h < hmin, propagate forcibly with warning.
!
!-----------------
!info = -9 (maybe path the discontinue points)
! = 0 (Running now)
! = 1 (x reach xbound)
!-----------------
!
implicit none
integer,intent(in)::N
real*16,intent(in)::xbound,tol
real*16,intent(inout)::x,h,y(1:N)
integer,intent(inout)::info
real*16,external::grk

integer::i,j,FLAG,key
real*16::R,delta,tx,Sy,err
real*16,allocatable::tmp(:),K(:,:)
real*16,parameter::hmin=1q-20,hmax=0.5q0
integer,parameter::s=6
real*16::a(1:s,1:s),b1(1:s),b2(1:s),c(1:s),Rc(1:s)

c(1:6)=(/0q0, 0.25q0, 0.375q0,&
0.9230769230769230769230769230769230769231q0, 1q0, 0.5q0/)
a(1:6,1:6)=0q0
a(1,1:6)=(/0q0, 0q0, 0q0, 0q0, 0q0, 0q0/)
a(2,1:6)=(/0.25q0, 0q0, 0q0, 0q0, 0q0, 0q0/)
a(3,1:6)=(/0.09375q0, 0.28125q0, 0q0, 0q0, 0q0, 0q0/)
a(4,1:6)=(/0.8793809740555302685480200273099681383705q0, &
-3.277196176604460628129267182521620391443q0, &
3.320892125625853436504324078288575329995q0, 0q0, 0q0, 0q0/)
a(5,1:6)=(/2.032407407407407407407407407407407407407q0,-8q0, &
7.173489278752436647173489278752436647173q0, &
-0.2058966861598440545808966861598440545809q0, 0q0, 0q0/)
a(6,1:6)=(/-0.2962962962962962962962962962962962962963q0,2q0, &
-1.381676413255360623781676413255360623782q0, &
0.4529727095516569200779727095516569200780q0,-0.275q0,0q0/)
b2(1:6)=(/0.1185185185185185185185185185185185185185q0, 0.q0,&
0.5189863547758284600389863547758284600390q0, &
0.5061314903420166578061314903420166578061q0, &
-0.18q0, 0.03636363636363636363636363636363636363636q0/)
b1(1:6)=(/0.1157407407407407407407407407407407407407q0, 0q0,&
0.5489278752436647173489278752436647173489q0, &
0.5353313840155945419103313840155945419103q0, -0.2q0, 0q0/)
Rc(1:6)=(/0.002777777777777777777777777777777777777778q0,0q0, &
-0.02994152046783625730994152046783625730994q0, &
-0.02919989367357788410419989367357788410420q0, 0.02q0, &
0.03636363636363636363636363636363636363636q0/)

if(info.eq.-1)info=0

key=0
allocate(tmp(1:N),K(1:s,1:N))
tmp=0q0; K=0q0

if(abs(h).ge.hmax)then
h=sign(1q0,h)*hmax
endif

if(h.ge.abs(xbound-x))h=xbound-x

FLAG=1
if(abs(x-xbound).le.hmin)then
info=1
FLAG=0
endif

do while(FLAG.eq.1)

tx=x
do j=1,s
tx=x+c(j)*h
tmp(1:N)=y(1:N)
do i=1,j-1
tmp(1:N)=tmp(1:N)+K(i,1:N)*a(j,i)
enddo
do i=1,N
K(j,i)=h*grk(N,tx,tmp,i)
enddo
enddo

!step 4
R=0q0
do i=1,N
R=R+(Rc(1)*K(1,i)+Rc(3)*K(3,i)+Rc(4)*K(4,i)+Rc(5)*K(5,i)+Rc(6)*K(6,i))**2q0
enddo
R=abs(sqrt(R)/h/(N*1q0))

Sy=0q0
do i=1,N
Sy=Sy+(y(i)*y(i))
enddo
Sy=sqrt(Sy)
if(Sy.ge.1q0)then
err=tol*Sy
else
err=tol
endif

!step 5
if(R.le.err.or.key.eq.1)then
x=x+h
y(1:N)=y(1:N)+b1(1)*K(1,1:N)+b1(3)*K(3,1:N)+b1(4)*K(4,1:N)+b1(5)*K(5,1:N)
FLAG=0
endif

!step 6
! Avoid zero deviding.
if(R.ge.1q-50)then
delta=(err/(2q0*R))**0.25q0
else
delta=4q0
endif

!step 7
if(delta.le.0.1q0)then
!function changes dramatically.
h=0.1q0*h
elseif(delta.ge.4q0)then
!function changes loosely.
h=4q0*h
else
!function changes moderately.
h=delta*h
endif

!step 8
if(abs(h).ge.hmax)then
h=sign(1q0,h)*hmax
elseif(abs(h).lt.hmin)then
h=sign(1q0,h)*hmin
key=1
endif

!step 9
if(abs(xbound-x).le.abs(h))then
h=xbound-x
if(abs(h).le.hmin)then
info=1
FLAG=0
endif
end if

if(h.le.0q0.and.xbound-x.ge.0q0)then
info=1
FLAG=0
elseif(h.ge.0q0.and.xbound-x.le.0q0)then
info=1
FLAG=0
endif
enddo

if(key.eq.1)then
write(6,'(A,f10.5,A,f10.5)')"Strange point between ",x-h," and ",x
info=-9
endif

deallocate(tmp,K)
return
end subroutine qrkf45

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

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

ニューメリカルレシピ・イン・シー日本語版―C言語による数値計算のレシピ

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

物理学

運動量演算子の導出

2015年4月13日 sikino コメントする

量子力学の運動量演算子$\hat{p}$が
$
\displaystyle \hat{p}=-i\hbar \frac{d}{d x}
$
であらわされることを導出します。
ここでの”導出”は運動量演算子として妥当なものを導出する、と言った方がいいかもしれません。

平行移動演算子$\hat{T}$を出発点とします。
この平行移動演算子はある波動関数$\psi_{(x)}$を平行移動させる演算子で、これを作用させると
$
\hat{T}_{(\Delta x)}\psi_{(x)}=\psi_{(x+\Delta x)}
$
となる演算子です。

よく説明される、シュレディンガー方程式の平面波解からの出発はしません。

この$\hat{T}$の持つであろう性質を考え、その特性から$\hat{T}$の具体的な形を推定し、運動量演算子を導きます。

平行移動演算子$\hat{T}$が満たすべき性質

当たり前と思われる性質から、平行移動演算子に関して以下の4つが言えます。

波動関数$\psi_{(x)}$が1に規格化済みであれば、平行移動した波動関数$\psi_{(x+\Delta x)}$もまた規格化されているはずである。
すなわち、
$
\begin{align}
\langle\psi_{(x)}|\psi_{(x)}\rangle &= \langle\psi_{(x)}|\hat{T}^{\dagger}_{(\Delta x)}\hat{T}_{(\Delta x)}|\psi_{(x)}\rangle \\
&= \langle\psi_{(x+\Delta x)}|\psi_{(x+\Delta x)}\rangle \\
&\rightarrow \hat{T}^{\dagger}_{(\Delta x)}\hat{T}_{(\Delta x)}=\mathbf{1}
\end{align}
$
となり、3本目の式から$\hat{T}$はユニタリー演算子であることがわかります。
$\Delta x$進めた後、$\Delta x^{\prime}$を進める操作は$\Delta x + \Delta x^{\prime}$進める操作に等しいはずである。
すなわち、
$
\hat{T}_{(\Delta x^{\prime})}\hat{T}_{(\Delta x)}=\hat{T}_{(\Delta x+\Delta x^{\prime})}
$
負の方向への平行移動$\hat{T}_{(-\Delta x)}$は正の方向への平行移動$\hat{T}_{(\Delta x)}$の逆のはずである。
$
\hat{T}_{(-\Delta x)}=\hat{T}^{-1}_{(\Delta x)}
$
$\Delta x$が0のとき、恒等操作なはずである。
$
\begin{align}
\hat{T}_{(0)}\psi_{(x)}&=\psi_{(x)} \\
&\rightarrow \hat{T}_{(0)}=\mathbf{1}
\end{align}
$

平行移動演算子の具体的な形

1,より、$\hat{T}_{(\Delta x)}$はユニタリー演算子である事が分かりました。ユニタリー演算子を作る1つの方法として、エルミート演算子$\hat{A}_{(\Delta x)}$を用いて、
$
\hat{T}_{(\Delta x)}=e^{i\hat{A}_{(\Delta x)}}
$
と置けば、$\hat{T}_{(\Delta x)}$はユニタリー演算子となります(エルミート演算子とユニタリー演算子の性質)。
この仮定の下、
$
\hat{T}_{(\Delta x)}\psi_{(x)}=\psi_{(x+\Delta x)}
$
に代入して、
$
e^{i\hat{A}_{(\Delta x)}}\psi_{(x)}=\psi_{(x+\Delta x)}
$
無限小移動を考えて$\Delta_{(x)}\rightarrow 0$と約束すれば、
$
(1+i\hat{A}_{(\Delta x)})\psi_{(x)}=\psi_{(x+\Delta x)}
$
変形して、
$
\psi_{(x+\Delta x)}-\psi_{(x)}=i\hat{A}_{(\Delta x)}\psi_{(x)}
$
辺々を$\Delta x$で割れば、
$
\begin{align}
\displaystyle \lim_{\Delta x\rightarrow 0}\frac{\psi_{(x+\Delta x)}-\psi_{(x)}}{\Delta x} = \lim_{\Delta x\rightarrow 0} i\frac{\hat{A}_{(\Delta x)}}{\Delta x}\psi_{(x)}
\end{align}…(i)
$
となります。左辺は位置$x$における波動関数の傾きです。
ここで右辺を見ます。$\hat{A}_{(\Delta x)}$は$\Delta x$のなんらかの関数であるはずですが、$\hat{A}_{(\Delta x)}$が分母の$\Delta x$と打ち消す形でなければ0、もしくは発散する事になります。
この場合、波動関数の傾きが到る点で0になるということは、全空間で定数ということになり非物理的な解です。また、発散してしまうのであれば波動関数が至る所で発散してしまい、これまた日物理的です。

具体的に、もしも$\hat{A}_{(\Delta x)}=(\Delta x)^2$だったら、
$
\displaystyle \lim_{\Delta x\rightarrow 0} \frac{(\Delta x)^2}{\Delta x}=\lim_{\Delta x\rightarrow 0} \Delta x =0
$
であるし、もしも$\hat{A}_{(\Delta x)}=(\Delta x)^{1/2}$だったら
$
\displaystyle \lim_{\Delta x\rightarrow 0} \frac{(\Delta x)^{1/2}}{\Delta x}=\lim_{\Delta x\rightarrow 0} \frac{1}{\sqrt{\Delta x}} =\infty
$
となります。
すなわち、言いたいことは、$\hat{A}_{(\Delta x)}$は$\Delta x$に関して1次であるべきで、
$
\hat{A}_{(\Delta x)}=\hat{B}\cdot \Delta x
$…(ii)
と書けるはずです。ここで$\hat{B}$は$\Delta x$に依らない演算子です(下に補足説明あり)。

よって、式(i)は$\Delta x\rightarrow 0$のとき、
$
\begin{align}
\displaystyle \lim_{\Delta x\rightarrow 0}\frac{\psi_{(x+\Delta x)}-\psi_{(x)}}{\Delta x} &= \lim_{\Delta x\rightarrow 0} i\frac{\hat{A}_{(\Delta x)}}{\Delta x}\psi_{(x)} \\
\frac{d \psi_{(x)}}{dx}=i\hat{B}\psi_{(x)}
\end{align}…(iii)
$
と書けるため、
$
\psi_{(x)}=e^{i\hat{B}\cdot x}\psi_{(0)}
$
となります。(この式で右辺の$\psi_{(0)}$が0である必要はなく、定数であればいいんです。)

平行移動演算子と運動量演算子との関係

次元について考えましょう！
今、$\hat{T}_{(\Delta x)}=e^{i\hat{B}\cdot \Delta x}$であり、$\Delta x$の次元は$\mathrm{[m]}$であるため、$\hat{B}$の次元は$\mathrm{[m^{-1}]}$でなければなりません。よって、$\hat{B}$は波数の次元を持つ演算子であると考えられるため、$\hat{B}$を$c\hat{k}$と記述することにします。ここで、$c$は無次元の定数です。
式(iii)より、
$
\displaystyle \frac{d \psi_{(x)}}{dx}=ic\hat{k}\psi_{(x)}
$
辺々に$-i\hbar$を掛けると、
$
\displaystyle -i\hbar \frac{d}{dx} \psi_{(x)}=c\hbar \hat{k}\psi_{(x)}
$…(iv)

!!!

$p=\hbar k$は運動量を表すものでした。よって,$\hat{p}=\hbar \hat{k}$と書けば、式(iv)より
$
\begin{align}
\hat{p}\psi_{(x)}&=\frac{1}{c}\left[-i\hbar\frac{d}{d x}\right] \psi_{(x)} \\
\end{align}
$
となります。よって運動量演算子は
$
\begin{align}
\hat{p}\propto -i\hbar\frac{d}{d x} \\
\end{align}
$
となります。この考えからは定数倍を含めてしまうので、比例する、という事しか言えません。

なぜ決まらないのか考えてみますと、時間の依存性に関する考え方が式のどこにも含んでいないからです。
時間の単位が決まらないと速度の定義が出来ないわけで、これ以上進めることが出来ないのです。

各成分について言えるので、多次元の場合では
$
\hat{p}_x\propto -i\hbar\frac{\partial}{\partial x}
$
としても同じことです。
申し訳ないですが、数学的に厳密か？は保証しません。

補足説明

果たして運動量演算子として考えた$e^{i\hat{B}\cdot \Delta x}$は$\hat{T}$の特性1.~4.を満たすでしょうか？

$
\left(e^{i\hat{B}\cdot \Delta x}\right)^{\dagger} \left(e^{i\hat{B}\cdot \Delta x}\right)=e^{-i\hat{B}^{\dagger}\cdot \Delta x}e^{i\hat{B}\cdot \Delta x}=\mathbf{1}
$
ここで$\hat{B}$はエルミート演算子を考えているので、$\hat{B}^{\dagger}=\hat{B}$です。
$
e^{-i\hat{B}^{\dagger}\cdot \Delta x}e^{i\hat{B}\cdot \Delta x}=\mathbf{1}
$
ok!
$
\hat{T}_{(\Delta x^{\prime})}\hat{T}_{(\Delta x)}=\hat{T}_{(\Delta x+\Delta x^{\prime})}
$か？
$
(\mbox{左辺})=e^{i\hat{B}\cdot \Delta x^{\prime}} e^{i\hat{B}\cdot \Delta x}=e^{i\hat{B}\cdot (\Delta x^{\prime}+\Delta x)}
$
ok!
$
\hat{T}_{(-\Delta x)}=\hat{T}^{-1}_{(\Delta x)}
$か？
両辺の左から$\hat{T}_(\Delta x)$を掛けて、
$
\begin{align}
&\hat{T}_{(\Delta x)}\hat{T}_{(-\Delta x)}=\mathbf{1} \\
&= e^{i\hat{B}\cdot \Delta x}e^{i\hat{B}\cdot (-\Delta x)}=e^{i\hat{B}\cdot 0}=\mathbf{1}
\end{align}
$
ok!
$e^{i\hat{B}\cdot 0}=1$
ok!

よって、4つすべてを満たすので、$\hat{T}_{(\Delta x)}=e^{i\hat{B}\cdot \Delta x}$は平行移動演算子として適当であると考えることができるのです。

プログラミングと数値計算

数値積分(等間隔)

2015年4月7日 sikino コメントする

数値計算での積分方法、特に等間隔の分点の場合であるニュートン・コーツ積分(とロンバーグ積分)に関する理論とのプログラムを載せます。

ニュートンコーツ型の公式
ルンゲ現象
ロンバーグ積分
サブルーチン”integral”のコードと使い方
- サブルーチン”integral”を用いた例題
サブルーチン”romberg”のコードと使い方
- サブルーチン”romberg”を用いた例題

等間隔の分点で積分を行う場合、よく使われる方法は台形則、もしくはシンプソン積分です。

台形則

台形則は、分点間を台形近似して面積を求める方法であり、下図のようなイメージで積分を近似する方法です。

積分$\displaystyle \int_a^b f(x) dx$を、
$
\displaystyle \int_a^b f(x) dx \sim \sum_{i=0}^{N-1} \frac{f(x_{i+1})+f(x_i)}{2}h
$
として近似します。この時、誤差のオーダーは$O(h^3)$となります。
言い換えると、
台形積分とは、隣り合う分点間を1次関数で近似して求める積分
言えます。

シンプソン積分

シンプソン積分は、隣り合う分点間を2次関数で近似して求める積分です。イメージは

な感じです。数式では
$
\displaystyle \int_a^b f(x) dx \sim \sum_{i=0}^{N/2} \frac{f(x_{2i})+4f(x_{2i+1})+f(x_{2i+2})}{3}h
$
となります。

高次≠高精度

これ以降もあります。分点間を3次、4次、…としていけば高次の積分をすることが可能になります。
この分点間の次数をどんどん上げていって積分をする方法をまとめてニュートンコーツの公式、と呼びます。
例えばの台形則はニュートンコーツの公式の1次に相当し、シンプソン則はニュートンコーツの公式の2次に相当しています。

じゃあ100次の公式作れば精度が凄い公式ができるんじゃないか？となりますが、これは違います。
高次であることと精度が高いことは違うのです。だから使うな、というわけではなくて知っておいてほしいことなんです。
ではなぜか？これはwikipediaのニュートンコーツの公式のページからとってきたものですが、この式に理由があります。

問題となるのが誤差項です。例えば台形近似を見ますと、誤差項は$h^3 f^{(2)}(\xi)$とあります。hは刻み幅なので、この項は刻み幅を小さくすれば小さくするほど精度がよくなることを言っています。
が、問題は$f^{(2)}(\xi)$の方です。この項は本当の関数が分からないと評価できません。
多項式だったら微分するごとに次数が1つづつ減少するのでこの項は高次になれば小さくなります。じゃあもしも微分するごとに値が増加していったら…？この場合、hを小さくしてもそれを上回るくらい$f^{(2)}(\xi)$が大きくなることがあります。これをルンゲ現象と呼びます。

ルンゲ現象

例えば関数
$
\displaystyle f(x)=\frac{1}{1+25x^2}
$
を考えます。wikipediaのページを参考にすれば、その導関数と値は、

となります。高次になればなるほど値が増加する関数なのです。
実際には$h^n f^{(m)}(\xi)$の積の値が小さいか大きいか？なので、一概には高次は悪い！とは言えません。高次は危ないくらいでしょう。

これを可能な限り回避するために、低次の積分を組み合わせて使う方法がとられます。
こうすることで高次の値を増加を抑えるのです。

そのほかの手法もいろいろあります。分点の位置に縛られなければガウス求積法が最高でしょう。素晴らしい方法です。ぜひ調べてみてください。

[adsense1]

ロンバーグ積分

さて、高次を追い求めるのもいいですが、一線を駕す…とまではいきませんが、ロンバーグ(Romberg)積分というものがあります。
このロンバーグ積分とは、台形則による計算を基本とし、無限小区間での積分結果を補外によって求めよう、というものです。
ここでいう補外の概念は以下のような感じです。

刻み幅1で台形則の結果から積分値S1が分かる
↓
刻み幅0.5で台形則による結果から積分値S2がわかる
↓
刻み幅0.25で台形則による結果から積分値S3がわかる
↓
刻み幅0.175で台形則による結果から積分値S4がわかる

とある程度計算していきます。そうすると刻み幅hを変化させるとそれに応じて積分値がS1→S2→S3→S4…と推移していくわけです。
この推移の仕方を計算すれば、刻み幅無限小の場合の積分値Sが分かる、こういう仕組みになっています。収束加速法と言ったり、リチャードソン補外、なんて言ったりします。
この補外による方法はかなり優れています。先ほどの高次がどうこう、といった問題が発生しないのです。台形則しか使ってないんですから。有限桁計算に適した方法である、とかwikipediaに載ってたりします。

ロンバーグ積分の詳しい計算方法や理論は
ニューメリカルレシピ・イン・シー日本語版―C言語による数値計算のレシピ

ニューメリカルレシピ・イン・シー日本語版―C言語による数値計算のレシピ

を見ると詳しく書かれています。

日本語では第6章数値積分と数値微分にロンバーグ積分の記述があります。
英語ではRomberg integralで調べると良いでしょう。参考までに、英語では
Romberg integration exampleや、Romberg Integrationなにかがいいと思います。

2016/03/01
上のRomberg Integrationを参考にして書いたプログラムはこちらです。これは、解析関数fをRomberg積分するプログラムです。
解析関数ではなく、サイズ$2^{n}+1$の配列に格納されている場合はこのページの下の方にあるプログラムを参照してください。

call romberg(a,b,s,pre)

で呼び出し、区間a~bの積分$s=\int_a^b f(x) dx$を精度preで求めるプログラムです。

▼クリックでこの場に展開されます

サブルーチン”integral”のコードと使い方

下のサブルーチン”integral”は1,2,3次元の短冊近似、台形則、シンプソン積分、シンプソン3/8積分、ブール則(ニュートン・コーツ型積分)をカバーします。
ただし、分点の個数に注意してください。分点の個数を気にせず使える近似方法は、短冊近似、台形近似のみです。

ロンバーグ積分では引数の関係上、ニュートン・コーツ型積分との間をうまく処理できなかったため、別のサブルーチンとして書きます。もしもロンバーグ積分を使いたい場合はもう少し先に進んでください。

まずはニュートン・コーツ型積分です。
下のサブルーチンを使ってください。
使い方は、

call integral(配列の大きさ,配列,刻み幅,積分結果,積分手法)

という感じです。
実際のプログラムでは、以下のように使用してください。配列yは複素数解列(ただし実軸上の積分)でもokです。

1次元

サブルーチンintegralの使い方
積分則	呼び方
短冊近似(長方形近似)	call integral(size(w,1),w,h,s,"box")	配列yの大きさに指定はない。
台形近似	call integral(size(y,1),y,h,s,"trapezoid")	配列yの大きさに指定はない。
シンプソン則	call integral(size(y,1),y,h,s,"simpson")	配列yの大きさが3,5,7,…個、$2n+1$個でないといけない
シンプソン3/8則	call integral(size(y,1),y,h,s,"simpson38")	配列yの大きさが4,7,10,…個、$3n+1$個でないといけない
ブール則	call integral(size(y,1),y,h,s,"boole")	配列yの大きさが5,9,13,…個、$4n+1$個でないといけない

2次元配列の場合は以下のように指定してください。

call integral(size(w,1),size(w,2),w,hx,hy,s,"simpson")

3次元の場合はこう指定してください。

call integral(size(w,1),size(w,2),size(w,3),w,hx,hy,hz,s,"simpson")

積分のためのモジュール

module integral_mod
!developer --> sikino
!date --> 2015/04/07
implicit none
interface integral
module procedure &
dintegral, &
dintegral2d, &
dintegral3d, &
cintegral, &
cintegral2d, &
cintegral3d
end interface integral
contains
subroutine dintegral(N,y,h,s,method)
integer,intent(in)::N
double precision,intent(in)::h,y(1:N)
character(*),intent(in)::method
double precision,intent(out)::s
integer::i
double precision::y0,y1,y2,y3

s=0.d0; y0=0.d0; y1=0.d0; y2=0.d0; y3=0.d0

if(trim(method).eq."box")then
s=h*sum(y(1:N-1))
elseif(trim(method).eq."trapezoid")then
y1=y(1)+y(N)
do i=2,N-1
y2=y2+y(i)
enddo
s=(y1+2.d0*y2)*h*0.5d0
elseif(trim(method).eq."simpson")then
if(mod(N,2).ne.1)then
write(6,*)"=====cannot calculation with simpson"
write(6,*)"=====program stop"
stop
endif

y1=y(1)+y(N)
do i=2,N-1,2
y2=y2+y(i)
enddo
do i=3,N-2,2
y3=y3+y(i)
enddo

s=(y1+4.d0*y2+2.d0*y3)*h/3.d0

elseif(trim(method).eq."simpson38")then

if(mod(N,3).ne.1)then
write(6,*)"=====cannot calculation with simpson38"
write(6,*)"=====program stop"
stop
endif

y0=y(1)+y(N)
do i=2,N-2,3
y1=y1+y(i)
enddo
do i=3,N-1,3
y2=y2+y(i)
enddo
do i=4,N-3,3
y3=y3+y(i)
enddo
s=(y0+3.d0*(y1+y2)+2.d0*y3)*3.d0*h/8.d0

elseif(trim(method).eq."boole")then
if(mod(N,4).ne.1)then
write(6,*)"=====cannot calculation with boole"
write(6,*)"=====program stop"
stop
endif

y0=y(1)+y(N)
do i=5,N-4,4
y1=y1+y(i)
enddo
do i=2,N-1,2
y2=y2+y(i)
enddo
do i=3,N-2,4
y3=y3+y(i)
enddo

s=(14.d0*y0+28.d0*y1+64.d0*y2+24.d0*y3)*h/45.d0
else
write(6,*)"=====cannot calculation in integral"
write(6,*)"=====program stop"
stop
end if

return
end subroutine dintegral

subroutine dintegral2d(Nx,Ny,z,hx,hy,s,method)
implicit none
integer,intent(in)::Nx,Ny
double precision,intent(in)::hx,hy,z(1:Nx,1:Ny)
character(*),intent(in)::method
double precision,intent(out)::s
integer::i
double precision::ty(1:Ny),r(1:Nx)

s=0.d0
ty(1:Ny)=0.d0
r(1:Nx)=0.d0
do i=1,Nx
ty(1:Ny)=z(i,1:Ny)
call integral(Ny,ty,hy,s,method)
r(i)=s
enddo
call integral(Nx,r,hx,s,method)

return
end subroutine dintegral2d

subroutine dintegral3d(Nx,Ny,Nz,w,hx,hy,hz,s,method)
implicit none
integer,intent(in)::Nx,Ny,Nz
double precision,intent(in)::hx,hy,hz,w(1:Nx,1:Ny,1:Nz)
character(*),intent(in)::method
double precision,intent(out)::s
integer::i
double precision::tyz(1:Ny,1:Nz),r(1:Nx)

s=0.d0
tyz(1:Ny,1:Nz)=0.d0
r(1:Nx)=0.d0
do i=1,Nx
tyz(1:Ny,1:Nz)=w(i,1:Ny,1:Nz)
call integral(Ny,Nz,tyz,hy,hz,s,method)
r(i)=s
enddo
call integral(Nx,r,hx,s,method)

return
end subroutine dintegral3d

subroutine cintegral(N,y,h,s,method)
integer,intent(in)::N
complex(kind(0d0)),intent(in)::y(1:N)
double precision,intent(in)::h
character(*),intent(in)::method
complex(kind(0d0)),intent(out)::s

double precision::res,ims

s=dcmplx(0d0,0d0); res=0.d0; ims=0.d0

call integral(N,dble(y),h,res,trim(method))
call integral(N,dimag(y),h,ims,trim(method))

s=dcmplx(res,ims)

return
end subroutine cintegral

subroutine cintegral2d(Nx,Ny,z,hx,hy,s,method)
integer,intent(in)::Nx,Ny
complex(kind(0d0)),intent(in)::z(1:Nx,1:Ny)
double precision,intent(in)::hx,hy
character(*),intent(in)::method
complex(kind(0d0)),intent(out)::s

double precision::res,ims

s=dcmplx(0d0,0d0); res=0.d0; ims=0.d0

call integral(Nx,Ny,dble(z),hx,hy,res,trim(method))
call integral(Nx,Ny,dimag(z),hx,hy,ims,trim(method))

s=dcmplx(res,ims)

return
end subroutine cintegral2d

subroutine cintegral3d(Nx,Ny,Nz,w,hx,hy,hz,s,method)
integer,intent(in)::Nx,Ny,Nz
complex(kind(0d0)),intent(in)::w(1:Nx,1:Ny,1:Nz)
double precision,intent(in)::hx,hy,hz
character(*),intent(in)::method
complex(kind(0d0)),intent(out)::s

double precision::res,ims

s=dcmplx(0d0,0d0); res=0.d0; ims=0.d0

call integral(Nx,Ny,Nz,dble(w),hx,hy,hz,res,trim(method))
call integral(Nx,Ny,Nz,dimag(w),hx,hy,hz,ims,trim(method))

s=dcmplx(res,ims)

return
end subroutine cintegral3d
end module integral_mod

サブルーチン”integral”を用いた例題

実際の使い方。参考にどうぞ。
下のコードは2次元の積分
$
\displaystyle \int\int xe^{-x^2}e^{-y}dxdy=\frac{e^{-x^2}}{2}e^{-y}+\mbox{const}
$
を数値積分します。積分範囲を$x=0\sim 3, \; y=0 \sim 5$にした場合、その解析解は
$
\begin{align}
\int_0^3 dx\int_0^5 dy xe^{-x^2}e^{-y}&=\frac{1}{2}(1-e^{-5}-e^{-9}+e^{-14})\\
&=0.4965697373627734784608751894356320535936864348993604\cdots
\end{align}
$
となります。実行すると

./a.out
simpson : 0.496569737366759E+000
romberg,3 : 0.496569737362773E+000

となります。コードは、

program main
use integral_mod
implicit none
integer::i,j,n,Nx,Ny
double precision::hx,hy,ans
double precision::xmin,xmax,ymin,ymax
double precision,allocatable::x(:),y(:)
double precision::s
double precision,allocatable::w(:,:)
double precision::f
external::f

n=10

xmin=0d0
xmax=3d0

ymin=0d0
ymax=5d0

Nx=2**(n)
Ny=2**(n)
allocate(x(0:Nx),y(0:Ny),w(0:Nx,0:Ny))
x=0d0; y=0d0; w=dcmplx(0d0,0d0)

hx=(xmax-xmin)/dble(Nx)
hy=(ymax-ymin)/dble(Ny)

do i=0,Nx
x(i)=xmin+hx*dble(i)
enddo
do i=0,Ny
y(i)=ymin+hy*dble(i)
enddo

do i=0,Nx
do j=0,Ny
w(i,j)=f(x(i),y(j))
end do
end do

call integral(size(w,1),size(w,2),w,hx,hy,s,"simpson")
write(6,'(A,e25.15e3)')"simpson : ",s

call integral(size(w,1),size(w,2),w,hx,hy,s,"romberg",3)
write(6,'(A,e25.15e3)')"romberg,3 : ",s

return
end program main

function f(x,y)
implicit none
double precision,intent(in)::x,y
double precision::f

f=x*exp(-x*x)*exp(-y)

return
end function f

3次元の計算の例題です。
下のコードは
$
\displaystyle \int\int\int x^4y^4z^4 dxdydz=\left(\frac{1}{5}\{0.7^5-(-1)^5\}\right)^3
$
を数値積分します。積分範囲を$x,y,z=-1\sim 0.7$にした場合、その値は
$
\displaystyle \int_{-1}^{0.7}\int_{-1}^{0.7}\int_{-1}^{0.7} x^4y^4z^4 dxdydz = 0.127496010896795\cdots
$
となります。実行すると、

$ gfortran integralmod.f90 main.f90
$ ./a.out
analysis : 0.127496010896795E-01
simpson : 0.10507613E-006 0.10507613E-006
romberg,3 : -0.17347235E-017 -0.17347235E-017

と得られます。中身はこうです。分点の数はx,y,z軸各々で違っていて構いません。

program main
use integral_mod
implicit none
integer::i,j,k,n,Nx,Ny,Nz
double precision::hx,hy,hz,ans
double precision,allocatable::x(:),y(:),z(:)
complex(kind(0d0))::s
complex(kind(0d0)),allocatable::w(:,:,:)

n=6
Nx=2**(n)
Ny=2**(n-1)
Nz=2**(n+1)

allocate(x(0:Nx),y(0:Ny),z(0:Nz),w(0:Nx,0:Ny,0:Nz))
x=0d0; y=0d0; z=0d0; w=dcmplx(0d0,0d0)

hx=1.7d0/dble(size(w,1)-1)
hy=1.7d0/dble(size(w,2)-1)
hz=1.7d0/dble(size(w,3)-1)

do i=0,Nx
x(i)=hx*dble(i)-1d0
enddo
do i=0,Ny
y(i)=hy*dble(i)-1d0
enddo
do i=0,Nz
z(i)=hz*dble(i)-1d0
enddo

do i=0,Nx
do j=0,Ny
do k=0,Nz
w(i,j,k)=dcmplx((x(i)*y(j)*z(k))**4,(x(i)*y(j)*z(k))**4)
end do
end do
end do

ans=(0.7d0**5.d0+1.d0)/5.d0
ans=ans**3.d0
write(6,'(A,e23.15e2)')"analysis : ",ans

call integral(size(w,1),size(w,2),size(w,3),w,hx,hy,hz,s,"simpson")
write(6,'(A,2e17.8e3)')"simpson : ",dble(s)-ans,dimag(s)-ans

call integral(size(w,1),size(w,2),size(w,3),w,hx,hy,hz,s,"romberg",3)
write(6,'(A,2e17.8e3)')"romberg,3 : ",dble(s)-ans,dimag(s)-ans

return
end program main

ロンバーグ積分のプログラム

ロンバーグ積分

call romberg(jx,x,y,s,pre)

配列yの大きさは2,3,5,10,…,$2^{jx}+1$個でないといけない

ロンバーグ積分では収束精度”pre”を指定すること。もしも収束精度に達しなかった場合、警告と共に、与えられた値での収束限界の積分結果を返す。積分精度を高めたければ刻み幅を小さくする操作(例えば分点数を増やす等)をしてください。”pre”は、積分結果の絶対値が1より大きくなる場合は相対誤差を、小さくなる場合は絶対誤差を取ります。この理由は以前の結果の補正を加え、収束をさせるためであり、非常に小さい誤差の場合いつまでたっても機械誤差のため収束しないからです。

module romberg_mod
!developer --> sikino
!date --> 2016/03/01
! 2016/03/03
!
! 1D case :
! romberg(jx,x,y,s,pre)
! | | | | +- (in) precision (e.g. 1d-8)
! | | | +--- (out) integration result
! | | +----- (in) y(1:2**jx+1) f(x)
! | +------- (in) x(1:2**jx+1) x
! +---------- (in) related to array size
!
! 2D case :
! romberg(jx,jy,x,y,z,s,pre)
! | | | | | | +-- (in) precision (e.g. 1d-8)
! | | | | | +---- (out) integration result
! | | | | +------ (in) z(1:2**jx+1,1:2**jy+1) f(x,y)
! | | | +-------- (in) y(1:2**jy+1) y
! | | +---------- (in) x(1:2**jx+1) x
! | +------------- (in) related to array size y
! +---------------- (in) related to array size x
!
! 3D case :
! romberg(jx,jy,jz,x,y,z,w,s,pre)
! | | | | | | | | +-- (in) precision (e.g. 1d-8)
! | | | | | | | +---- (out) integration result
! | | | | | | +------ (in) w(1:2**jx+1,1:2**jy+1,1:2**jz+1) f(x,y,z)
! | | | | | +-------- (in) z(1:2**jz+1) z
! | | | | +---------- (in) y(1:2**jy+1) y
! | | | +------------ (in) x(1:2**jx+1) x
! | | +--------------- (in) related to array size z
! | +------------------ (in) related to array size y
! +--------------------- (in) related to array size x
!
implicit none
interface romberg
module procedure &
dromberg, &
dromberg2d, &
dromberg3d
end interface romberg
contains
subroutine dromberg(jx,x,y,s,pre)
implicit none
integer,intent(in)::jx
double precision,intent(in)::x(1:2**jx+1),y(1:2**jx+1),pre
double precision,intent(out)::s

integer,parameter::jm=6 !--> precision: O(h^(2*jm))
integer,parameter::nm=2**jm
double precision,allocatable::tx(:),ty(:)
double precision::ts
integer::k

s=0d0
if(jx.ge.jm)then
allocate(tx(1:nm+1),ty(1:nm+1)); tx=0d0; ty=0d0
do k=0,2**(jx-jm)-1
ts=0d0
tx(1:nm+1)=x(k*nm+1:(k+1)*nm+1)
ty(1:nm+1)=y(k*nm+1:(k+1)*nm+1)

call romberg_sub(jm,tx,ty,ts,pre)
s=s+ts
enddo
deallocate(tx,ty)
else
call romberg_sub(jx,x,y,s,pre)
endif

return
end subroutine dromberg

subroutine romberg_sub(jx,x,y,s,pre)
! reference "http://www.math.usm.edu/lambers/mat460/fall09/lecture29.pdf"
implicit none
integer,intent(in)::jx
double precision,intent(in)::x(1:2**jx+1),y(1:2**jx+1),pre
double precision,intent(out)::s

integer::i,j,k,n,dn
double precision::h,ps,tmp
double precision::T(1:jx+1,1:jx+1)

n=2**jx+1

h=x(n)-x(1)
dn=(n-1)/2

T(1,1)=0.5d0*h*(y(1)+y(n))
s=T(1,1)
ps=s
h=0.5d0*h
do j=2,jx+1

! trapezoidal rule
tmp=0d0
do i=1,2**(j-2)
tmp=tmp+y(1+(2*i-1)*(dn))
enddo
T(j,1)=0.5d0*T(j-1,1)+h*tmp

do k=2,j
! Richardson extrapolation
T(j,k)=T(j,k-1)+(T(j,k-1)-T(j-1,k-1))/(dble(4**(k-1))-1d0)
enddo
s=T(j,j)

! precision check
if(abs(s).ge.1d0)then
if(abs((ps-s)/s).le.pre)exit
else
if(abs((ps-s)).le.pre)exit
endif
ps=s
h=0.5d0*h
dn=dn/2
enddo

if(j-1.eq.jx)then
write(6,'(A)')" -+-+- didn't converge at romberg integral -+-+- "
write(6,'(A)')" Please change stepsize h of array x "
endif

return
end subroutine romberg_sub

subroutine dromberg2d(jx,jy,x,y,z,s,pre)
implicit none
integer,intent(in)::jx,jy
double precision,intent(in)::pre,x(1:2**jx+1),y(1:2**jy+1)
double precision,intent(in)::z(1:2**jx+1,1:2**jy+1)
double precision,intent(out)::s

integer::i,nx,ny
double precision::ty(1:2**jy+1),r(1:2**jx+1)
nx=2**jx+1
ny=2**jy+1

s=0.d0
ty(1:ny)=0.d0
r(1:nx)=0.d0
do i=1,nx
ty(1:ny)=z(i,1:ny)
call romberg(jy,y,ty,s,pre)
r(i)=s
enddo
call romberg(jx,x,r,s,pre)

return
end subroutine dromberg2d

subroutine dromberg3d(jx,jy,jz,x,y,z,w,s,pre)
implicit none
integer,intent(in)::jx,jy,jz
double precision,intent(in)::pre,x(1:2**jx+1),y(1:2**jy+1),z(1:2**jz+1)
double precision,intent(in)::w(1:2**jx+1,1:2**jy+1,1:2**jz+1)
double precision,intent(out)::s
integer::i,nx,ny,nz
double precision::tyz(1:2**jy+1,1:2**jz+1),r(1:2**jx+1)

nx=2**jx+1; ny=2**jy+1; nz=2**jz+1
s=0.d0
tyz(1:ny,1:nz)=0.d0
r(1:nx)=0.d0
do i=1,nx
tyz(1:ny,1:nz)=w(i,1:ny,1:nz)
call romberg(jy,jz,y,z,tyz,s,pre)
r(i)=s
enddo
call romberg(jx,x,r,s,pre)

return
end subroutine dromberg3d
end module romberg_mod

ロンバーグ積分(romberg_mod)の例

必要なのは、上で紹介したモジュール “romberg_mod” と下のメインプログラムです。
以下のプログラムは1次元の定積分
$
\int_1^10 \frac{1}{x^2} dx = 0.9
$
を分点数$2^8+1(jx=8)$個でロンバーグ積分するものです。
精度は$O(h^{2\cdot min{jm,jx}})$です。
ここで$jm$はモジュールromberg_modの中にパラメータとして宣言されています。

コンパイルは

ifort romberg_mod.f90 main.f90

とでもすればいいでしょう。

program main
use romberg_mod
implicit none
integer::i,jx,n
double precision::h,a,b,s
double precision,allocatable::x(:),y(:)
double precision::f
external::f

a=1d0
b=10d0
jx=8
n=2**jx+1
allocate(x(1:n),y(1:n))
h=(b-a)/dble(n-1)
do i=1,n
x(i)=a+h*dble(i-1)
y(i)=f(x(i))
enddo

call romberg(jx,x,y,s,1d-8)
write(6,*)s,"romberg"

stop
end program main

function f(x)
implicit none
double precision::f,x

f=1d0/(x*x)

return
end function f

実行例(要求精度は$10^{-8}$)

>./a.out
0.900000000062642
>

2次元ではこう。

program main
use romberg_mod
implicit none
integer::i,j,Nx,Ny,jx,jy
double precision::hx,hy
double precision::xmin,xmax,ymin,ymax
double precision,allocatable::x(:),y(:)
double precision::s
double precision,allocatable::w(:,:)
double precision::f
external::f

jx=10
jy=8
xmin=0d0
xmax=3d0

ymin=0d0
ymax=5d0

Nx=2**(jx)
Ny=2**(jy)
allocate(x(0:Nx),y(0:Ny),w(0:Nx,0:Ny))
x=0d0; y=0d0; w=dcmplx(0d0,0d0)

hx=(xmax-xmin)/dble(Nx)
hy=(ymax-ymin)/dble(Ny)

do i=0,Nx
x(i)=xmin+hx*dble(i)
enddo
do i=0,Ny
y(i)=ymin+hy*dble(i)
enddo

do i=0,Nx
do j=0,Ny
w(i,j)=f(x(i),y(j))
end do
end do

call romberg(jx,jy,x,y,w,s,1d-8)
write(6,*)s

stop
end program main

function f(x,y)
implicit none
double precision,intent(in)::x,y
double precision::f

f=x*exp(-x*x)*exp(-y)

return
end function f

実行例

>./a.out
0.496569737374163
>

[adsense2]