「プログラミングと数値計算」カテゴリーアーカイブ

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 で最小二乗フィッティングする -ゴルディアスの涙目

ゼロ点を探す(二分法、挟み撃ち法、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プログラムはこんな感じになります。

ニュートン法(1、2次元、多次元)

ニュートン法に関するお話です。
”ニュートン法”と呼ばれる方法は
2種類(1.初期値近傍の極小値を求めるニュートン法、2. ゼロ点を求めるためのニュートン法)
があるようですが、ここでは2. 関数のゼロ点を求める方のニュートン法についてのお話です。

  1. まとめ
  2. 1次元のニュートン法
    1. 初期値の推定
    2. 導関数の近似
  3. 1次元ニュートン法のプログラム
  4. 2次元ニュートン法
  5. 2次元ニュートン法のプログラム
  6. 多次元ニュートン法
  7. 高次のニュートン法?
  8. 補)ニュートン法が使われる問題
  9. 参考文献

まとめ

1次元の場合

\(\displaystyle
x_{n+1}=x_{n}-\frac{f(x_{n})}{f'(x_{n})},~~n=0,1,2\cdots,
\)

ここで、\(x_0\)は解の近傍であること。

2次元の場合

\(
\begin{align}
\left(
\begin{array}{c}
x_{n+1} \\
y_{n+1}
\end{array}
\right)
=
\left(
\begin{array}{c}
x_n \\
y_n
\end{array}
\right)
-\left(
\begin{array}{cc}
f_x(x_n,y_n) & f_y(x_n,y_n) \\
g_x(x_n,y_n) & g_y(x_n,y_n)
\end{array}
\right)^{-1}
\left(
\begin{array}{c}
f(x_n,y_n) \\
g(x_n,y_n)
\end{array}
\right)
\end{align}
\)

\(
\begin{align}
&\left(
\begin{array}{cc}
f_x(x_n,y_n) & f_y(x_n,y_n) \\
g_x(x_n,y_n) & g_y(x_n,y_n)
\end{array}
\right)^{-1} \\
&~~~~=
\frac{1}{f_x(x_n,y_n) g_y(x_n,y_n)-f_y(x_n,y_n)g_x(x_n,y_n)}
\left(
\begin{array}{cc}
g_y(x_n,y_n) & -f_y(x_n,y_n) \\
-g_x(x_n,y_n) & f_x(x_n,y_n)
\end{array}
\right)
\end{align}
\)

ここで、
\(
\begin{align}
f_x(x_n,y_n) &= \left.\frac{\partial f(x,y)}{\partial x}\right|_{x=x_n,~y=y_n} ,~~& f_y(x_n,y_n) = \left.\frac{\partial f(x,y)}{\partial y}\right|_{x=x_n,~y=y_n} \\
g_x(x_n,y_n) &= \left.\frac{\partial g(x,y)}{\partial x}\right|_{x=x_n,~y=y_n} ,~~& g_y(x_n,y_n) = \left.\frac{\partial g(x,y)}{\partial y}\right|_{x=x_n,~y=y_n}
\end{align}
\)

多次元の場合

\(
\mathbf{x}_{n+1}=\mathbf{x}_{n} – \mathbf{J}_n^{-1}\mathbf{f}_n
\)

ここで、
\(
\begin{align}
\mathbf{x}_n=
\left(
\begin{array}{c}
x_1\\
x_2\\
\vdots \\
x_N
\end{array}
\right)_n
,~~
\mathbf{f}_n=
\left(
\begin{array}{c}
f_1(\mathbf{x}_n)\\
f_2(\mathbf{x}_n)\\
\vdots \\
f_N(\mathbf{x}_n)
\end{array}
\right)
,~~
\mathbf{J}_n=
\left.\left(
\begin{array}{cccc}
\frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_N} \\
\frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \cdots & \frac{\partial f_2}{\partial x_N} \\
\vdots&\vdots&\ddots&\vdots \\
\frac{\partial f_N}{\partial x_1} & \frac{\partial f_N}{\partial x_2} & \cdots &\frac{\partial f_N}{\partial x_1}
\end{array}
\right)\right|_{\mathbf{x}=\mathbf{x}_n}
\end{align}
\)

※\((~~)_n\)の添え字nはn回目の繰り返しにおけるベクトル\(\mathbf{x}\)を表しています。


スポンサーリンク


ここから本文

1次元のニュートン法


方程式
\(
f(x)=0
\)

を満たす\(x=a\)を求めることを考えます。

まず、この問題をグラフ
\(
y=f(x)
\)

のゼロ点(\(y=0\)となるような\(x=a\))を探す問題だ、という問いに置き換えます。

これを、解近傍で関数をテーラー展開で表し、その展開式のゼロ点を探すことで解を近似することで方程式を解きます。

関数\(f(x)\)が与えられたとき、任意の点\(x’\)周りにおけるテーラー展開は、
\(
f(x) = f(x’)+f'(x’)(x-x’) + O(\Delta^2),~~\Delta=|x-x’|
\)

と書くことが出来ます。もしも\(x\)と\(x’\)が近接していれば、\(\delta=|x-x’| \ll 1\)となるので、\(\Delta^2\)の項は無視できるほど小さくなると期待します。

さて、この式を導く際に用いた条件は、\(x\)と\(x’\)が近くにある、という条件のみです。なので、求めたい解\(a\)と\(a\)の近くの適当な点\(x_0\)を用いれば、
\(
\begin{align}
f(a) &= f(x_0)+f'(x_0)(a-x_0) +O(\Delta^2)~~~…(1)\\
f(x_0) &= f(a)+f'(a)(x_0-a) +O(\Delta^2) ~~~…(2)
\end{align}
\)

