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

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

スプライン補間（1, 2次元）

7件のコメント

3次スプライン補間に関する簡単な説明とfortran90による計算コードです。

概要

・与えられたデータ点すべてを通るように補間
・与えられたデータ点間を3次関数で近似する
・与えられたデータ点間の隣り合う補間関数（スプライン関数）の一階微分、二階微分が連続
・与えられたデータ点の両端は二階微分がゼロと仮定

1次元の場合

区間\([x_i, x_{i+1}]\)の間を3次関数
\(
S_i(x)=a_i(x-x_i)^3+b_i(x-x_i)^2+c_i(x-x_i)+d_i,~(i=0,1,\cdots N-1)
\)
で補間します。補間は与えられたデータ点間の隣り合う区間の補間関数の一階微分と二階微分が一致するように決めます。
すなわち
\(
\displaystyle \frac{d S_i(x_i)}{dx}=\frac{d S_{i-1}(x_i)}{dx} \\
\displaystyle \frac{d^2 S_i(x_i)}{dx^2}=\frac{d^2 S_{i-1}(x_i)}{dx^2} \\
\)
という式が成り立っているとします。与えられたデータ点の端、点\(x_0\)と\(x_N\)では自然境界条件と呼ばれる境界条件
\(
\displaystyle \frac{d^2 S_0(x_0)}{dx^2}=\frac{d^2 S_{N-1}(x_N)}{dx^2}=0
\)
を課します…この境界条件を課す理由は良く分かりません。[1]に端ではない領域で補間が滑らかになるようにこう決める、とあります。良く分かりません。きっと解く都合上こう決めると簡単に解けるからでしょう。実際そうですし…

詳細は[1]が詳しいのでそちらを参照してください。
係数は以下の通り与えられます。
\(
\begin{align}
a_i&=\frac{v_{i+1}-v_i}{6(x_{i+1}-x_i)} \\
b_i&=\frac{v_i}{2}\\
c_i&=\frac{y_{i+1}-y_i}{x_{i+1}-x_i}-\frac{1}{6}(x_{i+1}-x_i)(2v_i+v_{i+1}) \\
d_i&=y_i
\end{align}
\)
ここで、\(v_i\)は以下の方程式の解です。
\(
\scriptsize
\begin{eqnarray}
\left(
\begin{array}{ccccc}
2(h_0+h_1)& h_1 & 0 & \cdots & 0 \\
h_1 & 2(h_1+h_2) & h_2 & \cdots & 0\\
0 & h_2 & 2(h_2+h_3) & \cdots & 0\\
\ldots &&&& \\
0& \cdots & 0 & h_{N-2} & 2(h_{N-2}+h_{N-1})
\end{array}
\right)
\left(
\begin{array}{c}
v_1 \\
v_2 \\
v_3 \\
\ldots \\
v_{N-1}
\end{array}
\right)=\left(
\begin{array}{c}
w_1 \\
w_2 \\
w_3 \\
\ldots \\
w_{N-1}
\end{array}
\right)
\end{eqnarray}
\normalsize
\)
ここで、
\(
\begin{align}
h_i&=x_{i+1}-x_i,~ (i=0,1,\cdots,N-1)\\
w_i&=6\left(\frac{y_{i+1}-y_{i}}{h_i}-\frac{y_i-y_{i-1}}{h_{i-1}}\right),~ (i=1,2,\cdots,N-1)
\end{align}
\)
です。

計算コード

Fortranコードを以下に置いておきます。
注意しておきますが、[1]を元にしたコードではありません。[2]を元にしたコードです。

・配列 fdata(0:N),xdata(0:N) は与えられたデータで、既知の値です。
・ c3spline1p は任意の点 x の補間値 f 、一階微分 df 、二階微分 df2 を fdata(0:N),xdata(0:N) を元に計算します。
・ c3spline は任意の点配列 xa(0:M) の補間値 fa(0:M) を fdata(0:N),xdata(0:N) を元に計算します。
・ c3spline_integral はスプライン補間された関数の積分 s を計算します。

c3spline1p を用いて複数の点を計算することもできますが、c3spline に比べて遅いです。
計算時間は  c3spline1p : c3spline = 1 : 0.03です。

2次元の補間

2次元の補間を考えます…が、少し横着します。
「横着」は2次元のスプラインを調べるのが面倒で、思いついてしまったのでこれでいいだろう、という意味です。

格子状に与えられたデータ列\(f(x_i,y_j), i=0,\cdots,N_x, j=0,\cdots,N_y\)があり、点\((x,y)\)での補間された値が欲しいとします。
点\((x,y)\)は、それぞれ
\(
\begin{align}
& x_i\le x \lt x_{i+1} \\
& y_j\le y \lt y_{j+1}
\end{align}
\)
に存在するとします。
この時、考え方は以下の通りです。

まず、\(j\)を固定して\(x\)上の補間値を\(N_y\)個(点\(f(x,y_0),f(x,y_1),f(x,y_2),\cdots,f(x,y_{N_y})\))を得ます。
その後、補間された\(N_y\)個の補間値に対して補間を行い、点\((x,y)\)に補間値\(f(x,y)\)を得ることが出来ます。

ここで問題なのが、\(x\)方向、\(y\)方向のどちらを先に補間するのか?
そしてその2つの結果は変わるのか?ということです。

直感的に考えて良さそうなのは２つの平均を取ってしまう事でしょう。しかし、単純に考えて計算時間が2倍になります。
直ぐには分からなそうなので、プログラムではそれを指定できる形にする、にとどめましょう。

偏微分を求められます。最後に補間した方向の偏微分を得ることが出来ます。

c3spline2d1pは元となるデータ配列を入れると、点(x,y)の補間された値を返します。
c3spline2dは元となるデータ配列を入れると、配列で指定された格子状の点を返します。

上記コードを動かすと、以下の画像を作ることが出来ます。

赤い点は補間に用いたオリジナルのデータ列を表し、黒の線は補間された関数を表します。また、底面のカラーマップは本当の関数との相対誤差を表しています。xが大きいところの端で補間した結果の関数が本当の関数とあまり一致していません。

これは、元にするデータ点の両端において、元の関数の二階微分がゼロという仮定が満たされていないためです。

上の図を再現するgnuplotのコードは以下のものです。

set palette defined ( -1 '#ffffff', 0 '#000090',1 '#000fff',2 '#0090ff',3 '#0fffee',4 '#90ff70',5 '#ffee00',6 '#ff7000',7 '#ee0000',8 '#7f0000')

set ticslevel 0
set view 46,40,1,1

set xl "x"
set yl "y"
set cbr[0:1]
unset key
splot "interpolate_data1.d" u 1:2:6 w pm3d at b
replot "interpolate_data1.d" u 1:2:3 w l lw 1
replot "original_data.d" u 1:2:3 w p pt 7 ps 0.5

参考文献

[1]https://www.mk-mode.com/rails/docs/INTERPOLATION_SPLINE.pdf

[2]http://150.19.250.16/MULTIMEDIA/numeanal/node16.html

アイキャッチ画像のフォント：キルゴhttp://getsuren.com/killgoU.html

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

sin基底とquadpackによる1次元時間依存しないシュレーディンガー方程式

コメントする

確実に積分の計算が可能な、適応型の数値積分パッケージQUADPACKを用いて行列要素を計算します。
精度に影響するのは基底関数（sin関数）の数のみです。

固定端条件の下で時間依存しない1次元シュレーディンガー方程式
\(\displaystyle
\left[-\frac{1}{2}\frac{d^2}{dx^2}+V(x)\right]\psi(x)=E\psi(x)
\)
をsin基底で展開し、変分原理に基づいて対角化します。

sin基底：
\(\displaystyle
\varphi_n(x)=\sqrt{\frac{2}{b-a}}\sin\left(n\pi\frac{x-a}{b-a}\right)
\)

ここで紹介するコードは並列計算に対応しています。

例えば、有限深さの井戸型ポテンシャルのような不連続点があったとしても、
行列要素はquadpackを用いて確実に求める事が出来るため、後はどれだけ基底を積むか？にかかってきます。

必要なファイルは以下のものです。

tar

https://slpr.sakura.ne.jp/qp/supplement_data/tise1d_by_sinbasis_and_quadpack/tise1d_sinbasis_quadpack.tar.gz
もしくは

個別のf90ファイル

https://slpr.sakura.ne.jp/qp/supplement_data/tise1d_by_sinbasis_and_quadpack/main.f90
https://slpr.sakura.ne.jp/qp/supplement_data/tise1d_by_sinbasis_and_quadpack/quadpack.f90

対角化のパッケージを用いるため、mklやlapackと一緒にコンパイルしてください。例えば

ifort -mkl quadpack.f90 main.f90

等です。

デフォルトでは、時間依存しないシュレーディンガー方程式
\(\displaystyle
\left[-\frac{1}{2}\frac{d^2}{dx^2}+\frac{1}{2}x^2\right]\psi(x)=E\psi(x)
\)
を解きます。計算に用いているパラメータは、
区間\([-20:20]\),
基底の数\(80\),
で計算します。

プログラミングと数値計算, 物理学

水素様原子の波動関数のコード

コメントする

非相対論的な水素様原子の電子の波動関数を求めるfortranコードを書きます。

解く問題

電荷\(+Ze\)の原子核の周りに存在する、電荷\(-e\)を持つ電子が満たす時間依存しないシュレーディンガー方程式は、
\(
\displaystyle \left[-\frac{1}{2}\Delta-\frac{Z}{r}\right]\psi({\bf r})=E\psi({\bf r})
\)
です。
\(
\left. \psi({\bf r})\right|_{|{\bf r}|\to \infty}\to 0
\)
の境界条件を満たす正規直交化された解は、
\(
\begin{align}
\psi({\bf r})&=R_{n,l}(r)\cdot Y_{l,m}(\theta,\varphi)\\
&=\sqrt{\left(\frac{2}{n}\right)^3\frac{(n-l-1)!}{2n[(n+l)!]}} e^{-Zr/n} \cdot \left(2\frac{Zr}{n}\right)^l\cdot L_{n-l-1}^{(2l+1)}\left(2\frac{Zr}{n}\right)\cdot Y_{l,m}(\theta,\varphi)
\end{align}
\)
と表されます。ここで、

  • \(n,l,m\)  主量子数,方位量子数,磁気量子数
  • \(r,\theta,\varphi\)  動径方向,z軸方向の角度,xy平面の角度
  • \(L_n^{(\alpha)}(x)\)  ラゲール陪多項式[1]
  • \(Y_{l,m}(\theta,\varphi)\)  球面調和関数[2]

を意味します。

動径関数\(R_{n,l}(r)\)は微分方程式
\(
\displaystyle \left[-\frac{1}{2}\left(\frac{d^2}{dr^2}+\frac{2}{r}\frac{d}{dr}\right)+\frac{l(l+1)}{2r^2}-\frac{Z}{r}\right]R_{n,l}(r)=E_n R_{n,l}(r)
\)
を満たす正規直交化された関数で、固有値は
\(
\displaystyle E_n=-\frac{1}{2n^2}
\)
です。

波動関数を出力するFortranコードはこれ。

計算結果

3次元の波動関数を表示するには3つの座標と波動関数の値の計4つが必要です。
しかし、4次元空間を分かりやすく描画する方法はありません。

そこで、ここでは\(\mathrm{Re}{\psi(x,0,z,t)}\)のように、\(y=0\)にしてグラフを描画します。
\(z\)軸周りの構造がうまく表現できませんが、良しとします。

\(n=1\)

\(n=2\)

\(n=3\)

”color”と書いてあるのは、波動関数の形がはっきり映るようにカラーバーの範囲を調整したものです。
なので、色が飽和しているところはそのカラーバーの範囲で表現することはできない領域です。

\(n=20\)まで同じような画像を作ったのですが、画像の容量が大きくなりすぎました。
なので、見たい方は以下のリンクから元ファイルをダウンロードしてください。

\(n=1\)  Hatomc_n1.png（48KB）

\(n=2\)  Hatomc_n2.png（223KB）

\(n=3\)  Hatomc_n3.png（612KB）

\(n=4\)  Hatomc_n4.png（1MB）

\(n=5\)  Hatomc_n5.png（2MB）

\(n=6\)  Hatomc_n6.png（3MB）

\(n=7\)  Hatomc_n7.png（4MB）

\(n=8\)  Hatomc_n8.png（3MB）

\(n=9\)  Hatomc_n9.png（6MB）

\(n=10\)  Hatomc_n10.png（6MB）

\(n=11\)  Hatomc_n11.png（9MB）

\(n=12\)  Hatomc_n12.png（9MB）

\(n=13\)  Hatomc_n13.png（10MB）

\(n=14\)  Hatomc_n14.png（14MB）

\(n=15\)  Hatomc_n15.png（15MB）

\(n=16\)  Hatomc_n16.png（16MB）

\(n=17\)  Hatomc_n17.png（18MB）

\(n=18\)  Hatomc_n18.png（19MB）

\(n=19\)  Hatomc_n19.png（21MB）

\(n=20\)  Hatomc_n20.png（21MB）

n=1～n=20までのファイル hatom.zip（164MB）

例えば、n=13, n=20はこんな感じです（荒くしてます）。

特殊関数に関するプログラミングは特殊関数へ。

How to plot in gnuplot?
Just for note for myself.

You need 2 files which named : multi.plt and atommulti.plt

Before ploting, you must make directory “dat” and move every data files to “dat”.

and in gnuplot, type

call "multi.plt" 2

then, you can get figure.

multi.plt

atommulti.plt

参考

[1]ラゲール陪多項式\(L_m^{(\alpha)}(x)\)

  • 定義域

    \(\;\;\;\;
    0 \lt x \lt \infty
    \)

  • 微分方程式

    \(\;\;\;\;
    \displaystyle
    \left[x\frac{d^2}{dx^2}+(\alpha+1-x)\frac{d}{dx}+n\right]L_m^{(\alpha)}(x)=0
    \)

  • 具体例

    \(\;\;\;\;
    \begin{align}
    & L^{(0)}_0(x)=1 \\
    & L^{(0)}_1(x)=-x+1 \\
    & L^{(0)}_2(x)=\frac{1}{2}(x^2-4x+2) \\
    & L^{(0)}_3(x)=\frac{1}{6}(-x^3+9x^2-18x+6) \\
    & L^{(1)}_0(x)=1 \\
    & L^{(1)}_1(x)=-x+2 \\
    & L^{(1)}_2(x)=\frac{1}{2}\left(x^2-6x+6\right) \\
    & L^{(1)}_3(x)=\frac{1}{6}\left(-x^3+12x^2-36x+24\right)
    \end{align}
    \)

  • 漸化式

    \(\;\;\;\;
    \begin{align}
    L_n^{(\alpha)}(x)&=\binom{n+\alpha}{n}a_0(x) \\
    &\hspace{1em}a_{m-1}(x)=1-\frac{n-m+1}{m(\alpha+m)}\cdot x\cdot a_m(x) \\
    &\hspace{5em}(m=n,n-1,n-2,\cdots,1,\;\;a_n(x)=1)
    \end{align}
    \)
    [3]
    1階微分
    \(\;\;\;\;
    \begin{align}
    \frac{d}{dx}L_n^{(\alpha)}(x)=-L_{n-1}^{(\alpha+1)}(x)
    \end{align}
    \)
    2階微分
    \(\;\;\;\;
    \begin{align}
    \frac{d^2}{dx^2}L_n^{(\alpha)}(x)=L_{n-2}^{(\alpha+2)}(x)
    \end{align}
    \)

  • 直交性

    \(\;\;\;\;
    \displaystyle
    \int_{0}^{\infty} e^{-x}x^{\alpha}L^{(\alpha)}_n(x)L^{(\alpha)}_{n’}(x) dx =\frac{\Gamma(\alpha+n+1)}{n!}\delta_{n,n’},\;\;(\alpha>-1)
    \)

[2]球面調和関数\(Y_{l,m}(\theta,\varphi)\)
ルジャンドル多項式(\(P_l^m(x)\),[3])を用いて
\(
Y_{l,m}(\theta,\varphi)=(-1)^{(m+|m|)/2}\sqrt{\frac{2l+1}{4\pi}\frac{(l-|m|)!}{(l+|m|)!}}P_l^{|m|}(\cos\theta)e^{im\phi}
\)

