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

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

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


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

数値的に計算します。

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

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

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

を考えます。

数値計算方法


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

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

です。読み取れる情報は

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

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

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

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

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

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

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

ジョルダンの補題の利用


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

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

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

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

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

数値計算プログラム


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

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

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

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

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

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

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

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

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

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

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

  return
end function g

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

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

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

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

解法


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

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

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

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

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

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

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

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

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

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

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

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

固定端条件の場合


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

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

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

周期境界条件の場合


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

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

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


スポンサーリンク


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


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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

数値計算コード


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

スポンサーリンク

実例


自由粒子

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

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


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

gnuplotコードはこちら

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

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

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

gnuplotコードはこちら


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

参考文献


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

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

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

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

数値積分の種類

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

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

解析的に計算


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

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

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

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

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

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

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

です。

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


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


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

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

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

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

で実装できます。

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

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

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


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

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

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

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

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

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

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


スポンサーリンク


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


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

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

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

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

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

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

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

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


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

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

端点に特異性がある場合


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


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

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

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

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

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


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

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

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

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

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

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

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

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

二重積分を行いたい場合


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

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

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

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

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

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

表面積分を行いたい場合


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

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

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

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


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

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

3. Lebedev angular quadrature -Methods of Numerical Integration.

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

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

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


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

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

と求める事が出来ます。

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

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

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

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

スポンサーリンク

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


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

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

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

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

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

注釈

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

モンテカルロ積分

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

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

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

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


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

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

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

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

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

と求める事が出来ます。

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

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

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

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

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

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

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

乱数について


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

fortran90コード

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

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

下のコードを実行すると

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

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

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

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

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

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

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

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

gfortran mt19937.f90 main.f90

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

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

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

3次元積分


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

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

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

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

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

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

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

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

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

参考文献


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

スプライン補間(1, 2次元)

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つの平均を取ってしまう事でしょう。しかし、単純に考えて計算時間が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

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

個別のf90ファイル

http://slpr.sakura.ne.jp/qp/supplement_data/tise1d_by_sinbasis_and_quadpack/main.f90
http://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}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プログラム(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)
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プログラム


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

1次元実積分


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

1次元複素積分


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,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のコード
http://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のコード
http://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追)
http://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
スポンサーリンク

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


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

参考文献


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}\)の関係を用いて最後にまとめて定数倍しています。

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

という結果を得ます。

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