sikino のすべての投稿

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)\)は符号関数


実軸上の積分は複素平面の上から近づいた値と下から近づいた値の平均をとって定めるとします。すなわち、
\(
\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のフーリエ変換結果
\(
\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のフーリエ変換


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}
\)

で定義しています。

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


ヘヴィサイド関数\(\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}H(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}
\)

と求められるようです
図にするとこういう極限を取る、という意味です。

特に、実軸上\(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つの計算
(sinc関数の無限積分、デルタ関数のフーリエ変換、1/xのフーリエ変換、ヘヴィサイド関数のフーリエ変換)
が正しい場合、上記の解釈でなければ矛盾が出ます。

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

また、僕の立場として、コーシーの主値は全くの意味が無いものと考えています。

実際に現れる値は(実軸上の経路)+(特異点周りをまわる経路)であるからです。

反例がありましたらお知らせください。

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

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

の、Gauss-Lobatto求積法を利用した変分法よりも計算速度が遅いですが、確実に積分の計算が可能です。
ここで紹介するコードには並列計算は実装していません

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

必要なファイルは以下の2つです。
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

等です。
デフォルトでは、時間依存しないシュレーディンガー方程式
\(
\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プログラム
  6. ベンチマーク時のプログラム
  7. 参考文献

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


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

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

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

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

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

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


スポンサーリンク


ライセンスについて


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

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

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

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

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

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

NUMPAC

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

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

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

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

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

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

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

比較方法


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

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

結果


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

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

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

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


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


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


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


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

ちなみに


二重指数関数型数値積分が有効な例として
\(
\displaystyle \int_0^{1}\sin(1/\sqrt{x})/\sqrt{x}
\)

という例題が挙げられます。難しい事には変わりませんが二重指数関数型数値積分だけが出来るのではないことをここに注記しておきます。

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

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


まとめです。

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

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

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

QUADPACKプログラム


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

Quadpackのコード
http://slpr.sakura.ne.jp/qp/supplement_data/quadpack_dqag.f90

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

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

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

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

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

  return
end function g

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

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

コンパイルと実行は

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

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

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

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

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

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

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

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

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

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

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

  g=sqrt(z)
 
  return
end function g
スポンサーリンク

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


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

参考文献


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

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

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

フリーソフトの紹介

全てwindows10上で動きます。
Microsoft Powerpoint, Word以外は無料です。

  1. キャプチャ -Setuna2
  2. ランチャー -RocketDuck
  3. メール管理 -Thuderbird
  4. テキストエディタ -Sublime Text
  5. 文書作成 -Microsoft Word, -LaTeX
  6. グラフ作成 -gnuplot
  7. PDFビューワ- SumatraPDF
  8. 画像圧縮 -Caesium portable
  9. 画像編集 -Microsoft Powerpoint, -GIMP
  10. GIF画像作成、編集 -Giam
  11. 画像ビューワ -Leeyes
  12. 音楽再生 -VLCメディアプレイヤー
  13. 音楽ダウンロード -Download helper(firefoxアドオン)
  14. 動画再生 -GOMplayer
  15. 動画変換、音楽の抽出 -Craving Explorer
  16. mp3の音量調節 -MP3Gain
  17. 動画キャプチャ -oCam
  18. 動画編集 -Aviutl
  19. PCの情報取得 -Speccy
  20. PC間のデータ移動 -IP Messenger
  21. PCのデータ削除 -Bleachbit portable
  22. PCのウイルス除去 -Norton Power Eraser
  23. アプリケーション情報 -Process Explorer
  24. 広告除去 -Adblock Plus
  25. 仮想環境 -Vmware Player
  26. オフラインゲーム -Track Mania Forever

1. キャプチャ


Setuna2
http://www.clearunit.com/clearup/setuna2/

説明

画面上の任意の範囲を画像として保存するソフトウェアです。
画像保存だけでなく、キャプチャした領域は半透明で常に画面最前列に現れます。
そのため、アプリケーション間を行ったり来たりすることなく欲しい情報を表示させておけます。
現在サイトが閉鎖されている…ようです?


2. ランチャー


RocketDuck
https://rocketdock.com/

説明

ランチャーです。Macにある感じのやつです。画像の上に出てる枠を追加します。
良く使う実行ファイルやフォルダへのショートカットをデスクトップ上以外に表示させることが出来ます。


3. メール管理


Thuderbird
https://www.mozilla.org/ja/thunderbird/

説明

複数のメールアドレス(例えばyahoo,gmail,hotmailなど他社サービスであっても)をまとめて管理できます。
インターネットブラウザを立ち上げないでメールを見られます。


4. テキストエディタ


Sublime Text
https://www.sublimetext.com/3

説明

プログラムを書くためのエディタです。


5. 文書作成


Microsoft Word, LaTeX

説明

文章を書くためのソフトウェアです。
数式や記号のようなものを表現したければLaTeXが良いです。
LaTeXは導入までが難しいです。導入はEasyTeXが比較的簡単でしょう。
EasyTeX


6. グラフ作成


gnuplot
http://www.gnuplot.info/

説明

グラフを描く為のソフトウェアです。
windows用に作られたwxtターミナルはアンチエイリアスが聞いているので、滑らかなグラフが表示されます。
eps,gif画像も作ることが出来ます。


7. PDFビューワ


SumatraPDF
https://www.sumatrapdfreader.org/free-pdf-reader.html

説明

PDFファイルを開くためのソフトウェアです。
PDFファイルを見るだけしか出来ませんが、起動が高速で、動作中に固まったことはありません。
PDFビューワとして有名なソフトウェアAdobe Readerがありますが、起動までに数秒~数十秒かかり動作も重いです。


スポンサーリンク



8. 画像圧縮


Caesium portable
https://saerasoft.com/caesium/

説明

jpg,png等の画像を圧縮し,容量を減少させます。
インストールなしに使うことができ、手軽に使えます。


9. 画像編集


Microsoft Powerpoint
GIMP
https://www.gimp.org/

説明

・Powerpoint
本来プレゼンテーションを作るためのソフトです。
画像や写真に文字を加えたい時などに、ppt上で編集→図として保存を行うと比較的簡単に編集できます。

・GIMP
画像編集用ソフトです。
フォトショップ並みの編集が可能ですが、慣れ、知識、技術が必要ですので、初心者向けではありません。


10. GIF画像作成、編集


Giam
http://furumizo.net/tsu/giamd.htm

説明

GIF画像の編集ソフトです。
GIF画像を作ったり、編集を加えたり、コマ送りとして見ることが出来るソフトです。
ある特定の場面を画像として保存することも出来ます。


11. 画像ビューワ


Leeyes
http://www3.tokai.or.jp/boxes/leeyes/

説明

画像、または圧縮ファイルにおさめられた画像を解凍することなく見ることが出来る画像ビューワです。
プラグインを入れることによりrarなどの様々な画像形式に対応することが出来ます。


12. 音楽再生


VLCメディアプレイヤー
https://www.videolan.org/vlc/index.ja.html

説明

音楽、動画再生ソフトです。


13. 音楽,動画ダウンロード


Download helper(firefoxアドオン)
https://addons.mozilla.org/ja/firefox/addon/video-downloadhelper/

説明

firefoxのアドオンで、あらゆるサイトの動画のダウンロードが可能です。
3つの風船がトレードマークで、動画配信サイトに行きこの風船をクリックし、動画ダウンロードを行います。
動画のダウンロードは、著作権法に抵触しないようにしましょう。
※ブラウザに付随するコンソールでも、メディアファイルのダウンロードが可能です。


14. 動画再生


GOMplayer
http://www.gomplayer.jp/

説明

動画再生ソフトです。
コマ送り、スロー再生、逆再生、左右反転、上下反転が出来ます。
広告が邪魔ですが、これは環境設定→一般→スキンからクラシック[フォルダ]を選ぶと広告は消えます。


15. 動画変換、音楽の抽出


Craving Explorer
https://www.crav-ing.com/

説明

本来はインターネットブラウザです。
Craving explorerに付属する、オフライン状況下でも使える動画変換機能がシンプルで優秀です。
以前は動画、音楽ダウンロードを売りにしていたブラウザでしたが、著作権の問題からダウンロード機能が除去されていたようです。
2017年現在、ホームページを訪れると再び動画をダウンロードできる仕様になったようです。


16. mp3の音量調節


MP3Gain
http://mp3gain.sourceforge.net/

説明

mp3ファイルの音量を変えるソフトです。
複数ファイルの音量の均一化をすることも出来ます。
エンコードするわけではないので、操作は一瞬で完了します。


17. 動画キャプチャ


oCam
http://ohsoft.net/eng/

説明

画面を動画としてキャプチャします。
録画後に広告が出ますが、録画した動画に広告は出てきません。


18. 動画編集


Aviutl
http://spring-fragrance.mints.ne.jp/aviutl/

説明

avi動画の編集を行えます。
全く別の2つの動画から比較動画を作ったり、文字を入力したりできます。
初心者向けではありません。編集には慣れが必要です。
プラグインを用いれば、mp4を読み込ませ、編集し、mp4で出力させることが出来ます。
詳しい説明は
VIPで初心者がゲーム実況するには https://www18.atwiki.jp/live2ch/pages/165.html

ニコニコ動画まとめwiki http://nicowiki.com/AviUtl.html
をご覧ください。

スポンサーリンク


19. PCの情報取得


Speccy
https://www.piriform.com/speccy

説明

PCの構成情報を綺麗なウインドウで教えてくれるだけのソフトウェアです。
例えば、OS, CPU, メモリ(DDR3など)、モニターの型番、マウスの情報、IPアドレス等々です。
スマートです。


20. PC間のデータ移動


IP Messenger
https://ipmsg.org/

説明

ルーターを介してLANでつながったPC間のデータ移動を行います。
インターネット環境は必要ありません。ネット接続もしません。
高速、セキュリティともに良いと思います。開発から20年経過した現在(2017/09/08)でも更新がなされています。


21. PCのデータ削除


Bleachbit portable
https://www.bleachbit.org/

説明

PCに存在するゴミファイルを消去して空き容量を確保します。
ブラウザのキャッシュ、一時ファイル、データを揃えて空き容量を確保したりできます。
アプリケーションごとにどういったファイルを消すか指定することが出来ます。


22. PCのウイルス除去


Norton Power Eraser
https://security.symantec.com/nbrt/npe.aspx?lcid=1041

説明

感染したPCからウイルスを除去する強力なツールです。
徹底的にデータを消すため、無害なソフトを巻き込むこともあるようです。
システムの復元でも消えない、hao123の除去も出来ます。


23. アプリケーション情報


Process Explorer
https://docs.microsoft.com/en-us/sysinternals/downloads/process-explorer

説明

不審な動きをするアプリケーションを特定するために使えます。


24. 広告除去


Adblock Plus
https://addons.mozilla.org/ja/firefox/addon/adblock-plus/

Element Hiding helper for Adblock Plus
https://addons.mozilla.org/ja/firefox/addon/elemhidehelper/

説明

インターネットのサイトに表示させる広告を除去することが出来ます。
スマホ(Android)用のAdblockブラウザもありますが、動作が遅いです。
firefox、Chromeにはプラグインがあります。

Adblock Plusは画像広告に強く、
Element Hiding helper for Adblock Plusはテキスト広告に強いです。


25. 仮想環境


Vmware Player
https://www.vmware.com/products/player/playerpro-evaluation.html

説明

windows上にlinux環境を構築できます。
あくまでwindows上の仮想環境なので、linuxが壊れてしまった場合やlinuxを消したい場合も単にwindows上からフォルダを消すだけokです。


26. オフラインゲーム


Track Mania Forever
http://trackmaniaforever.com/
無料版ダウンロード http://www.4gamer.net/games/048/G004822/20080423014/

説明

オフラインで遊べるレーシングゲームです。
あるコースを時間内に走り抜けることが目的のゲームです。
コースエディタがあるので、自分で作ることが可能です。