sikino のすべての投稿

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

陰的ルンゲ=クッタ法の基本的な考えは
陰的ルンゲ=クッタ法
にて説明しました。
基本的には
\(
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. 計算方法
    1. 初期値の推定
    2. ニュートン法の停止
    3. 誤差判定
    4. 刻み幅制御
  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)を一緒にコンパイルしてください。

収束判定の余地


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

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

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

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

4倍精度とGECP


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

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

倍精度 d
4倍精度 q

プログラムの評価



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

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

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

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

結論


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

RADAU5の使い方


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

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

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

参考文献


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

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

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

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

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

グリーン関数による1次元の運動方程式の解法

グリーン関数を用いて1次元のニュートンの運動方程式を解きます。

  1. グリーン関数
  2. 重力下の質点
  3. 重力下のばね
    1. 1つ目の選び方
    2. 2つ目の選び方
  4. 注釈
  5. Thanks(フォントについて)

グリーン関数は非斉次微分方程式の特解を求める為に良く利用されます。
やっていることは、
微分方程式の満たす関数を、あたかも解いた風に書くこと
です。
なので実際、グリーン関数を使う、とは単なる式変形(微分方程式→積分方程式)をしているにすぎません。グリーン関数を本当に求め終わった時、はじめて解けた、と言えるのです。

物理現象では、グリーン関数が表現するのは「無限小時間に力が働く衝撃を受けた結果」を表現します。
どんな力も「複数の衝撃を積算したもの」と捉えることが出来ます。
衝撃に対する応答を考えてから、任意の力を衝撃の足し合わせと考え、問題を解きます。
すると、各々の力に対してその都度微分方程式を解いていくよりも
グリーン関数を求めた後の結果を足し合わせて任意の力に適応する、という考えの方が便利そうだと想像できます。

グリーン関数


微分方程式

をグリーン関数を用いて解くことを考えます。
ここで、左辺の\(\hat{D}\)は演算子を表し、例えば

のような形です。
グリーン関数で解くとは、解を(一般解)+(特解)で書くときの特解を与える手段です。
特解なので、何かしらの境界条件の下でグリーン関数を求める訳です。
物理の問題では、”因果律”と呼ばれる、結果は原因の後の時間に起こる、ということを仮定した境界条件を採用します。

全ての演算子を左辺に移してしまいましたが、右辺に時間に関する微分演算子が入っていても問題ありません。
式(1.1)の形式的な解は、グリーン関数\(G(t,t’)\)を用いて

と、式(1.3a),(1.3b)を満たす解を用いて書くことが出来ます。
ここで、\(\delta(x)\)はディラックのデルタ関数です。
実際に解(1.3)は式(1.1)を満たすことを確かめます。
式(1.3)に左から演算子\(\hat{D}\)を作用させます。
ここで、\(\hat{D}\)は\(t\)に関する演算子なので、積分と交換できます。よって

と変形でき、まさに解(1.3)は式(1.1)を満たす解として表されます。

重力下の質点


さて、質量\(m\)を持つ質点が、重力加速度\(g\)の重力下で運動する問題を考えます。運動方程式は

です。質量\(m\)で割って

を得ます。
もちろん、この問題の解は
\(
\displaystyle x(t)=-\frac{1}{2}gt^2+at+b
\)

です。この問題はグリーン関数を用いて解く必要は全くありませんが、解析解が簡単に求められる練習問題として良いでしょう。
グリーン関数は単なる形式ですが、一般的な解法であり、非常に強力です。
しかし、こんな簡単な問題も解けないのであればお話になりません。

さて、式(1.6)の形式的な解は

と書くことが出来ます。ここで、\(x_0(t), G(t,t’)\)は以下の微分方程式を満たします。

\(x_0(t)\)については簡単に解けてしまい、

です。ここで、\(\alpha, \beta\)は\(x(t)\)の初期条件によって決まる定数です。

続いてグリーン関数の微分方程式を解きましょう。
方程式の右辺を見ますと、\(t= t’\)だけで変なことが起こっています。
ですが、\(t\ne t’\)であれば右辺は完全にゼロです。
すなわち、\(G(t,t’)\)は\(t\ne t’\)で、微分方程式

を満たしています。これは簡単に解けてしまい、

です。ここで、\(A_0, A_1, B_0, B_1\)を適当な定数としておきました。
これらの係数間には、とある関係式が成り立っています。それを調べましょう。

デルタ関数を発散する点を含んで積分すると定数1になることを利用します。
式(1.18)の下の式の両辺を時刻\(t=t”\)周りで積分します。
すなわち、微小な時間\(\varepsilon\)を用いて、

を計算します。
もしも、たまたま\(t”=t’\)であれば、積分を実行して

を得ます。これは、時刻\(t=t’\)周りでは、グリーン関数の導関数は、右から極限と左から極限は一致しないことを言っており、不連続になると言っています。
\(t”\ne t’\)であれば、この極限は一致し、グリーン関数の導関数は連続であるということを言っています。

式(1.11)を代入すれば、係数間の関係式

を得ます。
式(1.7)のグリーン関数を含む項は特解ですので、不明な定数は存在しません。
今、4つの未知の定数\(A_0, A_1, B_0, B_1\)に対して式の数は2つなので、全ての係数を決めるには2つ条件式が足りません。
そこで、境界条件を加えます。この条件は式(1.7)の解の形を見て決められます。
運動方程式を解こうとする場合、通常は物理的な要請から、ある時刻の運動はその時刻より過去の運動のみから決定される、という因果律が課されます(※1)。
さて、因果律を仮定すると、事が起こる時刻\(t=t’\)より前には何も起こらない、すなわち\(G(t,t’)=0\)を意味しており、\(A_0=A_1=0\)となります。よって、

を得ます。よって、グリーン関数

を得ます。ここで、\(\theta(x)\)はヘヴィサイド関数を表します。
グリーン関数(t’=5に固定)を書いてみるとこんな感じ。

黒線がグリーン関数、赤線がその導関数です。導関数の右からと左からの極限の差がちょうどデルタ関数の積分、すなわち1になっていることが分かります。

これでグリーン関数が得られました。式(1.16)を使って式(1.7)の積分を計算すると、

を得ます。ここで、\(c\)は任意の値ではないことに注意してください。定まっている値です。

式(1.9)と式(1.17)の結果を代入すれば、運動方程式の解

を得ます。ここで、\(\beta+c\)を改めて\(\beta\)と書きました。
この結果は本当の解に一致します。

重力化のばね


続いて、鉛直方向に置かれたばねと重力の問題を考えましょう。
ばね定数を\(k\)と置けば、運動方程式は

です。質量\(m\)で割り、\(\omega=\sqrt{k/m}\)と置くと、

を得ます。
解は平衡点を中心とした振動運動です。平衡点\(x_0\)は力のつり合い
\(
-mg=kx_0~~~\to~~~x_0=-\frac{g}{\omega^2}
\)
から導くことが出来るので、
\(
\displaystyle x(t)=a\sin(\omega t +\varphi) -\frac{g}{\omega^2}
\)

です。これをグリーン関数を用いて導くことが出来れば目的達成です。

さて、グリーン関数の考え方として少なくとも2通りの方法が考えられます。
どちらの方法でも同じ解を与えるので、どちらをえらんでも構いません。
数学的に解きやすい方法を選ぶのが賢い方法です。

実際に試してみましょう。

1つ目の選び方

微分方程式を

として考えます。
解を

と考えるのです。\(x_0(t), G(t,t’)\)はそれぞれ

の微分方程式を満たします。式(1.9), 式(1.16)の結果を用いれば、容易に求められて、

です。式(1.22)の積分部分は

なので、解

