sikino のすべての投稿

スプライン補間(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

1/(x^n)のフーリエ変換

その他フーリエ変換
ヘヴィサイド関数のフーリエ変換
1のフーリエ変換
\(x\)のフーリエ変換
\(1/x\)のフーリエ変換
\(1/(x^n)\)のフーリエ変換 ← 今ここ

\(\frac{1}{x}\)のフーリエ変換から\(\frac{1}{x^n}\)のフーリエ変換を導けます。
\(
\begin{align}
\int_{-\infty}^{\infty}\frac{1}{x}e^{-ikx}dx =-\pi i ~\text{sgn}(k)
\end{align}
\)

(導出は\(\frac{1}{x}\)のフーリエ変換へ)

部分積分を用いると\(\frac{1}{x^2}\)のフーリエ変換
\(
\begin{align}
&\int_{-\infty}^{\infty}\frac{1}{x}\left(\frac{1}{-ik}e^{-ikx}\right)’dx \\
&=\left[\frac{1}{x}\cdot\frac{1}{-ik}e^{-ikx}\right]_{-\infty}^{\infty} -\frac{1}{ik}\int_{-\infty}^{\infty}\frac{1}{x^2}e^{-ikx} dx
\end{align}
\)
第一項はゼロになります。第二項と残りの項と\(1/x\)のフーリエ変換の右辺が等しいので、
\(
\begin{align}
\int_{-\infty}^{\infty}\frac{1}{x^2}e^{-ikx} dx =-\pi k ~\text{sgn}(k)
\end{align}
\)
が得られます。
これを繰り返すことで\(\frac{1}{x^n}\)のフーリエ変換
\(
\begin{align}
\int_{-\infty}^{\infty}\frac{1}{x^n}e^{-ikx}dx &=-\pi i \frac{(-ik)^{n-1}}{(n-1)!}~\text{sgn}(k)
\end{align}
\)

が導かれます。

1/xのフーリエ変換

\begin{align}
\int_{-\infty}^{\infty}\frac{1}{x}e^{-ikx}dx =-\pi i ~\text{sgn}(k)
\end{align}
※ただし、\(x=0\)での値は定義しない。また、\(\text{sgn}(x)\)は符号関数

数値計算による上式の確かめはexp(-ikx)/x の無限区間に渡る積分を数値的に計算したい。をご覧ください。

その他フーリエ変換
ヘヴィサイド関数のフーリエ変換
1のフーリエ変換
\(x\)のフーリエ変換
\(1/x\)のフーリエ変換 ← 今ここ
\(1/(x^n)\)のフーリエ変換


実軸上の積分は複素平面の上から近づいた値と下から近づいた値の平均をとって定めるとします。すなわち、
\(
\begin{align}
&\int_{-\infty}^{\infty}\frac{1}{x}e^{-ikx}dx \nonumber \\
&=\lim_{\varepsilon\to +0} \frac{1}{2}\left[\int_{-\infty}^{\infty}\frac{1}{x+i\varepsilon}e^{-ikx}dx +\int_{-\infty}^{\infty}\frac{1}{x-i\varepsilon}e^{-ikx}dx \right]
\end{align}
\)
で表されると仮定します。
そして、第1項、第2項をそれぞれ計算していきます。

\(
\displaystyle \lim_{\varepsilon\to +0}\int_{-\infty}^{\infty}\frac{1}{x\pm i\varepsilon}e^{-ikx}dx
\)

を計算するにあたり、複素x平面上の漸近で収束する半径が\(k\)の符号によって異なるので場合分けをしてます。

\(k\gt 0\)の時

まず、\(\varepsilon\)の符号が正である場合を考えます。
この実軸上の積分は実軸への近づき方を指定して
\(
\begin{align}
\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}\frac{1}{x+i\varepsilon}e^{-ikx}dx =
\left[\int_{-\infty}^{-\varepsilon’}+\int_{\Gamma_+}+\int_{\varepsilon’}^{\infty}\right]\frac{e^{ikx}}{x}dx
\end{align}
\)

です。ここで、\(\int_{\Gamma_+}\)は特異点周りを複素平面上半面を通る経路です。
この積分を計算するために閉路積分を考えます。すなわち、以下の図の左側のように回った積分を考えます。

