sikino のすべての投稿

量子コンピュータを使用した素因数分解アルゴリズム

量子コンピュータを使うと素因数分解が高速にできる。

が、実際はどうやって素因数分解をするんでしょう?
そして具体的に何が高速化するんでしょう?これを知るのが本稿の目的です。

  1. 素因数分解アルゴリズム(古典的)
  2. 実際のアルゴリズム(量子的)
    1. 初期状態の準備
    2. レジスタ1の状態に量子フーリエ変換を施す
    3. レジスタ2に、レジスタ1を元にした関数を作用させる
    4. レジスタ2を観測する
    5. レジスタ1を量子フーリエ変換する
    6. レジスタ1の観測
    7. 連分数展開によって周期\(r\)を見つける
  3. 最大公約数を求めるアルゴリズム(古典)

素因数分解アルゴリズム


量子的に高速に出来るのであれば、そのアイデア自体は古典的に考えられるはずです。

素因数分解する古典的アルゴリズムは以下のように記述されます。
\(m\)を素因数分解したい数とし、\(x\)をmと互いに素な整数(\(\text{gcd}(x,m)=1\)を満たす数)と表記します。

ここで、\(\text{gcd}(a,b)\)は、\(a\)と\(b\)の最大公約数(greatest common divisor)を出力する関数です。

  1. \(a(=0,1,\cdots )\)を引数として、関数\(f_m(a)=x^a~ \text{mod}~ m\)を計算
  2. \(f_m(a)\)の周期\(r\)を見つける。実は\(f_m(a)\)は周期\(r\)の周期関数(周期\(r\)は実際に計算してみるまで分からない。数式で表せば、\(f_n(a)=f_n(a+r)\))。
  3. 周期\(r\)が判明したら、\(p=\text{gcd}(x^{r/2}+1,m),~q=\text{gcd}(x^{r/2}-1,m)\)を計算
  4. \(p,q\)は\(m\)の素因数となっている

という風に素因数分解を行うことが出来ます。

具体的に\(m=21\)を素因数分解してみましょう。\(m\)と互いに素な整数として、\(x=11\)を選びます。
実際に横軸に\(a\)をとり、関数\(f_m(a)\)を計算すると以下のようになります。

この計算結果を見ますと、周期\(r\)が6だと分かります。

周期が\(6\)が分かったので、\(p=\text{gcd}(x^{r/2}+1,m),~q=\text{gcd}(x^{r/2}-1,m)\)を計算します。すると、
\(
\begin{align}
p&=\text{gcd}(x^{r/2}+1,m)=\text{gcd}(11^{6/2}+1,21)=3 \\
q&=\text{gcd}(x^{r/2}-1,m)=\text{gcd}(11^{6/2}-1,21)=7
\end{align}
\)

と分かります。\(21\)は\(3\times 7\)と書けるので、素因数分解を行うことが出来ました。

