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

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

非相対論的な水素様原子の電子の波動関数を求める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}P_l^m(x)\right]=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プログラム
  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)
http://slpr.sakura.ne.jp/qp/nc-integral/
非適応型 ロンバーグ積分
Sikino(de)
http://slpr.sakura.ne.jp/qp/de-fourmula/
非適応型 二重指数関数型数値積分
Sikino(rkf451)
http://slpr.sakura.ne.jp/qp/runge-kutta-ex/
適応型 ルンゲ-クッタ-フェールベルグ法による積分

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


スポンサーリンク


ライセンスについて


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

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問を選んで比較します。問題の抜粋は私が独断で選びました。

結果


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

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

を比較します。また、適応型の積分手法に対して、実際にどの点が評価されているかを知りたかったため、
実際に評価された点達に対し矩形波で畳み込みを行い、最大値を1にするように規格化し比較しました。
これが左から1枚目の画像です。
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プログラム


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

Quadpackのコード
http://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

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

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

Quadpackのコード
http://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,3d0)
  b=dcmplx(1d0,2d0)
  eps=1d-12

  call cqag_k2(g,a,b,eps,s,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=sqrt(z)
 
  return
end function g
スポンサーリンク

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


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

参考文献


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
シキノート
http://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”

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

有限差分法の正体

まとめ


微分を中心差分に変える
\(
\begin{align}
\frac{d}{dx}\psi(x)&\approx \frac{\psi(x_{i+1})-\psi(x_{i-1})}{2h} \\
\frac{d^2}{dx^2}\psi(x)&\approx \frac{\psi(x_{i+1})-2\psi(x_i)+\psi(x_{i-1})}{h^2}
\end{align}
\)

という操作によって求めたハミルトニアンの行列要素は、空間に局在化した正規直交で選点直交性を持つ基底関数
\(
\begin{eqnarray}
\varphi_i(x)=
\left\{
\begin{aligned}
& \frac{1}{\sqrt{h}} ~~~~~~(x_i-h/2 \lt x \lt x_i+h/2) \\
& \frac{1}{2\sqrt{h}} ~~~~~~(x = x_i\pm h/2) \\
& 0 ~~~~~~(上以外)
\end{aligned}
\right.
\end{eqnarray}
\)

で展開したものと同一です。すなわち、行列要素は
\(
\begin{align}
\langle \varphi_i|\hat{H}|\varphi_j \rangle
&=\int_{-\infty}^{\infty} \varphi_i(x)\left[-\frac{1}{2}\frac{d^2}{dx^2}+V(x)\right]\varphi_j(x)dx\\
&=-\frac{1}{2h^2}\left[\delta_{j,i+1}-2\delta_{j,i}+\delta_{j,i-1}\right]+V(x_j)\delta_{j,i}
\end{align}
\)

と求められます。また、関数\(f(x)\)の期待値は
\(
\begin{align}
&\int_{-\infty}^{\infty} \psi(x) f(x)\psi(x) dx\approx\sum_i f(x_i)
\end{align}
\)

と求められます。


有限差分法で固有値問題を数値計算で解きます。
例として、時間依存しないシュレーディンガー方程式を解いてみます。
同様の方法で波動方程式の固有振動数、固有モードを得ることが出来ます。

wikipediaをはじめとする解説書、解説ページでは、
微分方程式を解くために微分を有限差分近似(差分商)で置き換えて得られる差分方程式で近似する
する解法だ、とあります。

一番簡単な微分を近似する方法としてオイラー法があります。
オイラー法は
\(
\displaystyle \frac{df}{dx}\approx \frac{f(x+\Delta x)-f(x)}{\Delta x}
\)

の形で有限の間隔\(\Delta x\)で微分を差分に置き換える方法です。

この、微分を差分に置き換える近似によって固有値問題である時間依存しないシュレーディンガー方程式
\(
\displaystyle \left[-\frac{1}{2}\frac{d^2}{dx^2}+V(x)\right]\psi(x)=E\psi(x)
\)

を解きます。二階微分の差分近似式は、一階微分の近似式を2階用います。
すなわち、二階微分は一階微分\(\frac{df}{dx}\equiv g(x)\)の一階微分と表現できるので、前進差分を用いて
\(
\begin{align}
\frac{d^2f}{dx^2}&\approx\frac{d}{dx}\left[\frac{g(x+\Delta x)-g(x)}{\Delta x}\right] \\
\end{align}
\)

と表せます。各々に後退差分を適用すると、
\(
\begin{align}
\frac{d}{dx}g(x+\Delta x)&\approx\frac{f(x+\Delta x)-f(x)}{\Delta x} \\
\frac{d}{dx}g(x)&\approx\frac{f(x)-f(x-\Delta x)}{\Delta x}
\end{align}
\)

より、二階微分は
\(
\begin{align}
\frac{d^2f}{dx^2}\approx\frac{f(x+\Delta x)-2f(x)+f(x-\Delta x)}{\Delta^2}
\end{align}
\)

となります。これを原子単位系のシュレーディンガー方程式に代入すると
\(
\begin{align}
\displaystyle -\frac{1}{2}\frac{\psi(x+\Delta x)-2\psi(x)+\psi(x-\Delta x)}{\Delta x^2}+V(x)\psi(x)=E\psi(x)
\end{align}
\)

より
\(
\begin{align}
-\frac{1}{2\Delta x^2}\psi(x+\Delta x)+\left(\frac{1}{\Delta x^2}+V(x)\right)\psi(x)-\frac{1}{2\Delta x^2}\psi(x-\Delta x)=E\psi(x)
\end{align}
\)

という差分方程式が得られます。
固定端条件(\(\psi(x_0)=\psi(x_{N+1})=0\), 区間\([a\sim b]=[x_0\sim x_{N+1}]\)では)の下で、行列表現では、以下のように書き下せます。

周期境界条件の場合、区間\([a\sim b]=[x_1\sim x_N]\)では、

となります。

こののち、左辺を対角化することでエネルギー固有値と固有関数を得ることが出来ます。

差分で対角化→本当に固有値?


さて、以上のように求めた固有値はいったい何なのでしょうか。
何なのでしょうか?とは、
微分を差分で近似して、何か知らないけど対角化した。そしたらエネルギー固有値と固有ベクトルが出てきた。
これは何に基づいた計算なのでしょうか。この方法で求めた固有ベクトルは固有関数に何故一致するのでしょうか。

私はどうもここがすっきりせず、差分による固有値問題の解法が何となく使いたくありませんでした。

ここからは、変分法の原理に基づき、矩形関数による展開を行うと上記の差分による展開式と同じ行列表現が得られることを示していきます。

基底関数による考え


上記の通り、演算子が行列表現として表せるのであれば、それに応じた基底関数があってもおかしくありません。
上記のハミルトニアン行列を導き出す方法として、位置的に局在した長方形の関数を考えます。

この基底関数を考えることで上記の行列表現を導くことが出来ます。

空間の離散化


まず、波動関数\(\psi(x)\)をディラックのデルタ関数を用いて
\(
\displaystyle\psi(x)=\int \psi(x’)\delta(x-x’)dx
\)

と表します。これは明らかに成立していますね。右辺を計算すればまさに左辺です。
これは、

空間のある1点\(x\)における波動関数の値\(\psi(x)\)は、関数\(\psi(x’)\)にデルタ関数を掛けて積分したものだ。

と考えることが出来ます。別の言い方をすれば、滑らかな波動関数が無限個の点の集合で書かれているということです。

この、無限個の点の集合を有限個の点の集合と近似して、空間を離散化します。
離散化に物理的な意味は無く、数値計算の都合上と考えてよいです。


※ここで、\(\Delta x\)と\(h/latex]が混在していますが、どちらも同じものを指します。
いつも[latex]\Delta x=h\)です。

ディラックのデルタ関数と離散化した関数は積分が等しくなるように決めます。
図のように考えますと、結局
\(
\displaystyle \delta(x-x’)\xrightarrow{離散化}\frac{\delta_{i,j}}{\Delta x}
\)

とするのが良いと分かります。ここで、\(\delta_{i,j}\)はクロネッガーのデルタを表し、\(x\)を離散化した座標\(x_i\)がある座標\(x_j\)に等しい時だけ値\(1\)を持つことを意味します。それ以外はゼロです。

以上の事から、空間の離散化はクロネッガーのデルタを用いて表現できることが分かりました。
これはすなわち、波動関数
\(
\displaystyle\psi(x)=\int \psi(x’)\delta(x-x’)dx
\)

を離散化すると、
\(
\begin{align}
\psi(x_i)&=\sum _j \psi(x_j)\frac{\delta_{i,j}}{\Delta x}\cdot \Delta x \\
&=\sum _j \psi(x_j)\delta_{i,j}
\end{align}
\)

と表現できます。ここで、積分を区分求積に変えたため、\(\Delta x\)が掛けられています。

離散化した関数は、

区間\(\displaystyle \left[x_i-\frac{h}{2}, x_i+\frac{h}{2}\right]\)の波動関数\(\psi(x)\)の値を\(\psi(x_i)\)で代表させて近似する

ということです。

※この離散化の考えは自然ですが、非常に非物理的であることに注意してください。なぜならば、\(x_i\pm \frac{h}{2}\)で隣り合う区間の波動関数と、ほとんどの場合不連続になるためです。量子力学では波動関数は滑らかにつながっていなければなりません(カフス等を除いて)。

ハミルトニアンの行列表現


これから考えていく基底関数を用いてハミルトニアンを行列表示します。
ハミルトニアン\(\hat{H}\)は
\(
\begin{align}
\hat{H}=-\frac{1}{2}\frac{d^2}{dx^2}+V(x)
\end{align}
\)

です。


今回用いる基底関数は空間を離散化したもので、数式で表せば、
\(
\begin{eqnarray}
\varphi_i(x)=
\left\{
\begin{aligned}
& \frac{1}{\sqrt{h}} ~~~~~~(x_i-h/2 \lt x \lt x_i+h/2) \\
& \frac{1}{2\sqrt{h}} ~~~~~~(x = x_i\pm h/2) \\
& 0 ~~~~~~(上以外)
\end{aligned}
\right.
\end{eqnarray}
\)

となります。係数の値は正規直交化
\(
\displaystyle \int_{-\infty}^{\infty} \varphi_i(x)\varphi_j(x) dx =\delta_{i,j}
\)

を満たすように決定しています(実関数なので、複素共役は消しています)。

この基底関数は正規直交であるだけでなく、選点直交性と呼ばれる直交性をも持ちます。それは、
\(
\displaystyle \varphi_i(x_j)=\frac{1}{\sqrt{h}}\delta_{i,j}
\)

という直交性です。この性質は積分を計算する際に非常に役立ちます。

波動関数が
\(
\displaystyle \psi(x)=\sum_i c_i \varphi_i(x)
\)

で近似されることを仮定し、変分法の考えに基づいて係数\(c_i\)を求めます。

要素の計算


ハミルトニアン要素を計算していきます。

まず、ポテンシャルの項の方が簡単なのでこちらを計算してしまいます。
\(
\begin{align}
\int_{-\infty}^{\infty} \varphi_i(x)V(x)\varphi_j(x)dx
&=\frac{1}{\sqrt{h}}\int_{x_i-\frac{h}{2}}^{x_i+\frac{h}{2}} V(x)\varphi_j(x)dx \\
&=V(x_i)\delta_{i,j} \\
\end{align}
\)

と書けます。\(i,j\)が同じでなければ、基底関数がゼロになってしまうので値を持ちません。よって上が導けます。簡単ですね。
補足しますと、この基底関数系では任意の関数\(f(x)\)を基底関数で挟んだ値は
\(
\displaystyle \int_{-\infty}^{\infty} \varphi_i(x) f(x)\varphi_j(x) dx \approx f(x_i)\delta_{i,j}
\)

と計算できます。よって、任意の関数\(f(x)\)の系の期待値は
\(
\begin{align}
&\int_{-\infty}^{\infty} \psi(x) f(x)\psi(x) dx \\
&\approx \sum_i \sum_j \int_{-\infty}^{\infty} \varphi_i(x) f(x)\varphi_j(x) dx \\
&=\sum_i \sum_j f(x_i)\delta_{i,j} \\
&=\sum_i f(x_i)
\end{align}
\)

と計算されます。非常に簡単ですね。

では、運動エネルギー項を計算しましょう。
まず一階微分を含む積分は
\(
\begin{align}
\int_{-\infty}^{\infty} \varphi_i(x)\frac{d}{dx}\varphi_j(x)dx
&=\frac{1}{\sqrt{h}}\int_{x_i-\frac{h}{2}}^{x_i+\frac{h}{2}} \frac{d}{dx}\varphi_j(x)dx \\
&=\frac{1}{\sqrt{h}}\int_{\varphi(x_i-\frac{h}{2})}^{\varphi(x_i+\frac{h}{2})} d\varphi_j(x) \\
&=\frac{1}{\sqrt{h}}\left[\varphi_j(x_i+\frac{h}{2})-\varphi_j(x_i-\frac{h}{2})\right] \\
&=\frac{1}{\sqrt{h}}\left[\frac{1}{2\sqrt{h}}\left(\delta_{j,i+1}+\delta_{j,i}-\delta_{j,i-1}-\delta_{j,i}\right)\right] \\
&=\frac{1}{2h}\left(\delta_{j,i+1}-\delta_{j,i-1}\right)
\end{align}
\)

となります。これはまさに中央差分です。二階微分は、

\(
\begin{align}
&\int_{-\infty}^{\infty} \varphi_i(x)\frac{d^2}{dx^2}\varphi_j(x)dx \\
&=\frac{1}{\sqrt{h}}\left[\left.\frac{d\varphi_j}{dx}\right|_{x=x_i+h/2}-\left.\frac{d\varphi_j}{dx}\right|_{x=x_i-h/2}\right]\\
&~~~~~~~\left(\frac{d\varphi_j}{dx}=\frac{1}{\sqrt{h}}\left[\delta\left(x-(x_j-h/2)\right)-\delta(x-(x_j+h/2))\right]\mbox{より}\right) \\
&=\frac{1}{h}\left[\delta(x_i+h/2-(x_j-h/2))-\delta(x_i+h/2-(x_j+h/2)\right. \\
&~~~~~~~~~~~~ \left. -\delta(x_i-h/2-(x_j-h/2))+\delta(x_i+h/2-(x_j+h/2))\right] \\
&=\frac{1}{h}\left[\delta(x_i-x_j+h)-2\delta(x_i-x_j)+\delta(x_i-x_j-h)\right] \\
&~~~~~~~\delta(x_i-x_j)\rightarrow\frac{\delta_{i,j}}{h}\mbox{の関係を用いて} \\
&=\frac{1}{h^2}\left[\delta_{j,i+1}-2\delta_{j,i}+\delta_{j,i-1}\right] \\
\end{align}
\)

と計算できます。
最後の式はまさに微分を差分で近似したものと一緒です。

ということは、
微分を上記の差分に置き換えるということは、矩形型の基底関数で展開したものと同じ事が示せました。

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

まとめ


離散フーリエ変換(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}
\)

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


スポンサーリンク


目次


  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]です。

この2つの違いは周波数か、角周波数のどちらが本質であるか?という事だと思います。
物理では、周波数よりも角周波数の方が式が簡便になります。
ですが、人間の理解がしやすいのは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}
\)