すると、この閉経路に対して、
\(
\begin{align}
\oint \frac{e^{-ikx}}{x}dx&=
\lim_{R\to\infty}\left[\int_{-R}^{-\varepsilon’}+\int_{\Gamma_+}+\int_{\varepsilon’}^{R}+\int_{\Gamma’}\right]
\frac{e^{ikx}}{x}dx\\
&=\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}\frac{1}{x+i\varepsilon}e^{-ikx}dx+\lim_{R\to\infty}\int_{\Gamma’}\frac{e^{-ikx}}{x}dx
\end{align}
\)
が成立します。今、閉経路内に特異点が一つ存在し、負の方向に回っているので、左辺\(\oint \frac{e^{-ikx}}{x}dx=-2\pi i\)です。また、大きく外側を回る経路\(\Gamma’\)の積分はゼロです。よって、
\(
\begin{align}
\lim_{\varepsilon\to +0}\int_{-\infty}^{\infty}\frac{1}{x+i\varepsilon}e^{-ikx}dx =-2\pi i.
\end{align}
\)

が示されます。

また、特異点を下に回る積分(図の左の経路)では、閉経路内に特異点は存在しないので、\(\oint \frac{e^{-ikx}}{x}dx=0\)です。なので、
\(
\begin{align}
\lim_{\varepsilon\to +0}\int_{-\infty}^{\infty}\frac{1}{x-i\varepsilon}e^{-ikx}dx =0.
\end{align}
\)

と導けます。

\(k\lt 0\)の時

まず、\(\varepsilon\)の符号が正である場合を考えます。

kが正の時と異なる点は、漸近で収束する領域が複素x平面の上半面であるという点です。

図の左のような閉経路を考えます。

同様に計算を行うと
\(
\begin{align}
\oint \frac{e^{-ikx}}{x}dx&=
\lim_{R\to\infty}\left[\int_{-R}^{-\varepsilon’}+\int_{\Gamma_+}+\int_{\varepsilon’}^{R}+\int_{\Gamma’}\right]
\frac{e^{ikx}}{x}dx\\
&=\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}\frac{1}{x+i\varepsilon}e^{-ikx}dx+\lim_{R\to \infty}\int_{\Gamma’}\frac{e^{-ikx}}{x}dx
\end{align}
\)
今、閉経路内に特異点を含まないので、左辺\(\oint \frac{e^{-ikx}}{x}dx=0\)、また遠方ではゼロに近づくので、\(\int_{\Gamma’}\frac{e^{ikx}}{x}dx=0\)です。よって、
\(
\begin{align}
\lim_{\varepsilon\to +0}\int_{-\infty}^{\infty}\frac{1}{x+i\varepsilon}e^{-ikx}dx =0.
\end{align}
\)

また、特異点周りを正方向に回る場合、閉経路内に特異点を含むので、
\(
\begin{align}
\lim_{\varepsilon\to +0}\int_{-\infty}^{\infty}\frac{1}{x-i\varepsilon}e^{-ikx}dx =2\pi i.
\end{align}
\)

となります。

\(k= 0\)の時

\(k=0\)では被積分関数は漸近で収束しません。なので漸近で大きく回った時の値がゼロになりません。

ここからは誤解を恐れずに計算します。
今求めたい積分は
\(
\int_{-\infty}^{\infty}\frac{1}{x}e^{-i0x}dx =\int_{-\infty}^{\infty}\frac{1}{x}dx
\)

です。実軸上の積分がこのページで書いた
\(
\begin{align}
&\int_{-\infty}^{\infty}\frac{1}{x}e^{-ikx}dx \nonumber \\
&=\lim_{\varepsilon\to +0} \frac{1}{2}\left[\int_{-\infty}^{\infty}\frac{1}{x+i\varepsilon}e^{-ikx}dx +\int_{-\infty}^{\infty}\frac{1}{x-i\varepsilon}e^{-ikx}dx \right]
\end{align}
\)
で定義されているのだとすれば、
\(
\begin{align}
\int_{-\infty}^{\infty}\frac{1}{x}dx =\lim_{\varepsilon\to +0} \frac{1}{2}\left[\int_{-\infty}^{\infty}\frac{1}{x+i\varepsilon}dx +\int_{-\infty}^{\infty}\frac{1}{x-i\varepsilon}dx \right]
\end{align}
\)

です。簡単のため、特異点を含まない閉経路を取れば第1項は
\(
\begin{align}
\lim_{R\to \infty}
\left[\int_{-R}^{-\varepsilon’}+\int_{\Gamma_+}
+\int_{\varepsilon’}^{R}+\int_{\Gamma’} \right]\frac{1}{x}dx =\oint \frac{1}{x}dx=0
\end{align}
\)

なので、
\(
\begin{align}
\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}\frac{1}{x+i\varepsilon}dx +\int_{\Gamma’} \frac{1}{x}dx =0
\end{align}
\)