さて、量子的に素因数分解を行う場合、上のアルゴリズムはShorのアルゴリズムと呼ばれます。
上の計算は、

  • 最大公約数を求めるアルゴリズム(Euclidの互除法, 計算量\(O(\text{log} n)\))
  • 周期\(r\)を見つけるアルゴリズム(
    古典的:\(O(\exp[(\text{log}n)^{1/3}(\text{ln}\text{ln}n)^{2/3}])\),
    量子的:\(O((\text{log}n)^{2}(\text{log}\text{log}n)(\text{log}\text{log}\text{log}n))\)

です。最も時間が掛かる部分は、周期\(r\)を見つける部分で、これを見つけるために量子コンピュータを用いるのです。Shorのアルゴリズムは、この部分だけ量子コンピュータを使うのです。

量子的に高速に計算する、という部分は、
\(f_m(a)\)の周期\(r\)を見つける
という部分です。

実際のアルゴリズム


ここで紹介するモデルは、量子ビット(キュービット)を出さないで説明するものです。

キュービットの考えはShorのアルゴリズムの本質ではなく、実現できるか?の部分で重要です。
すなわち、ある状態を指定するために
・1つの系で、非常に高い量子番号に属する状態まで制御する(\(1[つの系]\times N[量子状態]\))
ではなく、
・複数の系で、基底状態と1つの励起状態を制御する(\(Q[つの系]\times 2[量子状態]\))
の方が現実で実現しやすいという考えです。

1. 初期状態の準備



2つの系(2つの記憶素子、レジスタ)を用意し、初期状態として各々の基底状態を表します。
この時点で、2つの系は独立に存在するもので、系同士に相互作用や、もつれ等は存在しません。

それぞれの系を添え字\(1,2\)で表すと、系全体の状態\(|\psi\rangle\)は(直積の記号を省略して)
\(
|\psi\rangle=|0\rangle_1|0\rangle_2
\)

と書き表せます。
ここで、\(|a\rangle\)は量子数\(a\)に属する固有状態を表しています。

図示すると、上の図のようになり、左はレジスタ1, 右はレジスタ2の量子状態を表していて、それぞれ、基底状態であることを表しています。

2. レジスタ1の状態に量子フーリエ変換を施す


レジスタ1に量子フーリエ変換を施すことで、レジスタ1の全状態の確率振幅を等確率にします。
量子フーリエ変換は量子状態に対して施す変換で、

\(
\begin{align}
\mathcal{F}|a\rangle&=\frac{1}{\sqrt{N}}\sum_{b=0}^{N-1}e^{-2\pi i a b/N}|b\rangle \\
\mathcal{F}^{-1}|a\rangle&=\frac{1}{\sqrt{N}}\sum_{a=0}^{N-1}e^{2\pi i a b/N}|a\rangle
\end{align}
\)

で定義します。基底状態\(|0\rangle\)に対して量子フーリエ変換を施すと、考えたい全励起状態を等確率にすることが出来ます。
\(
\begin{align}
\mathcal{F}|0\rangle&=\frac{1}{\sqrt{N}}\sum_{b=0}^{N-1}e^{-2\pi i 0 b/N}|b\rangle \\
&=\frac{1}{\sqrt{N}}\sum_{b=0}^{N-1}|b\rangle
\end{align}
\)

量子フーリエ変換を作用させると、上図のようになります。

3. レジスタ2に、レジスタ1を元にした関数を作用させる


レジスタ2に、レジスタ1を参照して計算した結果をレジスタ2に出力する、という演算を行います。
この操作を行うことで、レジスタ1と2をもつれさせます(もつれ、という言葉が正しいかは分かりません。ですが、言いたいことはレジスタ1と2を結びつける、という操作を行います)。

この演算を行う事が出来る装置があるかはわかりませんが、あると仮定します。この ”あるかどうか分からない演算を行う装置” はオラクル(神託装置)といいます。

このオラクルを用いて、レジスタ1の量子数\(a\)の値を引数として、レジスタ2の量子状態を\(|f_m(a)\rangle\)にしていきます。

この作用を\(\hat{U}_{1\to2}\)という演算子で表すと、レジスタ1の結果をレジスタ2に格納するので、

\(
\begin{align}
\hat{U}_{1\to2} \hat{\mathcal{F}}_1|0\rangle_1|0\rangle_2
&=\frac{1}{\sqrt{N}}\sum_{a=0}^{N-1}\hat{U}_{1\to2}|a\rangle_1|0\rangle_2 \\
&=\frac{1}{\sqrt{N}}\sum_{a=0}^{N-1}|a\rangle_1|x^a~ \text{mod}~ m\rangle_2
\end{align}
\)

4. レジスタ2を観測する


レジスタ2を観測し、レジスタ2の量子数を知ります。
この観測は以下の2つの影響を系に与えます。

・レジスタ2の量子状態を確定させる
・レジスタ1が取り得る量子状態に制限を加える

1回の”観測”の操作はただ1点だけを与えますが、数式上扱いづらいです。
なので、数式では何度も観測を行った時に得られる期待値を数式で表すと分かりやすくなります。
レジスタ2の状態を量子数\(k\)の状態に見出す確率を知るために、左から\(\langle k|_2\)を作用させます。
すると、以下のように式変形をすることが出来ます。

\(
\begin{align}
\langle k|\hat{U}_{1\to 2} \hat{\mathcal{F}}_1|0\rangle_1|0\rangle_2
&=\frac{1}{\sqrt{N}}\sum_{a=0}^{N-1}|a\rangle_1|x^a ~\text{mod}~ m\rangle_2 \\
&\approx\frac{1}{\sqrt{N}}\sum_{a=0}^{N-1}|a\rangle_1\langle k|x^a ~\text{mod}~ m\rangle_2 \\
&=\frac{1}{\sqrt{N}}\sum_{a=0}^{N-1}|a\rangle_1\left(\frac{1}{\sqrt{J}}\sum_{j=0}^J \delta_{a,a_j^{(k)}}\right)
\end{align}
\)
ここで、\(a_j^{(k)}\)は、
\(k=x^{a_{j}^{(k)}}~ \text{mod} ~m,~~(j=0,1,\cdots,J-1)\)を満たす\(J\)個の\(a\)です。

この観測により、レジスタ1の状態を特定の状態を\(|a_j^{(k)}\rangle,~~(j=0,1,\cdots, J-1)\)のみ存在させることになります。レジスタ1のとれる状態は周期的で周期\(r\)をもち、
\(
\displaystyle a_j^{(k)}=a_0^{(k)}+j\cdot r,~~(j=0,1,\cdots, J-1)
\)

という形で表すことが出来ます。

レジスタ1の状態を図示するとこのようになっているはずです。

係数をまとめると分かりやすくなります。
単なる式変形ですが、状態に適当な係数が掛かった状態の重ね合わせで表現されるので、なじみ深い形…だと私は思います。
\(
\begin{align}
\langle k|\hat{U}_{1\to 2} \hat{\mathcal{F}}_1|0\rangle_1|0\rangle_2
\approx=\sum_{a=0}^{N-1}\left(\frac{1}{\sqrt{NJ}}\frac{1}{\sqrt{J}}\sum_{j=0}^{J-1} \delta_{a,a_j^{(k)}}\right)|a\rangle_1
\end{align}
\)

5. レジスタ1を量子フーリエ変換する



レジスタ1の状態に対して量子フーリエ変換を行うことで、レジスタ1の持つ周期を知ることが出来ます。

量子フーリエ変換を行う前は量子の属する状態そのものに周期\(r\)の情報は含まれていませんが、量子フーリエ変換を行った後は量子数そのものに周期\(r\)に関係する値が含まれるようになります。

実際に観測を行うと、

\(
\begin{align}
\hat{\mathcal{F}}^{-1}_1\langle k|\hat{U}_{1\to 2} \hat{\mathcal{F}}_1|0\rangle_1|0\rangle_2
&\approx\sum_{a=0}^{N-1}\left(\frac{1}{\sqrt{NJ}}\frac{1}{\sqrt{J}}\sum_{j=0}^J \delta_{a,a_j^{(k)}}\right)\hat{\mathcal{F}}^{-1}_1|a\rangle_1 \\
&=\sum_{b=0}^{N-1}\left(\frac{1}{N\sqrt{J}}\sum_{j=0}^{J-1} e^{2\pi i a_j^{(k)}b/N}\right)|b\rangle_1 \\
&=\sum_{b=0}^{N-1}\left[\frac{1}{N\sqrt{J}}e^{2\pi i a_0^{(k)}b/N}\sum_{j=0}^{J-1} e^{2\pi i jrb/N}\right]|b\rangle_1
\end{align}
\)

レジスタ1の量子数\(b’\)に見出す確率は、左から\(\langle b’|_1\)を作用させて絶対値の2乗をとると分かります。

まず、左から\(\langle b’|_1\)を作用させると、
\(
\begin{align}
C'(b’)=\langle b’|_1\hat{\mathcal{F}}^{-1}_1\langle k|\hat{U}_{1\to 2} \hat{\mathcal{F}}_1|0\rangle_1|0\rangle_2
&\approx\frac{1}{N\sqrt{J}}e^{2\pi i a_0^{(k)}b’/N}\sum_{j=0}^{J-1} e^{2\pi i jrb’/N}
\end{align}
\)

となります。確率は上記\(C'(b)\)の値の絶対値の2乗で与えられるので(表記を\(b’\to b\)に変更)、
\(
\begin{align}
|C'(b)|^2&=\left|\frac{1}{N\sqrt{J}}e^{2\pi i a_0^{(k)}b/N}\sum_{j=0}^{J-1} e^{2\pi i jrb/N}\right|^2 \\
&=\frac{1}{N^2 J}\left|\sum_{j=0}^{J-1} e^{2\pi i jrb/N}\right|^2
\end{align}
\)
となります。確率密度は
\(
\displaystyle \left|\sum_{j=0}^{J-1} e^{i(2\pi jrb/N)}\right|
\)
に依存します。この形は離散フーリエ変換の場合に似ています。
つまり、この項が大きくなる時の状態がレジスタ1の量子状態として観測されることが期待され、
観測される状態とは、\(j\)が変化しても位相が一致しているとき、すなわち

\(2\pi rb/N=2\pi \times (整数)\)

を満たすような\(r\)の時に存在確率密度が大きくなる、ということになります。
整数をsと表記すると、

\(rb/N= s~\to~ b/N= s/r\)

と表されます。

6. レジスタ1の観測



確率密度分布は\(b=a_0^{(k)}+j\cdot r\)の量子状態を中心に分布しますが、あくまで、確率分布ですので、観測を行った場合は若干ずれた所に見出したりします。

実際にレジスタ1の観測を行った際に得られた量子数を\(A\)とします。すると、関係式

\(\frac{A}{N}\approx \frac{s}{r}\)

が近似的に成立していると考えることが出来ます。
今、左辺は完全に分かっている量ですが、右辺の分子・分母は整数を持つ、ということくらいしかわかりません。これを求めるために連分数展開を用います。

7. 連分数展開によって周期\(r\)を見つける


レジスタ1の観測によって得られた状態の量子数を\(A\)とします。今、
\(\displaystyle
\frac{A}{N}=\frac{s}{r}
\)

ここで、\(A\)は観測された量子数、
\(N\)はレジスタの全量子状態数、
\(s (\lt N)\)は任意の整数、
\(r\)は求めたい周期です。

これから行いたい操作は、\(\frac{A}{N}\)の分子分母を出来るだけ小さい数にする、という作業を行いたいのです。
これは単純に割るだけで済む話ではありません。なぜなら、\(\frac{A}{N}\)は綺麗に割り切れる数ではない可能性があるからです。

例えば、\(k=184, N=243\)という値であった場合、直感的に\(\frac{184}{243}\approx \frac{3}{4}\)だということが分かります。

これを機械的に行うためには連分数展開を利用します。
連分数展開とは、ある実数\(x\)を整数\(a_n\)を用いて以下のように展開することです。
\(\displaystyle
x=a_0+\frac{1}{a_1+\frac{1}{a_2+\frac{1}{\cdots}}}
\)

ここで、\(a_n~(n=0,1\cdots)\)は、
\(
\begin{align}
&a_0=\lfloor x \rfloor,~~b_0=\frac{1}{x-a_0} \\
&a_n=\lfloor b_{n-1} \rfloor,~~b_n=\frac{1}{b_{n-1}-a_n}
\end{align}
\)
と得られます。ここで\(\lfloor x \rfloor\)はガウスの記号で、実数の整数部分を表しています。
\(x\)が有理数であれば、連分数展開は有限の項数で終わります。

連分数展開後、\(x\)の分数による近似は
\(\displaystyle
x\approx\frac{d_n}{r_n}
\)

ここで、
\(
\begin{align}
& d_0=a_0,~d_1=1+a_0a_1,~d_n=a_n d_{n-1}+d_{n-2},\\
& r_0=1,~r_1=a_1,~r_n=a_n r_{n-1}+r_{n-2},
\end{align}
\)

として順々に求める事が出来ます。
また、連分数の打切りに関しては定理

「\(\frac{q}{p}\)を\(\left|\frac{q}{p}-x\right|\lt \frac{1}{2q^2}\)を満たす任意の有理数とする時、\(\frac{q}{p}\)は\(x\)の近似分数である。さらに、その近似分数は\(p,q\)の最大公約数は1である。」

を用いて、条件式\(\left|\frac{q}{p}-x\right|\lt \frac{1}{2q^2}\)から決めます。

連分数展開によって分子・分母を出力するプログラムはこんな感じで実装できます。

最大公約数を求めるアルゴリズム(古典)


最大公約数を求めるアルゴリズムはEuclidのアルゴリズムと呼ばれ、他のアルゴリズムと比べても早く終わります。そのプログラムは以下で実装できます。

function gcd(a,b)
  implicit none
  integer,intent(in)::a,b
  integer::gcd
  !compute the greatest common divisor
  integer::x,r
 
  x=a
  gcd=b
  r=mod(x,gcd)
  do while( r .gt. 0)
     x=gcd
     gcd=r
     r=mod(x,gcd)
  enddo
 
  return
end function gcd

参考文献

Elisa Bäumer, Jan-Grimo Sobez, Stefan Tessarini, Shor’s Algorithm
https://qudev.phys.ethz.ch/content/QSIT15/Shors%20Algorithm.pdf

C. P. ウィリアムズ、S. H. クリアウォータ著、西野哲郎、新井隆、渡邉昇訳「量子コンピューティング」(springer, 2000)

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

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

その他フーリエ変換
ヘヴィサイド関数のフーリエ変換
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}
\)
としてデルタ関数が表現されていることが分かります。