という定義と同一です。
ただし、これは積分をそのまま離散化したものではないため、積分の意味は薄れてしまいます。
その結果、畳み込みなど、フーリエ変換後同士の物を掛け算する時は違った結果になってしまいます。
このページではこの簡便さの為の離散フーリエ変換は使いません。

上の2つの定義が同じ理由は、順方向→逆方向を行い、元に戻った時に係数が\(\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のルーチンには起こっていません。どういうアルゴリズムか知りませんが、優秀なのでしょう。

スポンサーリンク


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

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求積法

クレンショウ-カーティス(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=10
  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

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]ランダウ=リフシッツ,「量子力学1」第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 (計算例)をご覧ください。


スポンサーリンク



解説


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

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


時間依存シュレーディンガー方程式は
\(
\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\)と置きました。

すなわち、2つの式
\(
\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)に対して2つの方向に対して変分を取ると、
\(
\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)間接法…変分パラメータとしてある正規直交基底関数系を採用する
の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を用いるのが良いでしょう。
対角化のプログラムに関する詳細は対角化へどうぞ。

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

スポンサーリンク

コンパイルと実行は

$ 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]ランダウ=リフシッツ,「量子力学1」第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くらいの時間で計算できるため、圧倒的に早いです。


スポンサーリンク


方針