となります。外側を回る経路\(\Gamma’\)を通った積分は
\(
\begin{align}
\int_{\Gamma’} \frac{1}{x}dx &=\int_0^\pi \frac{1}{re^{i\theta}}ire^{i\theta} d\theta \\
&=\pi i
\end{align}
\)

と計算できてしまうので、代入すれば
\(
\begin{align}
\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty} \frac{1}{x+i\varepsilon} dx=-\pi i
\end{align}
\)

となります。また、同様にして
\(
\begin{align}
\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty} \frac{1}{x-i\varepsilon} dx=\pi i
\end{align}
\)

を得ます。
よって、実軸上の積分の式に代入すれば
\(
\begin{align}
\int_{-\infty}^{\infty}\frac{1}{x}dx =\lim_{\varepsilon\to +0} \frac{1}{2}\left[(-\pi i) +\pi i \right] =0
\end{align}
\)

を得ます。すなわち、フーリエ変換後の\(k=0\)の値はゼロということです。

計算結果のまとめ


以上をまとめると、3つの関係式
\(
\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}
\)

\(
\begin{align}
\int_{-\infty}^{\infty}\frac{1}{x}e^{-i0x}dx=0
\end{align}
\)

を得ることが出来ました。
これを実軸上の積分の式に代入すると、
\(
\begin{align}
\int_{-\infty}^{\infty}\frac{1}{x}e^{-ikx}dx =-\pi i ~\text{sgn}(k)
\end{align}
\)

となります。ここで、\(\text{sgn}(x)\)は符号関数で、
\(
\begin{eqnarray}
\text{sgn}(x) =
\left\{
\begin{aligned}
-1~~~(x\lt 0)\\
0~~~(x=0) \\
1~~~(x\gt 0)
\end{aligned}
\right.
\end{eqnarray}
\)
で定義される関数です。

xのフーリエ変換

その他フーリエ変換
ヘヴィサイド関数のフーリエ変換
1のフーリエ変換
\(x\)のフーリエ変換 ← 今ここ
\(1/x\)のフーリエ変換
\(1/(x^n)\)のフーリエ変換

1のフーリエ変換結果
\(
\begin{align}
\int_{-\infty}^{\infty}e^{-ikx}dx =2\pi\delta(k)
\end{align}
\)

を利用します(1のフーリエ変換)。

両辺を\(k\)で微分すると\(x\)のフーリエ変換
\(
\begin{align}
\int_{-\infty}^{\infty}(-i)xe^{-ikx}dx &=2\pi\delta'(k) \\
\int_{-\infty}^{\infty}xe^{-ikx}dx &=i2\pi\delta'(k)
\end{align}
\)

が導けます。

これを繰り返せば、\(x^n\)のフーリエ変換
\(
\begin{align}
\int_{-\infty}^{\infty}x^ne^{-ikx}dx &=i^n 2\pi\delta^{(n)}(k)
\end{align}
\)

も導くことが出来ます。

1のフーリエ変換

\begin{align}
\int_{-\infty}^{\infty}e^{-ikx}dx &=-\frac{1}{i}\lim_{\varepsilon\to +0}\left[\frac{1}{k+i\varepsilon }-\frac{1}{k-i\varepsilon }\right] \\
&=-\frac{1}{i}\left[\frac{1}{k+i0}-\frac{1}{k-i0}\right] \\
&=2\pi\delta(k)
\end{align}

その他フーリエ変換
ヘヴィサイド関数のフーリエ変換
1のフーリエ変換 ← 今ここ
\(x\)のフーリエ変換
\(1/x\)のフーリエ変換
\(1/(x^n)\)のフーリエ変換

1のフーリエ変換


1のフーリエ変換を求めるためには、ヘヴィサイド関数の無限積分が
\(
\begin{align}
\int_{-\infty}^{\infty} \theta(x)e^{\pm ikx}dx &=\int_{0}^{\infty} e^{\pm ikx}dx\\
&=\lim_{\varepsilon\to +0}\mp\frac{1}{i}\frac{1}{k- i\varepsilon} \\
&=\mp\frac{1}{i}\frac{1}{k\pm i0}
\end{align}
\)
であることを用います(ヘヴィサイド関数のフーリエ変換)。

1のフーリエ変換
\(
\displaystyle \int_{-\infty}^{\infty} \theta(x)e^{-ikx}dx
\)

は、ヘヴィサイド関数の無限積分を用いて
\(
\begin{align}
&\int_{-\infty}^{\infty}e^{-ikx}dx \\
&=\int_{-\infty}^{0}e^{-ikx}dx+\int_{0}^{\infty}e^{-ikx}dx \nonumber \\
&=\int_{0}^{\infty}e^{ikx}dx+\int_{0}^{\infty}e^{-ikx}dx \nonumber \\
&=-\frac{1}{i}\frac{1}{k+i0} +\frac{1}{i}\frac{1}{k-i0} \nonumber \\
&=-\frac{1}{i}\left[\frac{1}{k+i0}-\frac{1}{k-i0}\right]
\end{align}
\)
と求められます。2行目から3行目の式変形は、第1項目に関して変数変換\(x\to -x\)を行っています。

ここで、計算結果
\(
\begin{align}
\int_{-\infty}^{\infty}e^{-ikx}dx=-\frac{1}{i}\left[\frac{1}{k+i0}-\frac{1}{k-i0}\right]
\end{align}
\)

の右辺が示す意味を考えます。
適当な関数に右辺を作用させたときを考えると、この右辺がデルタ関数の定数倍に等しいことが分かります。

\(
\begin{align}
g(x)=-\frac{1}{i}\left[\frac{1}{x+i0}-\frac{1}{x-i0}\right]
\end{align}
\)

とおいて、ある関数\(f(x)\)に掛けて\(x\)で積分すると

\(
\begin{align}
&\int_{-\infty}^{\infty}f(x)g(x)dx \nonumber\\
&=\int_{-\infty}^{\infty}f(x)\left[ -\frac{1}{i }\left(\frac{1}{x+i0}-\frac{1}{x-i0}\right)\right]dx \nonumber\\
&=\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}f(x)\left[ -\frac{1}{i }\left(\frac{1}{x+i\varepsilon}-\frac{1}{x-i\varepsilon}\right)\right]dx\nonumber \\
&=-\frac{1}{i}\left[\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}f(x) \frac{1}{x+i\varepsilon}dx-\lim_{\varepsilon\to +0}\int_{-\infty}^{\infty}f(x)\frac{1}{x-i\varepsilon}dx\right] \nonumber\\
&=-\frac{1}{i}\left[
\left\{ \mathcal{P}\int_{-\infty}^{\infty}\frac{f(x)}{x}dx-\pi i f(0)\right\}
-\left\{\mathcal{P}\int_{-\infty}^{\infty}\frac{f(x)}{x}dx+\pi i f(0)\right\}
\right]\nonumber \\
&=2\pi f(0) \nonumber\\
&=\int_{-\infty}^{\infty}f(x)\cdot 2\pi \delta(x)dx
\end{align}
\)