を得ます。
さて、あらわな形ではないものの、解が得られました。
式(1.27)を満たす\(x(t)\)を探さなければならないのですが、
正直なところ解を見つけられる気がしません。これは自己無撞着の方程式であり、
例えば散乱問題等のように、グリーン関数を得るのが難しい時などでこの形式にとどめておく方法がとられます。
あらわに求める事が出来ずとも、他の式へ代入が出来たり、摂動展開なんかに用いられるので、微分方程式の形式よりその後の計算が扱いやすいのです。
ただ、今の目的からそれてしまうので、ここではこれ以上は話しません。

2つ目の選び方

微分方程式を

のように選びます。すると、解は

であり、\(x_0(t), G(t,t’)\)はそれぞれ

の微分方程式を満たします。
\(x_0(t)\)はただのばねの問題なので、

という解を得ます。一方、グリーン関数も\(t=t’\)を境に場合分けすると

を得ます。時刻\(t=t”\)周りの微小時間で積分すれば、

を得て、たまたま\(t”=t’\)であれば、

の関係式を得ます。ここで、グリーン関数自体はいかなる時刻でも右からと左から極限が一致するという事実をもちいました。
因果律を境界条件として使用すれば、係数は

であるので、グリーン関数は

と求められます。ここで、\(n\)は三角関数の周期性からくる任意の整数です。具体的なグリーン関数(t’=5に固定)の格好はこんな感じ。

黒線がグリーン関数、赤線がその導関数です。導関数の右からと左からの極限の差がちょうどデルタ関数の積分、すなわち1になっていることが分かります。

式(1.29)の積分を計算すると

と求められます。三角関数を時刻無限に飛ばした時の値ですが、拡張リーマン・ルベーグの補題から得られる事実を用います(※2)。これによると、

という事実が得られるので、式(1.37)の最後の項はゼロです。よって、

を得ます。
したがって、解

を得ます。これは求めたかった平衡点を中心に振動する運動です。

さて、1つ目の選び方で得た解(1.27)に解(1.42)を代入してみましょう。解(1.27)が本当にばねの問題なら、式は満たされるはずです。
式(1.27)の積分は

と計算できます。ここで、\(d\)を任意の定数ではない定数と置きました。
式(1.27)に代入すると、(RHS:right-hand side, 右辺)

を得ますが、今、\(\alpha’, \beta’\)は自由に決めることが出来ます。
もしも\(\alpha’=0,~~\beta’+d=0\)となるように定数を決めれば、解(1.27)は式(1.42)で表される解と同じで、微分方程式(1.20)の解となっていることが分かります。

注釈

※1) 因果律によるグリーン関数の境界条件は、妥当っぽく思われますが、人間が勝手にそうだろう、として経験的に与えている条件にすぎません。導くことは出来ない境界条件です。数学的には特解が何でも1つ与えられればいいので、時刻が未来であろうが特定の時間であろうが何でも構わないのです。
 今回の場合は解きたい微分方程式の右辺が定数(重力のみ)という非常に特別な場合ですが、一般的にはグリーン関数を用いた特解の被積分関数には、左辺の解である\(x(t)\)が含まれます。
この時、因果律を仮定しないと、『ある時刻の運動を計算するために、その時刻より未来の運動を知らなければならない』という理解しがたい状況が生まれます。言い換えれば、『現時刻の運動は、未来の振る舞いが決定されているので、どんな力を新たに加えようとも決して運動を変えることは出来ない』となってしまい、非物理的に感じてしまいます。
 もう一度言いますが、これは人間が「時間は必ず過去から未来に進む」という経験的に得られた事実を使用しているのであって、この根拠以外に因果律を採用する理由は存在しません。数学的に解こうとする際にはどんな条件でも構わないので、因果律を境界条件としても特解として正しい答えが得られるので良いのです。
 位置に対する微分方程式である場合は、例えば位置が無限大でゼロである条件だったり、ある点で位置と微分が決まっていたり、漸近で特定の形を持っていなければならない、という条件が当てはめられます。

(※2)私は物理の人間です。数学の人間ではないので簡単に使えます。微分と積分の入れ替えも、エイヤってやっちゃいます。直感的にこの定理は理解できますが、詳細は分かりません。物理的に欲しい結果を与えてくれるので使っちゃいました。

Thanks

アイキャッチ画像のフォントはmagicringを使用しています。
配布ページ
http://inatsuka.com/extra/magicring/
配布元のサークル様
森の中の猫の小屋 http://inatsuka.com/

陰的ルンゲ=クッタ法

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


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

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

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

陽的解法と陰的解法


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

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

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

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

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

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

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

一階の微分方程式


問題

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

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

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

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

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

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

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

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

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

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

で表される行列です。

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

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

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

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

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

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

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

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

行列要素の計算

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

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

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

簡易Newton法

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

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

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

ということです。

複数の微分方程式


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

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

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

です。それぞれ

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

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

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

です。簡易Newton法では

として得られます。

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


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

2段4次

3段6次

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

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

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

プログラム


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

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

3段6次 陰的Runge-Kutta

2段4次 陰的Runge-Kutta

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

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

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

刻み幅制御


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

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

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

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


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

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

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

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

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

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

となります。


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

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

プログラムの名づけ方

倍精度 d
4倍精度 q

LAPACK使用 l
GECP使用 g

openMP使用 m
不使用 n

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

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

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

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

GECPとは

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

参考文献


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

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

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

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

その他参照したページ


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

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

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

-1のn乗の計算速度

fortran90で、
\(a=(-1)^n, ~(n=0,1,\cdots, M)\)を計算する早いアルゴリズムは

a=1
if(mod(n,2).eq.1)a=-1

です。


物理、数学をやっていると至る所で\((-1)^n\)を見受けます。
これを数値計算する際に早い計算方法は何でしょうか?

試してみるのは以下のものです。


\(
(-1)^n
\)


\(
(-1)^n
\)
※①を倍精度型で計算


\(
e^{i\pi n}
\)
※\(i\)は虚数単位、\(\pi\)は円周率