方針としては

  • 漸化式
  • 相反公式
  • 逆ガンマ関数の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
スポンサーリンク

計算速度の比較


大浦様のコードと計算速度を比較します。
比較方法は、
乱数を[-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
です。全桁一致しました。

スポンサーリンク



以下、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);
スポンサーリンク

参考文献


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

逆行列の数値計算

Lapackを用いて数値計算で逆行列を求めます。
基本的に逆行列の計算は桁落ちなどの問題があるため、基本的に求めるべきではありません。
そうはいっても使いたい時があります。そんな状況を考えましょう。

Lapack, もしくはインテルのマス・カーネル・ライブラリ(MKL)を用います。

用いるルーチンは

?getrf, ?getri

の2つです。

ここに用意したのは配列の大きさNと配列Aを入れると逆行列を出力するプログラムです。
倍精度、複素倍精度のルーチンを書きました。

ifort -mkl main.f90

でコンパイルしてください。

program main
  implicit none
  integer::i,j,N
  complex(kind(0d0)),allocatable::A(:,:),B(:,:)
 
  N=3
  allocate(A(1:N,1:N),B(1:N,1:N))
  A(1,1:3)=(/(2d0,0d0),(3d0,0d0),(0d0,0d0)/)
  A(2,1:3)=(/(2d0,0d0),(9d0,0d0),(1d0,3d0)/)
  A(3,1:3)=(/(3d0,4d0),(5d0,0d0),(3d0,1d0)/)

  B=A
 
  do i=1,N
     do j=1,N
        write(6,'(2f15.10,$)')A(i,j)
     enddo
     write(6,*)
  enddo
  write(6,*)

  call zinvmat(size(A,1),A)
 
  do i=1,N
     do j=1,N
        write(6,'(2f15.10,$)')A(i,j)
     enddo
     write(6,*)
  enddo
  write(6,*)

  A=matmul(B,A)
 
  do i=1,N
     do j=1,N
        write(6,'(2f15.10,$)')A(i,j)
     enddo
     write(6,*)
  enddo
  write(6,*)

  stop