の2通りの式を考えることが出来ます。(2)に含まれる\(f'(a)\)は本当の\(a\)が分からないと評価するのが難しいので、(1)を選びます。

(1)の左辺は条件よりゼロです。残りの右辺の\(a\)について解けば、

\(
\begin{align}
0 = f(x_0)+f'(x_0)(a-x_0)+O(\Delta^2) \\
\to a = x_0-\frac{f(x_0)}{f'(x_0)}+O(\Delta^2)
\end{align}
\)


ここで、\(f(x_0)=f'(x_0)=O(\Delta^0)\)を想定しているため、\(f'(x_0)\)の割り算を行っても誤差のオーダーは変化せず、\(\Delta\)に対して2乗のままとなります。

\(O(\Delta^2)\)を無視すると、得られる解は近似値となります。この近似値を\(x_1\)と置くと
\(\displaystyle
x_1 = x_0-\frac{f(x_0)}{f'(x_0)}
\)

得られる解\(x_1\)は\(x_0\)よりも解\(a\)に近いため、\(O(\Delta^2)\)はさらに小さくなります。なので、もう一度計算します。この時の近似値を\(x_2\)と置くと
\(\displaystyle
x_2 = x_1-\frac{f(x_1)}{f'(x_1)}
\)

ここで得られた\(x_2\)は更に\(a\)に近いため、・・・と繰り返していきます。

すなわち、解\(a\)に近い初期値\(x_0\)を与えた時、漸化式
\(\displaystyle
x_{n+1}=x_{n}-\frac{f(x_{n})}{f'(x_{n})},~~n=0,1,2\cdots,
\)

を繰り返すことで解\(a\)に近づいていく、ということです。
この式から分かる通り、解近傍で導関数がゼロになる場合はゼロ割発生するかもしれないので危険です。

もしもこのステップが無限回繰り返されれば、無限回目のステップの近似値\(x_\infty\)は
\(
x_\infty=a
\)

に収束します。

以上のように導関数を用いて方程式の解を求める方法は、ニュートン法(またはニュートン=ラフソン法)と呼ばれます。
この方法は、解に近いに所の値とその導関数が分かりさえすればいいので、この方法は1次元、多次元でも使用する事ができ、また複素関数でさえ使用することが出来ます。
さらに、誤差が\(O(\Delta^2)\)で減衰していくため、一度繰り返すごとに厳密解に一致する桁数が2倍になっていきます。

しかし問題があり、
1. 初期値を解の近傍に取らないと失敗する
2. 導関数が分からないと使えない
という大きな問題を抱えています。

また、解近傍でテーラー展開可能でなければならず、解近傍で導関数がゼロになる点があってはいけません。

  • 初期値の推定


    ニュートン法は解の近傍で、テーラー展開可能な領域から始めなければ用いることは出来ません。
    これは工夫をするしかありません。

    解近傍の初期値を推定する方法として、以下の2つの方法が良く取られます。
    1. 解の範囲を限定し、その近傍からニュートン法を行う(1次元のみ)。
    2. 適当な項を落とせば方程式が解けてしまう場合、その解を出発点とし、残りの項を摂動的に加えながら逐次解く。

    1. 解の範囲を限定する

      これは1次元の場合に使える方法です。
      単純な考え方で、例えば広い区間\([x_a,x_b]\)にゼロ点があることだけが分かっているとします。
      この区間を\(N+1\)分割し、
      \(
      x^{(0)}=x_a,x^{(1)},\cdots,x^{(k)},x^{(k+1)},\cdots, x^{(N-1)},x^{(N)}=x_b
      \)

      とします。もしも\(x^{(k)}\)と\(x^{(k+1)}\)の符号が違う場合、この区間でゼロを横切っているということになりますので、この点からニュートン法を始めれば良いのです。

      問題は、解の範囲を狭めるために関数を\(N+1\)回評価しなければならない、という点で計算コストがかかってしまう、という点です。

    2. 既知の解から摂動的に逐次解く

      これは若干特殊な方法で、あてずっぽうな上の考えとは違います。
      本来解きたい式を
      \(
      f(x)=f_0(x)+g(x)
      \)

      と分けます。ここで、
      \(
      f_0(x)=0
      \)

      の解は簡単に解けて、\(x=a_0\)が関数\(f_0(x)\)のゼロ点であるとします。この\(a_0\)を初期値とし、\(\lambda\)を十分小さい値にとり
      \(
      f_0(x)+\lambda_1 g(x)=0
      \)

      を解きます。ここで、\(\lambda_m\)は\(0\lt\lambda_1\lt\lambda_2\cdots \lt\lambda_{M-1} \lt \lambda_M=1\)を満たすとします。
      この方程式の解\(a_1\)が得られたら、この解を初期値とし、
      \(
      f_0(x)+\lambda_2 g(x)=0
      \)

      を解きに行き、これを\(\lambda_M=1\)まで続けます。
      こうすることにより、最後の答え\(a_M\)がもともと解きたかった式の答えとなっているのです。

      関数の一回の評価に時間がかかってしまう場合、この解き方の方が早く解けます。
      しかし、関数\(g(x)\)によってゼロ点の位置が増えたり、消えたりしてしまう時にこの方法は使うことが出来ず、連続的に解が動いてくれないと使うことが出来ない点に注意してください。

  • 導関数の近似


    導関数を求める簡単な方法は、差分を使うことです。すなわち、導関数\(f'(x)\)を差分に置き換えます。
    \(\displaystyle
    f'(x)\approx \frac{f(x+h)-f(x)}{h} + O(h)
    \)

    すると、漸化式は
    \(\displaystyle
    x_{n+1}=x_{n}-h\frac{f(x_{n})}{f(x_n+h)-f(x_n)},~~n=0,1,2\cdots,
    \)

    と得られます。
    この方法は導関数を必要とせず、実装が簡単であるという意味で有用な式です。

    導関数を求める際の前進差分の刻み幅\(h\)は、計算機の扱える丸め誤差\(\epsilon\)(単精度ならば約\(\epsilon=10^{-8}\)、倍精度ならば約\(\epsilon=10^{-16}\))の\(1/2\)乗の2倍、つまり
    \(
    \displaystyle h\sim 2\sqrt{\epsilon}
    \)

    とすると良いです。

    この刻み幅の導出は折りたたんでおきますが、以下のように導出することが出来ます。

参考文献[1][3][5]

1次元ニュートン法のプログラム


以下のプログラムは、1次元のニュートン法のプログラムで、
方程式
\(
x^2-4 = 0
\)

を初期値\(x_0=3\)から始めて解くプログラムです。

このコードを動かすと、ステップを繰り返すごとに急激に解に近づいていっていることが分かります。
(デフォルトではコメントアウトしてあります。)

$ gfortran main.f90
$ ./a.out
    1       2.1666666616021075
    2       2.0064102565311974
    3       2.0000102400381690
    4       2.0000000000262466
    5       2.0000000000000000
   2.0000000000000000       ---result
$

上には出力していませんが、初期値の評価もしています。
なので、ニュートン法は6回呼び出されたことになります。
ニュートン法1回当たり関数は2回呼び出されますので、12回関数が呼び出されたことになります。

例えば比較対象として、二分法を考えますと、1回目に関数は3回,そのほかで1回ずつ呼び出されます。
二分法の解は1ステップ当たり区間の半分になりますので、16桁一致するまでには50回ほど繰り返す必要があります。そのため、52回位関数が呼び出されます。

ニュートン法の方が圧倒的に早いことが分かります。

複素関数の1次元ニュートン法のプログラム
$ gfortran main.f90
$ ./a.out
                    1 ( -1.2933333304320405     , 0.72000000272341036     )
                    2 (-0.78207788228079722     , 0.60930726266833468     )
                    3 (-0.43844290526861285     , 0.73503785880633421     )
                    4 (-0.50851135934742864     , 0.89043164725052448     )
                    5 (-0.50009896196621662     , 0.86666395460253720     )
                    6 (-0.49999991066297655     , 0.86602581130981315     )
                    7 (-0.49999999999985084     , 0.86602540378453374     )
                    8 (-0.50000000000000000     , 0.86602540378443860     )
 (-0.50000000000000000     , 0.86602540378443860     )  ----result
$

ソースコードはこちら。

2次元のニュートン法


続いて方程式
\(
\begin{eqnarray}
\left\{~~
\begin{aligned}
f(x,y)&=0,\\
g(x,y)&=0
\end{aligned}
\right.
\end{eqnarray}
\)

を満たす\(x=a,y=b\)を見つけることを考えます。

1次元の時と同様に、任意の点\(x’,y’\)周りにおけるテーラー展開は、
\(
\begin{eqnarray}
\left\{~~
\begin{aligned}
f(x,y) &= f(x’,y’)+f_x(x’,y’)\cdot(x-x’)+f_y(x’,y’)\cdot(y-y’)+ O(\varepsilon^2)\\
g(x,y) &= g(x’,y’)+g_y(x’,y’)\cdot(x-x’)+g_y(x’,y’)\cdot(y-y’)+ O(\varepsilon^2)
\end{aligned}
\right.
\end{eqnarray}
\)

ここで、
\(\displaystyle \left. f_x(x’,y’)=\frac{\partial f(x,y)}{\partial x}\right|_{x=x’,y=y’}\)
\(\displaystyle \left. f_y(x’,y’)=\frac{\partial f(x,y)}{\partial y}\right|_{x=x’,y=y’}\)
\(\displaystyle \Delta x=|x-x’|=O(\varepsilon^1),~~\Delta y=|y-y’|=O(\varepsilon^1)\)
と書きました。

1次元の時と同様に初期値\((x_0,y_0)\)が解\((a,b)\)に近ければ、
\(
\begin{eqnarray}
\left\{
\begin{aligned}
0 &= f(x_0,y_0)+f_x(x_0,y_0)\cdot(a-x_0)+f_y(x_0,y_0)\cdot(b-y_0) +O(\varepsilon^2)\\
0 &= g(x_0,y_0)+g_x(x_0,y_0)\cdot(a-x_0)+g_y(x_0,y_0)\cdot(b-y_0) +O(\varepsilon^2)
\end{aligned}
\right.
\end{eqnarray}
\)

を同時に満たす\((a,b)\)を探せばよい、ということになります。

行列形式で書けば
\(
\begin{align}
\left(
\begin{array}{cc}
f_x(x_0,y_0) & f_y(x_0,y_0) \\
g_x(x_0,y_0) & g_y(x_0,y_0)
\end{array}
\right)
\left(
\begin{array}{c}
a-x_0 \\
b-y_0
\end{array}
\right)
=
-\left(
\begin{array}{c}
f(x_0,y_0) \\
g(x_0,y_0)
\end{array}
\right)+O(\varepsilon^2)
\end{align}
\)

と書けます。\(a,b\)について解けば
\(
\begin{align}
\left(
\begin{array}{c}
a \\
b
\end{array}
\right)
=
\left(
\begin{array}{c}
x_0 \\
y_0
\end{array}
\right)
-\left(
\begin{array}{cc}
f_x(x_0,y_0) & f_y(x_0,y_0) \\
g_x(x_0,y_0) & g_y(x_0,y_0)
\end{array}
\right)^{-1}
\left(
\begin{array}{c}
f(x_0,y_0) \\
g(x_0,y_0)
\end{array}
\right)+O(\varepsilon^2)
\end{align}
\)

となります。丁度一次元のニュートン法と似た形になりました。

右辺の
\(
\begin{align}
\mathbf{J}_0=\left(
\begin{array}{cc}
f_x(x_0,y_0) & f_y(x_0,y_0) \\
g_x(x_0,y_0) & g_y(x_0,y_0)
\end{array}
\right)
\end{align}
\)

は\(x=x_0,y=y_0\)のヤコビアンです。また、
\(
\mathbf{a}=\left(
\begin{array}{c}
a \\
b
\end{array}
\right),~~
\begin{align}
\mathbf{x}_0=\left(
\begin{array}{c}
x_0 \\
y_0
\end{array}
\right),~~
\mathbf{f}_0=\left(
\begin{array}{c}
f(x_0,y_0) \\
g(x_0,y_0)
\end{array}
\right)
\end{align}
\)

と表記すれば、2次元のニュートン法は
\(
\mathbf{a}=\mathbf{x}_0 – \mathbf{J}_0^{-1}\mathbf{f}_0+O(\varepsilon^2)
\)

と書くことが出来ます。
近似として、\(O(\varepsilon^2)\)の項を無視すれば
\(
\mathbf{x_1}=\mathbf{x}_0 – \mathbf{J}_0^{-1}\mathbf{f}_0
\)

と書けます。繰り返せば、漸化式
\(
\mathbf{x}_{n+1}=\mathbf{x}_{n} – \mathbf{J}_n^{-1}\mathbf{f}_n
\)

を解いていけば良い、ということになります。

ニュートン法を実行するにあたり、ヤコビアンの逆行列を解くことが一番の問題となります。

良く知られているように、2次元の逆行列は、計算コストもさほどかからずに解けてしまい、
\(
\begin{align}
\mathbf{J}_n^{-1}=
\frac{1}{f_x(x_n,y_n) g_y(x_n,y_n)-f_y(x_n,y_n)g_x(x_n,y_n)}
\left(
\begin{array}{cc}
g_y(x_n,y_n) & -f_y(x_n,y_n) \\
-g_x(x_n,y_n) & f_x(x_n,y_n)
\end{array}
\right)
\end{align}
\)

として求める事が出来ます。

参考文献[2],[6]

スポンサーリンク

2次元ニュートン法のプログラム


以下のプログラムは、2次元のニュートン法のプログラムで、
方程式
\(
\begin{eqnarray}
\left\{~~
\begin{aligned}
f_1(x,y)&=x^2+y^2-1 = 0 \\
f_2(x,y)&=y-x^3 = 0
\end{aligned}
\right.
\end{eqnarray}
\)

を初期値\(x_0=2,~y_0=1\)から始めて解くプログラムです。

このコードを動かすと、ステップを繰り返すごとに解に近づいていっていることが分かります。

$ gfortran main.f90
$ ./a.out
    1       1.3571428532359113       0.2857142813732352
    2       0.9844126813492445       0.4401111876910193
    3       0.8485699953676349       0.5590406123395349
    4       0.8265084709901402       0.5633730814110012
    5       0.8260315915076198       0.5636240770681390
    6       0.8260313576542442       0.5636241621612319
    7       0.8260313576541870       0.5636241621612585
  0.82603135765418700       0.56362416216125855       ---result
$

この問題の答えは\(f_1=0,~f_2=0\)の陰関数の交点です。
図示すれば、下図の紫色の経路を辿って解に収束していくことが分かります。

ひとつ、注目しておかなければならない点として、解に近いからと言って、その解に収束するとは限らないということに注意してください。
図中の水色の線は初期値こそ右上の解に近いですが、実際に計算してみるともう一つの解に収束していることが分かります。

適当な初期値を用意し、ニュートン法が収束する点を調べるという問題はニュートン写像と呼ばれる分野です。
力学系とかそういう言葉が出てきます。
私も過去に複素関数でやったことがありますので、リンクだけ載せておきます。
ニュートン写像に現れる綺麗な画像 -シキノート

2次元の複素関数のニュートン法

置いておきます。

多次元のニュートン法

さて、この形まで持ってくると多次元への拡張が簡単にできます。
多次元のニュートン法は
\(
\mathbf{x}_{n+1}=\mathbf{x}_n – \mathbf{J}_n^{-1}\mathbf{f}_n
\)

の形がそのまま使えます。
ここで、
\(
\begin{align}
\mathbf{x}_n=
\left(
\begin{array}{c}
x_1\\
x_2\\
\vdots \\
x_N
\end{array}
\right)_n
,~~
\mathbf{f}_n=
\left(
\begin{array}{c}
f_1(\mathbf{x}_n)\\
f_2(\mathbf{x}_n)\\
\vdots \\
f_N(\mathbf{x}_n)
\end{array}
\right)
,~~
\mathbf{J}_n=
\left.\left(
\begin{array}{cccc}
\frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_N} \\
\frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \cdots & \frac{\partial f_2}{\partial x_N} \\
\vdots&\vdots&\ddots&\vdots \\
\frac{\partial f_N}{\partial x_1} & \frac{\partial f_N}{\partial x_2} & \cdots &\frac{\partial f_N}{\partial x_1}
\end{array}
\right)\right|_{\mathbf{x}=\mathbf{x}_n}
\end{align}
\)

です。

問題となるのは(ヤコビアンの逆行列)×(f)の計算です。
この計算を数値計算的に何とかするべく改良が重ねられてきたのが準ニュートン法(本稿では書きません)と呼ばれるアルゴリズムです。

逆行列を数値的に計算することは宜しくないですので、
\(
\begin{align}
\mathbf{r}&=\mathbf{J}_n^{-1}\mathbf{f}_n \\
\to \mathbf{J}_n \mathbf{r} = \mathbf{f}_n
\end{align}
\)

という方程式を解く問題に帰着させます。これが解ければ、ニュートン法の次のステップを
\(
\mathbf{x}_{n+1}=\mathbf{x}_n – \mathbf{r}
\)

として求める事が出来るのです。

高次のニュートン法?


1次元のニュートン法は
\(\displaystyle
a = x_0-\frac{f(x_0)}{f'(x_0)}+O(\Delta^2)
\)

の形で求めていきますが、二階微分が分かっていればもっと収束が早くなりそうです。
実際、
\(\displaystyle
a = x_0-\frac{f^{\prime}(x_0)}{f^{\prime\prime}(x_0)}\pm\sqrt{{f’}^2(x_0)-2f(x_0)f^{\prime\prime}(x_0)}+O(\Delta^3)
\)

と解けます。しかし、2次関数ですので
解が見付からないかもしれない
という可能性があります。1次のニュートン法であれば、導関数がゼロでない限り、必ず\(y=0\)と交わる点が見付かります。しかし、2次で解に近づけると交点が見つからないことがあるかもしれません。

なので、1次のニュートン法を繰り返して使う方が(2次と比べたら)安定な解法なのです。

参考ページ

[1]1 Newton 法:壱変数の場合
https://www.astr.tohoku.ac.jp/~chinone/Planck/Planck-node6.html
[2]2 Newton 法:弐変数の場合
https://www.astr.tohoku.ac.jp/~chinone/Planck/Planck-node7.html
[3]Newton法
http://www.misojiro.t.u-tokyo.ac.jp/~murota/lect-suchi/newton130805.pdf
[4]アルゴリズムによる誤差
http://www.aoni.waseda.jp/ykagawa/chap2html/node9.html
[5]Newton 法による方程式の近似解法
http://www.math.u-ryukyu.ac.jp/~suga/C/2004/7/node9.html
[6]多次元のニュートン・ラフソン法
http://www.slis.tsukuba.ac.jp/~fujisawa.makoto.fu/cgi-bin/wiki/index.php?%C2%BF%BC%A1%B8%B5%A4%CE%A5%CB%A5%E5%A1%BC%A5%C8%A5%F3%A1%A6%A5%E9%A5%D5%A5%BD%A5%F3%CB%A1
[7]山本 有作, 行列計算における高速アルゴリズム
http://www.cms-initiative.jp/ja/events/0620yamamoto.pdf
[8]中田 和秀, 大規模線形方程式を解くためのクリロフ部分空間法の前処理
http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/1288-6.pdf

微分方程式を解く際の計算量の工夫

運動方程式を解いて、ある特定の時刻の値が欲しいとします。
特定の時刻がいつなのか、あらかじめ分かっている場合、与えられた初期条件から特定の時刻まで計算し、データを貯めておいて後で使えば良いです。

特定の時刻が分からない時に計算量の問題が生じます。
例えば、時刻\(t=0\)の初期条件が与えられた時、時刻\(t_1=100,t_2=101, t_3=102\)の時刻が知りたいとします。ここで、時刻\(t_2\)の値は実は\(t_1\)の関数の値を評価してから初めて決まるとします。

この場合の計算方法としては、

  1. 毎回、時刻\(t_0\)から初めて\(t_1,t_2,t_3\)まで計算
  2. あらかじめ微分方程式を一度細かく計算しておき、欲しい時刻\(t=t’\)を補間で求める
  3. あらかじめ微分方程式を一度細かく計算しておき、欲しい時刻\(t=t’\)に近い点から計算

だと思います。最も簡単に実装できるのは1の方法で、欲しい時刻が数点であれば計算時間は問題にならないほど少なくなるでしょう。2の方法は精度を気にする場合には使用できません。ここでは3番目の方法を考えます。

計算方法を図示したものが以下のものです。

初期条件からいちいち出発して\(t=t’\)に到達するのではなく、あらかじめ計算していた参照点から出発して計算する、ということです。
ここで問題なのが、もし横軸が時刻だった場合、時間を逆向きに遡って良いか?という問題です。

これは物理の問題で言えば、解こうとしている問題は時間可逆か?ということで、言い換えれば系のエネルギーは保存しているか?ということです。もしも抵抗力など、エネルギーが散逸する場合は、いくら参照点が近いからといって逆向きに進んではいけないということです。

ここでは、微分方程式を解く際に戻ることはせず、必ず求めたい時刻よりも前の参照点からスタートする、というプログラムを作ります。

以下のように実装できます。

シンプルQUADPACK

QUADPACKのコード(http://www.netlib.org/quadpack/)は優秀ですが、2000行近くあって色々いじる際に面倒です(最速・高精度の数値積分)。

コードを覗くと、求積法の点数だけが違うプログラムがあるので長いことに気が付きます。
なのでコード量を減らしましょう。

使用する積分方法を20-41点ガウス求積法(key=4)に限定しました。key=4を選んだ理由は、最速・高精度の数値積分の結果から推測しました。

key=2の低次の方法(10-21点)であると滑らかな関数でも評価回数が多くなる反面、不連続点があっても少なく済みます。
key=6の高次の方法(30-61点)であると滑らかな関数で評価回数が少なくなる反面、不連続点があると多くなります。

QUADPACKにはkey=1,2,3,4,5,6のどれかを選べますが、基本的に滑らかな関数の積分を行うことが多いので若干高次のkey=4を選びました。

350行程度にまで減らし、ルーチンを一つ無くました。変更箇所が怖い方は、オリジナルのものを使用してください。

コードの変更をした結果、計算結果が違くなるなど問題が起こる場合、原因は私にあります。
著作者とは関係ありませんので、その点はよろしくお願いします。

——————————————

実軸上、複素関数でkey=4の場合

——————————————

key=4, 複素平面上の積分

——————————————
Key=1の場合
(7-15点ガウス求積法)
・不連続性や、端点特異性がある時に有利です。

exp(-ikx)/x の無限区間に渡る積分を数値的に計算したい。

問題:以下の積分を数値的に計算せよ。
\( \displaystyle
-\frac{1}{2\pi i}\int_{-\infty}^{\infty} \frac{e^{-ikx}}{x}dx,~~(kは任意の実数)
\)


1/xのフーリエ変換である
\( \displaystyle
\int_{-\infty}^{\infty} \frac{e^{-ikx}}{x}dx,~~(kは任意の実数)
\)

数値的に計算します。

この積分値は一意には決まりません。
複素関数論によると、この積分値は\(x=0\)周りに存在する特異点を上周りに迂回して積分するか、下周りに迂回して積分するか、で答えが異なります。

特異点を上周りに回る場合、答えは
\(
\begin{eqnarray}
\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}\frac{e^{-ikx}}{x+i\varepsilon}dx =
\left\{
\begin{aligned}
-2\pi i~~~(k\gt 0)\\
0~~~(k\lt 0)
\end{aligned}
\right.
\end{eqnarray}
\)
であり、特異点を下周りに回る場合、
\(
\begin{eqnarray}
\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}\frac{e^{-ikx}}{x-i\varepsilon}dx =
\left\{
\begin{aligned}
0~~~(k\gt 0)\\
2\pi i~~~(k \lt 0)
\end{aligned}
\right.
\end{eqnarray}
\)
となります。
私が疑問に思うのは、積分すると本当にステップ状の関数が現れるのか?ということです。にわかには信じられません。

ここでは、上を回った時の計算にのみ注目します。
答えに出てくる係数を減らしたいため、被積分関数
\( \displaystyle
-\frac{1}{2\pi i}\int_{-\infty}^{\infty} \frac{e^{-ikx}}{x}dx
\)

を考えます。

数値計算方法


さて、問題を整理しましょう。今取り組む問題は、

問題:以下の積分を数値的に計算せよ。
\( \displaystyle
-\frac{1}{2\pi i}\int_{-\infty}^{\infty} \frac{e^{-ikx}}{x}dx,~~(kは任意の実数)
\)

です。読み取れる情報は

・始点は、実軸上のマイナス無限大
・終点は、実軸上のプラス無限大
・\(x=0\)に特異点が存在する

という情報です。これらは理由なく勝手に変えてはいけません。
被積分関数が収束するならば、xがどんなに大きくても無限遠まで振動していきますが、フレネル積分のように積分結果は収束することを期待しましょう。

さて、\(x=0\)に特異点があるために、実軸上のマイナス無限からプラス無限まで積分していこうとすると\(x=0\)を通過しようとする時に無限大になるため、オーバーフローします。
ということで、積分経路を特異点周りで迂回します。

始点と終点は決められていますが、積分経路の指定はありませんので何も問題ありません。
積分経路は以下の図のように取ります。

このようにとると、積分経路上に特異点はないですので計算することが出来ます。
始点と終点は十分大きい値、ここでは(±2000,0)という値にします。

横軸をkとして、積分を計算すると以下のようになります。

複素関数論から、特異点を複素平面上方を迂回した場合に得られる解は
\(
\begin{eqnarray}
\lim_{\varepsilon\to +0} -\frac{1}{2\pi i}\int_{-\infty}^{\infty}\frac{e^{-ikx}}{x+i\varepsilon}dx =
\left\{
\begin{aligned}
1~~~(k\gt 0)\\
0~~~(k\lt 0)
\end{aligned}
\right.
\end{eqnarray}
\)
です。数値計算によって得られた解は、\(k\lt 0\)で値0, \(k\gt 0\)で値1であるので、数値的にも正しそうです。
しかし、厳密にそれぞれの値になっているわけではないことに気が付きます。なぜならば、\(k=0\)周りでは値は収束していないように見えます。
この原因は始点と終点としてとった、大きい値(±2000,0)が十分ではない事を意味しています。kが小さいと、被積分関数の\(\exp(-ikx)\)による振動はゆっくりになるため、収束が遅くなっている、と考えられます。

ジョルダンの補題の利用


この振動を無くすためには被積分関数が漸近領域で収束するように積分経路を変えればいいのです。

始点と終点は変わりません。図は\(k\)が負の時の経路です。

経路の取り方は無数にあり、上の経路は一つ前の計算で選んだ経路を出来る限り変更しないように経路をとっているだけにすぎません。

ジョルダンの補題を利用すれば、被積分関数が収束する領域であれば積分の寄与はゼロなので、数値的に積分する場合の視点、終点を複素平面上へ移すことが出来ます。
この考えの下実際に計算すると積分結果はこうなります。

綺麗なエッジが作られました。これならば解析解と数値解が一致した、と言えるでしょう。複素平面上の数値計算は一筋縄ではいかないことが実感できました。

数値計算プログラム


プログラムは以下のものを用いています。

数値積分はquadpackのガウス=ルジャンドル求積法の自動積分コードを複素関数用に変更して用いています。
実際に使用しているquadpackコードは以下のもので、
http://slpr.sakura.ne.jp/qp/supplement_data/quadpack_cqag.f90
これと、下のはmainプログラムです。

program main
  implicit none
  integer::ier,Nt,i,j,sn
  double precision::eps
  complex(kind(0d0))::a,b,s,s0
  complex(kind(0d0))::xa(1:5),xb(1:5)
  complex(kind(0d0)),external::g
     
  eps=1d-10
  ier=1

  Nt=1000
  do i=1,Nt
     t=i*(1d0-(-1d0))/dble(Nt)-1d0
     s=dcmplx(0d0,0d0)
     
     sn=1
     if(t.gt.0d0)sn=-1
     if(t.eq.0)sn=0
     
     xa(1)=dcmplx(-2000d0,sn*2000d0)
     xb(1)=dcmplx(-1d0,0d0)

     xa(2)=dcmplx(-1d0,0d0)
     xb(2)=dcmplx(-1d0,1d0)

     xa(3)=dcmplx(-1d0,1d0)
     xb(3)=dcmplx( 1d0,1d0)

     xa(4)=dcmplx( 1d0,1d0)
     xb(4)=dcmplx( 1d0,0d0)

     xa(5)=dcmplx( 1d0,0d0)
     xb(5)=dcmplx( 2000d0,sn*2000d0)  

     do j=1,5
        s0=dcmplx(0d0,0d0)
        call cqag_sk(g,xa(j),xb(j),eps,s0,3,ier)
        s=s+s0
     enddo
     write(10,*)t,dble(s),dimag(s)
  enddo  
     
  stop
end program main

function g(z)
  implicit none
  complex(kind(0d0))::g
  complex(kind(0d0)),intent(in)::z

  double precision::pi=dacos(-1d0)
 
  g=exp(-dcmplx(0d0,1d0)*z*t)/(z)
  g=g*(-1d0/(2d0*pi*dcmplx(0d0,1d0)))

  return
end function g

オリジナルのquadpackコードは
http://www.netlib.org/quadpack/(パブリックドメイン)です。

時間依存シュレーディンガー方程式の数値解法 クランク=ニコルソン法

解きたい問題は以下の形の時間依存シュレーディンガー方程式です。

これをクランクニコルソン法という、微分を差分で近似する数値的な時間発展方法で求めていきます。

解法


※途中、式番号が抜けたり(例えば式15)、前後していますが仕様です。

時間を初期時刻\(t=0\)から\(t=\Delta t\)まで時間発展したいと考えます。
微分を差分で離散化しますが、時刻\(t=\frac{\Delta t}{2}\)だけ進んだ時刻からの前進差分をまず考えます。
すると、

と導けます。時間,空間に関して\(O(\Delta t),O(\Delta x^2)\)です。

続いて、時刻\(t=\Delta t\)から\(\frac{\Delta t}{2}\)との後進差分をまず考えますと

です。時間,空間に関して\(O(\Delta t),O(\Delta x^2)\)です。

式(1)と式(2)を足して2で割りますと、時間に関して中央差分を取っていることになりますので、次数が上がります。代入しますと、

となります。時間,空間に関して\(O(\Delta t^2),O(\Delta x^2)\)です。

左辺に時刻\(t+\Delta t\)の波動関数をまとめますと、

となります。これで離散化が出来ました。これがクランクニコルソン法です。

式が複雑なので、少し簡単にまとめます。
位置\(x\)を等間隔に離散化して、以下の通り離散化と記述を変えます。

この表記に倣いますと、離散化した時間依存シュレーディンガー方程式は簡潔に書けて

と書けます。
ここで、係数\(a,b,c,r\)は式(5)の各係数をまとめたものですが、境界条件によって多少異なります。
ここでは、固定端条件と周期境界条件の場合を特に考えてみましょう。

固定端条件の場合


固定端条件は両端での関数の値が時間によって変化しないとする条件です。数式で書けば

という条件が波動関数に課されます。この時、行列表式で式(8)を記述すれば、

と書き表すことが出来ます。ここで、\(\psi_0^{(+)},\psi_N^{(+)}\)に掛かる係数は境界条件を満たすように決定しています。行列の演算で知られているように、第1行目を\(a_1\)倍して第2行目を引き、第N行目を\(c_{N-1}\)倍して第N-1行目を引けばオレンジで四角く囲った枠内のみを計算すればよいことが分かります。
各々の係数の明らかな形は以下の通り得られます。

周期境界条件の場合


周期境界条件は両端での関数の値が等しく、その微分も等しいという条件です。数式で書けば

という条件が波動関数に課されます。この時、行列表式で式(8)を記述すれば、

と書き表すことが出来ます。ここで、\(\psi_N^{(+)}\)に掛かる係数は境界条件を満たすように決定しています。行列の演算で知られているように、第N行目は既に解けているのでこの要素を考える必要はありません。なので、オレンジで四角く囲った枠内のみを計算すればよいことが分かります。
各々の係数の明らかな形は以下の通り得られます。


スポンサーリンク


三重対角行列の方程式の解法


さて、あとは式(9),(13)を数値的に解けば良いのです。
固有値問題ではありませんので、連立方程式を解けば良いのです。

固定端条件の式(9)を解くためのアルゴリズムはThomasのアルゴリズムと呼ばれる方法が一般的です。計算量のオーダーは\(O(N)\)です。

周期境界条件を含んだ三重対角行列の解法はLU分解を利用した方法[2]や、Sherman-morrison formulaという公式を利用した方法[3]があります。ここでは、擬コードが載せられている[2]を元に作成していきます。これも計算量のオーダーは\(O(N)\)ですが、固定端条件の場合よりも約二倍ほど計算量が多くなります。

メモとして書いておきますが、周期境界条件を含む5重対角行列の解法もありまして、[1]に見ることが出来ます。

フォン・ノイマンの安定性解析


これは、クランク=ニコルソン法のような線形変微分方程式を有限差分法で解く際に、数値的に安定なのか発散するのかを調べる手法です。
注意しなければならないのは、数値的に安定なので答えも正しい、という事ではない事です。
この解析により、有限差分を行う際に最低限満たさなければいけない空間の刻み幅と時間の刻み幅の条件が決定されます。

\(V(x)=0\)を考えます。
ざっくりと説明しますと、偏微分方程式における独立な解、すなわち固有の波数\(k\)で特徴づけられる解が発散しない条件から導かれます。
ある時刻の解が1ステップの時間発展をした時、複素数の増幅因子(amplification factor) \(\xi\)を伴って

\(
\begin{align}
\psi_j &= \xi^0 e^{ikj\Delta x}\\
&\downarrow \\
\psi^{(+)}_j &= \xi^1 e^{ikj\Delta x} \\
\end{align}
\)

で時間発展する、とします。時間発展によって、複素数の整数乗が掛かるのです。
もしも増幅因子の絶対値、すなわち\(|\xi|\)が1より小さければ、時間発展しても独立な解は収束するため、数値的に安定だと言えます。
しかし、\(|\xi|\)が1より大きければ、独立な解の時間発展が指数関数で発散していく事を意味するため、数値的に不安定になります。

この\(\xi\)を求めるためには、波数\(k\)をもつ固有モードの式を有限差分で書き表した式(ここでは式(5))に代入すれば求められます。kはいかなる値を持ち得ると考えます。

計算の結果、\(\xi\)は以下の通り求められます[3]。
\(
\begin{align}
\xi=\frac{1-2\beta \sin^2\left(\frac{k\Delta x}{2}\right)}{1+2\beta \sin^2\left(\frac{k\Delta x}{2}\right)} \\
\beta=\frac{\alpha}{i}\frac{\Delta t}{(\Delta x)^2}
\end{align}
\)

※検算はしていません。[3]を参照しました。

この式から分かることは、位置の間隔\(\Delta x\)が与えられたとき、いかなる\(\Delta t\)であっても\(|\xi|\lt 1\)を取るので、数値的に安定である、ということが分かります。

これはクランク=ニコルソン法(陰解法)の特徴となっており、数値的には安定な計算方法です。

詳細は[3]や[4]をご覧ください。

数値計算コード


数値計算コードは以下のように書けます。
以下のプログラムは下の例2を計算するプログラムです。
実用上は、固定端、周期境界のどちらかだけを使えばいいので、適宜消してください。

スポンサーリンク

実例


自由粒子

初期状態:ガウシアン型, 運動量を持つ波束が時間とともにどう発展していくか。
\(
\begin{align}
i\frac{\partial }{\partial t}\psi(x,t)&=-\frac{1}{2}\frac{\partial^2 }{\partial x^2}\psi(x,t) \\
\psi(x,t=0)&=\left(\frac{1}{\pi a^2}\right)^{1/4}\exp\left[-\frac{x^2}{2a^2}+ik_0 x\right]
\end{align}
\)

\(
\begin{align}
\psi(x,t)&=\frac{1}{\pi^{1/4}}\left(\frac{1}{a+it/a}\right)^{1/2}
\exp\left[-\frac{1}{2}\frac{(x-k_0 t)^2}{a^2+it}+i k_0 x -i\frac{k_0^2}{2} t\right]\\
|\psi(x,t)|^2&=\frac{1}{\pi^{1/2}}\left(\frac{1}{a^2+(t/a)^2}\right)^{1/2} \exp\left[-\frac{(x-k_0 t)^2}{a^2+(t/a)^2}\right]
\end{align}
\)


上:存在確率密度\(|\psi(x,t)|^2\)。黒→解析解、青→固定端条件の数値解、赤→周期境界条件の数値解
下:解析解との相対誤差\(|(|\psi(x,t)|^2-|f(x,t)|^2)/|f(x,t)|^2|\)。

gnuplotコードはこちら

調和振動子の固有状態の重ね合わせ

初期状態:調和振動子の基底状態+第一励起状態
\(
\begin{align}
i\frac{\partial }{\partial t}\psi(x,t)&=\left[-\frac{1}{2}\frac{\partial^2 }{\partial x^2}+\frac{1}{2}x^2\right]\psi(x,t) \\
\psi(x,t=0)&=\frac{1}{\sqrt{2}}\left[\frac{2^{1/2}}{\pi^{1/4}}xe^{-x^2/2}+\frac{1}{\pi^{1/4}}e^{-x^2/2}\right]
\end{align}
\)

\(
\displaystyle \psi(x,t)=\frac{1}{\sqrt{2}}\left[\frac{2^{1/2}}{\pi^{1/4}}xe^{-x^2/2}e^{-i\frac{3}{2}t}+\frac{1}{\pi^{1/4}}e^{-x^2/2}e^{-it/2}\right]
\)

gnuplotコードはこちら


上:存在確率密度\(|\psi(x,t)|^2\)。黒→解析解、青→固定端条件の数値解、赤→周期境界条件の数値解
下:解析解との相対誤差\(|(|\psi(x,t)|^2-|f(x,t)|^2)/|f(x,t)|^2|\)。

参考文献


[1]Tomohiro Sogabe, New algorithms for solving periodic tridiagonal and periodic
pentadiagonal linear systems, Appl. Math. Comput. 202 (2008) 850-856

[2]A. A. Karawia, A computational algorithm for solving periodic pentadiagonal linear systems, Appl. Math. Comput. 174 (2006) 613–618

[3]Van Wijngaarden-Dekker-Brent法, William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery著, 丹慶勝市ら訳(1993) 『Numerical Recipes in C [日本語版]C言語による数値計算のレシピ』 技術評論社

[4]フォン・ノイマンの安定性解析 -wikipedia

数値積分の種類

ある関数の数値積分を考えます。
解析的に積分を行う事をまず考えましょう。1回の関数の評価で、全桁一致するからです。

  1. 解析的に計算
  2. どんな値でも関数値が得られる場合
  3. 特定の位置の値だけ関数値が得られる場合

解析的に計算


例えば,すぐには思いつかなそうな積分
\(
\displaystyle \int \sqrt{x}e^{-x^2} dx
\)

http://www.wolframalpha.com/input/?i=integral+sqrt(x)*exp(-x%5E2)
でも、wolfram alphaで計算させますと、不完全ガンマ関数を用いることで計算することが出来ることが分かります。

また、
森口繁一, 宇田川銈久, 一松信著『数学公式Ⅰ』、『数学公式Ⅱ』、『数学公式Ⅲ』岩波書店(1960)
にwolframでも答えてくれないような一般的な関数の積分も載っていますので、こちらも参照しましょう。

残念ながら解析的に積分出来ない場合、数値積分を行います。
ほとんどの数値積分法は、関数\(f(x)\)の積分を
\(
\displaystyle \int f(x) dx \approx \sum_{i} w_i f(x_i)
\)

と右辺の形で近似します。
ここで右辺に現れる、
あらかじめ決められた分点の位置\(x_i\)と
あらかじめ決められた重み\(w_i\)
をどのように決めるのか?で数値積分の良さが決まります。
上の例にあてはまらないものとして乱数を利用するモンテカルロ積分が挙げられます。

数値積分は大きく2つに分けられます。

1. どんな位置の値でも関数値が得られる場合
2. 特定の位置の値だけ関数値が得られる場合

です。

どんな位置の値でも関数値が得られる場合


とても簡単に実装、1-2桁程度合っていれば良い場合


モンテカルロ積分が簡単で、被積分関数が非連続性を持っていても、積分が存在すれば求めることが出来ます。
欠点としては、再現性が乱数に依存してしまう点でしょう。シード値を記録して置けばよいですが…。
良い乱数の発生方法として、メルセンヌツイスタと呼ばれる乱数の発生方法が良いそうです。

一様乱数を用いて1次元積分を行う場合、
\(
\displaystyle \int_a^b f(x)dx \approx \frac{b-a}{N}\sum_{i=1}^N f(x_i)+O(\frac{1}{\sqrt{N}})
\)

として計算できます。ここで\(x_i\)は位置\([a,b]\)に渡って発生される一様乱数です。
コードでは積分値をsとすると

s=0d0
do i=1,N
   x=([0,1)の範囲で発生する乱数)*(b-a)+a
   s=s+f(x)
enddo
s=s*(b-a)/dble(N)

で実装できます。

多重積分の場合、通常の積分法の収束が指数関数的に遅くなるので、モンテカルロ積分は良く使われます。
詳細はモンテカルロ積分をご覧ください。

モンテカルロ法 (Monte Carlo method)
確率密度関数からモンテカルロ積分まで

簡単に実装、2-4桁程度合っていれば良い場合


台形則による積分が良いでしょう。
区間\([a,b]\)で積分したい場合、区間を\(N+1\)点で分割すると
・等間隔の分点の場合
\(
\displaystyle \int_a^b f(x) dx \approx \frac{b-a}{N}\left(\frac{1}{2}[f(a)+f(b)]+\sum_{i=1}^{N-1} f(i\frac{b-a}{N}+a)\right)
\)

・非等間隔の分点の場合
\(
\displaystyle \int_a^b f(x) dx \approx
\frac{1}{2}\left[(x_1-x_0)f(x_0)+(x_N-x_{N-1})f(x_N)
+\sum_{i=1}^{N-1} (x_{i+1}-x_{i-1})f(x_i)\right]
\)

と求める事が出来ます。ここで、\(x_0=a,~x_N=b\)です。

s=0d0
do i=1,N
   c=1d0
   if(i.eq.1.or.i.eq.N)c=0.5d0
   x=(i-1)*(b-a)/dble(N-1)+a
   s=s+c*f(x)
enddo
s=s*(b-a)/(N-1)

3次元での台形則は以下のようなプログラムで実装できます。

s=0d0
do ix=1,Nx
   cx=1d0
   if(ix.eq.1.or.ix.eq.Nx)cx=0.5d0
   x=(ix-1)*(xb-xa)/dble(Nx-1)+xa
   do iy=1,Ny
      cy=1d0
      if(iy.eq.1.or.iy.eq.Ny)cy=0.5d0
      y=(iy-1)*(yb-ya)/dble(Ny-1)+ya
      do iz=1,Nz
         cz=1d0
         if(iz.eq.1.or.iz.eq.Nz)cz=0.5d0
         z=(iz-1)*(zb-za)/dble(Nz-1)+za
         s=s+cx*cy*cz*f(x,y,z)
      enddo
   enddo
enddo
s=s*((xb-xa)/dble(Nx-1))*((yb-ya)/dble(Ny-1))*((zb-za)/dble(Nz-1))

3,4桁以上の計算精度が欲しい時、台形則では間に合いません。
刻み幅を小さくすると精度は確実に上がりますが、刻み幅を1/10すると精度が1桁上がるだけなので、なかなか収束しません。
しかし、台形則も馬鹿にはできません。例えば、周期関数の1周期に渡る積分や、両端で定数に近づいていく積分の場合に、台形則はとんでもない精度をたたき出します(たしかガウス求積法と同じ精度まで上がると思いましたが、詳しい文献が見つけられないので参照はありません)。
この考えを元に考え出されたのが二重指数関数型数値積分です。端点特異性がある場合に特に強いため、学んでおいて損は無いでしょう。


スポンサーリンク


積分する関数に性質が見当たらないが、高い精度が欲しい場合


関数の形によって分点の取り方を自動的に決めてくれる適応型の自動積分が良いです。この種の積分で有名なのは、

NetlibにあるQUADPACK
http://www.netlib.org/quadpack/
最速・高精度の数値積分, QUADPACKプログラム(1/2次元, 実/複素積分) -シキノート)、
もしくは
名古屋大学大型計算機センター、Numpacにあるaqnn9d
http://www2.itc.nagoya-u.ac.jp/center/ja/numpac/index.html
を参照してください。

QUADPACKはガウス-ルジャンドル求積法とガウス-クロンロッド求積法による適応型自動積分であり、
aqnn9dはニュートン・コーツ9点則に基づく適応型自動積分です。
おおよそ単精度ではaqnn9dの方が関数の評価回数が少なくて済み、倍精度以上ではQUADPACKの方が少なくて済みます。
この単精度/倍精度の違いは積分アルゴリズムそのものの特性です。

自動積分について、頒布されているfortranコードの良し悪しを最速・高精度の数値積分で比較したので、詳しいことを知りたい方はそちらをご覧ください。

また、刻み幅制御されたルンゲクッタ法による積分も良いでしょう。
定積分
\(
\displaystyle g(x)= \int_a^x f(x’)dx’
\)

を行うには、初期条件\(g(a)=0\)の下で、
\(
\displaystyle \frac{d g(x’)}{dx’} = f(x’)
\)

を点\(x\)まで上の微分方程式を解けば良いのです。

滑らかな関数(※1)であることが分かっている場合


積分区間内に特異性(※2)は無いとします。
この場合はチェビシェフ多項式による展開に基づくクレンショウ-カーティス(Clenshaw–Curtis)積分が優秀です。
定積分を
\(
\begin{align}
\displaystyle \int_{-1}^1 f(x) dx &\sim \sum_{k=0}^{n} \omega_k f(x_k)\\
x_k&=\cos\left(\frac{k\pi}{N}\right),\;\;\;k=0,1,\cdots,N\\
\omega_0 &=\omega_N =\frac{1}{N^2 -1} \\
\omega_s &=\omega_{N-s} = \frac{4}{n^2}\sum^{n/2}_{j=0}{\prime\prime }\frac{1}{1-4j^2}\cos\left(\frac{2j\pi s}{n}\right),\;\;\;s=0,1,\cdots,\frac{N}{2}
\end{align}
\)
ここで、\(\sum”\)は和の最初と最後の項に\(\frac{1}{2}\)を掛けることを意味します。

クレーンショー・カーチス数値積分則 Ooura’s Mathematical Software Packages
または
Clenshaw–Curtis求積法 シキノート
最速・高精度の数値積分
をご覧下さい。

端点に特異性がある場合


計算区間内の途中に特異点は無く、計算区間の端点に特異性がある場合、変数変換型の数値積分である二重指数関数型数値積分が有効です。
具体的には例えば
\(
\displaystyle \int_0^1 \sqrt{x} dx
\)


\(
\displaystyle \int_{-1}^{1} x^{-0.8} dx
\)

の事です。これらの積分を行う際に二重指数関数型数値積分は非常に少ない分点数で非常に高精度の結果を返します。

重みと分点位置は以下の通り決められます。積分を
\(
\begin{align}
\int_{-1}^{1} f(x) dx &\approx \sum_{i=-N_-}^{N^+} w_i f(x_i),\\
x_i&=\tanh\left(\frac{\pi}{2}\sinh(ih)\right), \\
w_i&=\frac{h\frac{\pi}{2}\cosh(ih)}{\cosh^2\left(\frac{\pi}{2}\sinh(ih)\right)}
\end{align}
\)
として近似します。\(N−,N+\)は離散化誤差と打ち切り誤差が等しくなるように決められます。
\(N−,N+\)はプログラム上では、\(x_i\)がコンピュータの扱える桁数を超えず、\(w_i\)がアンダーフローを起こさない範囲で決められます。

詳しくは二重指数関数型数値積分 -シキノート
をご覧ください。
プログラム自体は大浦様による二重指数関数型数値積分公式の方が優秀なので、こちらを参照すると良いでしょう。

被積分関数が(特定の関数関数)×(多項式)の場合


被積分関数\(f(x)\)が(重み関数\(\omega(x)\))×(多項式\(g(x)\))で記述される特別な場合、ガウス求積法による数値積分が有効です。

多項式\(g(x)\)とは\(x^n\)(nは整数)の足し算で書ける関数です。例えば\(x^4+3x^2+7x\)とかのことです。\(x^2+2\sqrt{x}\)とか、\(x\)の0以上の整数乗ではない項があってはいけません。

重み関数\(\omega(x)\)とは、特定の形をした関数の事です。例えば

\(\omega(x)=1~\to\) ガウス-ルジャンドル求積法、ガウス-ロバート求積法、ガウス-radau求積法
\(\omega(x)=x^n e^{-x}~\to\) ガウス-ラゲール求積法
\(\omega(x)=e^{-x^2}~\to\) ガウス-エルミート求積法
\(\omega(x)=\frac{1}{\sqrt{1-x^2}}~\to\) ガウス-チェビシェフ(第一種)求積法
\(\omega(x)=\sqrt{1-x^2}~\to\) ガウス-チェビシェフ(第二種)求積法
\(\omega(x)=(1-x)^\alpha (1+x)^\beta~\to\) ガウス-ヤコビ求積法

の形を数値積分したい場合、ガウス求積法を考えてみてはどうでしょうか。
上では省略しましたが、ガウス求積法を用いる際には積分区間も決まっています。しかし、これは簡単な変数変換によって書き換えることが出来ます。

例えば、ガウス-ルジャンドル求積法は区間\([-1,1]\)出なければ使うことはできませんが、変数変換
\(
\displaystyle \int_a^b f(x)dx=\frac{b-a}{2}\int_{-1}^{1}f(\frac{b-a}{2}+\frac{a+b}{2}) dx
\)

の変換が出来るので、左辺を計算しようと思ったら右辺を計算しにいけばいいのです。

ガウス-ルジャンドル求積法 http://mathworld.wolfram.com/Legendre-GaussQuadrature.html
ガウス-radau求積法 http://mathworld.wolfram.com/RadauQuadrature.html
ガウス-ロバート求積法 http://mathworld.wolfram.com/LobattoQuadrature.html
ガウス-ラゲール求積法 http://mathworld.wolfram.com/Laguerre-GaussQuadrature.html
ガウス-エルミート求積法 http://mathworld.wolfram.com/Hermite-GaussQuadrature.html
ガウス-チェビシェフ求積法 http://mathworld.wolfram.com/ChebyshevQuadrature.html
ガウス-ヤコビ求積法 http://mathworld.wolfram.com/Jacobi-GaussQuadrature.html

二重積分を行いたい場合


一番理解しやすい方法は直積を考えて計算する方法です。
「二重積分」を「一重積分した結果を一重積分する」という問題に焼直します。
すなわち、
\(
\displaystyle \int dx \int dy f(x,y)
\)

という問題は、
\(
\begin{align}
\int dx \int dy f(x,y)=\int dy\left[\int f(x;y) dx\right]
\end{align}
\)

と考えるのです。ここで、\(f(x;y)\)のセミコロンの意味は、

『\(f\)は\(x\)と\(y\)の変数だよ。だけど、今変数として注目したいのは\(x\)についてだけ。だから\(y\)は一応書いておくけど、今あらわな変数として扱わないよ。』

ということを示しています。
さて、計算を進めましょう。\([~~]\)内は今、固定したある\(y\)について\(x\)のみの変数なので、求積法の形にして
\(
\begin{align}
\int dy\left[\int f(x;y) dx\right]&\approx \int dy\left[\sum_i w^{(x)}_i f(x_i,y)\right]\\
&= \sum_i w^{(x)}_i \int f(x_i,y)dy \\
&\approx \sum_i w^{(x)}_i \sum_j w^{(y)}_j f(x_i,y_j)
\end{align}
\)

となります。ここで、\(w^{(x)}_i\)は\(x\)方向の\(i\)番目の分点の重みを表します。
これは、例えば被積分関数が\(x\)方向については端点特異性を持っていて、\(y\)方向については簡単な多項式で描けている、といった場合に別々の求積法を使うことで計算量を減らし、高精度の結果が得られるということを示しています。

表面積分を行いたい場合


表面積分
\(
\displaystyle \int_0^\pi d\theta \int_0^{2\pi} d\varphi f(\theta,\varphi) \sin\theta \approx 4\pi \sum_{i} w_i f(\theta_i, \varphi_i)
\)

を行う場合、Lebedev求積法かGaussian-quadrature(ガウス-ルジャンドル求積法と台形則の組み合わせ)が良いです。
lebedev求積法は、球面調和関数の直交性を利用した求積法で、\(f(\theta,\varphi)\)が最大\(l=131\)までの球面調和関数で展開されるのであれば、厳密な値を返すという積分法です。

fortranのプログラムはhttp://www.ccl.net/cca/software/SOURCES/FORTRAN/Lebedev-Laikov-Grids/Lebedev-Laikov.F
にあります。このプログラムと

program main
  implicit none
  integer::i,N
  double precision,allocatable::x(:),y(:),z(:),w(:)
  double precision::s
  double precision,parameter::pi=dacos(-1d0)
  double precision,external::f
 
  N=146
  allocate(x(1:N),y(1:N),z(1:N),w(1:N))
  call LD0146(x,y,z,w,N)
  s=0d0
  do i=1,N
     s=s+w(i)*f(x(i),y(i),z(i))
  enddo
  write(6,*)4d0*pi*s
 
  stop
end program main


function f(x,y,z)
  implicit none
  double precision::x,y,z,f
 
  f=x*x+y*y+z*z
 
  return
end function f

を一緒にコンパイルしてください。プログラムでは\(\theta,\varphi\)での指定ではなく、\(x,y,z\)のデカルト座標系のものになっています。動かすと半径1の球面の立体角に関する積分、すなわち\(4\pi\cdot 1^2\)の値を返します。

3. Lebedev angular quadrature -Methods of Numerical Integration.

Gaussian-quadratureは、\(\theta\)方向についてはガウス-ルジャンドル求積法、\(\varphi\)方向については台形則により積分する、という方法です。良く使われるのでgaussian-quadratureと名前がついているようです。
\(\varphi\)方向で台形則が使われる理由は、\(\varphi\)方向について周期的であるため、台形則は非常に高精度の結果を与えるためです。
\(\theta\)方向でガウス-ルジャンドル求積法が用いられるのは、球面調和関数の\(m=0\)がルジャンドル多項式で展開されるためです。他の\(m\ne 0\)についてはルジャンドル多項式は完全系なので、その性質から計算される太郎という思惑なのだと思います。

lebedev求積法による計算は最大\(l=131\)の求積点と重みが知られていますが、これ以上の\(l\)に関してのlebedev求積点、重みに関する論文は見つけられませんでした。
もしも\(l=131\)以上が出てくることが分かっている場合、現状ではGaussian-quadratureによる求積が良いでしょう。

多次元(6,7次元以上)の積分を行いたい場合


この場合はモンテカルロ積分が良いです。

一様乱数によるD次元の積分は、
\(
\begin{align}
\displaystyle &\int_{x^{(1)}_{a}}^{x^{(1)}_{b}}\int_{x^{(2)}_{a}}^{x^{(2)}_{b}}\cdots\int_{x^{(D)}_{a}}^{x^{(D)}_{b}}
f(x^{(1)},x^{(2)},\cdots,x^{(D)})dx^{(D)}\cdots dx^{(2)} dx^{(1)}\\
&\approx
\frac{(x^{(1)}_{b}-x^{(1)}_{a})(x^{(2)}_{b}-x^{(2)}_{a})\cdots(x^{(N)}_{b}-x^{(N)}_{a})}{N}
\sum_{i=1}^N f(x^{(1)}_i,x^{(2)}_i,\cdots,x^{(D)}_i)
+O(\frac{1}{\sqrt{N}})
\end{align}
\)

と求める事が出来ます。

モンテカルロ積分の誤差はサンプル数にのみに依存します。
\(O(\frac{1}{\sqrt{N}})\)なので、サンプル数\(N\)を100倍に増やすと積分精度が\(\frac{1}{\sqrt{100}}=\frac{1}{10}\)、つまり1桁上がる事を意味しています。

一般的な数値積分法、例えば台形則を用いて1次元積分を\(N\)点で積分するとします。
すると、2次元の積分では\(N^2\)点、D次元では\(N^D\)点被積分関数を評価する必要があります。
台形則の組み合わせによる1次元当たりの積分の精度は\(N\cdot O(h^2)\sim O(\frac{1}{N})\)で決まるので、
1桁、積分精度を上げたいとすれば被積分関数の評価回数は\(10^{D}\)倍に増やさなければなりません。

しかしモンテカルロ積分では積分の精度は次元に依らずサンプル数にのみで決定されます。これは、
1桁、積分精度を上げたいとすれば被積分関数の評価回数は\(100\)倍に増やさなければならないことを意味しています。

台形則と比べるとおおよそ3重積分を解こうとする際に、モンテカルロ法の方が有利になる、ということが分かります。

スポンサーリンク

特定の\(x\)だけ関数値\(f(x)\)が得られる場合


”特定の\(x\)”が等間隔であったり、特に整然としていないと仮定します。
この場合は、

台形則による積分
もしくは
3次スプライン補間による積分

がどんな分点数、分点の間隔であっても、プログラムの変更をほとんどしないで計算できるので良いでしょう。

もしも分点が等間隔の場合で、データ配列の大きさが\(1,2,…,2N+1, (Nは整数)\)の場合はシンプソン則による計算でも良いでしょう。
また、分点が等間隔で、データの配列の大きさが\(1,2,…, 2^N+1(Nは整数)\)の場合はロンバーグ積分法が使えます。

3次スプライン補間による積分は
スプライン補間(1, 2次元) -シキノートに置いてあるc3spline_integralを用いると良いでしょう。

注釈

(※1)”滑らか”とは、
積分区間内の至る所で被積分関数の\(x\)のn階微分の右からの極限と左からの極限が一致すること、を意味します。
(※2)”特異性”とは、
積分区間のどこか一点でも被積分関数の\(x\)のn階微分の右からの極限と左からの極限が一致しない点がある、もしくは積分区間の端点で被積分関数の\(x\)のn階微分が無限大に発散する、ことを意味しています。

モンテカルロ積分

モンテカルロ積分は乱数を利用した積分方法です。

モンテカルロ積分の精度を1桁上げるためには評価点の数を100倍に増やす必要があります。
ですが、この性質は次元数に依らないため、多重積分(5~8次元以上など)を行うときに有利な性質を持ちます。

また、とにかく実装が簡単なので非常に使い易いです。

1次元のモンテカルロ積分


1次元のモンテカルロ積分は定積分を
\(
\displaystyle \int_a^b f(x)dx \approx \frac{1}{N}\sum_{i=1}^N\frac{f(x_i)}{p(x_i)}+O(\frac{1}{\sqrt{N}})
\)

として求めます。
ここで、\(x_i\)は\(p(x)\)に従って\(i\)番目に発生させた乱数で、
\(p(x)\)は区間\([a,b]\)内の乱数の分布を表す関数で、
\(
\displaystyle \int_a^b p(x) dx =1
\)

を満たします。一様乱数であれば
\(
\displaystyle p(x)=\frac{1}{b-a}
\)

です。モンテカルロ積分は一様乱数のみで出来るのではなく、任意の乱数分布に対して適用することが出来ます。

一様乱数によるでD次元の積分は、
\(
\begin{align}
&\int_{x^{(1)}_{a}}^{x^{(1)}_{b}}\int_{x^{(2)}_{a}}^{x^{(2)}_{b}}\cdots\int_{x^{(D)}_{a}}^{x^{(D)}_{b}}
f(x^{(1)},x^{(2)},\cdots,x^{(D)})dx^{(D)}\cdots dx^{(2)} dx^{(1)} \\
&\hspace{1em}\approx
\frac{(x^{(1)}_{b}-x^{(1)}_{a})(x^{(2)}_{b}-x^{(2)}_{a})\cdots(x^{(D)}_{b}-x^{(D)}_{a})}{N}
\sum_{i=1}^N f(x^{(1)}_i,x^{(2)}_i,\cdots,x^{(D)}_i)
+O(\frac{1}{\sqrt{N}})
\end{align}
\)

と求める事が出来ます。

モンテカルロ積分の誤差はサンプル数にのみに依存します。
モンテカルロ積分の誤差は\(O(\frac{1}{\sqrt{N}})\)なので、サンプル数\(N\)を100倍に増やすと積分精度が\(\frac{1}{\sqrt{100}}=\frac{1}{10}\)、つまり1桁上がる事を意味しています。

一般的な数値積分法、例えば台形則を用いて1次元積分を\(N\)点で積分するとします。
すると、2次元の積分では\(N^2\)点、D次元では\(N^D\)点被積分関数を評価する必要があります。
台形則を組み合わせて積分する場合、精度は\(N\cdot O(h^2)\sim O(\frac{1}{N})\)で決まるので、
1桁、積分精度を上げたいとすれば被積分関数の評価回数は\(10^{D}\)倍に増やさなければなりません。

しかしモンテカルロ積分では積分の精度は次元に依らずサンプル数にのみで決定されます。これは、
1桁、積分精度を上げたいとすれば被積分関数の評価回数は\(100\)倍に増やさなければならないことを意味しています。

台形則と比べると収束速度は、おおよそ3重積分を解こうとする際に、モンテカルロ法の方が有利になる、ということが分かります。

(注意)
3次元の定積分を考えます。
この時、台形則では100点で真値と1桁、モンテカルロ積分では10000点で真値と2桁合っているとします。
この状態から精度を1桁上げようとすると関数の評価回数は
台形 \(100\times 10^3= 100,000 \)点
モンテカルロ積分 \(10000\times 100\approx 1,000,000\)点
となり、この場合では台形の方が評価回数が少なくて済みます。ずっと繰り返せばモンテカルロ積分の方が少なくなりますが、現実的な評価回数ではなくなると思います。
なので、評価回数が増えていく割合が変わるのは3重積分あたりですが、最終的に評価回数が少なくなるのはもっと多次元の積分を実行する時かもしれません。

台形則と比較しましたが、台形則の代わりにガウス求積法を用いた場合、もう少し高次元で逆転します。

実用的には6~10次元以上でモンテカルロ積分が使われるようです。

乱数について


モンテカルロ積分は乱数を利用して積分を行うため、乱数の性質が重要になってきます。
私は乱数について詳しくありませんが、調べた限りメルセンヌツイスタという乱数発生法が良い乱数を発生するようです。
メルセンヌツイスタの詳細はMersenne Twister Home Page にあるのでご参照ください。
fortran90コードはhttps://gist.github.com/ykonishi/5569005にあります。

fortran90コード

1次元の積分
\(
\displaystyle \int_0^2 \sin x dx \approx 1.416146836547142\cdots
\)

を実際にメルセンヌツイスタを用いて、計算します。
関数の評価点数は10000点です。

下のコードを実行すると

> gfortran -fno-range-check mt19937.f90 main2.f90
> ./a.out
       10000   1.4167388988307803        1.4161468365471424  
>

という出力を得ます。1列目が関数の評価数、2列目がモンテカルロ積分による積分値、3列目が真値です。
1万点の評価で、ただのsin関数を計算しているのに4桁程度というのは経験的に非常に良くない精度です。
しかも、乱数を用いているので、4桁一致していると判断するのは危険です。
もっと乱数を発生させて、値が変わらないのかチェックを欠かすことはできません。

モンテカルロ積分は以下のプログラムで可能です。

program main
  use mt19937
  implicit none
  integer::i,j,N
  double precision::a,b,s,x,y,z,ans
  double precision,external::f

  call sgrnd(1)
 
  a=0d0
  b=2d0
  N=10000
  ans=cos(a)-cos(b)
     
  s=0d0
  do i=1,N
     x=grnd()*(b-a)+a
     s=s+f(x)
  enddo
  s=s*(b-a)/dble(N)

  write(6,*)N,s,ans
 
  stop
end program main

function f(x)
  implicit none
  double precision::f,x
 
  f=sin(x)
 
  return
end function f

コンパイルは上のファイルをmain.f90と名付けた場合

gfortran mt19937.f90 main.f90

とします。もしもエラーが出てしまう場合

gfortran -fno-range-check mt19937.f90 main.f90

とします。
このエラーはgfortranコンパイラ(ver4.8.4で確認)の問題です。
fortranの整数は
-2147483648から2147483647までで定義されているのにもかかわらず、
gfortran(ver4.8.4)は
-2147483647から2147483647の範囲でなければエラーを出すようです。
しかし、実際に整数に-2147483648を代入し、このエラーを無視する上のオプションを付けて実行すると、ちゃんと値が入ります。なので、コンパイラに問題があると判断しました。

3次元積分


3次元の定積分を2つ考えます。
\(
\displaystyle \int_0^1 dz \int_0^1 dy \int_0^1 dx \sin(x)\sin(y)\sin(z) \approx 0.097144222323873819\cdots
\)

\(
\displaystyle \int_0^1 dz \int_0^1 dy \int_0^1 dx 100 e^{x}\sin(30y)z^5 \approx 0.80735242505763782\cdots
\)

比較する積分法はモンテカルロ積分、台形則、quadpackによる自動積分です。
quadpackは15-31点ガウス求積法を用います。3次元関数のfortranによるquadpackコードはこちら

結果は以下の図の方になりました。

quadpackについては余りにも精度が違い過ぎるので図には重ねていません。
それぞれ、関数の評価回数と積分値を載せました。下線は真値と一致している桁です。

quadpack, key=3
3375points, 0.097 144 222 323 873 861
50625points,0.807 352 425 057 637 71

左の図で、台形則では一次元方向の分点数が10倍(3次元なのでグラフ上では1000倍)になると積分の精度が約1桁上昇し、モンテカルロ積分では点数が100倍増えると積分の精度が約1桁上昇しているのが見て取れます。
右の図ではあたかも台形則の方が傾きが急峻そうに見えますが、まだ安定している積分になってなく、傾きが緩やかになっていきそうなことが分かります。

QUADPACKは流石と言いますが、次元が増えると辛いですがまだまだモンテカルロ積分を使わなくてよさそうです。関数の評価点数が5桁なのにほぼ全桁一致しています。

用いたメインのプログラムは以下のものです。これと上で紹介した
メルセンヌツイスタのfortran90プログラム,3次元quadpackプログラムを一緒にコンパイルしてください。

参考文献


モンテカルロ法 (Monte Carlo method)
確率密度関数からモンテカルロ積分まで
メルセンヌツイスタMersenne Twister Home Page
メルセンヌツイスタのfortran90コードhttps://gist.github.com/ykonishi/5569005