計算前と比較すると
\(
\begin{align}
g(x)&=2\pi \delta(x) \\
&\to \delta(x)=-\frac{1}{2\pi i}\left[\frac{1}{x+i0}-\frac{1}{x-i0}\right]
\end{align}
\)
としてデルタ関数が表現されていることが分かります。

ヘヴィサイド関数のフーリエ変換

\begin{align}
\int_{-\infty}^{\infty} \theta(x)e^{-ikx}dx
=\lim_{\varepsilon\to +0}\frac{1}{i}\frac{1}{k- i\varepsilon}
=\frac{1}{i}\frac{1}{k- i0}
\end{align}
を導出します。
ここでは、アーベル総和の考えを用いません。使うのはフーリエ変換のユニタリー性です。

※いろいろ数学的に自信が持てない点があります。これは僕が納得した解釈であるので、完全に正しい導出方法ではないかもしれません。

また、ここではフーリエ変換を
\(
\begin{align}
& g(k)=\int_{-\infty}^{\infty}f(x)e^{-ikx} dx \\
& f(x)=\int_{-\infty}^{\infty}g(k)e^{ikx} \frac{dk}{2\pi}
\end{align}
\)

で定義しています。

その他フーリエ変換
ヘヴィサイド関数のフーリエ変換 ← 今ここ
1のフーリエ変換
\(x\)のフーリエ変換
\(1/x\)のフーリエ変換
\(1/(x^n)\)のフーリエ変換

ヘヴィサイド関数のフーリエ変換