[3]ルジャンドル陪多項式\(P_l^m(x)\)

  • 定義域

    \(\;\;\;\;
    1 \lt x \lt -1
    \)

  • 微分方程式

    \(\;\;\;\;
    \displaystyle (1-x^2)\frac{d^2}{dx^2}P_l^m(x)-2x\frac{d}{dx}P_l^m(x)+\left[l(l+1)-\frac{m^2}{1-x^2}\right]P_l^m(x)=0
    \)

  • 具体例

    \(\;\;\;\;
    \begin{align}
    & P_0^0(x)=1 \\
    & P_{1}^{-1}(x)=\frac{1}{2}(1-x^2)^{1/2} \\
    & P_{1}^{0}(x)=x \\
    & P_{1}^{1}(x)=-(1-x^2)^{1/2} \\
    & P_{2}^{-2}(x)=\frac{1}{8}(1-x^2) \\
    & P_{2}^{-1}(x)=\frac{1}{2}x(1-x^2)^{1/2} \\
    & P_{2}^{0}(x)=\frac{1}{2}(3x^2-1) \\
    & P_{2}^{1}(x)=-3(1-x^2)^{1/2} \\
    & P_{2}^{2}(x)=3(1-x^2)
    \end{align}
    \)

  • 漸化式

    用いている式
    \(\;\;\;\;
    \begin{align}
    & P_{|m|}^{|m|}(x)=(-1)^{|m|}(2|m|-1)!!(1-x^2)^{\frac{|m|}{2}} \\
    & P_{l}^{-m}(x)=(-1)^m\frac{(l-m)!}{(l+m)!}P_l^m(x) \\
    & P_{|m|}^{|m|+1}(x)=(2|m|+1)xP_{|m|}^{|m|}(x) \\
    & P_{|m|}^{|m|+q}(x)=\left(\frac{2|m|-1}{q}+2\right)xP_{|m|}^{|m|+q-1}(x)-\left(\frac{2|m|-1}{q}+1\right)P_{|m|}^{|m|+q-2}(x)
    \end{align}
    \)[1]
    1階微分
    \(\;\;\;\;
    \begin{align}
    m=\pm 1&\; \\
    &:\hspace{1em} \frac{d}{dx}P_l^m(x)=\frac{lx}{(x^2-1)}P_l^m(x)-\frac{l+m}{(x^2-1)} P_{l-1}^m(x) \\
    m\ne\pm 1\; \\
    &:\hspace{1em} \frac{d}{dx}P_l^m(x)=c_{2}P^{m+2}_{l-1}(x)+c_0 P^{m}_{l-1}(x)+c_{-2} P^{m-2}_{l-1}(x) \\
    &\hspace{1em} c_{2}=\frac{1}{4(m+1)}\\
    &\hspace{1em} c_{0}=\frac{l+m}{2}\left(1+\frac{l}{1-m^2}\right)\\
    &\hspace{1em} c_{-2}=-\frac{(l+m)(l+m-1)(l+m-2)(l-m+1)}{4(m-1)}
    \end{align}
    \)
         ※\(m=\pm 1\)の時、\(x=\pm 1\)で発散します。
    2階微分
    \(\;\;\;\;
    \begin{align}
    m=\pm 1,\pm 3&\; \\
    &:\hspace{1em} \frac{d^2}{dx^2}P_l^m(x)=\frac{1}{1-x^2}\left[2x\frac{dP_l^m(x)}{dx}+\frac{(l+1)(l+m)}{2l+1}\frac{dP_{l-1}^m(x)}{dx}-\frac{l(l-m+1)}{2l+1}\frac{dP_{l+1}^m(x)}{dx}\right] \\
    m\ne\pm 1,\pm 3&\; \\
    &:\hspace{1em} \frac{d^2}{dx^2}P_l^m(x)=c_{2}\frac{dP^{m+2}_{l-1}(x)}{dx}+c_0 \frac{dP^{m}_{l-1}(x)}{dx}+c_{-2} \frac{dP^{m-2}_{l-1}(x)}{dx}
    \end{align}
    \)
         ※\(m=\pm 1,\pm 3\)の時、\(x=\pm 1\)で発散します。

  • 直交性

    \(l\)に対する直交性
    \(\;\;\;\;
    \displaystyle
    \int_{-1}^{1} P_{m}^{l}(x)P_{m}^{l’}(x) dx =\frac{2}{2l+1}\cdot \frac{(l+m)!}{(l-m)!}\delta_{l,l’}
    \)

[3] Abramowitz & Stegun (1965), ISBN 978-0486612720 (Numerical Method, 22.18, Page 789)
[4] Abramowitz & Stegun (1965), ISBN 978-0486612720 (Numerical Method, 22.2, Page 774,775)

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

最速・高精度の数値積分

コメントする

無償で公開されている1次元数値積分コードで優秀なルーチンは何か？を見つけます。

何よりも早さが足りない人へ。
言語はFortranです。

まとめ

QUADPACK(ルーチン名dqag, key=2,6)
http://www.netlib.org/quadpack/
が広い場合で使えて優秀です。しかもパブリックドメイン！
直ぐ使えるように私が変更を加えたQUADPACKプログラムが下の方にあります。

もう少し問題に対して最適化をしたければ以下のチャートに従ってください。

もくじ

  1. 比較対象について(QUADPACK, NUMPAC, Ooura’s MSP)
  2. ライセンスについて
  3. 比較方法
  4. 結果
  5. QUADPACKプログラム(1/2次元の実/複素積分,3次元実積分)
  6. ベンチマーク時のプログラム
  7. 参考文献

比較対象(QUADPACK, NUMPAC, Ooura’s MSP)

どんな関数が入ってきても計算時間と精度を両立させる万能な数値積分手法は存在しません。一長一短があります。

この一長一短は非常に極端です。
被積分関数に特定の性質があるのであればその方法を使うべきです。
特定の積分範囲、被積分関数が(重み関数)×(多項式)で表現されるのであれば、ガウス求積法が最速で高精度です。
被積分関数が端点で特異性(例えば端点で\(x^{-1/2}\))を持つのであれば二重指数関数型数値積分が最速で高精度です。

ここでは、無償公開されている自動積分コードのみを比較します。
使っている手法が実行したい積分とマッチすれば自動積分は最善の方法です。

比較するコードは以下のものです。

比較コード
パッケージ
or 作成者（ルーチン名）		 積分戦略 積分の手法
NUMPAC(aqnn9d)
http://www2.itc.nagoya-u.ac.jp/center/ja/numpac/ 		適応型 ニュートン・コーツ9点則
NUMPAC(rombgd)
http://www2.itc.nagoya-u.ac.jp/center/ja/numpac/		 非適応型 ロンバーグ積分
Ooura(intcc)
http://www.kurims.kyoto-u.ac.jp/~ooura/intcc-j.html		 非適応型 クレンショウ-カーティス積分
Ooura(intde)
http://www.kurims.kyoto-u.ac.jp/~ooura/index-j.html		 非適応型 二重指数関数型数値積分
QUADPACK(dqag,key=2)
http://people.sc.fsu.edu/~jburkardt/f_src/quadpack/quadpack.html		 適応型 7-15点ガウス-ルジャンドル、
ガウス-クロンロッド求積法
QUADPACK(dqag,key=6)
http://people.sc.fsu.edu/~jburkardt/f_src/quadpack/quadpack.html		 適応型 30-61点ガウス-ルジャンドル、
ガウス-クロンロッド求積法
Sikino(romberg)
https://slpr.sakura.ne.jp/qp/nc-integral/		 非適応型 ロンバーグ積分
Sikino(de)
https://slpr.sakura.ne.jp/qp/de-fourmula/		 非適応型 二重指数関数型数値積分
Sikino(rkf451)
https://slpr.sakura.ne.jp/qp/runge-kutta-ex/		 適応型 ルンゲ-クッタ-フェールベルグ法による積分

　適応型というのは、被積分関数の形に依って分点を増やす減らすを決める自動積分法です。
非適応型というのは、被積分関数の形に依らず分点を増やす減らすを決める自動積分法です。

[adsense1]

ライセンスについて

各プログラムの利用条件に関して注記しておきます。
私の理解であるので、解釈に間違いがあるかもしれません。最終的な判断は個人の判断にお任せします。

NetlibにあるQUADPACK
http://www.netlib.org/quadpack/
パブリックドメイン。Fortran数値積分用のパッケージ。
著作権表示無しに、商用非商用問わず複製、利用、改変、配布して良いもの。
ただし、所有権や人格権などを侵してはならない。

私がNetlibにあるQUADPACKをfortran90用に改変したものを置いておきます。下のQUADPACKプログラムへ。

people.sc.fsu.eduにあるQUADPACK
GNU LGPLライセンス。Fortran数値積分用のパッケージ。

! Licensing:
!
! This code is distributed under the GNU LGPL license.

とあり、http://people.sc.fsu.edu/~jburkardt/f_src/quadpack/quadpack.html
にコードがあります。こちらの方が使うまでが簡単ですが、ライセンスがあることが面倒。

NUMPAC

NUMPACを利用した論文には，NUMPAC利用の記述を記載してください。
また，著者から修正プログラムの提供があったときには置き換えてください。
NUMPAC邦文マニュアル – XML

とあります。著作権表示や複製、利用、改変、配布に関する記述は無いので、これらをするべきではありません。

Ooura’s Mathematical Software Packages
プログラムのreadmeの中に

copyright
Copyright(C) 1996 Takuya OOURA (email: ooura@mmm.t.u-tokyo.ac.jp).
You may use, copy, modify this code for any purpose and
without fee. You may distribute this ORIGINAL package.

とあります。翻訳すれば、
”著作権はTakuya OOURAにあります。どんな目的に対しても無償でコードを使用、複製、改変して良いです。そしてオリジナルパッケージを配布して良いです。”
とあります。
例えば上の説明にあてはまらないため、許諾なしに改変したものを配布してはいけません

本サイト(シキノート)のプログラム
私のは最終的な出版物や何かに、私の著作権表示、利用した旨を記述してもらえれば商用非商用であっても複製、利用、改変、配布して良いです。
論文に使用した場合、表示は要りません。
改変したものを配布したい場合、著作権、責任は改変を行ったものにあります。また、オリジナルのコードはシキノートにあることを注記してください。オリジナルのコード、また改変したものについてもそれによって生じた責任は負いません。

すなわち本当に自由に使えるのはNetlibで公開されているQUADPACKだけです。

比較方法

Kahaner’s 21 test problemsより、問題を抜粋して比較します。
Kahaner’s 21 test problemsとは21個のベンチマーク用の問題で、全て解析解が存在するものです。
これが解けたら大した数値積分プログラムだ！ということです。
被積分関数に特異性や発散などを含んでいる問題です。

より、以下の10問を選んで比較します。問題の抜粋は私が独断で選びました。

結果

計算結果を並べます。比較は、

  • 積分プログラムに積分区間と精度を要求し、その精度に達したか？
  • 要求精度を実現するまでに被積分関数が評価された回数

を比較します。また、適応型の積分手法に対して、実際にどの点が評価されているかを知りたかったため、
実際に評価された点達に対し矩形波で畳み込みを行い、最大値を１にするように規格化し比較しました。
これが左から１枚目の画像です。
2枚目は、真値との相対誤差を要求精度\(\varepsilon\)の関数として、
3枚目は、要求精度が達成された場合の被積分関数の評価回数です。
3枚目が一番重要な情報です。

この問題(2),(3)は被積分関数に不連続性、特異性がある場合の問です。
(2)の有効な手法は適応型の数値積分法であることが分かります。シンプルな問題ですが、大域的に精度を決める自動積分では3-4桁程度にとどまり、全く計算できていません。
(3)は左の端点で一回微分が発散する関数です。最も有効な手法は変数変換型の数値積分、二重指数関数型数値積分法で、次いで適応型の積分法です。一回微分が発散が意味するのは、テーラー展開による誤差評価で、誤差項が発散することを意味しています。


(3)はラグランジュ補間が難しい関数です。ラグランジュ補間が難しい場合、チェビシェフ多項式による補間が有効です。
(9)は振動型の関数です。sin(x)のテーラー展開が難しいように、これもチェビシェフ多項式による補間が有効です。


(13)は早く振動する関数です。関数に節の数が多いため次数の高い手法が必要、もしくは分割するしかありません。
(14)は振動型の関数です。sin(x)のテーラー展開が難しいように、これもチェビシェフ多項式による補間が有効です。


(17)はx=0回りだけに値が存在します。適応型でなければ、局在部分を評価するために要らない部分を増やさなければならない無駄が生じてしまいます。
(18)は程々の振動型の関数です。チェビシェフ多項式による補間が有効です。


(20)はローレンツ型の関数です。ローレンツ型はラグランジュ補間ではうまく補間されません。
(21)は突如値が大きくなる地点が存在する関数です。うまく変化する点を見つけないと的外れな計算結果になります。

ちなみに

二重指数関数型数値積分が有効な例として
\(
\displaystyle \int_0^{1}\sin(1/\sqrt{x})/\sqrt{x}
\)
という例題が挙げられます。難しい事には変わりませんが二重指数関数型数値積分だけが出来るのではないことをここに注記しておきます。

二重指数関数型数値積分の凄いところは、適当に作った私のプログラムでさえ有名な積分プログラムと張り合えるという点でしょう。

非適応型の数値積分法の重みは被積分関数の形に依存しないで決まることを確認しておきましょう。いつでもこんな感じの重みになります。

まとめです。

問題番号 被積分関数の特徴 最も有効な積分法 次いで有効な積分法
(2) 不連続性 NUMPAC(aqnn9d) QUADPACK(dqag, key=2)
(3) 端点特異性 Ooura(intde) NUMPAC(aqnn9d)
(5) ローレンツ Ooura(intcc) QUADPACK(dqag, key=6)
(9) 振動関数 Ooura(intcc) QUADPACK(dqag, key=6)
(13) 高振動 Ooura(intcc) QUADPACK(dqag, key=6)
(14) 局在 NUMPAC(aqnn9d) Ooura(intcc)
(17) 局在、振動 Ooura(intcc) QUADPACK(dqag, key=6)
(18) 振動 Ooura(intcc) QUADPACK(dqag, key=6)
(20) ローレンツ Ooura(intcc) QUADPACK(dqag, key=6)
(21) 突如変化する関数 NUMPAC(aqnn9d) NUMPAC(rombgd)

ある精度で積分の問題を解きたい時に、最高のパフォーマンスを出したければ以下のチャートに従うと良いと思います。ライセンスの問題を含めて、QUADPACKが優秀でしょう。

QUADPACKのkey=2がよいか、key=6がよいかは良く分かりません。上の表にはkey=6が優秀に見えますが、\(\sqrt{x}\)の積分などでは圧倒的にkey=2が優秀だからです。

QUADPACKプログラム

私が手を加えているので、オリジナルと比較し相違点が生じるかもしれません。
また、使用するにあたっては、私の作成したものと同様、元にしたものが私の物である、と表示さえしていただければ、商用非商用問わず複製、利用、改変、配布してくださって結構です。ただし、責任は負いません。
問題ないと思いますが、もし公開しているプログラムに著作権的にダメな点があればお知らせください。判明次第、ページを非公開にし、早急に対応いたします。

1次元実積分

Quadpackのコード
https://slpr.sakura.ne.jp/qp/supplement_data/quadpack_dqag.f90
と
プログラムのサンプルコード

program main
  implicit none
  integer::ier
  double precision::a,b,s,eps,g
  external::g
 
  a=0d0
  b=1d0
  eps=1d-12

  call dqag_k2(g,a,b,eps,s,ier)
  write(6,*)s,ier
  ! ier = 0 : success
  ! ier > 0 : fail

  call dqag_k6(g,a,b,eps,s,ier)
  write(6,*)s,ier
 
  stop
end program main

function g(x)
  implicit none
  double precision::g,x
 
  g=sqrt(x)

  return
end function g

g　　：　被積分関数\(g(x)\)
a,b　：　積分区間\([a,b]\)
eps　：　要求する相対誤差
s　　：　\(g(x)\)の数値積分値
ier　：　ier=0 成功、収束した。
　　　　　ier>0 失敗、収束しなかった。

本来のquadpackには色々な引数がありますが、それを省略するためにdqag_k2、dqag_k6というサブルーチンを挟んでいます。これはquadpack_dqag.f90の中に入っているので参照してください。
dqag_k2は7-15点ガウス-ルジャンドル、ガウス-クロンロッド求積法を用いた適応型の数値積分法です。
dqag_k6は30-61点ガウス-ルジャンドル、ガウス-クロンロッド求積法を用いた適応型の数値積分法です。

コンパイルと実行は

$ gfortran quadpack_dqag.f90 main.f90
$ ./a.out
  0.66666666666666718                0
  0.66666666666666796                0
$