end program main
   
subroutine invmat(N,A)
  implicit none
  integer,intent(in)::N
  double precision,intent(inout)::A(1:N)

  integer::info,lwork,ipiv(1:N)
  double precision::tw(1:1)
  double precision,allocatable::work(:)
 
  info=0
  call dgetrf(N,N,A,N,ipiv,info)
  if(info.ne.0)then
     write(6,*)"Error at inverse matrix ?getrf, info",info
     stop
  endif
 
  !size query
  call dgetri(N,A,N,ipiv,tw,-1,info)
  lwork=nint(tw(1))
  allocate(work(1:lwork+1))
  work=0d0
 
  call dgetri(N,A,N,ipiv,work,lwork,info)
  if(info.ne.0)then
     write(6,*)"Error at inverse matrix ?dgetri, info",info
     stop
  endif
 
  return
end subroutine invmat


subroutine zinvmat(N,A)
  implicit none
  integer,intent(in)::N
  complex(kind(0d0)),intent(inout)::A(1:N)

  integer::info,lwork,ipiv(1:N)
  complex(kind(0d0))::tw(1:1)
  complex(kind(0d0)),allocatable::work(:)
 
  info=0
  call zgetrf(N,N,A,N,ipiv,info)
  if(info.ne.0)then
     write(6,*)"Error at inverse matrix ?getrf, info",info
     stop
  endif
 
  !size query
  call zgetri(N,A,N,ipiv,tw,-1,info)
  lwork=nint(dble(tw(1)))
  allocate(work(1:lwork+1))
  work=dcmplx(0d0,0d0)
 
  call zgetri(N,A,N,ipiv,work,lwork,info)
  if(info.ne.0)then
     write(6,*)"Error at inverse matrix ?dgetri, info",info
     stop
  endif
 
  return
end subroutine zinvmat