ヘヴィサイド関数\(\theta(x)\)は
\(
\begin{eqnarray}
\theta(x) =
\left\{
\begin{aligned}
0~~~(x\gt 0)\\
1~~~(x\lt 0)
\end{aligned}
\right.
\end{eqnarray}
\)

と定義されます。この関数のフーリエ変換、すなわち積分
\(
\displaystyle \int_{-\infty}^{\infty}\theta(x)e^{- ikx} dx
\)

を考えます。上の積分は\(x\lt 0\)ではヘヴィサイド関数はゼロなので、
\(
\begin{align}
\int_{-\infty}^{\infty}\theta(x)e^{- ikx} dx=\int_{0}^{\infty}e^{-ikx}dx
\end{align}
\)

となります。後はこの積分が計算出来れば良いです。

\(k\ne0\)の場合


もしも\(k\ne 0\)であれば、この積分は計算可能です。
被積分関数を複素x平面に解析接続して、以下の図のような閉経路を考えます。

今、閉経路内に特異点は無いので、一周回った時の周回積分の値はゼロです。また、大きく外側を回る経路\(\Gamma’\)は、被積分関数が漸近で収束する領域なので、無限遠方での線積分の値はゼロです。このことから、
\(
\begin{align}
\oint f(x) dx &= \lim_{R\to\infty} \left[\int_0^R f(x)dx + \int_{\Gamma’}f(x)dx+\int_{R}^0 f(re^{i\theta})re^{i\theta}dr\right]\nonumber \\
0&= \lim_{R\to\infty}\left[ \int_0^R f(x)dx +\int_{R}^0 f(re^{i\theta})re^{i\theta}dr\right] \nonumber\\
\lim_{R\to\infty} \int_0^R f(x)dx &=\lim_{R\to\infty}\int_{0}^R f(re^{i\theta})re^{i\theta}dr \nonumber
\end{align}
\)
が成立します。

\(e^{-ikx}\)で\(k\ne 0\)の場合、
\(\pi \lt \theta\lt 0 (k\gt 0),~~0 \lt \theta\lt \pi (k\lt 0) \)では
\(
\begin{align}
\int_0^{\infty}e^{-ikx}dx &= \lim_{R\to\infty}\int_0^{R e^{i\theta}}e^{-ikz}dz \nonumber \\
&= \lim_{R\to\infty}\left[\frac{e^{-ikz}}{-ik}\right]^{R e^{i\theta}}_{0} \nonumber \\
&=\frac{1}{i}\frac{1}{k}
\end{align}
\)
となります。
しかし、\(k=0\)は漸近で収束しないので同様の方法で計算することができません。

\(k=0\)の場合


\(k=0\)で、被積分関数は無限大に発散します。
なぜなら、
\(
\begin{align}
\int_{0}^{\infty}e^{-i0x}dx= \int_{0}^{\infty}dx =\infty
\end{align}
\)

だからです。ここで問題なのは、無限大への発散の仕方です。

これが示すことは、フーリエ変換後の空間で\(k=0\)は特異点ということです。

以上から
\(
\begin{eqnarray}
\int_{-\infty}^{\infty}\theta(x)e^{- ikx} dx =
\left\{
\begin{aligned}
\frac{1}{i}\frac{1}{k}~~~(k\ne 0)\\
\infty~~~(k= 0)
\end{aligned}
\right.
\end{eqnarray}
\)
というフーリエ変換結果を得ました。

特異点の周り方の指定

さて、前回の計算結果から、\(k=0\)に特異点がある事が分かりました。
この特異点の性質を知るために、逆フーリエ変換を考えます。

なぜなら、フーリエ変換はユニタリー変換であることを考えれば、フーリエ変換した関数を逆フーリエ変換した時には元に戻っていて欲しい、という条件が課されているからです。

そして、逆フーリエ変換時にはこれを避けて回らねばならりません。
特異点を回る方法として考えられるのは以下の3つです。

  1. 複素k平面上で特異点周りを上方向を回る
  2. 複素k平面上で特異点周りを下方向を回る
  3. 複素k平面上で特異点周りを上方向と下方向を回った結果の平均をとる

結果的に正しいのは下方向を回った時の2です。そしてこの特異点を回る方法をフーリエ変換後の結果に含めると、
\(
\displaystyle \int_{-\infty}^{\infty}\theta(x)e^{- ikx} dx = \lim_{\varepsilon \to +0 }\frac{1}{i}\frac{1}{k-i\varepsilon}
\)

という結論を得ることが出来ます。
ここで、\(-i\varepsilon\)が意味するのは、複素k平面を積分する時には下方向を回りなさい、ということを数学的に表しています。

実際にそれぞれ計算して確かめてみましょう。
被積分関数が漸近で収束する領域は\(e^{ikx}\)なので
\(x\gt 0\)の時、複素k平面の上半面、
\(x\lt 0\)の時、複素k平面の下半面になります。
一位の特異点を閉経路内に含み、正の方向に回る場合、コーシーの積分定理より値は\(2\pi i\)です。

1. 複素k平面上で特異点周りを上方向を回る場合

フーリエ変換後の関数\(f(x)\)は
\(
\begin{eqnarray}
f(x)=\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}\frac{-i}{k+i\varepsilon}e^{ikx}\frac{dk}{2\pi} =
\left\{
\begin{aligned}
1~~~(x\lt 0) \\
0~~~(0\lt x)
\end{aligned}
\right.
\end{eqnarray}
\)

なので、これは元のヘヴィサイド関数ではないので特異点を回る方法として欲しいものではありません。

2. 複素k平面上で特異点周りを下方向を回る

フーリエ変換後の関数\(f(x)\)は
\(
\begin{eqnarray}
f(x)=\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty}\frac{-i}{k-i\varepsilon}e^{ikx}\frac{dk}{2\pi} =
\left\{
\begin{aligned}
0~~~(x\lt 0) \\
1~~~(0\lt x)
\end{aligned}
\right.
\end{eqnarray}
\)

なので、これはフーリエ変換前のヘヴィサイド関数に一致します。よってこれが適切です。

3. 複素k平面上で特異点周りを上方向と下方向を回った結果の平均をとる

フーリエ変換後の関数\(f(x)\)は1. 2. の結果から、
\(
\begin{eqnarray}
f(x)=\left\{
\begin{aligned}
1/2~~~(x\lt 0) \\
1/2~~~(0\lt x)
\end{aligned}
\right.
\end{eqnarray}
\)

なので、これは元のヘヴィサイド関数ではないので特異点を回る方法として欲しいものではありません。

以上から、複素k平面上の\(k=0\)にある特異点の回り方を指定して、ヘヴィサイド関数
\(
\displaystyle \int_{-\infty}^{\infty}\theta(x)e^{- ikx} dx = \lim_{\varepsilon \to +0 }\frac{1}{i}\frac{1}{k-i\varepsilon}
\)

と導くことが出来ます。良く、\(\varepsilon\to 0\)の表記を省略して
\(
\displaystyle \int_{-\infty}^{\infty}\theta(x)e^{- ikx} dx = \frac{1}{i}\frac{1}{k-i0}
\)

と記述されます。

演算子としての表現


ヘヴィサイド関数のフーリエ変換後の結果
\(
\displaystyle \int_{-\infty}^{\infty}\theta(x)e^{- ikx} dx = \frac{1}{i}\frac{1}{k-i0}
\)

の右辺に現れる特徴的な量
\(
\displaystyle \frac{1}{k-i0}
\)

を演算子として捉えた場合を考えます。この時、\(k\)を\(x\)に入れ替えて、ある関数\(f(x)\)に作用させると
\(
\begin{align}
\int_{-\infty}^{\infty} \left(\frac{1}{x-i0}\right) f(x)dx&=
\lim_{\varepsilon\to +0} \int_{-\infty}^{\infty} \frac{f(x)}{x-i\varepsilon} dx\nonumber \\
&=\lim_{R\to \infty}
\left[\int_{-R}^{-\varepsilon’} \frac{f(x)}{x} dx+\int_{\Gamma_-} \frac{f(x)}{x}dx+\int_{\varepsilon’}^{R} \frac{f(x)}{x}dx\right]
\end{align}
\)
となります。

ここで、\(\Gamma_-\)は特異点を下に回る経路です。
注記しますが、コーシーの積分定理があるため、この段階で\(\varepsilon’\)をゼロに近づける必要は全くありません。

特に\(\lim_{\varepsilon’\to +0}\)を考えますと、コーシーの主値を用いて
\(
\begin{align}
&=\mathcal{P}\int_{-\infty}^{\infty}\frac{f(x)}{x} dx +\pi i f(0)\nonumber \\
&=\mathcal{P}\int_{-\infty}^{\infty}\frac{f(x)}{x} dx +\pi i \int_{-\infty}^{\infty}f(x)\delta(x)dx \nonumber \\
&=\int_{-\infty}^{\infty}\left(\mathcal{P} \frac{1}{x}+\pi i \delta(x)\right) f(x)dx \nonumber
\end{align}
\)
と書き表せます。以上から、右辺に現れる特徴的な量は
\(
\begin{align}
\frac{1}{x-i0}=\mathcal{P} \frac{1}{x}+\pi i \delta(x)
\end{align}
\)

と書き換えて解釈しても良い、となります。

実軸上の積分とコーシーの主値積分

結論


実軸上の積分とは、複素平面からの近づき方を指定して、
\(
\begin{align}
\int_{-\infty}^{\infty}f(x)dx=
\frac{1}{2}\left[\lim_{\varepsilon\to +0}\int_{-\infty+i\varepsilon}^{\infty+i\varepsilon}f(x)dx+
\lim_{\varepsilon\to +0}\int_{-\infty-i\varepsilon}^{\infty-i\varepsilon}f(x)dx\right]
\end{align}
\)

と求められるようです
このように表現されていると考えなければ、\(\int_{-\infty}^{\infty} e^{ikx}/x dx\)の解が得られません。
図にするとこういう極限を取る、という意味です。

特に、実軸上\(x=0\)に1位の特異点が1つ存在する場合、
\(
\begin{align}
\int_{-\infty}^{\infty}\frac{f(x)}{x}dx=
\int_{-\infty}^{-\varepsilon}\frac{f(x)}{x}dx+
\frac{1}{2}\left[
\int_{\Gamma_+}\frac{f(x)}{x}dx+\int_{\Gamma_-}\frac{f(x)}{x}dx
\right]+
\int_{\varepsilon}^{\infty}\frac{f(x)}{x}dx
\end{align}
\)

と書かれます。ここで、\(\varepsilon\to +0\)である必要はありません

\(\Gamma_+\)は実軸上の点\(x=-\varepsilon\)から複素平面で、特異点のを通って実軸上\(x=\varepsilon\)に戻る経路、
\(\Gamma_-\)は実軸上の点\(x=-\varepsilon\)から複素平面で、特異点のを通って実軸上\(x=\varepsilon\)に戻る経路
です。

図では

という経路です。
右図の点線の経路でも積分結果が変わりません。なぜなら、点線と実線で囲まれた領域内に特異点が無い状況を考えているです。

コーシーの主値とは、複素平面の迂回部分を消した値であり、さらに\(\varepsilon\to +0\)を課した実軸上のみの積分結果です。\(x=0\)に1位の特異点がある場合、
\(
\begin{align}
\mathcal{P}\int_{-\infty}^{\infty}\frac{f(x)}{x}dx=
\lim_{\varepsilon \to +0}\left[
\int_{-\infty}^{-\varepsilon}\frac{f(x)}{x}dx+
\int_{\varepsilon}^{\infty}\frac{f(x)}{x}dx
\right]
\end{align}
\)
と表現されるはずです。

コーシーの主値を用いて、本来の実軸上の積分を評価しようとすれば、
\(
\begin{align}
\int_{-\infty}^{\infty}\frac{f(x)}{x}dx=
\mathcal{P}\int_{-\infty}^{\infty}\frac{f(x)}{x}dx+
\frac{1}{2}\lim_{\varepsilon\to +0}\left[
\int_{\Gamma_+}\frac{f(x)}{x}dx+\int_{\Gamma_-}\frac{f(x)}{x}dx
\right]
\end{align}
\)

と書けるはずです。ここで、複素平面を回る経路は特異点を回る半径を無限小にした時の値でなければなりません。

実軸上に特異点が無い場合、積分の\(\Gamma_{+/-}\)からの寄与はゼロになります。なので、コーシーの主値積分と本来の実軸上の積分は一致します。すなわち、
\(
\displaystyle \int_{-\infty}^{\infty}f(x)dx=\mathcal{P}\int_{-\infty}^{\infty}f(x)dx
\)

が成立します。

特異点をずらす操作について


特異点が実軸上にある時、

  • 特異点を迂回する経路を通って積分
  • 特異点をずらして積分

の二通りが頻出します。結局は同じことをしています。重要なのは、特異点の上を通るか下を通るかという点のみです。
1位の特異点が実軸上に存在する時、特異点をずらす場合、ずらし方として
\(
\frac{f(x)}{x+i\varepsilon}, \frac{f(x)}{x-i\varepsilon}
\)
となります。
ここで、\(\lim \varepsilon\to 0\)を取らなければなりませんが、この操作の本質は、特異点の無い積分領域の端では実軸上にいなければならない、という要請から来るものです。

特異点周りを半径無限小で回らなければならない理由はありません。

なので、特異点をずらしてから実軸上に近づける操作は、上、もしくは下を迂回して積分することと同じです。

具体的な関数をずらした時はこんな感じ。\(\varepsilon=0.2\)にとって描画しています。

注意


これは私の解釈です。

以下の4つの計算が正しい場合、上記の解釈でなければ正しい結果が出ません。
ヘヴィサイド関数のフーリエ変換
1のフーリエ変換
\(x\)のフーリエ変換
\(1/x\)のフーリエ変換
\(1/(x^n)\)のフーリエ変換

しかし、上記の説明を行っている教科書、解説ページを見受けることはできませんでした。
なので、あくまで私が正しいと思い、この考えでないと導出が出来ず、また反例が今のところ見つかっていないだけです。

コーシーの主値に意味を見出すことができません。
実際に現れる値は(実軸上の経路)+(特異点周りをまわる経路)であるからです。

コーシーの主値でなければならない理由がありましたら、是非ともお知らせいただけると幸いです。
特異点の周りを半径+0で回らなくていいことは、具体的に複素平面上の数値積分を行って確かめています。

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

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

固定端条件の下で時間依存しない1次元シュレーディンガー方程式
\(\displaystyle
\left[-\frac{1}{2}\frac{d^2}{dx^2}+V(x)\right]\psi(x)=E\psi(x)
\)

をsin基底で展開し、変分原理に基づいて対角化します。

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

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

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

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

tar

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

個別のf90ファイル

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

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

ifort -mkl quadpack.f90 main.f90

等です。

デフォルトでは、時間依存しないシュレーディンガー方程式
\(\displaystyle
\left[-\frac{1}{2}\frac{d^2}{dx^2}+\frac{1}{2}x^2\right]\psi(x)=E\psi(x)
\)

を解きます。計算に用いているパラメータは、
区間\([-20:20]\),
基底の数\(80\),
で計算します。

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

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

解く問題

電荷\(+Ze\)の原子核の周りに存在する、電荷\(-e\)を持つ電子が満たす時間依存しないシュレーディンガー方程式は、
\(
\displaystyle \left[-\frac{1}{2}\Delta-\frac{Z}{r}\right]\psi({\bf r})=E\psi({\bf r})
\)

です。
\(
\left. \psi({\bf r})\right|_{|{\bf r}|\to \infty}\to 0
\)

の境界条件を満たす正規直交化された解は、
\(
\begin{align}
\psi({\bf r})&=R_{n,l}(r)\cdot Y_{l,m}(\theta,\varphi)\\
&=\sqrt{\left(\frac{2}{n}\right)^3\frac{(n-l-1)!}{2n[(n+l)!]}} e^{-Zr/n} \cdot \left(2\frac{Zr}{n}\right)^l\cdot L_{n-l-1}^{(2l+1)}\left(2\frac{Zr}{n}\right)\cdot Y_{l,m}(\theta,\varphi)
\end{align}
\)

と表されます。ここで、

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

を意味します。

動径関数\(R_{n,l}(r)\)は微分方程式
\(
\displaystyle \left[-\frac{1}{2}\left(\frac{d^2}{dr^2}+\frac{2}{r}\frac{d}{dr}\right)+\frac{l(l+1)}{2r^2}-\frac{Z}{r}\right]R_{n,l}(r)=E_n R_{n,l}(r)
\)

を満たす正規直交化された関数で、固有値は
\(
\displaystyle E_n=-\frac{1}{2n^2}
\)
です。

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

計算結果


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

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

\(n=1\)

\(n=2\)

\(n=3\)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

n=1~n=20までのファイル hatom.zip(164MB)

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

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


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

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

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

and in gnuplot, type

call "multi.plt" 2

then, you can get figure.

multi.plt

atommulti.plt

参考


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

  • 定義域

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

  • 微分方程式

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

  • 具体例

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

  • 漸化式

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

  • 直交性

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

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

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

  • 定義域

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

  • 微分方程式

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

  • 具体例

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

  • 漸化式

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

  • 直交性

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

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

最速・高精度の数値積分

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

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

まとめ


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

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

もくじ

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

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


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

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

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

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

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

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

[adsense1]

ライセンスについて


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

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

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

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

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

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

NUMPAC

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

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

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

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

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

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

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

比較方法


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

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

結果


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

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

を比較します。また、適応型の積分手法に対して、実際にどの点が評価されているかを知りたかったため、
実際に評価された点達に対し矩形波で畳み込みを行い、最大値を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のコード
https://slpr.sakura.ne.jp/qp/supplement_data/quadpack_dqag.f90

プログラムのサンプルコード

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

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

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

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

  return
end function g

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

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

コンパイルと実行は

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

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

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

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

1次元複素積分


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

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

Quadpackのコード
https://slpr.sakura.ne.jp/qp/supplement_data/quadpack_cqag.f90

プログラムのサンプルコード

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

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

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

  g=1d0/sqrt(z)

  return
end function g

2次元実積分


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

Quadpackのコード
https://slpr.sakura.ne.jp/qp/supplement_data/quadpack_dqag2d.f90

プログラムのサンプルコード

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

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

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

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

2次元複素積分


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

Quadpackのコード
https://slpr.sakura.ne.jp/qp/supplement_data/quadpack_cqag2d.f90

プログラムのサンプルコード

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

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

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

  return
end function g

3次元実積分


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

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

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

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

  write(6,*)s
 
  stop
end program main

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

[adsense2]

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


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

参考文献


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

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

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