でokです。積分値が0の時、相対誤差の条件は満たされないためier=2が返ってきます。
問題に応じて無視するかどうか決めてください。

もしも値が存在するはずなのにier=1が帰ってくるようであれば、サブルーチンdqag_k2、dqag_k6の中のlimitの値を増やしてください。こんなふうに。

  ier=1
  epsabs=-1d0  
  !limit=200
  limit=1000
  lenw=limit*4

1次元複素積分

2017/11/02 追記)
複素平面上の点\(a=(a_x,a_y)\)から、\(b=(b_x,b_y)\)に至る経路を直線で結ぶ、複素関数\(g(z)\)の線積分は以下のコードです。

ただし、多価関数は不可能です。偏角の範囲はfortran90では\([-\pi\lt \theta \le \pi]\)出なければなりません。

Quadpackのコード
https://slpr.sakura.ne.jp/qp/supplement_data/quadpack_cqag.f90
と
プログラムのサンプルコード

program main
  implicit none
  integer::ier
  double precision::eps
  complex(kind(0d0))::a,b,s
  complex(kind(0d0)),external::g
 
  a=dcmplx(0d0,0d0)
  b=dcmplx(1d0,1d0)
  eps=1d-12

  call cqag_sk(g,a,b,eps,s,2,ier)
  write(6,*)s,ier
 
  stop
end program main

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

  g=1d0/sqrt(z)

  return
end function g

2次元実積分

2018/01/03追）
wolfram alphaによる解析解
http://www.wolframalpha.com/input/?i=%5Cint_%7B0%7D%5E2+%5Cint_%7B0%7D%5E1+sin(x)cos(10y%5E2)+dxdy+10digits
\(
\displaystyle \int_{0}^2 dy \int_{0}^1 dx sin(x)cos(10y^2) \approx 0.09975138562
\)

Quadpackのコード
https://slpr.sakura.ne.jp/qp/supplement_data/quadpack_dqag2d.f90
と
プログラムのサンプルコード

program main
  implicit none
  integer::ier,key
  double precision::eps
  double precision::xa,xb,ya,yb,s
  double precision,external::g,h
 
  xa=0d0
  xb=1d0
  ya=0d0
  yb=2d0
  eps=1d-8
  key=2

  call dqag_sk2d(h,xa,xb,ya,yb,eps,s,ier,key)
  write(6,*)s,ier
 
  stop
end program main

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

  h=sin(x)*cos(y*y*10d0)
 
  return
end function h

2次元複素積分

2018/01/03追）
wolfram alphaによる解析解
http://www.wolframalpha.com/input/?i=int+from+1%2Bsqrt(-1)+to+3+int+from+-2%2Bsqrt(-1)+to+5+of+exp(-sqrt(x)*y)+dx+dy
\(
\displaystyle \int_{1+i}^3 dy \int_{-2+i}^5 dx e^{-\sqrt{x}y} \approx -3.41684-2.47923i
\)

Quadpackのコード
https://slpr.sakura.ne.jp/qp/supplement_data/quadpack_cqag2d.f90
と
プログラムのサンプルコード

program main
  implicit none
  integer::ier,key
  double precision::eps
  complex(kind(0d0))::xa,xb,ya,yb,s
  complex(kind(0d0)),external::g
 
  xa=dcmplx(-2d0,1d0)
  xb=dcmplx(5d0,0d0)
  ya=dcmplx(1d0,1d0)
  yb=dcmplx(3d0,0d0)
  eps=1d-6

  key=6
  call cqag_sk2d(g,xa,xb,ya,yb,eps,s,key,ier)
  write(6,*)s,ier
 
  stop
end program main

function g(x,y)
  implicit none
  complex(kind(0d0))::g
  complex(kind(0d0)),intent(in)::x,y
 
  g=exp(-sqrt(x)*y)

  return
end function g

3次元実積分

2018/01/30追）
https://slpr.sakura.ne.jp/qp/supplement_data/quadpack_dqag3d.f90

integrate{sin(x)*cos(y*z)*z,{x,1,3},{y,2,4},{z,0,1}} -wolfram alpha
\(
\displaystyle \int_1^3 \int_2^4 \int_0^1 \sin(x) \cos(y z) z dz dy dx \approx -0.450 920 512 214 481 28…
\)

program main
  implicit none
  integer::ier,key
  double precision::xa,xb,ya,yb,za,zb,s,eps
  double precision,external::f

  xa=1d0
  xb=3d0
  ya=2d0
  yb=4d0
  za=0d0
  zb=1d0
  key=3  
  eps=1d-8
  call dqag_sk3d(f,xa,xb,ya,yb,za,zb,eps,s,ier,key)

  write(6,*)s
 
  stop
end program main

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

[adsense2]

ベンチマーク用プログラム

自分の持っているプログラムのベンチマークを行いたい場合のプログラムを公開します。

参考文献

QUADPACK Netlibhttp://www.netlib.org/quadpack/
QUADPACK_DOUBLE Numerical Integration http://people.sc.fsu.edu/~jburkardt/f_src/quadpack_double/quadpack_double.html
NUMPAC http://www2.itc.nagoya-u.ac.jp/center/ja/numpac/
Ooura’s Mathematical Software Packages http://www.kurims.kyoto-u.ac.jp/~ooura/i
シキノートhttps://slpr.sakura.ne.jp/qp/nc-integral/

kahaner’s 21 test problemに関して
Takemitsu Hasegawa, Susumu Hibino, Yohsuke Hosoda and Ichizo Ninomiya,
“An extended doubly-adaptive quadrature method based on the combination of the Ninomiya and the FLR schemes”

幸谷智紀、鈴木千里 “補外法に基づく数値積分法の実装と性能評価”

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

離散フーリエ変換と畳み込みの計算

コメントする

こちらのページは後程消去いたします。

以下のページに統合しますので、ご参照ください。
https://slpr.sakura.ne.jp/qp/dft-for-numerical-calculation/

本ページにはミスもあり、上の方が正しいですのでご参照ください。

まとめ