\(
\begin{eqnarray}
\left\{
\begin{aligned}
& 1 ,~~(n=0,\mbox{偶数})\\
& -1 ,~~(n=\mbox{奇数})\\
\end{aligned}
\right.
\end{eqnarray}
\)


\(
a=-a
\)で逐次計算

使ったプログラムは以下のものです。

2次元時間依存しないシュレーディンガー方程式の数値解法

2次元の時間依存しないシュレーディンガー方程式を変分原理に基づいて解きます。

対象とする問題は以下の時間依存しないシュレーディンガー方程式です。

ここで、\(H_0\)は数値計算で用いる基底関数のハミルトニアン、\(V(\mathbf{r})\)はその他のポテンシャルです。

この系の固有関数\(\phi(\mathbf{r})\)を

として\(H_0\)の固有関数で展開します。ここで、\(\varphi_{\mathbf{r}}\)は

の固有値問題の解である、固有値\(E’_n\)に属する固有関数です。また、\(\varphi_{\mathbf{r}}\)は

として規格直交化されています。

式(2)を式(1)に代入した後、左から\(\varphi^*_{\mathbf{r}}\)を掛けて全空間で積分すれば

が得られます。ここで、\(V_{m,n}\)を

と置きました。式(7)を別の表現をすれば

と書くことが出来ます。なので、左辺のハミルトニアンの行列形式を対角化すれば、求めたい系の固有値、固有ベクトルが得られるわけです。固有ベクトルが得られたら式(2)に従って波動関数を構成すれば良いのです。

これが一般的な変分原理による解法です。
以降は2次元の問題に限り考えていきます。

2次元の場合


2次元の場合かつ\(H_0\)はxだけ含む部分とyだけ含む部分とで分離することが出来る場合を考えます。

この場合、固有関数は変数分離することが出来るので、2次元の基底関数を1次元系の基底関数の直積をとして考える方法を取ることが出来ます。すなわち、x,y方向の量子数\(n_x,n_y\)を用いて量子数\(n\)と以下のように関係して構成される、と考えるのです。

ここで、\(n\)と\(n_x,n_y\)には1対1の関係があれば順番は関係ないことに注意してください。それぞれの基底関数は

を満たします。もちろん、1次元系の基底関数がそれぞれ正規直交化されていれば二次元の場合でもその関係は保たれ、

が成立します。この基底関数で式(7)を書き換えてやれば

となります。

基底関数の対応


この問題はある長径、短径で指定される楕円内に存在する格子点の個数を数える問題と同じです。

nの順番を記録しておきさえすればどんな順番でも構いません。適当に求めればよいです。
ただし、ここではシュレーディンガー方程式を解いた結果、エネルギーの低い状態から順に欲しいので、\((n_x, n_y)\)の組み合わせでエネルギーが低い順に決定していきます。

今、エネルギーは\(x,y\)方向ともに単調に増加します。
\((n_x,n_y)\)の考えられる組み合わせで、
一番小さいエネルギーは\((1,1)\)なので、まず\(n=1\)は\((1,1)\)に対応させます。
次に小さいエネルギーで考えられるのは\((1,2),(2,1)\)のどちらかです。
ここで重要なのは現在より前に出てきた\((n_x,n_y)\)の組に1を足した組み合わせしか候補にあがりません。
なので、探索するのはこの条件を満たすものだけで充分です。

これをプログラムしたものは、下のtar.gz内に”qnumber”という名前でサブルーチンに入っています。

プログラム


こちらです。
http://slpr.sakura.ne.jp/qp/supplement_data/highprecision_tise2d.tar.gz

解凍すると、中にintegrate.f90 main.f90というファイルがはいっています。
Lapackを使うのでそれにリンクしてコンパイルしてください。intel®fortranならば、

ifort -mkl integrate.f90 main.f90

で良いと思います。
実行すると、

$ ./a.out
 Solving TISE...
     10 /    100
     20 /    100
     30 /    100
     40 /    100
     50 /    100
     60 /    100
     70 /    100
     80 /    100
     90 /    100
    100 /    100
     1.807[CPU sec]
                    1   1.0000086267300243    
                    2   2.0001412369101184    
                    3   2.0001696393476047    
                    4   3.0004659233630600    
                    5   3.0006690578629835    
                    6   3.0013176072105336    
                    7   4.0029644525394392    
                    8   4.0032488578583107    
                    9   4.0091516049507794
$

という結果が得られるでしょう。ここで、** / 100という数字はハミルトニアン行列要素の計算の経過を表しています。100は行列が\(100\times 100\)の行列であるという意味です。
1 1.0000086267300243
は、対角化後の固有エネルギーの一番低い状態から順に出力しています。
デフォルトでは二次元の調和振動子ですので、
\(
E_n=(n_x+1/2)+(n_y+1/2)
\)

の順になります。一番低いエネルギーから順に\(1,2,2,3,3,3,4,4…\)が解析解となります。

具体的な計算例


三角形型のポテンシャル

\(
\begin{equation}
V(x,y)=\left\{
\begin{array}{c}
0,~~~(y\lt x \cap x\lt a \cap y\gt -a) \\
10,~~~(otherwise)
\end{array}
\right.
\end{equation}
\)
ここで、\(a=3\).

2つの四角形のポテンシャル

\(
\begin{equation}
V(x,y)=\left\{
\begin{array}{c}
0,~~~(y\lt (x+b) \cap y \gt -(x+b) \cap y\lt -(x+b)+a \cap y \gt (x+b)-a) \\
0,~~~(y\gt (x-b) \cap y \lt -(x-b) \cap y\gt -(x-b)+a \cap y \lt (x-b)+a) \\
10,~~~(otherwise)
\end{array}
\right.
\end{equation}
\)
ここで、\(a=3\)

1次元時間依存シュレーディンガー方程式を解かねばならぬ

1次元時間依存シュレーディンガー方程式を確実に解くことを目的とします。

本稿の方針は
計算時間が掛かってもいいので、高精度に解く
ことを主眼においています。対象とする問題は、
ポテンシャルに空間的、時間的に不連続な点があっても高精度に解く
事を目的とします。

ポテンシャルが滑らかで、変なことが起こらないのならば、空間や時間を離散化して解く
クランク=ニコルソン法による時間発展が簡単で優秀な方法です。

本稿では導出する過程では\(N\)次元の問題として扱っていますが、プログラムは1次元の問題のみです。

方針


空間に関する積分を適応型数値積分、時間に関する積分を適応刻み幅ルンゲ=クッタ法で行います。
要は空間も時間も離散化しないで解く、ということです。

時間依存しないシュレーディンガー方程式は適応型数値積分であるQUADPACKを使うことで求められます。
詳しくはこちらへ。

解法


原子単位系における時間依存シュレーディンガー方程式は

として書けます。
波動関数を

のように展開します。ここで、基底関数\(\varphi_n(\mathbf{r})\)はハミルトニアン\(H_0\)の固有関数で、以下の固有値問題

の固有値\(E_n\)に属する固有関数です。
この固有関数は以下の通り規格化されているとします。

関数で展開するのはハミルトニアンに含まれる空間に対する二階微分を解析的に行い、消す為です。
式(1)に代入すると、

となり、左から\(\varphi^*_m(\mathbf{r})\)を掛けて全空間で積分すると

という式が得られます。ここで、\(U_{m,n}(t)\)を

と書きました。

この\(U_{m,n}(t)\)は

の性質を満たします。すなわち、ポテンシャル\(V(\mathbf{r},t)\)が実数値関数であれば、\(U_{m,n}(t)\)はエルミート行列となります。行列形式で書くとすれば、

と表せられます(補遺1)。

ここでは行列形式で書いても見やすくなるだけで、今回は使わないので、式(6)を高精度に解くことに集中します。

高精度に解くために

この問題を解くために高精度に求めなければならない重要な部分は

1. 空間に対する積分
\(\displaystyle
\int \varphi^*_m(\mathbf{r})f(\mathbf{r})\varphi_n(\mathbf{r})d\mathbf{r}
\)


2. 時間に対する積分(連立微分方程式)
\(\displaystyle
\frac{d c_m(t;c_1,~\cdots,c_{N})}{dt}=F_m(t;c_1,~\cdots,c_{N})
\)

という箇所を、不連続箇所があっても高精度に求める手法を使えばよい訳です。

空間に関してはQUADPACKによる適応型の数値積分, 時間に関しては刻み幅制御のルンゲ=クッタ法を用いることで対処します。
各々の詳細は以下のページをご覧ください。
最速・高精度の数値積分
ルンゲクッタ法の説明と刻み幅制御

ここまで確実に解けると仮定すると、精度の不安は基底関数の数だけに依存します。
基底関数の数を積めば積むほど精度が高い計算が出来ますが、計算時間が莫大になります。
ここは計算機の性能と折り合いをつけるしかありません。

ここから考える数値計算上の基底関数\(\varphi_n(x)\)は計算区間の端でゼロになるsin基底関数で、区間[x_a,x_b]で定義され、
\(\displaystyle
\varphi_n(x)=\sqrt{\frac{2}{x_b-x_a}}\sin\left(n\pi\frac{x-x_a}{x_b-x_a}\right)
\)

です。これは、固有値問題
\(\displaystyle
-\frac{1}{2}\frac{d^2}{dx^2}\varphi_n(x)=E_n\varphi_n(x),~~E_n = \frac{1}{2}\left(\frac{n\pi}{x_b-x_a}\right)^2
\)

境界条件
\(\varphi_n(x_a)=\varphi_n(x_b)=0\)
の元での解となっています。

さて、定式化は終わりましたので、あとはプログラムして解くだけです。

莫大な計算時間


問題は計算時間が膨大にかかりすぎる事です。

計算時間を評価してみます。

・次元数を\(D\)
・1次元辺りの基底関数の数は\(N\)
・合計\(N^D\)元1次連立微分方程式を解く
・1つの連立微分方程式の右辺を計算するのに\(N^D\)個の積分
・適応刻み幅ルンゲクッタ法を\(N_{t}\)回繰り返す
・適応刻み幅ルンゲクッタ法は時間1ステップ辺り\(5\)回評価
・積分は空間の分点数\(N_{\mathbf{r}}^D\)回に比例
・ポテンシャルが対称行列であることを利用すれば\(1/2\)倍
・CPU個数を\(N_{\mathrm{cpu}}\)とし、完璧に並列計算が行われれば\(1/N_{\mathrm{cpu}}\)倍
・関数を1回呼び出すのに必要な時間を\(t_f\)

であるので、全計算時間(関数が評価される回数に比例する、と考える)は
\(\frac{5}{2N_{\mathrm{cpu}}}\times N^D_{\mathbf{r}}N_{t}N^{2D}t_f\)と求められます。

現実的な量で計算時間を見積もってみます。
10[原子単位]秒の時間発展を区間20[原子単位]メートルで計算することを考えます。

\(D=1,~ N_{\mathbf{r}}=1000,~N_{t}=2000,~N=50,~N_{\mathrm{cpu}}=1\)とすると、\(1.25\times 10^{10}\)回関数を呼び出さなければなりません。

手持ちのノートパソコン(1.9GHz、1コア使用)で、\(2.5\times 10^8\)回の関数\(f(x)=\sin(x)\)の呼び出しに掛かるCPU時間を測ると約7[cpu秒]かかりました。
ということは、関数を1回呼び出すのにかかる時間\(t_f\)は約\(3\times 10^{-8}\)秒です。

よって、\(1.25\times 10^{10}\)回の関数の呼び出しだけで約400秒、すなわち8分かかります。

現実には他の処理も入るので、8分は絶対にかかる、という事です。
仮にCPUが10個あれば40秒、CPU速度が倍程度の物を使えば20秒位になるでしょうか。
まぁまぁ受け入れられる時間になります。

2次元を計算しようと思ったら基底関数が二乗、かつ空間の分点も二乗になるので、CPUが10個,4GHz近くのCPUを使えても\(20\times 1000\times 50=1000000\)秒程度、すなわち11日です。

11日の計算は現実的ではありません。
二乗で効いてくる基底関数の数を減らすか、空間の積分の精度を犠牲にする、CPUを増やす等の対策をしなければなりません。

また、ここで示した時間は最低かかる時間ですので、実際に動かすのであればこれ以上は確実にかかるでしょう。関数の中の関数を呼び出す等の操作があれば倍々に増えていくのです。

計算量を減らす工夫


ポテンシャルが変わらない場合、\(U_{m,n}(t)\)の積分は、前の時刻と値が変わらないのでいちいち計算するのは無駄です。

数値的に”ポテンシャルが変わっていない”という情報を取り出すことが出来れば、この無駄を排除することが出来ます。

”ポテンシャル関数が変わっていない”というのはどうすれば評価できるのでしょうか?
私が考えたのは、Quadpackの自動積分を行う時に評価された引数の点における関数の値が全て一致すれば関数が変わっていない、と判断するという評価方法です。

実際にこれはうまく行きまして、ポテンシャルが一度だけ切り替わる、という問題に対して計算時間が
1085[CPU sec] → 25[CPU sec]
にまで減少しました。約1/50の計算時間になりました。

先ほどの11日で~と言っていたのはこの工夫をし無ければ、ですので、ポテンシャルの変化が殆どない問題に限れば1時間程度で計算が終わる、ということです。

もちろん、ポテンシャルが変化し続ける場合も試しましたが、変化はほぼありません。
83[CPU sec] → 84[CPU sec]
程度になりました。
この場合は2次元の例で出した通り、11日間掛かる推測です。

初期状態の準備


式(6)を解こうとしても初期条件が与えられなければ数値的に解くことは出来ません。
初期条件で良く使われるのは、1)解析的な解を与える場合と 2)系の固有状態を与える場合の2種類があります。

  1. 解析的な解を与える場合

    この場合、時間発展を考える場合、上記の時間発展方法ではsin関数基底に射影して係数\(c_n(t)\)を求めなければなりません。
    時刻\(t=t_a\)で初期状態\(\psi(x,t_a)\)が解析的に与えられた場合、係数\(c_n(t)\)は、

    として与えられます。

  2. 系の固有状態を与える場合

    この場合は時刻\(t_a\)のポテンシャル中の時間依存しないシュレーディンガー方程式
    \(\displaystyle
    \left[H_0 + V(\mathbf{r},t_c)\right]\phi_n(\mathbf{r})=E’_n\phi_n(\mathbf{r})
    \)

    を実際に解かなければなりません。
    そのためには、時間発展に使う基底と同じ基底\(\varphi_n(\mathbf{r})\)で展開し、対角化するのが賢い方法です。
    すなわち、

    として展開し、ハミルトニアンを行列表示にし、対角化して固有値、固有ベクトルを得ればいいのです。

計算の解析


計算結果を解析するにあたり、ある時刻の固有関数系の基底にどの位存在確率があるかという情報が良く使われます。

初期状態で利用した固有状態の存在確率を調べましょう。
今、求めたいのは時刻\(t_c\)のポテンシャル中の固有関数で展開した時の係数\(a_m(t)\)を求める事です。

ここで、\(\phi_n(\mathbf{r})\)は時刻\(t_c\)のポテンシャル中の固有関数であり、

を満たしています。また、\(\varphi_n(\mathbf{r})\)は数値計算を行う際の基底関数であり、

を満たしています。
任意の時刻\(t\)における波動関数を

として表したいので、\(a_n(t)\)について求めていけば、

という関係式を導くことが出来ます。

数値計算プログラム


プログラムは以下のリンクにあります。
http://slpr.sakura.ne.jp/qp/supplement_data/highprecision_1Dtdsetise.tar.gz
そのままのプログラムでは、以下の設定となっています。
プログラムの利用に関してはhttp://slpr.sakura.ne.jp/qp/about/をお読みください。

解凍すると、4つのファイル(main.f90, integrate.f90,anime.plt,anime_project.plt)が入っています。計算パラメータの設定は、main.f90に入っているので各自調整してください。

繰り返しますが、全て原子単位系です。デフォルトでは、
————————————————–
計算時間の範囲\([ta,tb]\) … \([-5, 20]\)秒
出力時間枚数\(Nt\) … \(100\)枚
1ステップ当りの時間積分の精度\(\mathrm{RKeps}\) … \(10^{-5}\)
1ステップ当りの時間の最小刻み幅\(\mathrm{hmin}\) … \(10^{-4}\)
※理想は最小刻み幅は倍精度ならば\(10^{-14}\)程度が良いと思いますが、時間変化するポテンシャルを計算しようとする時、いつまでも最小刻み幅で進む時間がたまにあります。こういう時間はさっさと抜け出たいので、それを考慮して設定できるようにしています。

位置空間の範囲\([xa,xb]\) … \([-10, 10]\)
出力時の分点数\(Nx\) … 200点
時間積分の精度\(\mathrm{QPeps}\) … \(10^{-6}\)

初期状態を時間依存しないシュレーディンガー方程式を解くか、解析解で与えるかiname … “tise”
inameが”tise”の場合、どの時刻\(t_{ini}\)の固有状態を使うか … \(t_a\)

基底関数の数\(N\) … \(30\)個
計算に用いるスレッド数\(N{\mathrm{cpu}}\) … 1つ
空間積分の分点数を強制的に増やすか否かQPforce … 0

初期状態波動関数… 時間依存しないシュレーディンガー方程式の基底状態
ポテンシャル\(
\begin{eqnarray}
V(x,t)=\left\{
\begin{aligned}
0.5 x^2 (t\lt 0)\\
0.1 x^2 (t\ge 0)
\end{aligned}
\right.
\end{eqnarray}\)
————————————————–
で解くプログラムとなっています。

lapackと、必要があればopemMPを使います。

ifort -mkl -openmp integrate.f90 main.f90
./a.out

もしくは、gfortranでlapackにリンクしてintegrate.f90 main.f90を一緒にコンパイル、実行してください。
実行後、以下のファイルを生成します。

  • wf_***.d
    波動関数の各時刻におけるデータです。時間順にwf_0.d(初期状態), wf_1.d, wf_2.d,…という風に指定された等間隔の時刻で出力します。アウトプットのファイル数はNtで指定した数です。
  • project_ratio.d
    各時刻でのsin関数基底への射影した際の式(2)の係数\(c_n(t)\)の絶対値二乗を出力します。最も高い状態が振幅を持ってしまう場合、基底数が足りていないことを示しています。\( |c_n(t)|^2\)の事です。
  • unitarity.d
    ユニタリー性を確認します。要は数値計算の発散とかで全計算区間内の存在確率が変化していないか見るためのものです。\(\sum_{n=1}^N |c_n(t)|^2\)の事です。
  • timerange.d
    時間グリッドを確認します。物理的な意味は無く、どの時間で計算負荷が大きいかを見るためのものです。
  • timeindex_correspondence.d
    時間インデックス(wf.***.dの***の値と、原子単位系の時刻の変換)を出力します。

もしも初期状態を時間依存しないシュレーディンガー方程式を解いて準備していた場合、以下のファイルを追加で出力します。

  • initial_eigenvalue.d
    初期状態のエネルギー固有値を出力します。式(14)の\(E’_n\)です。
  • project_tdse_to_tisebasis.d
    各時刻における初期状態の固有状態に射影した時の係数の絶対値の二乗を出力します。式(13)の係数\(a_n(t)\)の絶対値二乗である\(|a_n(t)|^2\)を出力します。

波動関数をgnuplot上でアニメーションとして見たい場合、同封してあるファイル”anime.plt”をgnuplot上でloadすれば良いです。

gnuplot> load "anime.plt"

また、gifアニメが欲しい場合、

gnuplot> load "anime_project.plt"

とすれば、波動関数と時間依存しないシュレーディンガー方程式を解いた時の基底関数への係数の絶対値二乗を共に出力します。

実行例


デフォルトの設定で計算を行います。
実行すると

$ ./a.out
0 / 100        1.0000000000
10 / 100        1.0000001064
20 / 100        1.0000001082
30 / 100        1.0000012466
40 / 100        1.0000079960
50 / 100        1.0000082264
60 / 100        1.0000114047
70 / 100        1.0000154242
80 / 100        1.0000154993
90 / 100        1.0000193449
     1.948[CPU sec]
$

という文が出力されます。
左の***/100とは、今計算が100枚中***枚まで終わった、ということを示しており、枚数は変数Ntに対応しています。
右の1に近い数字は、全空間の存在確率を表しており、計算が発散しないかどうかを確認できます。これはルンゲクッタ法の要求精度\(RKeps\)に大きく依存します。
最後のCPUsecは計算に要したCPU時間であり、1コアではおおよそ実時間と同一です。なので、この計算は約2秒で終わる、ということです。ちなみにCPUは1.9GHzの物です。

anime_project.pltによってgifアニメを生成すると、以下のgifアニメが得られます。

その他の出力されるデータは以下の通り。

計算例


ポテンシャルによる反射


※この上図の縦軸は数値計算上の基底である\(\varphi_n(x)\)への射影した時の係数の絶対値二乗です。

ラビ周波数


基底状態と第一励起状態の間の周波数でポテンシャルを揺らしています。

鬼畜なポテンシャル


こんな変なポテンシャルでも計算できます、というデモです。

補遺1

もし仮に、\(U_{m,n}(t)\)が時間によって変化しないのであれば、
\(\displaystyle
\frac{\partial \mathbf{x}}{\partial t}=\mathbf{A}\mathbf{x}
\)

の形の微分方程式であるので、指数関数の解を仮定して\(\mathbf{A}\)を対角化、その固有値固有ベクトルを求めて、線形結合で表現すればokです(線形代数II/連立線形微分方程式等にあるのでそちらを参照してください)。

gnuplot上で行う任意関数のフィッティング

gnuplotには、”fit”というコマンドがあります。
これを利用すると、データ列に対して、任意の自由度を含んだ関数\(f(x)\)でフィッティングを行うことが出来ます。
手法はMarquardt-Levenberg法に基づいた非線形最小二乗法によるフィッティングです。

gnuplot ver4.6の場合

データ列”data.d”を\(f(x)=ax^2+bx+c\)でフィッティングする\(a,b,c\)を求めたい場合、gnuplot上で

f(x)=ax**2+bx+c
fit f(x) "data.d" u 1:2 via a,b,c

※gnuplot ver5.0以降の場合、エラーバーで重みづけしたフィッティングが出来るようです。
http://gnuplot.sourceforge.net/demo_5.0/fit.html

説明


以下のデータ列(ファイル名”data.d”)に対してgnuplot上でフィッティングを行います。

データは、関数
\(
\displaystyle f(x)=\frac{1}{1+e^{x-15}}
\)

に乱数を与えて生成したデータです。

このデータ列に対して、
\(
\displaystyle f(x)=\frac{a}{1+e^{(x-b)/c}}
\)

でフィッティングを行い、未知の定数\(a,b,c\)を決めたいと思います。

gnuplot上で

f(x)=a/(1+exp((x-b)/c))
a=1e0
b=10e0
c=0.5e0
fit f(x) "data.d" u 1:2 via a,b,c

と打つと、\(a,b,c\)にフィッティングした後の定数が入ります。
ここで、フィッティング前に\(a,b,c\)の値を入れているのは初期値です。
初期値があまりにも違うとうまくフィッティングが失敗するので、出来るだけ近い値を入れましょう。

実際に動かすと

> fit f(x) "data.d" u 1:2 via a,b,c
...

Final set of parameters            Asymptotic Standard Error
=======================            ==========================

a               = 0.995618         +/- 0.005679     (0.5704%)
b               = 15.3268          +/- 0.051        (0.3327%)
c               = 0.970465         +/- 0.04439      (4.574%)

...
>

という文章が出力されます。

フィッティングの結果、
\(a=0.995618,b=15.3268,c=0.970465\)と求められました。
乱数を与える前のデータは\(a=1,b=15,c=1\)ですので、元の関数に近いことが分かります。
fitコマンドを動かした後は既に\(a,b,c\)に値が代入されているので、そのまま

plot f(x)

とすれば、フィッティング結果をプロットすることが出来ます。

データと共に載せれば以下のようになります。

フィッティングする範囲を制限したければ、

fit [10:20] f(x) "data.d" u 1:2 via a,b,c

とすると、\(x\)の範囲を\([10,20]\)に制限することが出来ます。

2変数関数のフィッティング


2変数関数であろうとフィッティングは可能です。
元のデータ
(http://slpr.sakura.ne.jp/qp/supplement_data/data2.d)
は、関数
\(
g(x,y)=x\sin(xy+0.5)
\)

に適当な乱数を足して作られたデータです。グラフにすれば

となります。

このデータに対して、gnuplot上で関数
\(
g(x,y)=ax\sin(bxy+c)
\)

によるフィッティングを行います。コマンドは

g(x,y)=a*x*sin(b*x*y+c)
a=1.3e0
b=0.7e0
c=0.3e0
fit g(x,y) "data2.d" u 1:2:3 via a,b,c

です。
実行すると

>fit g(x,y) "data2.d" u 1:2:3 via a,b,c

...
Final set of parameters            Asymptotic Standard Error
=======================            ==========================

a               = 1.16313          +/- 0.004072     (0.3501%)
b               = 1.00073          +/- 0.001672     (0.1671%)
c               = 0.495134         +/- 0.003432     (0.6932%)
...

>

と得られますので、フィッティング結果と共にプロットすれば

となります。

初期値が答えに割と近いですが、これは失敗を防ぐためです。全ての初期値のパラメータを0で行った所失敗しました。

参考サイト


Gnuplotでの関数fitting
gnuplot demo script: fit.dem -ver5.0
gnuplot demo script: fit.dem -ver4.6
Fit (manual)
gnuplot で最小二乗フィッティングする -ゴルディアスの涙目

PV数の推移

シキノートは2014年12月に開始したので、現在(2018年8月)までに3年半経過しました。
おかげ様です。

自己満足のためにシキノートの閲覧データを公開します。

PV(ページビュー)数の推移



灰色の二つの線は線形補間と二次関数による補間です。
線形補間では負の切片を持ってしまったため、
恐らく二次関数的にPV数が増えている、というのが正しいでしょう。嬉しい限りです。

年ごとのPV数の推移は(1年は1月~12月)
2014年:64
2015年:23,993
2016年:67,647
2017年:125,982
2018年:141,141(8月17日まで)
となっています。

開始から43ヶ月
合計PV数:358,827
記事数:154
なので、1ヶ月当たり3.6記事です。

Analyticsのデータ


現在までの人気記事


タイトル 総PV数
gnuplotのカラーマップ 24,016
球体の空気抵抗と係数 17,803
Nexus5の購入とセルスタンバイ問題の解決 15,352
ルンゲクッタ法の説明と刻み幅制御 15,111
4次ルンゲ・クッタ法 10,836
レジスト系によるダメージ減少量 10,373
弾道計算(BB弾)の理論 10,280
gnuplotとアニメーション 9,260
離散フーリエ変換と高速フーリエ変換(fortran90) 8,299
三角波、のこぎり波、矩形波、その他の数式 7,768
弾道計算(BB弾)の結果 7,541
fortranコンパイラのコマンド 7,455
演算子の種類と説明 6,508
latexで数式全体を小さくする 6,109
乱数 5,849
正多角形とスピログラフの数式 5,834
fortranでのエラーメモ 5,543

二体衝突後に物体が持ち得る速度

2体が弾性衝突(エネルギー保存)を起こす時を考えます。
目的は、衝突後に2体が持つ速度ベクトルを知ることです。

2体を
物体1… 重さ\(m_1\)
物体2… 重さ\(m_2\)
衝突前の速度ベクトルを\(\mathbf{v}^{(i)}_1, \mathbf{v}^{(i)}_2\)、
衝突後の速度ベクトルを\(\mathbf{v}^{(f)}_1, \mathbf{v}^{(f)}_2\)
と書くと、衝突後に考えられる速度ベクトルは、
\(
\begin{align}
\mathbf{v}^{(f)}_1&=R(\theta)(\mathbf{v}^{(i)}_1-\mathbf{V})+\mathbf{V}\\
\mathbf{v}^{(f)}_2&=-R(\theta)(\mathbf{v}^{(i)}_2-\mathbf{V})+\mathbf{V} \\
\end{align}
\)

と表せられます。ここで、\(\mathbf{V}\)は重心の速度ベクトルで
\(
\begin{align}
\mathbf{V}=\frac{m_1\mathbf{v}^{(i)}_1+m_2\mathbf{v}^{(i)}_2}{m_1+m_2}
\end{align}
\)

を表し、\(R(\theta)\)はユニタリーな回転行列を表します。衝突が起こる際の角度はこの式からは求められません、あくまで、衝突によって起こり得る角度です。
ここで示しているのは、いろんな方向に衝突が起こることを考えた時の表記で、衝突によって変化することが出来る速度ベクトルを表しています。

特に二次元上での衝突を考える場合、回転行列は
\(
\begin{equation}
R(\theta)=
\left(
\begin{array}{rr}
\cos\theta& -\sin\theta\\
\sin\theta& \cos\theta\\
\end{array}
\right)
\end{equation}
\)

と書くことが出来ます。

二体の衝突


ある2つの物体が衝突することを考えます。
下添え字を物体を区別する番号、上添え字を衝突前\(i\) (initial), 衝突後\(f\) (final)として記述すると、
物体1… 重さ\(m_1\), 位置ベクトル\(\mathbf{r}_1\), 速度ベクトル\(\mathbf{v}_1\)
物体2… 重さ\(m_2\), 位置ベクトル\(\mathbf{r}_2\), 速度ベクトル\(\mathbf{v}_2\)
と書けます。但し、ここでは衝突前、後を表す上添え字\(i,f\)は省いています。

また、重心の位置,速度ベクトル\(\mathbf{R},\mathbf{V}\)は
\(
\begin{align}
\mathbf{R}&=\frac{m_1\mathbf{r}_1+m_2\mathbf{r}_2}{m_1+m_2} \\
\mathbf{V}&=\frac{m_1\mathbf{v}_1+m_2\mathbf{v}_2}{m_1+m_2}
\end{align}
\)

と書けます。

続いて重心から見た物体1,2の位置、速度ベクトルを\(‘\)を付けて表しますと、

物体1(重心系)… 重さ\(m_1\), 位置ベクトル\(\mathbf{r}’_1\), 速度ベクトル\(\mathbf{v}’_1\)
物体2(重心系)… 重さ\(m_2\), 位置ベクトル\(\mathbf{r}’_2\), 速度ベクトル\(\mathbf{v}’_2\)
と書くことが出来ます。
ここで、重心系における速度ベクトルは明示的に
\(
\begin{align}
\mathbf{v}’_1&=\mathbf{v}_1-\mathbf{V} \\
\mathbf{v}’_2&=\mathbf{v}_2-\mathbf{V}
\end{align}
\)

と表すことが出来ます。

重心系から見ると重心の位置、速度ベクトルはもちろんゼロとなります。
\(
\begin{align}
\mathbf{R}’&=\frac{m_1\mathbf{r}’_1+m_2\mathbf{r}’_2}{m_1+m_2}=\mathbf{0}\\
\mathbf{V}’&=\frac{m_1\mathbf{v}’_1+m_2\mathbf{v}’_2}{m_1+m_2}=\mathbf{0}
\end{align}
\)

重心系から見た系の全運動量\(\mathbf{p}’\)もゼロにならなければなりません。
このことから、
\(
\mathbf{p}’=m_1\mathbf{v}_1+m_2\mathbf{v}_2=\mathbf{0}
\)

という関係式が成り立ちます。
この式が言っているのは、重心系を考える限り、物体1,2の速度ベクトルは\(\mathbf{v}’_1,\mathbf{v}’_2\)は衝突が起こる起こらないにかかわらず、上式が満たさなければならない、ということです。
繰り返しますが、衝突が起こる起こらないにかかわらず、です。

衝突が起こり、どの角度に飛んでいくか?を与えられた情報からでは決めることは出来ません。
実験結果や、物体の形状を表す別の条件式が必要となります。

ここで出来ることは、起こりうる角度を知ることだけです。
衝突後に起こりうる角度を求めてみましょう。
重心系では系の合計の運動量\(\mathbf{p}’\)がゼロベクトルにならなければなりません。
このことは衝突前後でノルムが変わる衝突は許されないことを示しています。
これは、重心系においてユニタリーな変換が許されることを意味しています。

衝突後に重心系で取りうる角度は重心系で考えられる角度方向全てなので、
ユニタリーな回転行列\(R(\theta)\)を用いて、衝突後に重心系で速度ベクトルが
\(
\begin{align}
\mathbf{p}’_1&\to R(\theta) \mathbf{p}’_1\\
\mathbf{p}’_2&\to -R(\theta) \mathbf{p}’_2
\end{align}
\)

を満たす変化のみ許される、ということです。マイナスの符号は、\(\mathbf{p}’\)がゼロでなければならないという条件からきているものです

実験系での表現に戻せば、運動量ベクトルは
\(
\begin{align}
\mathbf{p}^{f}_1&=m_1R(\theta) (\mathbf{v}^{(i)}_1-\mathbf{V})+m_1\mathbf{V}\\
\mathbf{p}^{f}_2&=-m_2R(\theta) (\mathbf{v}^{(i)}_2-\mathbf{V})+m_2\mathbf{V}
\end{align}
\)

であり、速度ベクトルは
\(
\begin{align}
\mathbf{v}^{f}_1&=R(\theta) (\mathbf{v}^{(i)}_1-\mathbf{V})+\mathbf{V}\\
\mathbf{v}^{f}_2&=-R(\theta) (\mathbf{v}^{(i)}_2-\mathbf{V})+\mathbf{V}
\end{align}
\)

と書き表すことが出来ます。
特に二次元の場合、デカルト座標において回転行列は
\(
\begin{equation}
R(\theta)=
\left(
\begin{array}{rr}
\cos\theta& -\sin\theta\\
\sin\theta& \cos\theta\\
\end{array}
\right)
\end{equation}
\)

と書くことが出来るので、
\(
\begin{equation}
\left(
\begin{array}{r}
v^{(f)}_x\\
v^{(f)}_y
\end{array}
\right)
=
\left(
\begin{array}{rr}
\cos\theta& -\sin\theta\\
\sin\theta& \cos\theta\\
\end{array}
\right)
\left(
\begin{array}{r}
v^{(i)}_x-V_x\\
v^{(i)}_y-V_y
\end{array}
\right)
+
\left(
\begin{array}{r}
V_x\\
V_y
\end{array}
\right)
\end{equation}
\)

と表現することが出来ます。

衝突角度について


衝突角度を簡単に推測してみましょう。
衝突に際し重要なのは、物体と物体が接する点における法線ベクトル\(\mathbf{n}\)だと推測できます。
恐らく、実験系において法線ベクトルを軸として、反転させた入射ベクトル\(-\mathbf{v}^{(i)}_{1,2}\)を\(\pi\)回転させる角度になるでしょう、と予測できます。
これは私の直感なので間違っているかもしれません。

ともあれ、必要な関係式は\(\mathbf{v}^{(i)}_1,\mathbf{v}^{(i)’}_1 \)を軸とした時の、実験系の角度\(\theta_1\)と重心系の角度\(\theta’_1\)との間の関係式です。

衝突後の速度ベクトルは
\(
\begin{align}
\mathbf{v}^{(f)}_1&=R(\theta) (\mathbf{v}^{(i)}_1-\mathbf{V})+\mathbf{V}\\
\mathbf{v}^{(f)}_2&=-R(\theta) (\mathbf{v}^{(i)}_2-\mathbf{V})+\mathbf{V}
\end{align}
\)

で表されていますので、実験系での角度、すなわち\(\mathbf{v}^{(i)}_1\)と\(\mathbf{v}^{(f)}_1\)のなす角度\(\theta_1\)は
\(
\begin{align}
\theta=\text{acos}
\left(
\frac{\mathbf{v}^{(i)}_1\cdot \mathbf{v}^{(f)}_1}{|\mathbf{v}^{(i)}_1||\mathbf{v}^{(f)}_1|}
\right)
\end{align}
\)

となります。
今、弾性散乱を考えているので、必要であれば\(|\mathbf{v}^{(i)}_1|=|\mathbf{v}^{(f)}_1|\)が成り立っています。
なので、
\(
\begin{align}
\theta=\text{acos}
\left(
\frac{\mathbf{v}^{(i)}\cdot [R(\theta) (\mathbf{v}^{(i)}_1-\mathbf{V})+\mathbf{V}]}{|\mathbf{v}^{(i)}_1||\mathbf{v}^{(f)}_1|}
\right)
\end{align}
\)

です。
変形すれば、
\(
|\mathbf{v}^{(i)}_1||\mathbf{v}^{(f)}_1|\cos\theta_1-\mathbf{v}^{(i)}_1\mathbf{V}=\mathbf{v}^{(i)}_1R(\theta’)(\mathbf{v}_1^{i}-\mathbf{V})
\)
という関係式が導けます。
この式は\(\theta\)と\(\theta’\)を結びつける式になっていますので、実験系での角度が分かれば、重心系での角度に焼直せるため、衝突した時の方向も決められるのです(簡単には求められなさそうですが・・・)。

ゼロ点を探す(二分法、挟み撃ち法、Brent法、Newton法、Steffensen法)

ここでは、非線形方程式の解を数値的に求める方法である、
二分法、
挟み撃ち(False position)法、
Brent法、
Newton法、
Steffensen法
について紹介します。

  1. 問題設定
  2. まとめ
  3. 数値解法
  4. 複数のゼロ点を見つけたい場合

問題設定


区間\([a,b]\)で定義された実関数\(f(x)\)に
\(
f(x)=0
\)

を満たす\(x\)が1つある。\(x\)を求めよ。但し\(f(a),f(b)\ne 0\)。

まとめ


1変数の場合:Brent法(二分法と逆2次補間、挟み撃ち法の組み合わせ)
多変数の場合:Newton法
を使うのが良いでしょう。

数値解法


1変数の問題で、
教育上良く使われる方法は二分法
実用的な手法はBrent法
です。
またNewton法は多変数の場合にも複素数の場合でも容易に拡張できるので、その意味で実用的です。

区間\([a,b]\)にゼロ点が1つ存在する事を数学的に言えば、\(f(a)f(b)\lt 0\)ということです。

二分法


要点

・解が確実に見つかります。
・関数の振る舞いが非連続であっても、たちが悪い関数でも確実ですが、収束は遅いです。

計算方法

二分法は位置\([a,b]\)における関数の符号の情報だけを見てゼロ点を推測します。
なので、区間内にゼロ点がある事が分かっていれば、そこを境として符号が違うので確実に解を囲い込んでいくことが出来ます。
教育的に分かりやすい方法であり、根を探すにあたり失敗が無い方法ですが、
滑らかな関数であっても計算時間がかかるため、あまり実用的ではありません。

収束速度について

2分法では、1回の計算あたり範囲が\(1/2\)になるため、漸化式
\(
\displaystyle \varepsilon_n=\varepsilon_{n-1}\cdot 2^{-1}
\)

が成立するため、
\(
\varepsilon_n=\varepsilon_{0}\cdot 2^{-n}
\)

が導けます。ここで\(\varepsilon_{0}\)は初期のxの範囲(\(\varepsilon_0=|b-a|\))を表します。
収束精度と計算回数の関係を決めるためにはnについて解いて、
\(
\begin{align}
\varepsilon_n &=2^{-n} \varepsilon_{0} \\
2^n &=\frac{\varepsilon_0}{\varepsilon_n} \\
n &=\log_2\left(\frac{\varepsilon_0}{\varepsilon_n}\right) \\
\mbox{または}& \\
n &=\frac{\ln\left(\varepsilon_0/\varepsilon_n\right)}{\ln{2}}
\end{align}
\)
より、例えば初期の解の範囲(\(\varepsilon_0\))を1, 解を\(10^{-14}\)まで収束させようとすれば、大体\(n=47\)となり、47回、2分法を繰り返す必要があります。
逆に言えば、どんな状況でも47回繰り返せば相対誤差を14桁収束させることが出来るということです。

収束の判定は
要求精度\(\epsilon_{\text{tol}}\)、機械精度\(\epsilon_{\text{mac}}\)、解の存在範囲\(\delta\)、新たな解の推定値\(x\)、
とすると、
\(
\delta \lt 4\epsilon_{\text{mac}}|x|+\epsilon_{\text{tol}}
\)

もしくは
\(
f(x)=0
\)

で判定します。
この意味は、解の存在範囲を\(\pm 2\epsilon_{\text{mac}}\)の範囲まで特定するか、解の存在範囲が要求精度\(\epsilon_{\text{tol}}\)まで達するかのどちらかが満たされるかで判定します。


Fortranプログラムはこちら。

二分法でやると以下の通り20回必要になります。

$ gfortran bisection.f90
$ ./a.out
      0.5000000000000000     10.0000000000000000
      0.5000000000000000      5.2500000000000000
      2.8750000000000000      5.2500000000000000
      2.8750000000000000      4.0625000000000000
      2.8750000000000000      3.4687500000000000
      2.8750000000000000      3.1718750000000000
      2.8750000000000000      3.0234375000000000
      2.9492187500000000      3.0234375000000000
      2.9863281250000000      3.0234375000000000
      2.9863281250000000      3.0048828125000000
      2.9956054687500000      3.0048828125000000
      2.9956054687500000      3.0002441406250000
      2.9979248046875000      3.0002441406250000
      2.9990844726562500      3.0002441406250000
      2.9996643066406250      3.0002441406250000
      2.9999542236328125      3.0002441406250000
      2.9999542236328125      3.0000991821289062
      2.9999542236328125      3.0000267028808594
      2.9999904632568359      3.0000267028808594
      2.9999904632568359      3.0000085830688477
      2.9999995231628418      3.0000085830688477
      2.9999995231628418      3.0000040531158447
      2.9999995231628418      3.0000017881393433
      2.9999995231628418      3.0000006556510925
   3.0000000894069672    
$

6桁収束するために24回も評価回数が必要です。
計算時間を気にしないのならば良い方法ですが、やはり遅いです。

挟み撃ち法(False position法)



要点

・解が確実に見つかります。
・勢いよくゼロ点を横切っていない場合、収束が遅くなります。

計算方法

二分法との違いはゼロ点の推測に用いる情報量です。
二分法では位置\([a,b]\)における関数の符号だけを見てゼロ点を推測しますが、
挟み撃ち法では位置\([a,b]\)における関数のも見てゼロ点を推測しています。
この挟み撃ち法の推測は、点\((a,f(a)),(b,f(b))\)を結ぶ直線が通るゼロ点して値を推測するのです。

一つ、挟み撃ち法の説明を見ていた時に、解の囲い込み方法が2通りあるようです。
オリジナルのものは片方の点を絶対に変えない方法のようですが、二分法と同じように決めている物もあります。

ここでは、二分法と同じように決めていく後者の方法を載せておきます。
二分法の新しい点を決める所だけ変えればプログラムは完了です。

要求精度を6桁にすると、実行結果は

$ gfortran false.f90
$ ./a.out
      0.5000000000000000     10.0000000000000000
      0.5000000000000000      4.8581311942172727
      2.6310989812608572      4.8581311942172727
      2.6310989812608572      3.1997221817009311
      2.9967206978247356      3.1997221817009311
      2.9967206978247356      3.0000354381074144
      2.9999999998952847      3.0000354381074144
   3.0000000000000000    
$

となり、初めの1回は無視すると6回で収束に至っていることが分かります。早いですね。

Brent法



要点

・解が確実に見つかります。
・基本的には逆二次補間(inverse quadratic interpolation)で根を探索します。
・逆二次補間が原理的に出来ない場合(二つ以上の関数値が同じ場合など)では挟み撃ち法を行います。
・”探索が失敗”した時は二分法に切り替えて根を探します。
・1変数の場合、実用的な方法として推奨されています。

計算方法

逆二次補間は3点を結ぶ逆二次関数(上下に凸ではなく、左右に凸をもつ関数)としてラグランジュ補間をして、その根をゼロ点の位置として用いる方法です。
しかし、逆二次補間を行う際に必要な3点の内、2点以上が同じ関数値を持つと解が発散する、という問題があるため、その時は挟み撃ち法に切り替えて計算を行います。

二分法は、逆二次補間や挟み撃ち法を用いた為に解の収束が遅くなる、と判定されたときに使われます。

実際に動かしてみると逆二次補間よりも挟み撃ち法が使われる頻度が多いようです。

詳細は(Van Wijngaarden-Dekker-Brent法, William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery著, 丹慶勝市ら訳(1993) 『Numerical Recipes in C [日本語版]C言語による数値計算のレシピ』 技術評論社 261pp.)に書いてありますので、こちらを参照してください。

プログラム自体はパブリックドメインのNetlibにあるもの(zeroin.f)を基本としましょう。ただし、Fortran77で掛かれているので並列計算を考えた時に厄介です。文番号を消したものを置いておきます。

要求精度を6桁にすると、実行結果は

$ gfortran main.f90
$ ./a.out
      0.5000000000000000     10.0000000000000000
      0.5000000000000000      4.8581311942172727
      2.6310989812608572      4.8581311942172727
      2.6310989812608572      3.5172292241144740
      2.6310989812608572      3.0069570368988154
      2.9997489497326733      3.0069570368988154
      2.9997489497326733      3.0000000032534229
      2.9999995032534215      3.0000000032534229
   3.0000000032534229    
$

となり、初めの1回は無視すると7回で収束に至っていることが分かります。挟み撃ち法とだいたい同じです。

Netlibがpublicdomainであることの記述は以下の本にも見ることが出来ます。
https://books.google.de/books?id=2GPNBQAAQBAJ&pg=PA109&lpg=PA109&dq=netlib+public+domain&source=bl&ots=d6Uw8D5JMv&sig=Dw3wXTAQThFKLVrL7pm_qnCaq9s&hl=de&sa=X&ved=0ahUKEwjO05Xxt_XKAhXCwxQKHUBlC5sQ6AEIXjAI#v=onepage&q=netlib%20public%20domain&f=false

ちなみにですが、wikipediaにあるBrent法の疑似コードを実装したのですが、うまく働きませんでした。
Brent’s method -wikipedia en
ブレント法 -wikipedia ja

Newton法



要点

・解近傍のただ1点から収束させるので必要な情報が少なく済みます。
・解近傍のただ1点上における関数の微分を用いて収束させ、一回の計算で一致桁数が倍になります。
・本当に解の近くから始めないと、失敗します。
・多変数への拡張が容易ですが、関数の微分を数値微分で代用すると効率が落ちます。

計算方法

ニュートン法について詳しくは
ニュートン法(1、2次元、多次元)
に書いたので、こちらを参照してください。


プログラムは以下のもので、

実際に実行すると

$ gfortran newton.f90
$ ./a.out
      4.0000000000000000
      2.4339000410871483
      3.0980975267667068
      2.9994762813770324
      3.0000000000829825
   3.0000000000000000
$

という結果を得ます。

番外編:Steffensen法(Steffensen’s method)


Steffensen法は自己無撞着方程式を解くために使われる方法で、Newton法に似ています。
特に良く使われるのは、方程式がたまたま
\(
x=g(x)
\)

と書けている特定の方程式を満たす\(x\)を探す場合に使われます。

ですが、
\(
f(x)=x-g(x)
\)

を考えれば
\(
f(x)=0
\)

を考える問題と同じ、ということになります。
ここでは、方程式
\(
x=x^2/2
\)

を満たす\(x\)を探す問題を考えます。

この問題は
\(
x-x^2/2=0
\)

と同じです。

計算手法は
htt://park.itc.u-tokyo.ac.jp/kato-yusuke-lab/nagai/note_070420_self.pdf
でまとめられている通り、

\(
\begin{align}
a&=g(x_n), \\
b&=g(a),\\
x_{n+1}&=x_n-\frac{(a-x_n)^2}{b-2a+x_n}
\end{align}
\)
の漸化式で求められます。

プログラムではこちら。

複数のゼロ点を見つけたい場合



複数のゼロ点を見つけたい場合、一つの方法として、

解のある区間の特定→収束

を繰り返すことで得られる、と考えられるでしょう。

具体的には計算区間を荒く区切り、符号が変わる範囲を特定、その符号が変わる範囲に対してのみゼロ点を探すのです。

上の図の縦線は等間隔に荒く探した時の、関数が正になる時に青、負の時に赤色に表示しています。
収束した結果が緑の点です。

fortranプログラムはこんな感じになります。