離散フーリエ変換(Discrete Fourier Transform,DFT)を
\(
\begin{eqnarray}
\left\{
\begin{aligned}
f(\tilde{k}_n)&=\mathcal{F}[f(x)](\tilde{k}_n)=\sum_{m=0}^{N-1}f(x_m)e^{-i2\pi\tilde{k}_nx_m}\Delta x \\
f(x_m)&=\mathcal{F}^{-1}[f(\tilde{k})](x_m)=\sum_{n=0}^{N-1}f(\tilde{k}_n)e^{i2\pi\tilde{k}_nx_m}\Delta \tilde{k}
\end{aligned}
\right.
\end{eqnarray}
\)
と定義し、畳み込み(Convolution)を
\(
\displaystyle (f\ast g)(x_l)= \sum_{m} f(x_m) g(x_l-x_m) \Delta x
\)
と定義します。この時、畳み込みは上で定義した離散フーリエ変換を用いて以下の通り書けます。
\(
\displaystyle (f\ast g)(x_l)=\mathcal{F}^{-1}\left[\mathcal{F}[f(x)](\tilde{k})\cdot \mathcal{F}[g(x)](\tilde{k})\right](x_l)
\)

intel®マス・カーネル・ライブラリ(MKL)に実装されているフーリエ変換は
\(
\begin{eqnarray}
\left\{
\begin{aligned}
(-1)^n f(\tilde{k}_n)&=\sum_{m=0}^{N-1}f(x_m)e^{-i2\pi nm/N} \\
f(x_m)&=\sum_{n=0}^{N-1}(-1)^n f(\tilde{k}_n)e^{i2\pi nm/N}
\end{aligned}
\right.
\end{eqnarray}
\)
の形で定義されています。ここで、\(f(\tilde{k}_n)\)は、正しいフーリエ変換後の関数。計算したもの。この時、上で定義されている場合、畳み込みは、
\(
\displaystyle (f\ast g)(x_m)=\frac{\Delta x}{N}\mathcal{F}_{MKL}^{-1}\left[\mathcal{F}_{MKL}[f(x)]\cdot\mathcal{F}_{MKL}[g(x)] (-1)^n\right](x_m)
\)
で計算できます。ここで、\(N\)は離散フーリエ変換に用いた分点数、\(n\)は順方向フーリエ変換後の空間を離散化した時のインデックスです。

intel®MKLの離散フーリエ変換ルーチンを用いた畳み込みは、プログラムでは\(h(x)=(f\ast g)(x)\)とした時、

\(
\begin{align}
& Status = DftiComputeForward(hand,f) \\
& Status = DftiComputeForward(hand,g) \\
& h(\tilde{k}_n)=z(\tilde{k}_n)\times g(\tilde{k}_n)\times (-1)^n ~~~~~(\mbox{for all } k_n) \\
& Status = DftiComputebackward(hand,h) \\
& h(x_m)=h(x_m)\times \frac{\Delta x}{N} ~~~~~(\mbox{for all } x_m)\\
\end{align}
\)

とする必要があります。
※離散フーリエ変換はフーリエ変換の積分を短冊近似したものです。

[adsense1]

目次

  1. フーリエ変換の定義
  2. 離散フーリエ変換の定義
  3. たたみ込みの定義
  4. 離散たたみ込みと離散フーリエ変換
  5. 数値計算で使われている離散フーリエ変換アルゴリズム (例: intel®MKL)

フーリエ変換の定義

フーリエ変換の定義方法には複数の慣例があり、
数値計算分野、化学、音の解析等においては
\(
\begin{align}
f(\tilde{k})&=\int_{-\infty}^{\infty} f(x)e^{-i2\pi\tilde{k}x}dx \\
f(x)&=\int_{-\infty}^{\infty} f(\tilde{k})e^{i2\pi\tilde{k}x}d\tilde{k}
\end{align}
\)
の形で良く定義されます。
ここで、例えば\(x\)は位置[m]であれば\(\tilde{k}\)は波数[1/m]です。
また、\(x\)が時間[s]であれば\(\tilde{k}\)は周波数[1/s]です。

このページでは扱いませんが、物理科学の世界では
\(
\begin{align}
f(k)&=\int_{-\infty}^{\infty} f(x)e^{-ikx}dx \\
f(x)&=\int_{-\infty}^{\infty} f(k)e^{ikx}\frac{dk}{2\pi}
\end{align}
\)
の形で良く定義されます。
ここで、例えば\(x\)は位置[m]であれば\(k\)は角波数[1/m]です。
また、\(x\)が時間[s]であれば\(k\)は角周波数[1/s]です。

この２つの違いは周波数か、角周波数のどちらが本質であるか？という事だと思います。
物理では、周波数よりも角周波数の方が式が簡便になります。
ですが、人間の理解がしやすいのは1秒当たり何回振動するか？という周波数だろう、と思うため工学よりの分野では周波数が使われるのでしょう。もしかしたら、単に慣例だけかもしれません。

どちらも定数倍の変数変換の違いだけなので、本質は変わりません。
数値計算依りの話をしていきたいので、ここからは前者の位置⇔波数で考えていきます。

また、

順方向フーリエ変換：\(e^{-i2\pi\tilde{k}x}\)の因子が掛かっている式、空間→波数
逆方向フーリエ変換：\(e^{+i2\pi\tilde{k}x}\)の因子が掛かっている式、波数→空間

と呼び、
”関数\(f(x)\)の順方向フーリエ変換して波数\(\tilde{k}\)の空間に移る”
という操作を華文字\(\mathcal{F}\)を用いて
\(
\mathcal{F}[f(x)](\tilde{k})
\)
と表します。また、
”関数\(f(k)\)の逆方向フーリエ変換して位置\(x\)の空間に移る”
という操作を
\(
\mathcal{F}^{-1}[f(\tilde{k})](x)
\)
と表します。

離散フーリエ変換の定義

このページでは離散フーリエ変換を
\(
\begin{align}
f(\tilde{k}_n)&=\sum_{m=0}^{N-1}f(x_m)e^{-i2\pi\tilde{k}_nx_m}\Delta x \\
f(x_m)&=\sum_{n=0}^{N-1}f(\tilde{k}_n)e^{i2\pi\tilde{k}_nx_m}\Delta \tilde{k}
\end{align}
\)
と定義します。
\(\Delta x\), \(\Delta \tilde{k}\)は、それぞれ位置、波数空間の刻み幅を表します。

この上の定義は簡便さのために良く使われる
\(
\begin{align}
f(\tilde{k}_n)&=\frac{1}{N}\sum_{m=0}^{N-1}f(x_m)e^{-i2\pi\tilde{k}_nx_m} \\
f(x_m)&=\hspace{1.5em}\sum_{n=0}^{N-1}f(\tilde{k}_n)e^{i2\pi\tilde{k}_nx_m}
\end{align}
\)
という定義と同一です。
ただし、これは積分をそのまま離散化したものではないため、積分の意味は薄れてしまいます。
その結果、畳み込みなど、フーリエ変換後同士の物を掛け算する時は違った結果になってしまいます。
このページではこの簡便さの為の離散フーリエ変換は使いません。

上の２つの定義が同じ理由は、順方向→逆方向を行い、元に戻った時に係数が\(\Delta x\Delta \tilde{k}=\frac{1}{N}\)になるためです。

\(x_m\) \(\tilde{k}_n\)
添え字 \(m=0,1,2,\cdots,N-1\) \(n=0,1,2,\cdots,N-1\)
区間 \(\displaystyle [a,b]\) \(\displaystyle \left[-\frac{1}{2\Delta x},\frac{1}{2\Delta x}\right]\)
刻み幅 \(\displaystyle \Delta x=\frac{b-a}{N}\) \(\displaystyle \Delta \tilde{k}=\frac{1}{b-a}\)
分点の位置 \(\displaystyle x_m=m\Delta x +a\) \(\displaystyle \tilde{k}_n=n\Delta \tilde{k} -\frac{1}{2\Delta x}\)

ただし、\(x_m, \tilde{k}_n\)は以下の関係式を見たします。
\(
\begin{align}
& x_m+N\Delta x=x_m\\
& \tilde{k}_n+N\Delta \tilde{k}=\tilde{k}_n
\end{align}
\)

ナイキスト(Nyquist)周波数について

ナイキスト周波数とはサンプリングする間隔でどの周波数を持つ関数までを表現できるか？という下限と上限を与えます。
今回の場合、上の表の通り\(\tilde{k}\)空間で表現できる上限と下限は\(\pm\frac{1}{2\Delta x}\)です。この絶対値、すなわち\(\frac{1}{2\Delta x}\)をナイキスト周波数と呼びます。

ナイキスト周波数が意味しているのは、高周波の成分を取り出したければサンプリングする間隔を小さくしなければならない。ということです。”サンプリングの間隔”は周波数と同じ意味ですので、サンプルする周波数を高くしなければならない。と言い換えることもできます。

具体例を上げましょう。

上の図は関数\(\displaystyle f(x)=e^{i2\pi \tilde{k}x}\)をサンプル間隔\(\Delta x = 0.5\)で離散的な関数で表現したものです。
この場合、ナイキスト周波数は\(\frac{1}{2\times 0.5}\)より、\(1\)です。
すなわち、このサンプル間隔では\(\tilde{k}\)が\([-1,1]\)の範囲を持つものしか表現することが出来ません。
では、実際にサンプル出来ない関数\(\displaystyle f(x)=e^{i2\pi 1.8 x},~~(\tilde{k}=1.8)\)を考えた時に何が起こるのでしょうか。
それを示したのが下から2番目の関数です。
本来は早く振動する関数であるのに、あたかも非常にゆっくりとした関数だ、と捉えられてしまいます。しかも位相が反転しています。

これは、元々は\(e^{i\theta}\)の周期性によるものです。簡単な計算の通り、
\(
\begin{align}
e^{i(\theta+2\pi)}=e^{i\theta}e^{i2\pi}=e^{i\theta}
\end{align}
\)
となります。この周期性のために、\(\tilde{k}:[-1,1]\)の範囲でしか表現できない空間で、
その範囲外の\(\tilde{k}\)を表現しようとした時には、\(\tilde{k}\)空間の大きさ（…この場合は\(2\)）を足したり引いたりして\(\tilde{k}:[-1,1]\)の範囲に強制的に押し込まれます。
\(\tilde{k}=1.8\)の場合は、\(\tilde{k}=1.8-2=-0.2\)となります。
言葉で言えば、
\(\tilde{k}=1.8\)を持つ波形は\(\tilde{k}=-0.2\)を持つ波形としてとして表現される
ということです。

本来の関数を薄く描いて重ねたものが以下のものになります。

フーリエ変換が元に戻ることの証明

上記のフーリエ変換、離散フーリエ変換はある制限を課しています。
それは、
ある関数に順方向フーリエ変換を行い、続いて逆方向フーリエ変換を施した場合、元と同じ関数に戻っていて欲しい。
という要望です。
それは、以下のように証明できます。

フーリエ変換
\(
\begin{align}
f(x)&=\int_{-\infty}^{\infty}d\tilde{k} f(\tilde{k})e^{i2\pi\tilde{k}x} \\
&=\int_{-\infty}^{\infty}d\tilde{k} \left[\int_{-\infty}^{\infty}dx’ f(x’)e^{-i2\pi\tilde{k}x’}\right]e^{i2\pi\tilde{k}x} \\
&=\int_{-\infty}^{\infty}dx’ f(x’)\int_{-\infty}^{\infty}d\tilde{k}
e^{i2\pi\tilde{k}(x-x’)} \\
&=\int_{-\infty}^{\infty}dx’ f(x’)\delta(x-x’) \\
&=f(x)
\end{align}
\)

離散フーリエ変換

\(
\begin{align}
f(x_m)&=\sum_{n=0}^{N-1}f(\tilde{k}_n)e^{i2\pi\tilde{k}_nx_m}\Delta \tilde{k}\\
&=\sum_{n=0}^{N-1}\left[\sum_{m’=0}^{N-1}f(x_{m’})e^{-i2\pi\tilde{k}_nx_{m’}}\Delta x\right]e^{i2\pi\tilde{k}_nx_m}\Delta \tilde{k} \\
&=\Delta \tilde{k}\Delta x\sum_{m’=0}^{N-1}f(x_{m’})\sum_{n=0}^{N-1}e^{i2\pi\tilde{k}_n(x_m’-x_m)} \\
&=\frac{1}{N}\sum_{m’=0}^{N-1}f(x_{m’})\sum_{n=0}^{N-1}e^{i2\pi n (m’-m)/N +i2\pi(\theta_{m}-\theta_{m’})} \\
&=\frac{1}{N}\sum_{m’=0}^{N-1}f(x_{m’})e^{i2\pi(\theta_{m}-\theta_{m’})}N\delta_{m,m’} \\
&=f(x_m)
\end{align}
\)
ここで、
\(
\begin{align}
\tilde{k}_n x_{m}=\frac{nm}{N}+\frac{an}{N\Delta x}+\theta_{m} \\
\theta_{m}=-\frac{a}{2\Delta x}-\frac{m}{2}
\end{align}
\)
そして、等比級数の考えから、
\(
\begin{eqnarray}
\sum_{n=0}^{N-1}e^{i2\pi n a}=
\left\{
\begin{aligned}
&\frac{e^{i2\pi aN}-1}{e^{i2\pi a}-1}, \hspace{1em} ( a\ne 0) \\
& N, \hspace{4em} (a=0) \\
\end{aligned}
\right.
\end{eqnarray}
\)
という関係式を用いました。

畳み込み(convolution)の定義

このページでは、畳み込みを
\(
\displaystyle (f\ast g)(x)=\int f(\tau) g(x-\tau) d\tau
\)
という積分であると定義し、離散的な畳み込みを
\(
\displaystyle (f\ast g)(x_l)= \sum_{m} f(x_m) g(x_l-x_m) \Delta x
\)
と定義します。

畳み込みの定義には\(\Delta x\)が掛けられていないものが存在しますが、積分の意味が薄れてしまうため、離散フーリエ変換時と同様、このページでは\(\Delta x\)を掛けたものを畳み込み、と定義します。

区間\([a,b]\)で定義された離散的な畳み込みの場合、関数\(f(x),g(x)\)は区間が\([a,b]\)で決められている場合がほとんどです。
この場合、引数\(x_m-x_l\)が定義域内に収まらない場合が現れます。

定義域からはみ出てしまった場合は推定するしかありません。
推定の仕方によって
循環畳み込みと呼ばれる方法と直線畳み込みと呼ばれる方法が良くつかわれています。

循環畳み込みは領域に対して周期境界条件を課したもの、すなわち、上限bと下限aが繋がっていると考えます。
直線畳み込みは領域に対して固定端条件を課したもの、すなわち、[a,b]の範囲外の関数値はゼロと考えます。

問題に大きく依るので、どちらがいいかはありません。

離散畳み込みと離散フーリエ変換

離散畳み込みは定義通り
\(
\displaystyle (f\ast g)(x_l)= \sum_{m} f(x_m) g(x_l-x_m) \Delta x
\)
で計算することが出来ます。計算量は\(O(n^2)\)です。

しかし、実用上は離散フーリエ変換を利用すると計算量を\(O(n\log n)\)に抑えられます。
計算方法はたたみ込み定理を利用して、
\(
\displaystyle \mathcal{F}[(f\ast g)(x)](k)=\mathcal{F}[f(x)](k)\cdot \mathcal{F}[g(x)](k)
\)
より、\((f\ast g)(x)\)は
\(
\displaystyle (f\ast g)(x)=\mathcal{F}^{-1}\left[\mathcal{F}[f(x)](k)\cdot \mathcal{F}[g(x)](k)\right]
\)
と書けるため、フーリエ変換を介することで計算量を減らせます。
離散フーリエ変換を介して畳み込みを計算する場合、必ず循環畳み込みで計算することになります。

畳み込み定理の証明

連続畳み込み

\(
\displaystyle \mathcal{F}[(f\ast g)(x)](k)=\mathcal{F}[f(x)](k)\cdot \mathcal{F}[g(x)](k)
\)

左辺）
\(
\begin{align}
\mathcal{F}[(f\ast g)(x)](\tilde{k})
&=\int \left[\int f(\tau)g(x-\tau) d\tau\right] e^{-i2\pi \tilde{k} x} dx \\
& x-\tau=y\mbox{の変数変換を行って} \nonumber \\
&=\int dy\int d\tau f(\tau)g(y) e^{-i2\pi \tilde{k} (\tau+y)} \\
&=\int dx\int dx’ f(x)g(x’) e^{-i2\pi \tilde{k} (x+x’)}
\end{align}
\)

右辺）
\(
\begin{align}
\mathcal{F}[f(x)](\tilde{k})\cdot \mathcal{F}[g(x)](\tilde{k})
&=\int f(x) e^{-i2\pi \tilde{k} x} dx \cdot \int g(x’) e^{-i2\pi \tilde{k} x’} dx’ \\
&=\int dx\int dx’ f(x)g(x’) e^{-i2\pi \tilde{k} (x+x’)}
\end{align}
\)
より、証明終了。また、逆変換により戻ることを示します。
\(
\begin{align}
&\mathcal{F}^{-1}\left[\mathcal{F}[f(x)](\tilde{k})\cdot \mathcal{F}[g(x)](\tilde{k})\right] \\
&=\int \left[\int dx\int dx’ f(x)g(x’) e^{-i2\pi \tilde{k} (x+x’)}\right] e^{i2\pi \tilde{k} y} dk\\
&=\int dx\int dx’ f(x)g(x’) \int e^{-i2\pi \tilde{k} (x+x’-y)} dk\\
&=\int dx\int dx’ f(x)g(x’) \delta(x’-(y-x)) \\
&=\int dx f(x)g(y-x) \\
&=\int d\tau f(\tau)g(x-\tau) \\
\end{align}
\)
となり、フーリエ変換を介してたたみ込みが計算できることが示せました。

離散畳み込み

\(
\displaystyle \mathcal{F}[(f\ast g)](\tilde{k})=\mathcal{F}[f(x)](\tilde{k})\cdot \mathcal{F}[g(x)](\tilde{k})
\)

左辺）
\(
\begin{align}
\mathcal{F}[(f\ast g)(x)](\tilde{k}_n)
&=\mathcal{F}\left[\sum_{m=0}^{N-1}f(x_m)g(x-x_m)\Delta x\right] \\
&=\sum_{l=0}^{N-1}\left[\sum_{m=0}^{N-1}f(x_m)g(x_l-x_m)\Delta x\right]e^{i2\pi \tilde{k}_n x_l}\Delta x \\
&=(\Delta x)^2\sum_{l=0}^{N-1}\sum_{m=0}^{N-1}f(x_m)g(x_l-x_m) e^{i2\pi \tilde{k}_n x_l} \\
&\mbox{ここで、循環たたみ込みを想定し、}x_l-x_m=x_{m’}\mbox{となるように}x_{m’}\mbox{を決めると}\nonumber \\
&=(\Delta x)^2\sum_{m=0}^{N-1}\sum_{m’=0}^{N-1}f(x_m)g(x_{m’}) e^{i2\pi \tilde{k}_n (x_m+x_{m’})} \\
\end{align}
\)

右辺）
\(
\begin{align}
\mathcal{F}[(f(x)](k_n)\cdot \mathcal{F}[(g(x)](k_n)
&=\sum_{m=0}^{N-1}f(x_m) e^{i2\pi \tilde{k}_n (x_m)}\Delta x\cdot \sum_{m’=0}^{N-1}f(x_{m’}) e^{i2\pi \tilde{k}_n (x_{m’})}\Delta x \\
&=(\Delta x)^2\sum_{m=0}^{N-1}\sum_{m’=0}^{N-1}f(x_m)g(x_{m’}) e^{i2\pi \tilde{k}_n (x_m+x_{m’})}
\end{align}
\)
よって右辺と左辺が正しくなります。

逆変換によって離散畳み込みになることを示します。

\(
\begin{align}
\mathcal{F}^{-1}[\mathcal{F}[(f\ast g)(x)](k)](x_l)
&=\mathcal{F}^{-1}\left[(\Delta x)^2\sum_{m=0}^{N-1}\sum_{m’=0}^{N-1}f(x_m)g(x_{m’}) e^{i2\pi \tilde{k}_n (x_m+x_{m’})}\right]\\
&=\sum_{n=0}^{N-1}\left[(\Delta x)^2\sum_{m=0}^{N-1}\sum_{m’=0}^{N-1}f(x_m)g(x_{m’}) e^{i2\pi \tilde{k}_n (x_m+x_{m’})}\right]e^{-i2\pi \tilde{k}_n x_l}\Delta \tilde{k}\\
&=(\Delta x)^2 \Delta \tilde{k} \sum_{m=0}^{N-1}\sum_{m’=0}^{N-1} f(x_m)g(x_{m’})\sum_{n=0}^{N-1}e^{i2\pi \tilde{k}_n (x_m+x_{m’}-x_l)} \\
&=\frac{\Delta x}{N} \sum_{m=0}^{N-1}\sum_{m’=0}^{N-1} f(x_m)g(x_{m’})N\delta_{x_{m’},x_l-x_m} \\
&=\sum_{m=0}^{N-1}f(x_m)g(x_l-x_m)\Delta x \\
&=(f\ast g)(x_l)
\end{align}
\)

よって畳み込みの定理通りに戻るため、離散フーリエ変換を介して畳み込みの計算ができます。

※注意
\(\Delta x, \Delta \tilde{k}\)が掛けられていない方の離散フーリエ変換と畳み込みでは元に戻りません。畳み込みの計算時に係数倍違って出てきます。

数値計算で使われている離散フーリエ変換アルゴリズム (例: intel®MKL)

数値計算では、上記の離散フーリエ変換の形、指数の肩に\(\tilde{k},x\)を含む物はあまり使われません。
そして、数値計算上での離散フーリエ変換の定義方法は複数存在します。
マニュアルを見てください。

ここでは、intel®マス・カーネル・ライブラリ(MKL)に使われている離散フーリエ変換アルゴリズムについて言及します。

MKLの変換はどうやら
\(
\begin{align}
(-1)^n f(\tilde{k}_n)&=\sum_{m=0}^{N-1}f(x_m)e^{-i2\pi nm/N} \\
f(x_m)&=\sum_{m=0}^{N-1}(-1)^n f(\tilde{k}_n)e^{i2\pi nm/N}
\end{align}
\)
という定義に従い変換しているようです。規格化因子\(\frac{1}{N}\)はユーザー任せです。
ここで、
\(f(\tilde{k}_n)\)は本来のフーリエ変換後の関数、
\(f(x)\)は\(x\)空間での関数
を表します。

”しているようです”と書いたのは、マニュアルと異なっているからです。マニュアルでは\((-1)^n\)など存在しませんが、実際試してみると、こういう変換をしているからです。

MKLではどういう\(\tilde{k}\)を想定しているかを考えたところ、以下の状況のようです。

\(x_m\) \(\tilde{k}_n\)
添え字 \(m=0,1,2,\cdots,N-1\) \(n=0,1,2,\cdots,N-1\)
区間 \(\displaystyle [-a,a]\) \(\displaystyle \left[0,N\Delta \tilde{k}\right]\)
刻み幅 \(\displaystyle \Delta x=\frac{2a}{N}\) \(\displaystyle \Delta \tilde{k}=\frac{1}{2a}\)
分点の位置 \(
\begin{eqnarray}
x_m=\left\{
\begin{aligned}
& m\Delta x ~~~~(m=0,\cdots,\frac{N}{2}-1) \\
& m\Delta x-\frac{1}{\Delta \tilde{k}} (m=\frac{N}{2},\cdots,N-1)
\end{aligned}
\right.
\end{eqnarray}
\)		 \(\displaystyle \tilde{k}_n=n\Delta \tilde{k} ~~~~(n=0,\cdots,N-1)\)

ただし、\(x_m, \tilde{k}_n\)は以下の周期境界条件を見たします。
\(
\begin{align}
& x_m+N\Delta x=x_m\\
& \tilde{k}_n+N\Delta \tilde{k}=\tilde{k}_n
\end{align}
\)

区間[-a,a]で定義された\(x_m\)の順番がおかしく感じますが、これは周期境界条件を用いて、
区間[0,2a]で定義された\(x_m=m\Delta x ~~~~(m=0,\cdots,N-1)\)
と見ても良いです。
なので、実用上どうやって関数\(f(x)\)を入れればよいか？に困ることは無いでしょう。
この上記離散フーリエ変換を考えると、上で提示したMKLの離散フーリエ変換が導くことが出来、プログラム上でも一致します。

MKLの離散フーリエ変換の導出）

離散フーリエ変換時の注意点

上に記した、余分な係数\((-1)^n\)のためにMKLを利用して求められた離散フーリエ変換の値は、本来のフーリエ変換後の関数と異なります。
具体的に例を示しましょう。

\(e^{-\alpha x^2}\)のフーリエ変換は
\(
\displaystyle \int_{-\infty}^{\infty} e^{-\alpha x^2}e^{-i2\pi\tilde{k} x}dx=\sqrt{\frac{\pi}{\alpha}}e^{-\pi^2\tilde{k}^2/\alpha}
\)
です。離散フーリエ変換を定義通り計算したものと、MKLを用いて離散フーリエ変換すると以下のようになります。
左から\(f(x), \mathcal{F}[f(x)](\tilde{k}),\mathcal{F}^{-1}[\mathcal{F}[f(x)]](x)\)になっています。

ここで、因子\((-1)^n\)のために交互に符号が変わっていることが分かるでしょう。
後でやりますが、このせいで畳み込みの計算をする際には気を付けなければなりません。

計算コードはこちら。

MKLのルーチンは優秀です。誤差が余りたまらないようです。
定義通り、愚直に作ったプログラムでは計算速度が遅いだけでなく、計算精度も悪くなります。
下の画像は、離散フーリエ変換→逆変換を複数回繰り返した時、元の関数とどれだけずれるかを表しています。

定義通りの離散フーリエ変換を\(10^n\)回繰り返すと\(n\)に対し線形にあがっていきます。この誤差は丸め誤差に起因しているのだと思います。しかし、MKLのルーチンには起こっていません。どういうアルゴリズムか知りませんが、優秀なのでしょう。

[adsense2]

フーリエ変換を介する畳み込み時の注意点

MKLの離散フーリエ変換ルーチンを使って畳み込みを計算します。
例として、以下の関数を考えましょう。
\(
\begin{align}
h(x)&=\int \tau^5 e^{-\tau^2}\cdot e^{-4(x-\tau)^2}d\tau \\
&=\frac{1}{3125}\sqrt{\frac{\pi}{5}}x(1024x^4+1600x^2+375)e^{-\frac{4}{5}x^2}
\end{align}
\)
ここで、
\(
\begin{align}
f(\tau)=\tau^5 e^{-\tau^2}
g(\tau)=e^{-4\tau^2}
\end{align}
\)
です。

このページで紹介した定義の畳み込みでは証明した通り、
\(
h(x_m)=\mathcal{F}^{-1}\left[\mathcal{F}[f(x)](\tilde{k})\cdot\mathcal{F}[g(x)](\tilde{k})\right](x_m)
\)
で計算できます。
しかし、MKLの離散フーリエ変換(\(\mbox{DFT}_{MKL}\))の場合は簡単にはできません。計算は、
\(
\displaystyle (f\ast g)(x_m)=\frac{\Delta x}{N}\mathcal{F}_{MKL}^{-1}[\mathcal{F}_{MKL}[f(x)]\cdot\mathcal{F}_{MKL}[g(x)] (-1)^n](x_m)
\)
としなければなりません。\((-1)^n\)は順方向フーリエ変換が2回あり、MKLの逆方向離散フーリエ変換に必要なこの係数が落ちてしまうため入れなければなりません。また、MKLのルーチンでは係数\(\Delta x, \Delta\tilde{k}\)を含んでいないので計算しておく必要があります。

intel®MKLを用いる際、
順方向離散フーリエ変換にはDftiComputeForward、
逆方向離散フーリエ変換にはDftiComputeBackward
を使うと、擬コードはこのように書けるかと思います。

\(
\begin{align}
& Status = DftiComputeForward(hand,f) \\
& Status = DftiComputeForward(hand,g) \\
& h(\tilde{k}_n)=z(\tilde{k}_n)\times g(\tilde{k}_n)\times (-1)^n ~~~~~(\mbox{for all } k_n) \\
& Status = DftiComputebackward(hand,h) \\
& h(x_m)=h(x_m)\times \frac{\Delta x}{N} ~~~~~(\mbox{for all } x_m)\\
\end{align}
\)

\(\displaystyle \Delta x \Delta\tilde{k}=\frac{1}{N}\)の関係を用いて最後にまとめて定数倍しています。

上の計算方法に従うと正しく畳み込みが計算出来ます。
実際にやってみますと、

という結果を得ます。

畳み込みを計算するコードはこちら

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

Clenshaw–Curtis求積法

5件のコメント

クレンショウ-カーティス(Clenshaw–Curtis)型の数値積分公式です。

[2]の論文”IS GAUSS QUADRATURE BETTER THAN CLENSHAW–CURTIS?”の結論だけかいつまみますと、

・クレンショウ-カーティスの求積法は、基本的にはガウス-ルジャンドル求積法と同じ性能
・クレンショウ-カーティスの求積法は実装が簡単

ということです。
本気で計算精度を求めようとするならばガウス-ルジャンドル求積法、簡単に、でもかなりの精度で求めたければクレンショウ-カーティスの求積法なのだと思います。
ガウス-ルジャンドル求積法は分点の位置と重みの計算が簡単ではないのです。

積分
\(
\displaystyle \int_{-1}^1 f(x) dx \sim \sum_{k=0}^{n} \omega_k f(x_k)
\)
Nが偶数の時、分点\(x_k\)と重み\(\omega_k\)は
\(
\displaystyle x_k=\cos\left(\frac{k\pi}{N}\right),\;\;\;k=0,1,\cdots,N
\)
\(
\begin{align}
\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{\prime\prime }\)は和の最初と最後の項に\(\frac{1}{2}\)を乗じることを意味します。

N点の分点に対し、誤差のオーダーは\(O(N+1)\)です。
\(x_k\)はチェビシェフ多項式\(T_n(x)\)の極値点になっています。

利点

・重みの計算に掛かる時間が\(O(N log N)\).
・2N点での積分値を求める際にN点で計算した評価点を用いることが出来る。

ガウス求積法と比べて

・ガウス型の重みの計算に掛かる時間は\(O(N^2)\).ただし、分点と重みは一度求めてしまえばずっと使いまわせるので、さほど重要な問題ではありません。
・ガウス-ルジャンドル求積法の分点を\(N\)点利用し、それに\(N\)点追加することでガウス-クロンロッド求積法が使えます。ガウス-ルジャンドルはN点の求積で\(O(x^{2N})\)、ガウス-クロンロッドは2N点の求積で\(O(x^{3N})\)の多項式まで正確に表すことが出来ます。

結論

本気で計算精度が欲しいならばガウス-ルジャンドル求積法。
簡単に、でもかなりの精度で求めたければクレンショウ-カーティスの求積法なのだと思います。

プログラム

最適化は特にしていませんが、fortran90ではこんな感じになります。

program main
  implicit none
  integer::i,j,N
  double precision::a,b,s
  double precision,allocatable::x(:),w(:)
  double precision,parameter::pi=acos(-1d0)
  double precision::f
  external::f
 
  a=0d0
  b=pi
   
  N=11
  allocate(x(1:N),w(1:N))
  call clenshaw_curtis_xw(N,x,w)
  w=w*(b-a)*0.5d0
  x=(x*(b-a)+(b+a))*0.5d0

  s=0d0
  do i=1,N
     s=s+f(x(i))*w(i)
  enddo
  write(6,*)s

  stop
end program main
   
function f(x)
  implicit none
  double precision::f,x
  double precision::pi=dacos(-1d0)
 
  f=sin(x)*x
 
  return
end function f 

subroutine clenshaw_curtis_xw(N,x,w)
  implicit none
  integer,intent(in)::N
  double precision,intent(out)::x(1:*),w(1:*)
  double precision,parameter::pi=3.14159265358979323846264338327950288d0
  integer::i,j,M,k
  double precision::c
 
  if(mod(N,2).eq.0)then
     write(6,*)"N must be odd"
     stop
  endif
 
  M=N-1
 
  do i=0,M
     x(i+1)=cos(i*pi/dble(M))
  enddo
  w(1)=1d0/((M-1d0)*(M+1d0))
  w(M+1)=w(1)

  do i=1,M/2
     k=i+1
     w(k)=0.5d0*(1d0+cos(pi*i)/((1d0-M)*(1d0+M)))
     c=2d0*pi*i/M
     do j=1,M/2-1
        w(k)=w(k)+cos(c*j)/((1d0-2d0*j)*(1d0+2d0*j))
     enddo
     w(k)=w(k)*4d0/M
     
     w(M-i+1)=w(k)
  enddo

  return
end subroutine clenshaw_curtis_xw

参考文献

[1]P. J. Davis, P. Rabinowitz著、森正武訳 『計算機による数値積分法』p. 131-132,初版 (1980)
[2]LLOYD N. TREFETHEN, IS GAUSS QUADRATURE BETTER THAN CLENSHAW–CURTIS?

プログラミングと数値計算, 物理学

時間依存しないシュレーディンガー方程式と変分原理 2/2 （計算例）

コメントする

計算例

具体的に適当なポテンシャルの形を考えて解いてみましょう。

プログラムは
時間依存しないシュレーディンガー方程式と変分原理 1/2
にあります。

おさらいですが、間接法による変分原理は、

  • 固有値の上限を与える（変分原理によるエネルギー固有値は必ず厳密なエネルギー固有値よりも大きくなる）
  • 基底関数の数を増やすほど精度は上がる(完全系の場合)

ということを想定しています。
そして、時間依存しないシュレーディンガー方程式と変分原理 1/2 （計算手法の説明）
では、これから束縛状態のみを考えるので、sin関数系の場合の変分原理について計算を行っていきました。
上に想定した変分原理の想定は、

正しい境界条件の下で誤差なく行列要素が求められ、誤差なく対角化が行われた

場合にのみ成り立ちます。sin関数系の場合の境界条件は固定端条件、つまり区間[a,b]で
\(
\psi(x=a)=\psi(x=b)=0
\)
ということです。

ここでは、解析的な解が存在する

  1. 井戸型ポテンシャル
  2. 三角量子井戸
  3. 調和ポテンシャル
  4. F.モースによるポテンシャル[3]
  5. \(1/(\cosh^2)\)型ポテンシャル[3]
  6. クーロンポテンシャル(動径方向, \(l=0\))

について比較していきます。

シュレーディンガー方程式は全て原子単位系(\(m=1, \hbar=1\))で考え、
\(
\displaystyle \left[-\frac{1}{2}\frac{d^2}{dx^2}+V(x)\right]\psi(x)=E\psi(x)
\)
の形を考えます。

無限に深い井戸型ポテンシャル

ポテンシャル

\(
V(x)=\left\{
\begin{eqnarray}
0 \hspace{2em}&(& 0\lt x \lt L,)\\
\infty \hspace{2em}&(& x\le 0 \ or\ L\le x)
\end{eqnarray}
\right.
\)
※井戸の区間が\([a,b]\)であっても変数変換\(y=x-a,\;\;L=b-a\)によって上記ポテンシャルの形に変換されます。

固有値

\(
\displaystyle E_n=\frac{(n+1)^2\pi^2}{2 L^2},\;\;(n=0,1,2,\cdots)
\)

固有関数

\(
\psi_n(x)=\left\{
\begin{eqnarray}
\sqrt{\frac{2}{L}}\sin\left(\frac{(n+1)\pi}{L}(x+\frac{L}{2})\right) &(& 0\lt x \lt L,)\\
0 \hspace{6em}&(& x\le 0 \ or\ L\le x)
\end{eqnarray}
\right.
\)

解法

区間\(0\lt x \lt L,\)ではシュレーディンガー方程式の\(V(x)=0\)の解、sin,cosの形をしてなければなりません。
また、\(\psi(x=0)=0,\;\psi(x=L)=0\)では境界条件より、波動関数の値はゼロです。
よって、解はsin,cosの形で、そのゼロ点がちょうど\(x=0,x=L\)にくるような関数となります。

数値計算解と解析解との比較

数値計算の設定:
Gauss-Lobatto求積法の次数(Ngl):12
計算区間(a,b):[0,\(\pi\)]
基底関数の数(N)：200

エネルギー固有値

量子数\(n\) 数値計算解 解析解
0 0.5000000000000000 0.5000000000000000
1 2.0000000000000000 2.0000000000000000
2 4.5000000000000000 4.5000000000000000
3 8.0000000000000000 8.0000000000000000
4 12.5000000000000000 12.5000000000000000
…
10 60.500000000000000 60.500000000000000
20 220.50000000000000 220.50000000000000
30 480.50000000000000 480.50000000000000
40 840.50000000000000 840.50000000000000
50 1300.5000000000000 1300.5000000000000
考察

行列は元々対角です。なので、対角化する前から対角行列です。
誤差は入りようがありません。

三角量子井戸

ポテンシャル

\(
V(x)=\left\{
\begin{eqnarray}
\infty \hspace{2em}&(& x\lt 0)\\
\alpha x \hspace{2em}&(& 0 \lt x)
\end{eqnarray}
\right.
\)

固有値

\(
\displaystyle E_n=-\left(\frac{\alpha^2}{2}\right)^{1/3}a_n
\)
ここで、\(a_n\)はAiry関数[4]のゼロ点です。

固有関数

\(
\displaystyle \psi_n(x)=C\cdot \text{Ai}\left(\left(\frac{2}{\alpha^2}\right)^{1/3}(\alpha x-E_n)\right)
\)
※積分区間がエアリー関数のゼロ点から無限までであり、これは解析的に実行できません。なので規格化定数をあらわに書くことはできません。

解法

区間\(0\lt x \lt \infty,\)ではシュレーディンガー方程式の\(V(x)=\alpha x\)の解、Airy関数の形をしてなければなりません。
また、\(\psi(x=0)=0, \psi(x\to \infty)=0\)の境界条件を満たします。
よって、解はAiry関数の形で、\(x\)が漸近で減衰していく方の解(\(Ai(x),Bi(x)\)のAiの方)で、そのゼロ点がちょうど\(x=0\)にくるような関数となります。

数値計算解と解析解との比較

数値計算の設定:
Gauss-Lobatto求積法の次数(Ngl):12
計算区間(a,b):[0,60]
基底関数の数(N)：200
\(\alpha\)：1

エネルギー固有値

量子数\(n\) 数値計算解 解析解
0 1.8557572021449626 1.8557570814892386
1 3.2446077578572123 3.2446076240031596
2 4.3816713853293896 4.3816712392861303
3 5.3866139387958034 5.3866137807905003
4 6.3052631766282223 6.3052630065857747
…
10 10.8669422991860358 10.8669420487522821
20 16.8461591354108187 16.8461586901918032
30 21.8969186983512394 21.8969179157570224
40 26.4182318860569758 26.4182304793452509
50 30.5803380635310234 30.5803354172938597

考察

上記計算では基底状態から励起状態までほとんど精度が変わっていません。
また、変分原理が想定したとおり、全ての数値計算解は、厳密な解よりも高い結果を与えています。
厳密な解よりも高い結果を与えていることから、
行列要素は高精度で計算されているが、基底関数が足らない、もしくは計算区間が足らない
事が考えられます。

精度がなぜ良くないのかを考えましょう。
\(x=0\)の境界条件は\(\sin\)関数で良く表現されます。が、\(x\to\infty\)側は指数関数で減少していきます。
(一次の)指数関数で減衰していくことを\(\sin\)関数系ではうまく表現できないのではないか、と思います。
実際、Airy関数の指数関数での減衰は漸近で
\(
\displaystyle \frac{e^{-\frac{2}{3}x^{3/2}}}{x^{1/4}}
\)
と表されます。
この後、調和ポテンシャルとクーロンポテンシャルの場合も考えますが、
調和ポテンシャル(漸近で\(\exp(-x^2)\)、基底状態の精度14,15桁)
↓
三角量子井戸(漸近でおおよそ\(\exp(-x^{2/3})\)、基底状態の精度7,8桁)
↓
クーロンポテンシャル(漸近で\(\exp(-x)\)、基底状態の精度2,3桁)
の順に精度が悪くなっていきます。
指数関数的な減衰が早い≒無限に深い量子井戸
を意味するはずなので、恐らくこうなのではと思います。

調和ポテンシャル


\(
\displaystyle V(x)=\frac{1}{2}x^2
\)

固有値

\(
\displaystyle E_n=n+\frac{1}{2},\;\;(n=0,1,2,\cdots)
\)

固有関数

\(
\displaystyle \psi_n(x)=\frac{\pi^{-1/4}}{\sqrt{2^n n!}}e^{-\frac{x^2}{2}} H_n\left(x\right)
\)
ここで\(H_n(x)\)はエルミート多項式を表します。

解法

1次元調和振動子の直接解法

数値計算解と解析解との比較

数値計算の設定:
Gauss-Lobatto求積法の次数(Ngl):12
計算区間(a,b):[-30,30]
基底関数の数(N)：200

量子数\(n\) 数値計算解 解析解
0 0.5000000000000038 0.5
1 1.5000000000000340 1.5
2 2.4999999999998268 2.5
3 3.4999999999996776 3.5
4 4.4999999999999165
…
10 10.4999999999999680 10.5
20 20.4999999999997691 20.5
30 30.5000000000065015 30.5
40 40.5000064794500432 40.5
50 50.5426649449273313 50.5
考察

変分原理の典型的な傾向を見せています。というのは高い励起状態になればなるほど精度が悪くなっているためです。
高励起状態になるにつれて精度が悪くなっていくのは基底関数が足らない為でしょう。計算区間に関してはポテンシャルは境界で高いまま(\(V(a)=V(b)=450\))であり、求められたエネルギー固有値の値\(E_{50}\approx 50\)と比べても十分大きいです。なので、計算区間は足りています。

また、変分原理の固有値の上限を与えるに関しては必ずしも満たしているわけでは無いことに気が付きます。これは行列要素を求める際の誤差や丸め誤差に起因するものであると考えられます。なので、原理が破られているのではなく、数値計算上の問題に起因しています。\(n=30,40,50\)では厳密解よりも大きな値を与えています。これは、数値計算上の誤差よりも変分原理の原理的な値の方が大きくなっているためです。

F.モースによるポテンシャル[3]


\(
\displaystyle V(x)=A\left(e^{-2\alpha x}-2e^{-\alpha x}\right)
\)

固有値

\(
\displaystyle E_n=-A\left[1-\frac{\alpha}{\sqrt{2A}}\left(n+\frac{1}{2}\right)\right]^2
\)
ここで、\(n\)は正の整数で、ゼロから始まり、
\(
\displaystyle\frac{\sqrt{2A}}{\alpha}\gt n+\frac{1}{2}
\)
を満足する最大値\(n_{\text{max}}\)まで。

固有関数

\(
\begin{align}
\psi(x)&=e^{-\xi/2} \xi^s w(\xi) \\
& w(\xi)=F(-n,2s+1,\xi) \\
& \xi=\frac{2\sqrt{2A}}{\alpha}e^{-\alpha x} \\
& s=\frac{\sqrt{-2E}}{\alpha}
\end{align}
\)
ここで\(F(a,b,x)\)は合流型超幾何関数です。

この問題では離散スペクトルは有限個だけ存在します。もしも
\(
\displaystyle\frac{\sqrt{2A}}{\alpha}\gt n+\frac{1}{2}
\)
であれば、離散スペクトルは一般に存在しません。

解法

[3]を参照してください

数値計算解と解析解との比較

数値計算の設定:
Gauss-Lobatto求積法の次数(NGL):12
計算区間(a,b):[-3,80]
基底関数の数(N)：200
\(A\)：5
\(\alpha\)：1

エネルギー固有値

量子数\(n\) 数値計算解 解析解
0 –3.5437173567132851 -3.5438611699158100
1 –1.3814252908344455 -1.3815835097474309
2 –0.2185797147047042 -0.2193058495790518
3 -0.0049226123802256 –
4 0.0081675032299395 –
…
10 0.0527783252864115 –
20 0.2666761210195411 –
30 0.6404251854103465 –
40 1.1703593933291645 –
50 1.8543552927405405 –
考察

厳密解は3つしか存在しません。この3つの状態に対しては数値計算解の方が厳密解よりも大きくなり、想定通りです。
しかし精度があまり良くありません。これはポテンシャルが\(V(x\to\infty)=0\)であり、波動関数が速やかにゼロに向かわないためだろうと考えられます。
また、基底状態でもあまり良い結果が得られていません。上記ポテンシャルの閉じ込めが強くなく、多くの基底関数が必要となり、数値計算で用いた200個では十分ではないのでしょう。

4つ目の解に注目しましょう。厳密解は存在しないのに数値計算ではあたかも解が存在するように見えます(\(E\le 0\))。また、変分原理は上限を与えるため、4つ目の解が存在すると考えてしまうかもしれません。しかし、これは違います。明らかにするために4つ目の解の波動関数の形(青色の線)を拡大してみますと、

となり、物理的に意味を成していないことが分かります。なぜなら、波動関数は\(x\to\infty\)に向かうにつれて徐々に小さくなっていかなければなりません。もしもこれが正しければ、\(x=80\)で節を持っている解ということになり、これは束縛状態の解としては適しません。

エネルギー固有値だけを考えては求められないものだったということですね。

さらに、束縛状態以上ではすべて連続状態です。それにもかかわらず\(E\gt 0\)でも離散的に固有値が得られているのは基底関数に固定端の境界条件を課しているためです。これらに物理的な意味は無いので評価する際には気を付けましょう。

\(1/(\cosh^2)\)型ポテンシャル[3]


\(
\displaystyle V(x)=-\frac{V_0}{\cosh^2 \alpha x}
\)

固有値

\(
\displaystyle E_n=-\frac{\alpha^2}{8}\left[-(1+2n)+\sqrt{1+\frac{8V_0}{\alpha^2}}\right]^2
\)
で、\(n\)は
\(
\displaystyle n\lt s = \frac{1}{2}\left(-1+\sqrt{1+\frac{8V_0}{\alpha^2}}\right)
\)
から有限個の準位が決まります。

固有関数

\(
\begin{align}
\psi(x)&=(1-\xi^2)^{\xi/2}F\left[\varepsilon-s, \varepsilon+s+1, \varepsilon+1, \frac{1-\xi}{2}\right]\\
& \varepsilon=\frac{\sqrt{-2E}}{\alpha}\\
& \xi=\tanh \alpha x
\end{align}
\)
ここでの\(F(a,b,c,x)\)は超幾何関数です。

解法

[3]を参照してください

数値計算解と解析解との比較

数値計算の設定:
Gauss-Lobatto求積法の次数(NGL):12
計算区間(a,b):[-30,30]
基底関数の数(N)：200

エネルギー固有値

量子数\(n\) 数値計算解 解析解
0 –3.8286053565505553 -3.885009803959261
1 –1.9750293441823026 -1.9750294118777854
2 –0.6459145040771800 -0.7050490197963089
3 –0.0750649726143545 -0.0750686277148326
4 0.0076704714372625 –
…
10 0.1139873400292487 –
20 0.4849841979196881 –
30 1.1689850012563472 –
40 2.1388129318397651 –
50 3.3888397726738235 –
考察

エネルギー固有値の上限を確かに与えています。そのほかはおおよそF.モースによるポテンシャルと同じような結果を与えています。

精度に関して、エネルギー固有値の精度は\(n=0,2\)の時、1~2桁程度の精度ですが、\(n=1,3\)の時、は6~7桁一致と著しく精度が上がっています。計算区間を変えたりしたのですが、この傾向は変わりません。
なぜなのか…詳しくは分かりません。偶関数、奇関数の問題だと思います。すなわち、奇関数の場合はなにか\(x\gt 0\)の領域に誤差があっても同じ量の誤差で\(x\lt 0\)で打ち消し合ってくれるのですが、偶関数の場合は打ち消し合いが起きないのではないか、ということです。
ただ、これが行列要素の精度からくるものなのかは詳しく見ていないので分かりません…

クーロンポテンシャル(動径方向, \(l=0\))

ポテンシャル

\(
\displaystyle V(x)=-\frac{1}{x}
\)
※ここでは、一階微分の存在しない形にした時のシュレーディンガー方程式を考えています。

固有値

\(
\displaystyle E_n=-\frac{1}{2n^2},\;\;(n=1,2,\cdots,)
\)
です。

固有関数

\(
\begin{align}
\psi(x)=C \cdot e^{-x/n} x^l L_{n-1}^{(1)}(x)
\end{align}
\)
ここで,\(C\)は規格化定数, \(L_n^k(x)\)はラゲール倍多項式[5]で
\(
\begin{align}
L_0^k(x)&=1 \\
L_1^k(x)&=-x+k+1 \\
L_2^k(x)&=\frac{1}{2}[x^2-2(k+2)x+(k+1)(k+2)] \\
L_3^k(x)&=\frac{1}{6}[-x^3+3(k+3)x^2-3(k+2)(k+3)x+(k+1)(k+2)(k+3)]
\end{align}
\)
です。もちろん規格化定数をあらわに決めることが出来ますが、通常、この問題は3次元のシュレーディンガー方程式を解いた時の動径方向として出てくるので、その時の規格化は
\(
\displaystyle \int_0^\infty x^2 \psi(x)^2 dx =1
\)
で規格化されます。しかし、本稿の規格化は\(\displaystyle \int_0^\infty \psi(x)^2 dx =1\)
で規格化しているので、エネルギー固有値のみの比較を行います。

解法

\(x=0, x\to infty\)の漸近形を考えて、本当の解を
(\(x=0\)の漸近形)×(\(x=\infty\)の漸近形)×(未知関数)
と仮定します。これをシュレーディンガー方程式に代入すると未知関数がラゲール陪多項式だと分かります。

数値計算解と解析解との比較

数値計算の設定:
Gauss-Lobatto求積法の次数(NGL):12
計算区間(a,b):[0,80]
基底関数の数(N)：200

エネルギー固有値

量子数\(n\) 数値計算解 解析解
1 –0.4970775696080274 -0.5000000000000000
2 -0.1246297940063756 -0.1250000000000000
3 –0.0554455545682885 -0.0555555555555556
4 –0.0312035408553899 -0.0312500000000000
5 –0.0199704293764349 -0.0200000000000000
6 –0.0134298208377274 -0.0138888888888889
7 –0.0068151982214300 -0.0102040816326531
8 0.0018616612074302 -0.0078125000000000
…
10 0.0251247065596211 -0.0050000000000000
20 0.2439421094228725 -0.0012500000000000
考察

クーロンポテンシャルの場合、\(x=\infty\)では波動関数は指数関数で減衰していきます。
このため、減衰がゆっくりになってしまい、sin関数でうまく表現できないために精度が落ちているのだと推測できます。
また、計算区間を十分に大きく取ると励起状態の表現はうまくいきますが、反対に基底状態の表現ができなくなります。なぜなら、多くの区間で波動関数はゼロですが、\(x=0\)の近傍だけに大きなピークを持つため、sin関数ではうまく表現できません。その兆候は上の数値計算でも見えてます。基底状態は2桁程度ですが、量子数の増加と共に1桁くらい精度が上がっています。
クーロンの場合、ラゲール関数を基底関数としてとる方が良いでしょう。

参考文献

[1]P. J. Davis, P. Rabinowitz著, 森正武訳, 「計算機による数値積分法」, Gauss-Lobatto求積法 p.88
[2]P. J. Davis, P. Rabinowitz著, 森正武訳, 「計算機による数値積分法」, ゼロ点に挟まれる区間の積分 p.131
[3]ランダウ＝リフシッツ,「量子力学１」第3刷, p81-84
[4]Abramowitz and Stegun, HANDBOOK OF MATHEMATICAL FUNCTIONS, p.446,478 http://people.math.sfu.ca/~cbm/aands/page_446.htm

Airy関数のゼロ点は
エアリー関数 Ai, Bi（ゼロ点）-Keisan
より22桁の精度で求めたものを用いています。
ちなみにそれらの値(Airy関数のゼロ点)は、

0 -2.338107410459767038489
1 -4.087949444130970616637
2 -5.52055982809555105913
3 -6.78670809007175899878
4 -7.944133587120853123138

10 -13.6914890352107179283
20 -21.2248299436420969552
30 -27.58838780988244481195
40 -33.28488468190140187962
50 -38.52880830509424882263

です。ゼロ点を0,1,2,…と、今回の計算の都合上0からインデックスを付けています。
もしかすると最後の1~2桁はあってないかもしれませんが、本稿では16桁まであっていればいいので、問題ありません。

[5]
Associated Laguerre Polynomial -wolfram mathworld式(22)-(25)

プログラミングと数値計算, 物理学

時間依存しないシュレーディンガー方程式と変分原理 1/2 （計算手法の説明）

コメントする

時間依存しないシュレーディンガー方程式を数値的に解くための、変分原理について述べていきます。
変分原理はエネルギー固有値の上限を与えます(厳密な解よりも変分原理で求めた解は必ず大きくなる)。
基底関数を適切に選べば高精度の結果が与えられます。ですが、基底関数をどのように選べばよいか？の指標は無く、経験則に頼らざるを得ない状況が出てきます。

ここでは、変分原理によってなぜエネルギー固有値が求まるのかを述べて行きます。

もしも変分原理より正しい解を得たければ、shooting methodと呼ばれる方法も推奨します。
恐らくこちらの方が色々なポテンシャルに対して有効ですが、
・計算速度が遅い
・発散の可能性がある
・変分原理よりも経験値が必要
です。

まとめ

区間[a,b]で定義された時間依存しないシュレーディンガー方程式
\(
\displaystyle \left[-\frac{1}{2}\frac{d^2}{dx^2}+V(x)\right]\psi(x)=E\psi(x)
\)
の解である固有関数を
\(
\displaystyle \psi(x)=\sum_n c_n \varphi_n(x) \ \ …(6)
\)
の形で固定端条件(\(\psi(a)=\psi(b)=0\))の下で探します。
ここで、\(\varphi_n(x)\)は区間\([a,b]\)で正規直交化されたsin関数系
\(
\displaystyle \varphi_n(x)=\sqrt{\frac{2}{b-a}}\sin\left(\frac{n\pi}{b-a}(x-a)\right),\;\;\;(n=1,2,\cdots,N)
\)
です。この時、行列要素\(H_{ij}=\langle\varphi_i|-\frac{1}{2}\frac{d^2}{dx^2}+V(x)|\varphi_j\rangle\)は
\(
\begin{align}
H_{ij}&=K_{ij}+V_{ij} \\
&\hspace{1em}K_{ij}=\frac{1}{2}\left(\frac{j\pi}{b-a}\right)^2\delta_{ij} \\
&\hspace{1em}V_{ij}=\frac{1}{j}\sum_{m=2}^{M-1}\omega_m^{(GL)} \cos\left(\frac{\pi}{2}x_m^{(GL)}\right)
\sum_{n=0}^{j-1}(-1)^n \sin\left(\frac{i\pi}{2j}(x_m^{GL}+2n+1)\right)V(X_n^{(m)}) \\
&\hspace{2em}X_{n}^{(m)}=\frac{h}{2}(x_m^{(GL)}+2n+1)+a,\;\;h=\frac{b-a}{j}
\end{align}
\)
と求められます(\(j\ge i\)を推奨)。ここで、\(\omega_m^{(GL)},\ x_m^{(GL)}\)はM点のGauss-Lobatto求積法の区間[-1,1]における重みと求積点で、
\(
\displaystyle \omega_m=\frac{2}{M(M-1)[P_{M-1}(x_m)]^2},\ \ x_m\ne \pm 1
\)
であり、求積点は関数
\(
\displaystyle \frac{P_{M-1}(x)+P_{M}(x)}{x-1}
\)
のゼロ点です。\(P_{n}(x)\)は\(n\)次Legendre多項式を表します。

  1. 解説
  2. 時間依存から時間依存しないシュレーディンガー方程式へ
  3. 間接法
    1. 行列要素の計算
    2. K_{nm}の計算
    3. V_{nm}の計算
    4. Gauss-Lobatto求積法による積分の計算
    5. ゼロ点に挟まれる区間の積分
  4. 数値計算プログラム
  5. 正規直交系を用いない場合
  6. 研究分野での数値計算法について
  7. 参考文献

計算の具体例は
時間依存しないシュレーディンガー方程式と変分原理 2/2 （計算例）をご覧ください。

[adsense1]

解説

対象とする問題を導出します。何はともあれ、時間依存シュレーディンガー方程式から出発します。

時間依存から時間依存しないシュレーディンガー方程式へ

時間依存シュレーディンガー方程式は
\(
\displaystyle i\frac{\partial \Psi(x,t)}{\partial t}=\hat{H}(x)\Psi(x,t)
\)
と記述されます。ここで、簡単ために原子単位系(\(\hbar=1, m=1\))を考えます。
変数分離の形で解を仮定します。
\(
\displaystyle \Psi(x,t)=\psi(x)U(t)
\)
代入して移行すると以下の式を得ます。

\(
\begin{align}
i\psi(x)\frac{\partial U(t)}{\partial t}=U(t)\hat{H}\psi(x) \\
i\frac{1}{U(t)}\frac{\partial U(t)}{\partial t}=\frac{1}{\psi(x)}\hat{H}\psi(x) = E
\end{align}
\)

2番目の式の右辺\(E\)は何らかの定数です。
左辺は時間のみを含む式、真ん中の式は位置のみを含む式なので、いつどんな時でもこの関係式が成り立っていないといけないという条件からこれらはなんらかの定数でなければなりません。その定数を\(E\)と置きました。

すなわち、２つの式
\(
\begin{align}
i\frac{\partial U(t)}{\partial t}=E U(t)\;\;\; \cdots(1)\\
\hat{H}\psi(x) = E\psi(x)\;\;\; \cdots(2)
\end{align}
\)
が導けます。

\(t\)に関する式は\(\hat{H}(x)\)に依らず解けてしまい、
\(U(t)=e^{-i E t}\)(1)
です。位置に関する式は
\(
\hat{H}\psi(x) = E\psi(x)
\)(2)
を満たさなければなりません。この(2)が時間依存しないシュレーディンガー方程式と呼ばれているものです。
これを数値計算で解いていくことを考えます。

左から\(\psi^*\)を掛けて全空間で積分しますと、
\(
\begin{align}
\int \psi^* \hat{H} \psi dx – E \int \psi^* \psi dx = 0\;\;\; \cdots(3)
\end{align}
\)
という式を得ます。この問題は
条件\(\int \psi^* \psi dx=1\)の下で、\(\int \psi^* \hat{H} \psi dx\)の極値を求めよ、
という問いになります。
上記問題を考えた時、(3)と同一の式になります。

すなわち、原理的にはどんな\(\psi(x)\)を与えても\(E(=E[\psi])\)の値を計算できるわけです。
その無数に考えられる\(\psi\)の中で、特に(3)の解が停留する\(\psi\)、すなわち
\(
\displaystyle \delta \int \psi^* (\hat{H}-E) \psi dx = 0
\)
となる\(\psi\)を探せ、という問題になります。

↑がもしも停留する\(\psi\)が見付かれば、それはシュレーディンガー方程式と同一になるということは以下の通りに示せます。

\(
\displaystyle \delta \int \psi^* (\hat{H}-E) \psi dx =0 \;\;\;\cdots (4)
\)
厳密な\(\psi\)の変分に対して停留、すなわち\(\psi\to \psi+\delta \psi\)、もしくは複素平面で\(\frac{\pi}{2}\)ずれた方向\(\psi\to \psi+i\delta \psi\)に対して0でなければなりません。
(4)に対して２つの方向に対して変分を取ると、
\(
\begin{align}
\int (\delta \psi)^* (\hat{H}-E)\psi dx &+ \int \psi^* (\hat{H}-E)(\delta\psi) dx =0 \\
-i\int (\delta \psi)^* (\hat{H}-E)\psi dx &+ i\int \psi^* (\hat{H}-E)(\delta\psi) dx =0
\end{align}
\)
と2式が導けます。2式は変分をとった方向にかかわらず同時に満たされていなければなりません。同時に成立するためには
\(
(\hat{H}-E)\psi=0,\;\; \text{or}\;\; \psi^*(\hat{H}^{\dagger}-E^*)=0
\)
を満たしている必要があります。これは時間依存しないシュレーディンガー方程式と一致します。

以上をまとめますと、

厳密な\(\psi\)の変分はシュレーディンガー方程式に一致する
↓
近似的に\(\psi\)を変分するとシュレーディンガー方程式の近似的な解に帰着されるだろう
↓
近似的に\(\psi\)を構成し、変分問題を解く

代表的な方法として、
手法1)直接法…\(\psi\)の具体的な形を変分パラメータと共に与える
手法2)間接法…変分パラメータとしてある正規直交基底関数系を採用する
の２つがありますが、数値計算で使われる間接法を説明します。

間接法

間接法による変分原理は、

  • 固有値の上限を与える（変分原理によるエネルギー固有値は必ず厳密なエネルギー固有値よりも大きくなる）
  • 基底関数の数を増やすほど精度は上がる(完全系の場合)

という特性があります。

さて、今解きたい問題は
条件
\(
\displaystyle \int \psi^* \psi dx =1
\)
の下で
\(
\displaystyle \int \psi^* \hat{H} \psi dx
\)
の極値を与える\(\psi\)を探すこと、です。

すなわち
\(
\displaystyle \int \psi^* \hat{H} \psi dx – E \int \psi^* \psi dx =0 \ \ …(5)
\)
を考えて変分します。求めたい波動関数は
\(
\displaystyle \psi(x)=\sum_n c_n \varphi_n(x) \ \ …(6)
\)
と展開できる、と考えます。ここで、\(c_n\)は変分パラメータ、\(\varphi_n(x)\)は区間[a,b]で定義されている正規直交関数系です。正規直交関数系\(\{\varphi_n(x)\}\)は以下の性質を持ちます。
\(
\displaystyle \int_a^b \varphi_n(x) \varphi_m(x)dx =\delta_{nm}
\)
この制約が無くても変分は可能ですが、以降の式が簡略になるので、正規直交関数だけを考えます。

(6)を(5)に代入すると、
\(
\begin{align}
\sum_{n}\sum_{m} c_n^* c_m \left[H_{nm}-E \delta_{nm}\right]=0 \\
\end{align}\ \ …(7)
\)
なります。ここで、
\(
\begin{align}
H_{nm}&=\int_a^b \varphi_n(x) \hat{H} \varphi_m(x)dx
\end{align}
\)
とおきました。(7)式を満たす係数\(\{c_n\}\)の内、自明な解(\(c_n=0\ (for\ all\ n)\))はすぐ見つかります。しかし、この解は\(\psi(x)=0\)であり、条件\(
\displaystyle \int \psi^* \psi dx =1\)を満たしません。なので解としては不適切です。

では非自明な解を見つけましょう。(7)式の解を探すため、極値を探します。今、係数の組を変分パラメータとみなしているので、(変分パラメータで微分したもの)=0を計算します。すなわち
\(
\displaystyle \frac{\partial}{\partial c_n^*} (for\ all\ n)
\)
を(7)式に作用させます。すると、
\(
\begin{align}
\sum_{n}\sum_{m} \left[H_{nm}-E\delta_{nm}\right]c_m=0
\end{align}
\)
となります。\(n,m\)を行列の行、列ととらえれば、

\(
\begin{align}
\left(
\begin{array}{cccc}
H_{11} & H_{12} & \cdots & H_{1N}\\
H_{21} & H_{22} & \cdots & H_{2N}\\
\ddots & \ddots & \ldots & \ddots \\
H_{N1} & H_{N2} & \cdots & H_{NN}
\end{array}
\right)
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
\ddots \\
c_{N}
\end{array}
\right)
=
E
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
\ddots \\
c_{N}
\end{array}
\right)
\end{align}
\)

となるわけです。ということは、左辺の正方行列を対角化すれば右辺のE,すなわち固有値が求まる、という事になります。

行列要素の計算

行列を対角化するためには行列要素\(H_{nm}\)が分かっていないとなりません。
これからこの行列要素を計算していきます。

ハミルトニアンが
\(
\displaystyle \hat{H}=-\frac{1}{2}\frac{d^2 }{d x^2} + V(x)
\)
として表現されている場合、行列要素は
\(
\begin{align}
H_{nm}&=K_{nm}+V_{nm} \\
&\hspace{3em} K_{nm}=-\frac{1}{2}\int_a^b \varphi_n(x)\frac{d^2}{dx^2} \varphi_m(x)dx \\
&\hspace{3em} V_{nm}=\int_a^b \varphi_n(x)V(x) \varphi_m(x)dx
\end{align}
\)
と書き表され、これらを数値計算で求めます。

\(K_{nm}\)の計算

まず、微分を含んだ
\(
\displaystyle K_{nm}=-\frac{1}{2}\int_a^b \varphi_n(x)\frac{d^2}{dx^2} \varphi_m(x)dx
\)
の項を考えます。微分は数値計算で行いません。
というのも正規直交関数系はこちらで用意する知っている関数系なので、求めることが出来てしまうからです。
多くの正規直交関数系には漸化式と微分方程式が存在します。一番簡単な基底としてsin関数系を考えてみましょう。
区間\([a,b]\)で正規直交化されたsin関数系は
\(
\displaystyle \varphi_n(x)=\sqrt{\frac{2}{b-a}}\sin\left(\frac{n\pi}{b-a}(x-a)\right),\;\;\;(n=1,2,\cdots,N)
\)
と書き表せます。sin関数なので、2階微分はもちろん
\(
\displaystyle \frac{d^2}{dx^2}\varphi_n(x)=-\left(\frac{n\pi}{b-a}\right)^2\varphi_n(x)
\)
を満たします。これを積分の中に代入して、
\(
\begin{align}
K_{nm}&=-\frac{1}{2}\int_a^b \varphi_n(x)\frac{d^2}{dx^2} \varphi_m(x)dx \\
&=\frac{1}{2}\left(\frac{n\pi}{b-a}\right)^2\int_a^b \varphi_n(x) \varphi_m(x)dx \\
&=\frac{1}{2}\left(\frac{n\pi}{b-a}\right)^2\delta_nm
\end{align}
\)
となります。積分を計算することなく求める事が出来ました。

より一般には一回微分の存在しない形、もしくは漸化式によって手計算で微分を含む積分を片付けてしまったり、部分積分を用いて何とか計算します。
正規直交系であれば、うまく整理することによって誤差が全く入らず、かつ高速に計算できます。

\(V_{nm}\)の計算

さて、残る問題は
\(
\displaystyle V_{nm}=\int_a^b \varphi_n(x)V(x) \varphi_m(x)dx
\)
の計算です。基底関数の次数が低い場合、この計算は何も問題がありません。
調和ポテンシャル\(V(x)=x^2\)の時、\(n=1,m=2\)の時の被積分関数は

と非常になだらかですが、\(n=5,m=99\)の時の被積分関数は

となります。なので細かい分点をとることになり、非常に積分が難しいです。

離散フーリエ変換によって求められそうな気もしますが、周期境界条件が満たされていないことは明らかなのであきらめましょう。
有限区間で早く振動する関数の積分は難しいです。二重指数関数型数値積分も振動する関数は苦手です。

こういった場合の計算方法として、
・ゼロ点の位置が分かっている場合、そこで区間を分割して積分する
という直観的な方法があります[2]。
ゼロ点の位置が分かっていれば、その間の被積分関数の振る舞いはシンプルな形になるので滑らかな関数として考えられます。

また、端点での値がゼロだと分かっている場合、Gauss-Lobatto求積法が非常に有用です。

Gauss-Lobatto求積法による積分の計算

Gauss-Lobatto求積法はGauss-Legendre求積法と同じく被積分関数を多項式で近似します[1]。
Gauss-Legendre求積法と異なる点は端点で分点を持つということです。端点を明示的に書けば、積分は
\(
\displaystyle \int_{-1}^1 f(x)dx\sim \frac{2}{N(N-1)}(f(-1)+f(1))+\sum_{i=2}^{N-1}\omega_i f(x_i)
\)
と書けます。ここで、\(\omega_i,x_i\)はGauss-Lobatto求積法における重み、求積点です。
Gauss-Lobatto求積法の重みは
\(
\displaystyle \omega_i=\frac{2}{N(N-1)[P_{N-1}(x_i)]^2},\ \ x_i\ne \pm 1
\)
であり、求積点は
\(
\displaystyle \frac{P_{N-1}(x)+P_{N}(x)}{x-1}
\)
のゼロ点です。\(P_{n}(x)\)は\(n\)次Legendre多項式を表します。

Gauss-Lobatto求積法は
\(f(\pm 1)\)が既知の値(特に0)を持つ場合
に便利なことが分かります。
関数\(f(\pm 1)\)がゼロであれば、その点を計算する必要性がなくなりますし、見かけ上の特異点であってもこちらで値を指定すれば参照されることはなく、計算できないという状況に陥りません。

ゼロ点に挟まれる区間の積分

ゼロ点に挟まれる区間の積分は、結局以下の式になります[2]。
\(
\displaystyle \varphi_j(x)=\sqrt{\frac{2}{b-a}}\sin\left(\frac{j\pi}{b-a}\right)
\)
とすると、積分は
\(
\begin{align}
&\int_a^b \varphi_i(x) V(x) \varphi_j(x) dx \\
&=\sum_{n=0}^{j-1} \int_{x_n}^{x_{n+1}} \varphi_i(x) V(x) \varphi_j(x)dx \\
&\approx \sum_{n=0}^{j-1}\sum_{m=2}^{M-1}\frac{h}{2} \omega_m^{(GL)}\varphi_i(X_n^{(m)}) V(X_n^{(m)}) \varphi_j(X_n^{(m)}) \\
&=\frac{h}{2}\sum_{m=2}^{M-1}\omega_m^{(GL)}\varphi_j(X_0^{(m)})S_j^{(m)} \\
\end{align}
\)
ここで、
\(
\begin{align}
&x_n=nh+a,\ h=\frac{b-a}{j},\ (n=0,1,\cdots,j) \\
&X_n^{(m)}=\frac{h}{2}x_m^{(GL)}+\frac{x_{n+1}+x_n}{2} \\
&X_{n+1}^{(m)}=X_n^{(m)}+h \\
&S_j^{(m)}=\sum_{n=0}^{j-1}(-1)^n \varphi_i(X_n^{(m)})V(X_n^{(m)})
\end{align}
\)
と置きました。
\(x_m^{(GL)},\ \omega_m^{(GL)}\)は区間\([-1,1]\)でのM点Gauss-Lobatto求積法の求積点と重みです。
計算を可能な限り減らすため、\(\varphi_j(x)\)の周期性を用いています。

数値計算プログラム

行列要素を組み立てた後は適当なルーチンを用いて対角化します。
今回の場合、基底関数をsin関数に選んだため、強制的に固定端条件が課されます。
なので、この基底関数系を用いるといかなる問題も端点\(a,b\)に無限に高いポテンシャルの壁がある、と読み替えることが出来ます。
よって行列は対称行列、エルミートです。故に固有値は必ず実数です。

この問題を解く場合、実対称行列を対角化するルーチン、lapackのdsyevdを用いるのが良いでしょう。
対角化のプログラムに関する詳細は対角化へどうぞ。

以下のプログラムは行列要素を上記手法で作り、対角化し、得られた波動関数を規格化し、出力するものです。

[adsense2]

コンパイルと実行は

$ ifort -mkl main.f90
$ ./a.out

などとしてください。lapackを用いますので、それと一緒にコンパイルしてください。

ラグランジュの未定乗数について

ラグランジュの未定乗数は特に物理的な意味を持たない魔法の変数、として置きますが、今回の場合は未定乗数がエネルギー固有値というちゃんとした意味を持つのです。

正規直交系を用いない場合

\(
\begin{align}
\sum_{n}\sum_{m} c_n^* c_m \left[H_{nm}-E I_{nm}\right]=0 \\
\sum_{n}\sum_{m} c_n^* c_m \left[J_{nm} H_{nm}-E \right]=0
\end{align}\ \ …(7)
\)
となります。ここで、
\(
\begin{align}
H_{nm}&=\int_a^b \varphi_n(x) \hat{H} \varphi_m(x)dx\\
J_{nm}&=I_{nm}^{-1}\\
I_{nm}&=\int_a^b \varphi_n(x) \varphi_m(x)dx
\end{align}
\)
と置きました。行列の形では、
\(
\begin{align}
\left(
\begin{array}{cccc}
J_{11}H_{11} & J_{12}H_{12} & \cdots & J_{1N}H_{1N}\\
J_{21}H_{21} & J_{22}H_{22} & \cdots & J_{2N}H_{2N}\\
\ddots & \ddots & \ldots & \ddots \\
J_{N1}H_{N1} & J_{N2}H_{N2} & \cdots & J_{NN}H_{NN}
\end{array}
\right)
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
\ddots \\
c_{N}
\end{array}
\right)
=
E
\left(
\begin{array}{cccc}
1 & 0 & \cdots & 0\\
0 & 1 & \cdots & 0\\
\ddots & \ddots & \ldots & \ddots \\
0 & 0 & \cdots & 1
\end{array}
\right)
\left(
\begin{array}{c}
c_{1} \\
c_{2} \\
\ddots \\
c_{N}
\end{array}
\right)
\end{align}
\)
と書けます。

\(
\begin{align}
\sum_{n}\sum_{m} c_n^* c_m \left[H_{nm}-E I_{nm}\right]=0 \\
\sum_{n}\sum_{m} c_n^* c_m \left[I_{nm}^{-1} H_{nm}-E \right]=0
\end{align}
\)

研究分野での数値計算法について

実際の研究の場では対角化はしますが、上記方法によって行列要素を作ることはしません。遅いからです。
特別な基底関数を用いると行列要素\(H_{nm}\)が1~2項、しかも非常に簡単な式に置き換えることが出来ます。精度も速度も上記方法以上です。
この計算方法に関してはこれ以上は述べることはしませんが、こういった方法も存在する、とだけ述べておきます。

参考文献

[1]P. J. Davis, P. Rabinowitz著, 森正武訳, 「計算機による数値積分法」, Gauss-Lobatto求積法 p.88
[2]P. J. Davis, P. Rabinowitz著, 森正武訳, 「計算機による数値積分法」, ゼロ点に挟まれる区間の積分 p.131
[3]ランダウ＝リフシッツ,「量子力学１」第3刷, p81-84

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

ガンマ関数の数値計算

コメントする

複素平面で有効なガンマ関数を数値計算します。
目指すは大浦様が公開してくださっている複素ガンマ関数cdgamma,cqgammaの計算速度を上回ることです。

結果として、本稿の複素ガンマ関数の計算コードは大浦様のコードの0.85~0.9倍ほどの時間で計算できるコードになりました。

大浦様のコードは、Lanczos近似という方法に基づいています。
Lanczos近似はガンマ関数を

\(
\begin{align}
\Gamma(z)=\frac{(z+g-0.5)^{z-0.5}}{e^{z+g-0.5}}L_g(z) \\
L_g(z)=C_0+\sum_{k=1}^{N-1}\frac{C_N}{z+k-1}
\end{align}
\)
の形で展開します[6]。係数\(C_k\)は定数\(g\)の値に依存します。
定数\(g\)は魔法の数のようなもので、値によって収束具合が異なります。このあたりの理論に関しては[5-6]が詳しいのでそちらをご参考ください。
おおよそ10-20項程度で収束し、非常にシンプルなコードになっているのが特徴です。

同じ考えでは超えられるはずがないので、ガンマ関数の特性をフルに使って計算していきます。

実数のみのガンマ関数が欲しい場合、gfortranでは組み込み関数dgamma, gammaを最適化オプション 　-Oまたは-O3　とともにコンパイルすると500分の1くらいの時間で計算できるため、圧倒的に早いです。

[adsense1]

方針

方針としては

  • 漸化式
  • 相反公式
  • 逆ガンマ関数のTaylor展開
  • 対数ガンマ関数の漸近展開

を用いて計算していきます。

詳細は以下に折りたたんでおきます。

領域ごとに使っている関係式をまとめると、このような感じです。

※範囲変更(cgamma,円の大きさ8→9)

考え方は以上です。
実際にFortran90でプログラムしてみます。

fortran90によるプログラム

用意してあるのは、

複素倍精度用のガンマ関数cgamma
(本当のガンマ関数値との相対誤差 < 2e-14 [-10,10]の範囲, 8e-14 [-40,40]の範囲)
と

4倍精度複素ガンマ関数zgamma
(本当のガンマ関数値との相対誤差 < 3e-32 [-10,10]の範囲, 8e-32 [-40,40]の範囲)
です。

2017/01/10 逆ガンマ関数の項数の最適化
2017/10/27 正の整数の時、1ずれることの修正

function cgamma(z)
  ! sikinote (http://slpr.sakura.ne.jp/qp)
  ! author : sikino
  ! date   : 2017/01/09
  !          2017/01/10 optimize the number of a(M)
  !          2017/10/27 fix for positive integer
 
  implicit none
  complex(kind(0d0)),intent(in)::z
  complex(kind(0d0))::cgamma
 
  integer::i,n,M
  double precision,parameter::pi=3.14159265358979324d0
  double precision,parameter::e1=0.3678794411714423216d0
  double precision,parameter::ln2pi2=0.91893853320467274d0
  double precision,parameter::sq2pi=2.5066282746310005d0
  double precision::a(1:30),AB(1:14),dz,iz
  complex(kind(0d0))::q,s,w,r,q1,q2

  dz=dble(z)
  iz=imag(z)

  if(dz.eq.nint(dz).and.iz.eq.0d0)then
     if(dz.gt.0d0)then
        s=dcmplx(1d0,0d0)
        n=nint(dz)
        do i=2,n-1
           s=s*i
        enddo
     else
        s=1d+300
     endif
     cgamma=s
  else
     q=z
     if(dz.lt.0d0)q=1d0-z

     if(abs(iz).lt.1.45d0)then
        ! For x=arb, |y|\lt 1.45
        n=int(dble(q))
        s=1d0
        do i=1,n
           s=s*(q-i)
        enddo
        w=q-dble(n)

        a(1) = 1d0
        a(2) = 0.57721566490153286d0
        a(3) =-0.65587807152025388d0
        a(4) =-0.42002635034095236d-1
        a(5) = 0.16653861138229149d0
        a(6) =-0.42197734555544337d-1
        a(7) =-0.96219715278769736d-2
        a(8) = 0.72189432466630995d-2
        a(9) =-0.11651675918590651d-2
        a(10)=-0.21524167411495097d-3
        a(11)= 0.12805028238811619d-3
        a(12)=-0.20134854780788239d-4
        a(13)=-0.12504934821426707d-5
        a(14)= 0.11330272319816959d-5
        a(15)=-0.20563384169776071d-6
        a(16)= 0.61160951044814158d-8
        a(17)= 0.50020076444692229d-8
        a(18)=-0.11812745704870201d-8
        a(19)= 0.10434267116911005d-9
        a(20)= 0.77822634399050713d-11
        a(21)=-0.36968056186422057d-11
        a(22)= 0.51003702874544760d-12
        a(23)=-0.20583260535665068d-13
        a(24)=-0.53481225394230180d-14
        a(25)= 0.12267786282382608d-14
        a(26)=-0.11812593016974588d-15
        a(27)= 0.11866922547516003d-17
        a(28)= 0.14123806553180318d-17
        a(29)=-0.22987456844353702d-18
        a(30)= 0.17144063219273374d-19

        M=int(11.3*abs(w)+13)
        if(M.gt.30)M=30
        r=a(M)
        do i=M-1,1,-1
           r=r*w+a(i)
        enddo
       
        cgamma=s/(r*w)
     else
        ! For x=arb, |y|\gt 1.45
        s=1d0
        if(abs(q).lt.9d0)then
           do i=0,7
              s=s*(q+i)
           enddo
           s=1d0/s
           q=q+8d0
        endif

        AB(1) = 0.83333333333333333d-1
        AB(2) =-0.27777777777777778d-2
        AB(3) = 0.79365079365079365d-3
        AB(4) =-0.59523809523809524d-3
        AB(5) = 0.84175084175084175d-3
        AB(6) =-0.19175269175269175d-2
        AB(7) = 0.64102564102564103d-2
        AB(8) =-0.29550653594771242d-1

        q1=1d0/q
        q2=q1*q1
       
        r=AB(8)
        do i=7,1,-1
           r=r*q2+AB(i)
        enddo
        cgamma=s*exp((q-0.5d0)*log(q)-q+ln2pi2+r*q1)
     endif
     if(dz.lt.0d0)cgamma=pi/(cgamma*sin(pi*z))
  endif

  return
end function cgamma

function zgamma(z)
  ! sikinote (http://slpr.sakura.ne.jp/qp)
  ! author : sikino
  ! date   : 2017/01/09
  !          2017/01/10
  !          2017/10/27 fix for positive integer
  implicit none
  complex*32,intent(in)::z
  complex*32::zgamma
 
  integer::i,n,M
  real*16,parameter::    pi=3.14159265358979323846264338327950288q0
  real*16,parameter::    e1=0.36787944117144232159552377016146087q0
  real*16,parameter::ln2pi2=0.91893853320467274178032973640561764q0
  real*16,parameter:: sq2pi=2.50662827463100050241576528481104525q0
  real*16::a(1:50),AB(1:30),dz,iz,aw
  complex*32::q,s,w,r,q1,q2

  dz=real(z); iz=imag(z)
  if(dz.eq.nint(dz).and.iz.eq.0q0)then
     if(dz.gt.0q0)then
        s=dcmplx(1q0,0q0)
        n=nint(dz)
        do i=2,n-1
           s=s*i
        enddo        
     else
        s=1q+3000
     endif
     zgamma=s
  else
     q=z
     if(dz.lt.0q0)q=1q0-z  

     if(abs(imag(q)).lt.1.35q0)then
        ! For x=arb, |y|\lt 1.35
        n=int(real(q))
        s=1q0
        do i=1,n
           s=s*(q-i)
        enddo
        w=q-real(n)
       
        a(1) = 1q0
        a(2) = 0.577215664901532860606512090082402431q0
        a(3) =-0.655878071520253881077019515145390481q0
        a(4) =-0.420026350340952355290039348754298187q-1
        a(5) = 0.166538611382291489501700795102105236q0
        a(6) =-0.421977345555443367482083012891873913q-1
        a(7) =-0.962197152787697356211492167234819898q-2
        a(8) = 0.721894324666309954239501034044657271q-2
        a(9) =-0.116516759185906511211397108401838867q-2
        a(10)=-0.215241674114950972815729963053647806q-3
        a(11)= 0.128050282388116186153198626328164323q-3
        a(12)=-0.201348547807882386556893914210218184q-4
        a(13)=-0.125049348214267065734535947383309224q-5
        a(14)= 0.113302723198169588237412962033074494q-5
        a(15)=-2.056338416977607103450154130020572837q-7
        a(16)= 6.116095104481415817862498682855342867q-9
        a(17)= 5.002007644469222930055665048059991303q-9
        a(18)=-1.181274570487020144588126565436505578q-9
        a(19)= 1.043426711691100510491540332312250191q-10
        a(20)= 7.782263439905071254049937311360777226q-12
        a(21)=-3.696805618642205708187815878085766237q-12
        a(22)= 5.100370287454475979015481322863231803q-13
        a(23)=-2.058326053566506783222429544855237420q-14
        a(24)=-5.348122539423017982370017318727939949q-15
        a(25)= 1.226778628238260790158893846622422428q-15
        a(26)=-1.181259301697458769513764586842297831q-16
        a(27)= 1.186692254751600332579777242928674071q-18
        a(28)= 1.412380655318031781555803947566709037q-18
        a(29)=-2.298745684435370206592478580633699260q-19
        a(30)= 1.714406321927337433383963370267257067q-20
        a(31)= 1.337351730493693114864781395122268023q-22
        a(32)=-2.054233551766672789325025351355733797q-22
        a(33)= 2.736030048607999844831509904330982015q-23
        a(34)=-1.732356445910516639057428451564779799q-24
        a(35)=-2.360619024499287287343450735427531031q-26
        a(36)= 1.864982941717294430718413161878666894q-26
        a(37)=-2.218095624207197204399716913626860384q-27
        a(38)= 1.297781974947993668824414486330595180q-28
        a(39)= 1.180697474966528406222745415509988346q-30
        a(40)=-1.124584349277088090293654674261653627q-30
        a(41)= 1.277085175140866203990206677750563204q-31
        a(42)=-7.391451169615140823461289330001078362q-33
        a(43)= 1.134750257554215760954165252401478595q-35
        a(44)= 4.639134641058722029944804894433875948q-35
        a(45)=-5.347336818439198875077418068570325956q-36
        a(46)= 3.207995923613352622861238213439296964q-37
        a(47)=-4.445829736550756882101576639628564322q-39
        a(48)=-1.311174518881988712901096702293607192q-39
        a(49)= 1.647033352543813886818056435068597372q-40
        a(50)=-1.056233178503581218582902184722885523q-41

        aw=abs(w)
        if(aw.le.0.42q0)then
           M=int(35q0*aw+20q0)+1
        else
           M=int(15.5q0*aw+28.2q0)+1
        endif
        if(M.gt.50)M=50
        r=a(M)
        do i=M-1,1,-1
           r=r*w+a(i)
        enddo

        zgamma=s/(r*w)
     else
        ! For x=arb, |y|\ge 1.35
        s=1q0
        if(abs(q).lt.18q0)then
           do i=0,16
              s=s*(q+i)
           enddo
           s=1q0/s
           q=q+17q0
        endif

        AB(1) = 0.83333333333333333333333333333333333q-1
        AB(2) =-0.27777777777777777777777777777777778q-2
        AB(3) = 0.79365079365079365079365079365079365q-3
        AB(4) =-0.59523809523809523809523809523809524q-3
        AB(5) = 0.84175084175084175084175084175084175q-3
        AB(6) =-0.19175269175269175269175269175269175q-2
        AB(7) = 0.64102564102564102564102564102564103q-2
        AB(8) =-0.29550653594771241830065359477124183q-1
        AB(9) = 0.17964437236883057316493849001588940q0
        AB(10)=-1.3924322169059011164274322169059011q0
        AB(11)= 13.402864044168391994478951000690131q0
        AB(12)=-156.84828462600201730636513245208897q0
        AB(13)= 2193.1033333333333333333333333333333q0
        AB(14)=-36108.771253724989357173265219242231q0
        AB(15)= 691472.268851313067108395250775673468q0
        AB(16)=-1.52382215394074161922833649588867805q7

        q1=1q0/q; q2=q1*q1
       
        r=AB(16)
        do i=15,1,-1
           r=r*q2+AB(i)
        enddo

        zgamma=s*exp((q-0.5q0)*log(q)-q+ln2pi2+r*q1)  
     endif
     if(dz.lt.0q0)zgamma=pi/(zgamma*sin(pi*z))
  end if

  return
end function zgamma

[adsense2]

計算速度の比較

大浦様のコードと計算速度を比較します。
比較方法は、
乱数を[-20,20]の範囲で40000点振り、それを10回繰り返した時の平均cpu時間を計測します。
上の計算を20回繰り返し、その平均を取ることによって計算時間を比較します。
その結果は以下のようになりました。
大浦様のコードによる結果は黒、私のコードによる結果は赤です。

大浦様のコードと比べて、私のプログラムの方が計算時間が0.85~0.9倍になっていることが分かります。
嬉しい結果となりました。満足です。

コードの利用に関しては連絡先とサイトについてをご覧下さい。

ガンマ関数の概形

最後に、ガンマ関数の概形です。

参考文献

[1] Gamma Function -mathworld
[2]Reciprocal gamma function -wikipedia
[3]Gamma function – OeisWiki
[4] Wrench, J.W. (1968). Concerning two series for the gamma function. Mathematics of Computation, 22, 617–626
[5]Wrench, J.W. (1973). Erratum: Concerning two series for the gamma function. Mathematics of Computation, 27, 681–682
※[4]の係数表には間違いがあり、同著者によって[4]で訂正されています。
[6]The Lanczos Approximation -Boost
[7]Lanczos_approximation -wikipedia
[8] Abramowitz & Stegun (1965), ISBN 978-0486612720 (Numerical Method, 6.1.40, Page 257)
[8] Abramowitz & Stegun (1965), ISBN 978-0486612720 (Numerical Method, 6.1.34, Page 256)

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

Javaによる任意精度演算

コメントする

javaによる任意精度演算のプログラムです。
fortranでやろうとしましたがサポートされてないようです。

私はjavaに精通していません。とりあえず何のプログラムでもいいので任意精度演算を使おうとした際、javaが一番手軽だっただけです。

ここでは、ベルヌーイ数を計算することを考えます。
擬コードとして、wikipediaより、

for m in 0, 1, …, n do
     A[m] ← 1 / (m + 1)
     for j in m, m − 1, …, 1 do
         A[j − 1] ← j × (A[j − 1] - A[j])
     end for
end for
return A[0]    # which is Bn

でベルヌーイ数が得られます。
倍精度演算(10進数で16桁)では桁落ちのため\(B_{14}\)で6桁しか一致しません。

計算過程で100桁まで用いて(割り算で100桁目で四捨五入)、\(B_{50}\)を計算してみましょう。

import java.math.BigDecimal;
 
public class bernoulli {

  final static int divideN=100;
 
  public static void main(String[] args) {
    int m,j;
    BigDecimal b  = new BigDecimal("0");
    BigDecimal bm = new BigDecimal("0");
    BigDecimal a[]= new BigDecimal[51];

    a[0]=new BigDecimal("0");

    for(m=0;m<=50;m++){
       bm = BigDecimal.valueOf(m);
         a[m]=(new BigDecimal("1")).divide(bm.add(new BigDecimal("1")), divideN, BigDecimal.ROUND_HALF_UP);
       for(j=m;j>=1;j--){
             a[j-1]=(BigDecimal.valueOf(j)).multiply(a[j-1].subtract(a[j]));
         }
       b=a[0];
       String text=String.format("%.36e",b);
       System.out.println(text);
      }
    }
}

コンパイルと実行は

javac bernoulli.java
java bernoulli

です。

出力は36桁としています。上記プログラムによる結果は、
7.500866746076964366855720075757575758e+24
となり、wolframによる値Table[bernoulli b(i), {i,0,50}] with 36 digits
7.50086674607696436685572007575757576e+24
です。全桁一致しました。
[adsense1]

以下、javaの任意精度演算に関するメモです。

メモ

変数の定義方法
  • 整数型
    int i,j;
  • 倍精度実数型
    double x;
  • 任意精度型
    BigDecimal a,b,c;
配列の定義方法
BigDecimal r[]= new BigDecimal[51];

↑で配列r[0]~r[50]が確保される。

配列の初期化

for文でまわしましょう。

for(i=0;i<=N;i++){
    r[i]=new BigDecimal("0");
}
変数へ値の代入

aに任意精度型で\(1.2000000000\ldots\)を代入する。

a=new BigDecimal("1.2");
足し算
\(
a=b+1.2
\)
a=b.add(new BigDecimal("1.2"))

(整数型の型変換)
\(
a=b+i
\)
\(i\)は整数型

a=b.add(BigDecimal.valueOf(i));

引き算
\(
a=b-1.2
\)
a=b.subtract(new BigDecimal("1.2"));

掛け算
\(
a=b\times c
\)
a=b.multiply(c);

割り算

\(
a=b/c
\)
丸めの桁を100桁にする時

a=b.divide(c, 100, BigDecimal.ROUND_HALF_UP);

足し算と割り算

\(
a=b/(c+9.92)
\)
丸めの桁を100桁にする時

a=b.divide(c.add(new BigDecimal("9.92")), 100, BigDecimal.ROUND_HALF_UP);

累乗
\(
a=b^c
\)
a=b.pow(c);

ループ

\(
\displaystyle s=\sum_{i=0}^{N-1}a_i
\)

b=new BigDecimal("0");
for(i=0;i<=N-1;i++){
    b=b.add(r[i]);
}

端末への指数表記,有効桁数36桁で出力
String text=String.format("%.36e",a);
System.out.println(text);

[adsense2]

参考文献

【Java】任意精度整数(多倍長整数)で演算する -ヽ｜∵｜ゝ(Fantom) の 開発blog？
Java BigDecimalの四則演算とフォーマット処理 -ITのおもちゃ箱
intとBigDecimalの型変換 -Javaプログラミング講座