「数学」カテゴリーアーカイブ

ニュートン法(1、2次元、多次元)

ニュートン法に関するお話です。
”ニュートン法”と呼ばれる方法は
2種類(1.初期値近傍の極小値を求めるニュートン法、2. ゼロ点を求めるためのニュートン法)
があるようですが、ここでは2. 関数のゼロ点を求める方のニュートン法についてのお話です。

  1. まとめ
  2. 1次元のニュートン法
    1. 初期値の推定
    2. 導関数の近似
  3. 1次元ニュートン法のプログラム
  4. 2次元ニュートン法
  5. 2次元ニュートン法のプログラム
  6. 多次元ニュートン法
  7. 高次のニュートン法?
  8. 補)ニュートン法が使われる問題
  9. 参考文献

まとめ

1次元の場合

\(\displaystyle
x_{n+1}=x_{n}-\frac{f(x_{n})}{f'(x_{n})},~~n=0,1,2\cdots,
\)

ここで、\(x_0\)は解の近傍であること。

2次元の場合

\(
\begin{align}
\left(
\begin{array}{c}
x_{n+1} \\
y_{n+1}
\end{array}
\right)
=
\left(
\begin{array}{c}
x_n \\
y_n
\end{array}
\right)
-\left(
\begin{array}{cc}
f_x(x_n,y_n) & f_y(x_n,y_n) \\
g_x(x_n,y_n) & g_y(x_n,y_n)
\end{array}
\right)^{-1}
\left(
\begin{array}{c}
f(x_n,y_n) \\
g(x_n,y_n)
\end{array}
\right)
\end{align}
\)

\(
\begin{align}
&\left(
\begin{array}{cc}
f_x(x_n,y_n) & f_y(x_n,y_n) \\
g_x(x_n,y_n) & g_y(x_n,y_n)
\end{array}
\right)^{-1} \\
&~~~~=
\frac{1}{f_x(x_n,y_n) g_y(x_n,y_n)-f_y(x_n,y_n)g_x(x_n,y_n)}
\left(
\begin{array}{cc}
g_y(x_n,y_n) & -f_y(x_n,y_n) \\
-g_x(x_n,y_n) & f_x(x_n,y_n)
\end{array}
\right)
\end{align}
\)

ここで、
\(
\begin{align}
f_x(x_n,y_n) &= \left.\frac{\partial f(x,y)}{\partial x}\right|_{x=x_n,~y=y_n} ,~~& f_y(x_n,y_n) = \left.\frac{\partial f(x,y)}{\partial y}\right|_{x=x_n,~y=y_n} \\
g_x(x_n,y_n) &= \left.\frac{\partial g(x,y)}{\partial x}\right|_{x=x_n,~y=y_n} ,~~& g_y(x_n,y_n) = \left.\frac{\partial g(x,y)}{\partial y}\right|_{x=x_n,~y=y_n}
\end{align}
\)

多次元の場合

\(
\mathbf{x}_{n+1}=\mathbf{x}_{n} – \mathbf{J}_n^{-1}\mathbf{f}_n
\)

ここで、
\(
\begin{align}
\mathbf{x}_n=
\left(
\begin{array}{c}
x_1\\
x_2\\
\vdots \\
x_N
\end{array}
\right)_n
,~~
\mathbf{f}_n=
\left(
\begin{array}{c}
f_1(\mathbf{x}_n)\\
f_2(\mathbf{x}_n)\\
\vdots \\
f_N(\mathbf{x}_n)
\end{array}
\right)
,~~
\mathbf{J}_n=
\left.\left(
\begin{array}{cccc}
\frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_N} \\
\frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \cdots & \frac{\partial f_2}{\partial x_N} \\
\vdots&\vdots&\ddots&\vdots \\
\frac{\partial f_N}{\partial x_1} & \frac{\partial f_N}{\partial x_2} & \cdots &\frac{\partial f_N}{\partial x_1}
\end{array}
\right)\right|_{\mathbf{x}=\mathbf{x}_n}
\end{align}
\)

※\((~~)_n\)の添え字nはn回目の繰り返しにおけるベクトル\(\mathbf{x}\)を表しています。


スポンサーリンク


ここから本文

1次元のニュートン法


方程式
\(
f(x)=0
\)

を満たす\(x=a\)を求めることを考えます。

まず、この問題をグラフ
\(
y=f(x)
\)

のゼロ点(\(y=0\)となるような\(x=a\))を探す問題だ、という問いに置き換えます。

これを、解近傍で関数をテーラー展開で表し、その展開式のゼロ点を探すことで解を近似することで方程式を解きます。

関数\(f(x)\)が与えられたとき、任意の点\(x’\)周りにおけるテーラー展開は、
\(
f(x) = f(x’)+f'(x’)(x-x’) + O(\Delta^2),~~\Delta=|x-x’|
\)

と書くことが出来ます。もしも\(x\)と\(x’\)が近接していれば、\(\delta=|x-x’| \ll 1\)となるので、\(\Delta^2\)の項は無視できるほど小さくなると期待します。

さて、この式を導く際に用いた条件は、\(x\)と\(x’\)が近くにある、という条件のみです。なので、求めたい解\(a\)と\(a\)の近くの適当な点\(x_0\)を用いれば、
\(
\begin{align}
f(a) &= f(x_0)+f'(x_0)(a-x_0) +O(\Delta^2)~~~…(1)\\
f(x_0) &= f(a)+f'(a)(x_0-a) +O(\Delta^2) ~~~…(2)
\end{align}
\)

の2通りの式を考えることが出来ます。(2)に含まれる\(f'(a)\)は本当の\(a\)が分からないと評価するのが難しいので、(1)を選びます。

(1)の左辺は条件よりゼロです。残りの右辺の\(a\)について解けば、

\(
\begin{align}
0 = f(x_0)+f'(x_0)(a-x_0)+O(\Delta^2) \\
\to a = x_0-\frac{f(x_0)}{f'(x_0)}+O(\Delta^2)
\end{align}
\)


ここで、\(f(x_0)=f'(x_0)=O(\Delta^0)\)を想定しているため、\(f'(x_0)\)の割り算を行っても誤差のオーダーは変化せず、\(\Delta\)に対して2乗のままとなります。

\(O(\Delta^2)\)を無視すると、得られる解は近似値となります。この近似値を\(x_1\)と置くと
\(\displaystyle
x_1 = x_0-\frac{f(x_0)}{f'(x_0)}
\)

得られる解\(x_1\)は\(x_0\)よりも解\(a\)に近いため、\(O(\Delta^2)\)はさらに小さくなります。なので、もう一度計算します。この時の近似値を\(x_2\)と置くと
\(\displaystyle
x_2 = x_1-\frac{f(x_1)}{f'(x_1)}
\)

ここで得られた\(x_2\)は更に\(a\)に近いため、・・・と繰り返していきます。

すなわち、解\(a\)に近い初期値\(x_0\)を与えた時、漸化式
\(\displaystyle
x_{n+1}=x_{n}-\frac{f(x_{n})}{f'(x_{n})},~~n=0,1,2\cdots,
\)

を繰り返すことで解\(a\)に近づいていく、ということです。
この式から分かる通り、解近傍で導関数がゼロになる場合はゼロ割発生するかもしれないので危険です。

もしもこのステップが無限回繰り返されれば、無限回目のステップの近似値\(x_\infty\)は
\(
x_\infty=a
\)

に収束します。

以上のように導関数を用いて方程式の解を求める方法は、ニュートン法(またはニュートン=ラフソン法)と呼ばれます。
この方法は、解に近いに所の値とその導関数が分かりさえすればいいので、この方法は1次元、多次元でも使用する事ができ、また複素関数でさえ使用することが出来ます。
さらに、誤差が\(O(\Delta^2)\)で減衰していくため、一度繰り返すごとに厳密解に一致する桁数が2倍になっていきます。

しかし問題があり、
1. 初期値を解の近傍に取らないと失敗する
2. 導関数が分からないと使えない
という大きな問題を抱えています。

また、解近傍でテーラー展開可能でなければならず、解近傍で導関数がゼロになる点があってはいけません。

  • 初期値の推定


    ニュートン法は解の近傍で、テーラー展開可能な領域から始めなければ用いることは出来ません。
    これは工夫をするしかありません。

    解近傍の初期値を推定する方法として、以下の2つの方法が良く取られます。
    1. 解の範囲を限定し、その近傍からニュートン法を行う(1次元のみ)。
    2. 適当な項を落とせば方程式が解けてしまう場合、その解を出発点とし、残りの項を摂動的に加えながら逐次解く。

    1. 解の範囲を限定する

      これは1次元の場合に使える方法です。
      単純な考え方で、例えば広い区間\([x_a,x_b]\)にゼロ点があることだけが分かっているとします。
      この区間を\(N+1\)分割し、
      \(
      x^{(0)}=x_a,x^{(1)},\cdots,x^{(k)},x^{(k+1)},\cdots, x^{(N-1)},x^{(N)}=x_b
      \)

      とします。もしも\(x^{(k)}\)と\(x^{(k+1)}\)の符号が違う場合、この区間でゼロを横切っているということになりますので、この点からニュートン法を始めれば良いのです。

      問題は、解の範囲を狭めるために関数を\(N+1\)回評価しなければならない、という点で計算コストがかかってしまう、という点です。

    2. 既知の解から摂動的に逐次解く

      これは若干特殊な方法で、あてずっぽうな上の考えとは違います。
      本来解きたい式を
      \(
      f(x)=f_0(x)+g(x)
      \)

      と分けます。ここで、
      \(
      f_0(x)=0
      \)

      の解は簡単に解けて、\(x=a_0\)が関数\(f_0(x)\)のゼロ点であるとします。この\(a_0\)を初期値とし、\(\lambda\)を十分小さい値にとり
      \(
      f_0(x)+\lambda_1 g(x)=0
      \)

      を解きます。ここで、\(\lambda_m\)は\(0\lt\lambda_1\lt\lambda_2\cdots \lt\lambda_{M-1} \lt \lambda_M=1\)を満たすとします。
      この方程式の解\(a_1\)が得られたら、この解を初期値とし、
      \(
      f_0(x)+\lambda_2 g(x)=0
      \)

      を解きに行き、これを\(\lambda_M=1\)まで続けます。
      こうすることにより、最後の答え\(a_M\)がもともと解きたかった式の答えとなっているのです。

      関数の一回の評価に時間がかかってしまう場合、この解き方の方が早く解けます。
      しかし、関数\(g(x)\)によってゼロ点の位置が増えたり、消えたりしてしまう時にこの方法は使うことが出来ず、連続的に解が動いてくれないと使うことが出来ない点に注意してください。

  • 導関数の近似


    導関数を求める簡単な方法は、差分を使うことです。すなわち、導関数\(f'(x)\)を差分に置き換えます。
    \(\displaystyle
    f'(x)\approx \frac{f(x+h)-f(x)}{h} + O(h)
    \)

    すると、漸化式は
    \(\displaystyle
    x_{n+1}=x_{n}-h\frac{f(x_{n})}{f(x_n+h)-f(x_n)},~~n=0,1,2\cdots,
    \)

    と得られます。
    この方法は導関数を必要とせず、実装が簡単であるという意味で有用な式です。

    導関数を求める際の前進差分の刻み幅\(h\)は、計算機の扱える丸め誤差\(\epsilon\)(単精度ならば約\(\epsilon=10^{-8}\)、倍精度ならば約\(\epsilon=10^{-16}\))の\(1/2\)乗の2倍、つまり
    \(
    \displaystyle h\sim 2\sqrt{\epsilon}
    \)

    とすると良いです。

    この刻み幅の導出は折りたたんでおきますが、以下のように導出することが出来ます。

参考文献[1][3][5]

1次元ニュートン法のプログラム


以下のプログラムは、1次元のニュートン法のプログラムで、
方程式
\(
x^2-4 = 0
\)

を初期値\(x_0=3\)から始めて解くプログラムです。

このコードを動かすと、ステップを繰り返すごとに急激に解に近づいていっていることが分かります。
(デフォルトではコメントアウトしてあります。)

$ gfortran main.f90
$ ./a.out
    1       2.1666666616021075
    2       2.0064102565311974
    3       2.0000102400381690
    4       2.0000000000262466
    5       2.0000000000000000
   2.0000000000000000       ---result
$

上には出力していませんが、初期値の評価もしています。
なので、ニュートン法は6回呼び出されたことになります。
ニュートン法1回当たり関数は2回呼び出されますので、12回関数が呼び出されたことになります。

例えば比較対象として、二分法を考えますと、1回目に関数は3回,そのほかで1回ずつ呼び出されます。
二分法の解は1ステップ当たり区間の半分になりますので、16桁一致するまでには50回ほど繰り返す必要があります。そのため、52回位関数が呼び出されます。

ニュートン法の方が圧倒的に早いことが分かります。

複素関数の1次元ニュートン法のプログラム
$ gfortran main.f90
$ ./a.out
                    1 ( -1.2933333304320405     , 0.72000000272341036     )
                    2 (-0.78207788228079722     , 0.60930726266833468     )
                    3 (-0.43844290526861285     , 0.73503785880633421     )
                    4 (-0.50851135934742864     , 0.89043164725052448     )
                    5 (-0.50009896196621662     , 0.86666395460253720     )
                    6 (-0.49999991066297655     , 0.86602581130981315     )
                    7 (-0.49999999999985084     , 0.86602540378453374     )
                    8 (-0.50000000000000000     , 0.86602540378443860     )
 (-0.50000000000000000     , 0.86602540378443860     )  ----result
$

ソースコードはこちら。

2次元のニュートン法


続いて方程式
\(
\begin{eqnarray}
\left\{~~
\begin{aligned}
f(x,y)&=0,\\
g(x,y)&=0
\end{aligned}
\right.
\end{eqnarray}
\)

を満たす\(x=a,y=b\)を見つけることを考えます。

1次元の時と同様に、任意の点\(x’,y’\)周りにおけるテーラー展開は、
\(
\begin{eqnarray}
\left\{~~
\begin{aligned}
f(x,y) &= f(x’,y’)+f_x(x’,y’)\cdot(x-x’)+f_y(x’,y’)\cdot(y-y’)+ O(\varepsilon^2)\\
g(x,y) &= g(x’,y’)+g_y(x’,y’)\cdot(x-x’)+g_y(x’,y’)\cdot(y-y’)+ O(\varepsilon^2)
\end{aligned}
\right.
\end{eqnarray}
\)

ここで、
\(\displaystyle \left. f_x(x’,y’)=\frac{\partial f(x,y)}{\partial x}\right|_{x=x’,y=y’}\)
\(\displaystyle \left. f_y(x’,y’)=\frac{\partial f(x,y)}{\partial y}\right|_{x=x’,y=y’}\)
\(\displaystyle \Delta x=|x-x’|=O(\varepsilon^1),~~\Delta y=|y-y’|=O(\varepsilon^1)\)
と書きました。

1次元の時と同様に初期値\((x_0,y_0)\)が解\((a,b)\)に近ければ、
\(
\begin{eqnarray}
\left\{
\begin{aligned}
0 &= f(x_0,y_0)+f_x(x_0,y_0)\cdot(a-x_0)+f_y(x_0,y_0)\cdot(b-y_0) +O(\varepsilon^2)\\
0 &= g(x_0,y_0)+g_x(x_0,y_0)\cdot(a-x_0)+g_y(x_0,y_0)\cdot(b-y_0) +O(\varepsilon^2)
\end{aligned}
\right.
\end{eqnarray}
\)

を同時に満たす\((a,b)\)を探せばよい、ということになります。

行列形式で書けば
\(
\begin{align}
\left(
\begin{array}{cc}
f_x(x_0,y_0) & f_y(x_0,y_0) \\
g_x(x_0,y_0) & g_y(x_0,y_0)
\end{array}
\right)
\left(
\begin{array}{c}
a-x_0 \\
b-y_0
\end{array}
\right)
=
-\left(
\begin{array}{c}
f(x_0,y_0) \\
g(x_0,y_0)
\end{array}
\right)+O(\varepsilon^2)
\end{align}
\)

と書けます。\(a,b\)について解けば
\(
\begin{align}
\left(
\begin{array}{c}
a \\
b
\end{array}
\right)
=
\left(
\begin{array}{c}
x_0 \\
y_0
\end{array}
\right)
-\left(
\begin{array}{cc}
f_x(x_0,y_0) & f_y(x_0,y_0) \\
g_x(x_0,y_0) & g_y(x_0,y_0)
\end{array}
\right)^{-1}
\left(
\begin{array}{c}
f(x_0,y_0) \\
g(x_0,y_0)
\end{array}
\right)+O(\varepsilon^2)
\end{align}
\)

となります。丁度一次元のニュートン法と似た形になりました。

右辺の
\(
\begin{align}
\mathbf{J}_0=\left(
\begin{array}{cc}
f_x(x_0,y_0) & f_y(x_0,y_0) \\
g_x(x_0,y_0) & g_y(x_0,y_0)
\end{array}
\right)
\end{align}
\)

は\(x=x_0,y=y_0\)のヤコビアンです。また、
\(
\mathbf{a}=\left(
\begin{array}{c}
a \\
b
\end{array}
\right),~~
\begin{align}
\mathbf{x}_0=\left(
\begin{array}{c}
x_0 \\
y_0
\end{array}
\right),~~
\mathbf{f}_0=\left(
\begin{array}{c}
f(x_0,y_0) \\
g(x_0,y_0)
\end{array}
\right)
\end{align}
\)

と表記すれば、2次元のニュートン法は
\(
\mathbf{a}=\mathbf{x}_0 – \mathbf{J}_0^{-1}\mathbf{f}_0+O(\varepsilon^2)
\)

と書くことが出来ます。
近似として、\(O(\varepsilon^2)\)の項を無視すれば
\(
\mathbf{x_1}=\mathbf{x}_0 – \mathbf{J}_0^{-1}\mathbf{f}_0
\)

と書けます。繰り返せば、漸化式
\(
\mathbf{x}_{n+1}=\mathbf{x}_{n} – \mathbf{J}_n^{-1}\mathbf{f}_n
\)

を解いていけば良い、ということになります。

ニュートン法を実行するにあたり、ヤコビアンの逆行列を解くことが一番の問題となります。

良く知られているように、2次元の逆行列は、計算コストもさほどかからずに解けてしまい、
\(
\begin{align}
\mathbf{J}_n^{-1}=
\frac{1}{f_x(x_n,y_n) g_y(x_n,y_n)-f_y(x_n,y_n)g_x(x_n,y_n)}
\left(
\begin{array}{cc}
g_y(x_n,y_n) & -f_y(x_n,y_n) \\
-g_x(x_n,y_n) & f_x(x_n,y_n)
\end{array}
\right)
\end{align}
\)

として求める事が出来ます。

参考文献[2],[6]

スポンサーリンク

2次元ニュートン法のプログラム


以下のプログラムは、2次元のニュートン法のプログラムで、
方程式
\(
\begin{eqnarray}
\left\{~~
\begin{aligned}
f_1(x,y)&=x^2+y^2-1 = 0 \\
f_2(x,y)&=y-x^3 = 0
\end{aligned}
\right.
\end{eqnarray}
\)

を初期値\(x_0=2,~y_0=1\)から始めて解くプログラムです。

このコードを動かすと、ステップを繰り返すごとに解に近づいていっていることが分かります。

$ gfortran main.f90
$ ./a.out
    1       1.3571428532359113       0.2857142813732352
    2       0.9844126813492445       0.4401111876910193
    3       0.8485699953676349       0.5590406123395349
    4       0.8265084709901402       0.5633730814110012
    5       0.8260315915076198       0.5636240770681390
    6       0.8260313576542442       0.5636241621612319
    7       0.8260313576541870       0.5636241621612585
  0.82603135765418700       0.56362416216125855       ---result
$

この問題の答えは\(f_1=0,~f_2=0\)の陰関数の交点です。
図示すれば、下図の紫色の経路を辿って解に収束していくことが分かります。

ひとつ、注目しておかなければならない点として、解に近いからと言って、その解に収束するとは限らないということに注意してください。
図中の水色の線は初期値こそ右上の解に近いですが、実際に計算してみるともう一つの解に収束していることが分かります。

適当な初期値を用意し、ニュートン法が収束する点を調べるという問題はニュートン写像と呼ばれる分野です。
力学系とかそういう言葉が出てきます。
私も過去に複素関数でやったことがありますので、リンクだけ載せておきます。
ニュートン写像に現れる綺麗な画像 -シキノート

2次元の複素関数のニュートン法

置いておきます。

多次元のニュートン法

さて、この形まで持ってくると多次元への拡張が簡単にできます。
多次元のニュートン法は
\(
\mathbf{x}_{n+1}=\mathbf{x}_n – \mathbf{J}_n^{-1}\mathbf{f}_n
\)

の形がそのまま使えます。
ここで、
\(
\begin{align}
\mathbf{x}_n=
\left(
\begin{array}{c}
x_1\\
x_2\\
\vdots \\
x_N
\end{array}
\right)_n
,~~
\mathbf{f}_n=
\left(
\begin{array}{c}
f_1(\mathbf{x}_n)\\
f_2(\mathbf{x}_n)\\
\vdots \\
f_N(\mathbf{x}_n)
\end{array}
\right)
,~~
\mathbf{J}_n=
\left.\left(
\begin{array}{cccc}
\frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_N} \\
\frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \cdots & \frac{\partial f_2}{\partial x_N} \\
\vdots&\vdots&\ddots&\vdots \\
\frac{\partial f_N}{\partial x_1} & \frac{\partial f_N}{\partial x_2} & \cdots &\frac{\partial f_N}{\partial x_1}
\end{array}
\right)\right|_{\mathbf{x}=\mathbf{x}_n}
\end{align}
\)

です。

問題となるのは(ヤコビアンの逆行列)×(f)の計算です。
この計算を数値計算的に何とかするべく改良が重ねられてきたのが準ニュートン法(本稿では書きません)と呼ばれるアルゴリズムです。

逆行列を数値的に計算することは宜しくないですので、
\(
\begin{align}
\mathbf{r}&=\mathbf{J}_n^{-1}\mathbf{f}_n \\
\to \mathbf{J}_n \mathbf{r} = \mathbf{f}_n
\end{align}
\)

という方程式を解く問題に帰着させます。これが解ければ、ニュートン法の次のステップを
\(
\mathbf{x}_{n+1}=\mathbf{x}_n – \mathbf{r}
\)

として求める事が出来るのです。

高次のニュートン法?


1次元のニュートン法は
\(\displaystyle
a = x_0-\frac{f(x_0)}{f'(x_0)}+O(\Delta^2)
\)

の形で求めていきますが、二階微分が分かっていればもっと収束が早くなりそうです。
実際、
\(\displaystyle
a = x_0-\frac{f^{\prime}(x_0)}{f^{\prime\prime}(x_0)}\pm\sqrt{{f’}^2(x_0)-2f(x_0)f^{\prime\prime}(x_0)}+O(\Delta^3)
\)

と解けます。しかし、2次関数ですので
解が見付からないかもしれない
という可能性があります。1次のニュートン法であれば、導関数がゼロでない限り、必ず\(y=0\)と交わる点が見付かります。しかし、2次で解に近づけると交点が見つからないことがあるかもしれません。

なので、1次のニュートン法を繰り返して使う方が(2次と比べたら)安定な解法なのです。

参考ページ

[1]1 Newton 法:壱変数の場合
https://www.astr.tohoku.ac.jp/~chinone/Planck/Planck-node6.html
[2]2 Newton 法:弐変数の場合
https://www.astr.tohoku.ac.jp/~chinone/Planck/Planck-node7.html
[3]Newton法
http://www.misojiro.t.u-tokyo.ac.jp/~murota/lect-suchi/newton130805.pdf
[4]アルゴリズムによる誤差
http://www.aoni.waseda.jp/ykagawa/chap2html/node9.html
[5]Newton 法による方程式の近似解法
http://www.math.u-ryukyu.ac.jp/~suga/C/2004/7/node9.html
[6]多次元のニュートン・ラフソン法
http://www.slis.tsukuba.ac.jp/~fujisawa.makoto.fu/cgi-bin/wiki/index.php?%C2%BF%BC%A1%B8%B5%A4%CE%A5%CB%A5%E5%A1%BC%A5%C8%A5%F3%A1%A6%A5%E9%A5%D5%A5%BD%A5%F3%CB%A1
[7]山本 有作, 行列計算における高速アルゴリズム
http://www.cms-initiative.jp/ja/events/0620yamamoto.pdf
[8]中田 和秀, 大規模線形方程式を解くためのクリロフ部分空間法の前処理
http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/1288-6.pdf

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/(パブリックドメイン)です。

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}
\)
としてデルタ関数が表現されていることが分かります。

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

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

と求められるようです
このように表現されていると考えなければ、\(\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で回らなくていいことは、具体的に複素平面上の数値積分を行って確かめています。

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

まとめ


離散フーリエ変換(Discrete Fourier Transform,DFT)を
\(
\begin{eqnarray}
\left\{
\begin{aligned}
f(\tilde{k}_n)&=\mathcal{F}[f(x)](\tilde{k}_n)=\sum_{m=0}^{N-1}f(x_m)e^{-i2\pi\tilde{k}_nx_m}\Delta x \\
f(x_m)&=\mathcal{F}^{-1}[f(\tilde{k})](x_m)=\sum_{n=0}^{N-1}f(\tilde{k}_n)e^{i2\pi\tilde{k}_nx_m}\Delta \tilde{k}
\end{aligned}
\right.
\end{eqnarray}
\)

と定義し、畳み込み(Convolution)を
\(
\displaystyle (f\ast g)(x_l)= \sum_{m} f(x_m) g(x_l-x_m) \Delta x
\)

と定義します。この時、畳み込みは上で定義した離散フーリエ変換を用いて以下の通り書けます。
\(
\displaystyle (f\ast g)(x_l)=\mathcal{F}^{-1}\left[\mathcal{F}[f(x)](\tilde{k})\cdot \mathcal{F}[g(x)](\tilde{k})\right](x_l)
\)

intel®マス・カーネル・ライブラリ(MKL)に実装されているフーリエ変換は
\(
\begin{eqnarray}
\left\{
\begin{aligned}
(-1)^n f(\tilde{k}_n)&=\sum_{m=0}^{N-1}f(x_m)e^{-i2\pi nm/N} \\
f(x_m)&=\sum_{n=0}^{N-1}(-1)^n f(\tilde{k}_n)e^{i2\pi nm/N}
\end{aligned}
\right.
\end{eqnarray}
\)

の形で定義されています。ここで、\(f(\tilde{k}_n)\)は、正しいフーリエ変換後の関数。計算したもの。この時、上で定義されている場合、畳み込みは、
\(
\displaystyle (f\ast g)(x_m)=\frac{\Delta x}{N}\mathcal{F}_{MKL}^{-1}\left[\mathcal{F}_{MKL}[f(x)]\cdot\mathcal{F}_{MKL}[g(x)] (-1)^n\right](x_m)
\)

で計算できます。ここで、\(N\)は離散フーリエ変換に用いた分点数、\(n\)は順方向フーリエ変換後の空間を離散化した時のインデックスです。

intel®MKLの離散フーリエ変換ルーチンを用いた畳み込みは、プログラムでは\(h(x)=(f\ast g)(x)\)とした時、

\(
\begin{align}
& Status = DftiComputeForward(hand,f) \\
& Status = DftiComputeForward(hand,g) \\
& h(\tilde{k}_n)=z(\tilde{k}_n)\times g(\tilde{k}_n)\times (-1)^n ~~~~~(\mbox{for all } k_n) \\
& Status = DftiComputebackward(hand,h) \\
& h(x_m)=h(x_m)\times \frac{\Delta x}{N} ~~~~~(\mbox{for all } x_m)\\
\end{align}
\)

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


スポンサーリンク


目次


  1. フーリエ変換の定義
  2. 離散フーリエ変換の定義
  3. たたみ込みの定義
  4. 離散たたみ込みと離散フーリエ変換
  5. 数値計算で使われている離散フーリエ変換アルゴリズム (例: intel®MKL)

フーリエ変換の定義


フーリエ変換の定義方法には複数の慣例があり、
数値計算分野、化学、音の解析等においては
\(
\begin{align}
f(\tilde{k})&=\int_{-\infty}^{\infty} f(x)e^{-i2\pi\tilde{k}x}dx \\
f(x)&=\int_{-\infty}^{\infty} f(\tilde{k})e^{i2\pi\tilde{k}x}d\tilde{k}
\end{align}
\)

の形で良く定義されます。
ここで、例えば\(x\)は位置[m]であれば\(\tilde{k}\)は波数[1/m]です。
また、\(x\)が時間[s]であれば\(\tilde{k}\)は周波数[1/s]です。

このページでは扱いませんが、物理科学の世界では
\(
\begin{align}
f(k)&=\int_{-\infty}^{\infty} f(x)e^{-ikx}dx \\
f(x)&=\int_{-\infty}^{\infty} f(k)e^{ikx}\frac{dk}{2\pi}
\end{align}
\)

の形で良く定義されます。
ここで、例えば\(x\)は位置[m]であれば\(k\)は角波数[1/m]です。
また、\(x\)が時間[s]であれば\(k\)は角周波数[1/s]です。

この2つの違いは周波数か、角周波数のどちらが本質であるか?という事だと思います。
物理では、周波数よりも角周波数の方が式が簡便になります。
ですが、人間の理解がしやすいのは1秒当たり何回振動するか?という周波数だろう、と思うため工学よりの分野では周波数が使われるのでしょう。もしかしたら、単に慣例だけかもしれません。

どちらも定数倍の変数変換の違いだけなので、本質は変わりません。
数値計算依りの話をしていきたいので、ここからは前者の位置⇔波数で考えていきます。

また、

順方向フーリエ変換:\(e^{-i2\pi\tilde{k}x}\)の因子が掛かっている式、空間→波数
逆方向フーリエ変換:\(e^{+i2\pi\tilde{k}x}\)の因子が掛かっている式、波数→空間

と呼び、
”関数\(f(x)\)の順方向フーリエ変換して波数\(\tilde{k}\)の空間に移る”
という操作を華文字\(\mathcal{F}\)を用いて
\(
\mathcal{F}[f(x)](\tilde{k})
\)

と表します。また、
”関数\(f(k)\)の逆方向フーリエ変換して位置\(x\)の空間に移る”
という操作を
\(
\mathcal{F}^{-1}[f(\tilde{k})](x)
\)

と表します。

離散フーリエ変換の定義


このページでは離散フーリエ変換を
\(
\begin{align}
f(\tilde{k}_n)&=\sum_{m=0}^{N-1}f(x_m)e^{-i2\pi\tilde{k}_nx_m}\Delta x \\
f(x_m)&=\sum_{n=0}^{N-1}f(\tilde{k}_n)e^{i2\pi\tilde{k}_nx_m}\Delta \tilde{k}
\end{align}
\)

定義します
\(\Delta x\), \(\Delta \tilde{k}\)は、それぞれ位置、波数空間の刻み幅を表します。

この上の定義は簡便さのために良く使われる
\(
\begin{align}
f(\tilde{k}_n)&=\frac{1}{N}\sum_{m=0}^{N-1}f(x_m)e^{-i2\pi\tilde{k}_nx_m} \\
f(x_m)&=\hspace{1.5em}\sum_{n=0}^{N-1}f(\tilde{k}_n)e^{i2\pi\tilde{k}_nx_m}
\end{align}
\)

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

上の2つの定義が同じ理由は、順方向→逆方向を行い、元に戻った時に係数が\(\Delta x\Delta \tilde{k}=\frac{1}{N}\)になるためです。

\(x_m\) \(\tilde{k}_n\)
添え字 \(m=0,1,2,\cdots,N-1\) \(n=0,1,2,\cdots,N-1\)
区間 \(\displaystyle [a,b]\) \(\displaystyle \left[-\frac{1}{2\Delta x},\frac{1}{2\Delta x}\right]\)
刻み幅 \(\displaystyle \Delta x=\frac{b-a}{N}\) \(\displaystyle \Delta \tilde{k}=\frac{1}{b-a}\)
分点の位置 \(\displaystyle x_m=m\Delta x +a\) \(\displaystyle \tilde{k}_n=n\Delta \tilde{k} -\frac{1}{2\Delta x}\)

ただし、\(x_m, \tilde{k}_n\)は以下の関係式を見たします。
\(
\begin{align}
& x_m+N\Delta x=x_m\\
& \tilde{k}_n+N\Delta \tilde{k}=\tilde{k}_n
\end{align}
\)

ナイキスト(Nyquist)周波数について


ナイキスト周波数とはサンプリングする間隔でどの周波数を持つ関数までを表現できるか?という下限と上限を与えます。
今回の場合、上の表の通り\(\tilde{k}\)空間で表現できる上限と下限は\(\pm\frac{1}{2\Delta x}\)です。この絶対値、すなわち\(\frac{1}{2\Delta x}\)をナイキスト周波数と呼びます。

ナイキスト周波数が意味しているのは、高周波の成分を取り出したければサンプリングする間隔を小さくしなければならない。ということです。”サンプリングの間隔”は周波数と同じ意味ですので、サンプルする周波数を高くしなければならない。と言い換えることもできます。

具体例を上げましょう。

上の図は関数\(\displaystyle f(x)=e^{i2\pi \tilde{k}x}\)をサンプル間隔\(\Delta x = 0.5\)で離散的な関数で表現したものです。
この場合、ナイキスト周波数は\(\frac{1}{2\times 0.5}\)より、\(1\)です。
すなわち、このサンプル間隔では\(\tilde{k}\)が\([-1,1]\)の範囲を持つものしか表現することが出来ません。
では、実際にサンプル出来ない関数\(\displaystyle f(x)=e^{i2\pi 1.8 x},~~(\tilde{k}=1.8)\)を考えた時に何が起こるのでしょうか。
それを示したのが下から2番目の関数です。
本来は早く振動する関数であるのに、あたかも非常にゆっくりとした関数だ、と捉えられてしまいます。しかも位相が反転しています。

これは、元々は\(e^{i\theta}\)の周期性によるものです。簡単な計算の通り、
\(
\begin{align}
e^{i(\theta+2\pi)}=e^{i\theta}e^{i2\pi}=e^{i\theta}
\end{align}
\)

となります。この周期性のために、\(\tilde{k}:[-1,1]\)の範囲でしか表現できない空間で、
その範囲外の\(\tilde{k}\)を表現しようとした時には、\(\tilde{k}\)空間の大きさ(…この場合は\(2\))を足したり引いたりして\(\tilde{k}:[-1,1]\)の範囲に強制的に押し込まれます。
\(\tilde{k}=1.8\)の場合は、\(\tilde{k}=1.8-2=-0.2\)となります。
言葉で言えば、
\(\tilde{k}=1.8\)を持つ波形は\(\tilde{k}=-0.2\)を持つ波形としてとして表現される
ということです。

本来の関数を薄く描いて重ねたものが以下のものになります。

フーリエ変換が元に戻ることの証明


上記のフーリエ変換、離散フーリエ変換はある制限を課しています。
それは、
ある関数に順方向フーリエ変換を行い、続いて逆方向フーリエ変換を施した場合、元と同じ関数に戻っていて欲しい。
という要望です。
それは、以下のように証明できます。

フーリエ変換
\(
\begin{align}
f(x)&=\int_{-\infty}^{\infty}d\tilde{k} f(\tilde{k})e^{i2\pi\tilde{k}x} \\
&=\int_{-\infty}^{\infty}d\tilde{k} \left[\int_{-\infty}^{\infty}dx’ f(x’)e^{-i2\pi\tilde{k}x’}\right]e^{i2\pi\tilde{k}x} \\
&=\int_{-\infty}^{\infty}dx’ f(x’)\int_{-\infty}^{\infty}d\tilde{k}
e^{i2\pi\tilde{k}(x-x’)} \\
&=\int_{-\infty}^{\infty}dx’ f(x’)\delta(x-x’) \\
&=f(x)
\end{align}
\)

離散フーリエ変換

\(
\begin{align}
f(x_m)&=\sum_{n=0}^{N-1}f(\tilde{k}_n)e^{i2\pi\tilde{k}_nx_m}\Delta \tilde{k}\\
&=\sum_{n=0}^{N-1}\left[\sum_{m’=0}^{N-1}f(x_{m’})e^{-i2\pi\tilde{k}_nx_{m’}}\Delta x\right]e^{i2\pi\tilde{k}_nx_m}\Delta \tilde{k} \\
&=\Delta \tilde{k}\Delta x\sum_{m’=0}^{N-1}f(x_{m’})\sum_{n=0}^{N-1}e^{i2\pi\tilde{k}_n(x_m’-x_m)} \\
&=\frac{1}{N}\sum_{m’=0}^{N-1}f(x_{m’})\sum_{n=0}^{N-1}e^{i2\pi n (m’-m)/N +i2\pi(\theta_{m}-\theta_{m’})} \\
&=\frac{1}{N}\sum_{m’=0}^{N-1}f(x_{m’})e^{i2\pi(\theta_{m}-\theta_{m’})}N\delta_{m,m’} \\
&=f(x_m)
\end{align}
\)
ここで、
\(
\begin{align}
\tilde{k}_n x_{m}=\frac{nm}{N}+\frac{an}{N\Delta x}+\theta_{m} \\
\theta_{m}=-\frac{a}{2\Delta x}-\frac{m}{2}
\end{align}
\)

そして、等比級数の考えから、
\(
\begin{eqnarray}
\sum_{n=0}^{N-1}e^{i2\pi n a}=
\left\{
\begin{aligned}
&\frac{e^{i2\pi aN}-1}{e^{i2\pi a}-1}, \hspace{1em} ( a\ne 0) \\
& N, \hspace{4em} (a=0) \\
\end{aligned}
\right.
\end{eqnarray}
\)

という関係式を用いました。

畳み込み(convolution)の定義


このページでは、畳み込みを
\(
\displaystyle (f\ast g)(x)=\int f(\tau) g(x-\tau) d\tau
\)

という積分であると定義し、離散的な畳み込みを
\(
\displaystyle (f\ast g)(x_l)= \sum_{m} f(x_m) g(x_l-x_m) \Delta x
\)

定義します。

畳み込みの定義には\(\Delta x\)が掛けられていないものが存在しますが、積分の意味が薄れてしまうため、離散フーリエ変換時と同様、このページでは\(\Delta x\)を掛けたものを畳み込み、と定義します。

区間\([a,b]\)で定義された離散的な畳み込みの場合、関数\(f(x),g(x)\)は区間が\([a,b]\)で決められている場合がほとんどです。
この場合、引数\(x_m-x_l\)が定義域内に収まらない場合が現れます。

定義域からはみ出てしまった場合は推定するしかありません。
推定の仕方によって
循環畳み込みと呼ばれる方法と直線畳み込みと呼ばれる方法が良くつかわれています。

循環畳み込みは領域に対して周期境界条件を課したもの、すなわち、上限bと下限aが繋がっていると考えます。
直線畳み込みは領域に対して固定端条件を課したもの、すなわち、[a,b]の範囲外の関数値はゼロと考えます。

問題に大きく依るので、どちらがいいかはありません。

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


離散畳み込みは定義通り
\(
\displaystyle (f\ast g)(x_l)= \sum_{m} f(x_m) g(x_l-x_m) \Delta x
\)

で計算することが出来ます。計算量は\(O(n^2)\)です。

しかし、実用上は離散フーリエ変換を利用すると計算量を\(O(n\log n)\)に抑えられます。
計算方法はたたみ込み定理を利用して、
\(
\displaystyle \mathcal{F}[(f\ast g)(x)](k)=\mathcal{F}[f(x)](k)\cdot \mathcal{F}[g(x)](k)
\)

より、\((f\ast g)(x)\)は
\(
\displaystyle (f\ast g)(x)=\mathcal{F}^{-1}\left[\mathcal{F}[f(x)](k)\cdot \mathcal{F}[g(x)](k)\right]
\)

と書けるため、フーリエ変換を介することで計算量を減らせます。
離散フーリエ変換を介して畳み込みを計算する場合、必ず循環畳み込みで計算することになります。

畳み込み定理の証明

連続畳み込み

\(
\displaystyle \mathcal{F}[(f\ast g)(x)](k)=\mathcal{F}[f(x)](k)\cdot \mathcal{F}[g(x)](k)
\)

左辺)
\(
\begin{align}
\mathcal{F}[(f\ast g)(x)](\tilde{k})
&=\int \left[\int f(\tau)g(x-\tau) d\tau\right] e^{-i2\pi \tilde{k} x} dx \\
& x-\tau=y\mbox{の変数変換を行って} \nonumber \\
&=\int dy\int d\tau f(\tau)g(y) e^{-i2\pi \tilde{k} (\tau+y)} \\
&=\int dx\int dx’ f(x)g(x’) e^{-i2\pi \tilde{k} (x+x’)}
\end{align}
\)

右辺)
\(
\begin{align}
\mathcal{F}[f(x)](\tilde{k})\cdot \mathcal{F}[g(x)](\tilde{k})
&=\int f(x) e^{-i2\pi \tilde{k} x} dx \cdot \int g(x’) e^{-i2\pi \tilde{k} x’} dx’ \\
&=\int dx\int dx’ f(x)g(x’) e^{-i2\pi \tilde{k} (x+x’)}
\end{align}
\)

より、証明終了。また、逆変換により戻ることを示します。
\(
\begin{align}
&\mathcal{F}^{-1}\left[\mathcal{F}[f(x)](\tilde{k})\cdot \mathcal{F}[g(x)](\tilde{k})\right] \\
&=\int \left[\int dx\int dx’ f(x)g(x’) e^{-i2\pi \tilde{k} (x+x’)}\right] e^{i2\pi \tilde{k} y} dk\\
&=\int dx\int dx’ f(x)g(x’) \int e^{-i2\pi \tilde{k} (x+x’-y)} dk\\
&=\int dx\int dx’ f(x)g(x’) \delta(x’-(y-x)) \\
&=\int dx f(x)g(y-x) \\
&=\int d\tau f(\tau)g(x-\tau) \\
\end{align}
\)
となり、フーリエ変換を介してたたみ込みが計算できることが示せました。

離散畳み込み

\(
\displaystyle \mathcal{F}[(f\ast g)](\tilde{k})=\mathcal{F}[f(x)](\tilde{k})\cdot \mathcal{F}[g(x)](\tilde{k})
\)

左辺)
\(
\begin{align}
\mathcal{F}[(f\ast g)(x)](\tilde{k}_n)
&=\mathcal{F}\left[\sum_{m=0}^{N-1}f(x_m)g(x-x_m)\Delta x\right] \\
&=\sum_{l=0}^{N-1}\left[\sum_{m=0}^{N-1}f(x_m)g(x_l-x_m)\Delta x\right]e^{i2\pi \tilde{k}_n x_l}\Delta x \\
&=(\Delta x)^2\sum_{l=0}^{N-1}\sum_{m=0}^{N-1}f(x_m)g(x_l-x_m) e^{i2\pi \tilde{k}_n x_l} \\
&\mbox{ここで、循環たたみ込みを想定し、}x_l-x_m=x_{m’}\mbox{となるように}x_{m’}\mbox{を決めると}\nonumber \\
&=(\Delta x)^2\sum_{m=0}^{N-1}\sum_{m’=0}^{N-1}f(x_m)g(x_{m’}) e^{i2\pi \tilde{k}_n (x_m+x_{m’})} \\
\end{align}
\)

右辺)
\(
\begin{align}
\mathcal{F}[(f(x)](k_n)\cdot \mathcal{F}[(g(x)](k_n)
&=\sum_{m=0}^{N-1}f(x_m) e^{i2\pi \tilde{k}_n (x_m)}\Delta x\cdot \sum_{m’=0}^{N-1}f(x_{m’}) e^{i2\pi \tilde{k}_n (x_{m’})}\Delta x \\
&=(\Delta x)^2\sum_{m=0}^{N-1}\sum_{m’=0}^{N-1}f(x_m)g(x_{m’}) e^{i2\pi \tilde{k}_n (x_m+x_{m’})}
\end{align}
\)
よって右辺と左辺が正しくなります。

逆変換によって離散畳み込みになることを示します。

\(
\begin{align}
\mathcal{F}^{-1}[\mathcal{F}[(f\ast g)(x)](k)](x_l)
&=\mathcal{F}^{-1}\left[(\Delta x)^2\sum_{m=0}^{N-1}\sum_{m’=0}^{N-1}f(x_m)g(x_{m’}) e^{i2\pi \tilde{k}_n (x_m+x_{m’})}\right]\\
&=\sum_{n=0}^{N-1}\left[(\Delta x)^2\sum_{m=0}^{N-1}\sum_{m’=0}^{N-1}f(x_m)g(x_{m’}) e^{i2\pi \tilde{k}_n (x_m+x_{m’})}\right]e^{-i2\pi \tilde{k}_n x_l}\Delta \tilde{k}\\
&=(\Delta x)^2 \Delta \tilde{k} \sum_{m=0}^{N-1}\sum_{m’=0}^{N-1} f(x_m)g(x_{m’})\sum_{n=0}^{N-1}e^{i2\pi \tilde{k}_n (x_m+x_{m’}-x_l)} \\
&=\frac{\Delta x}{N} \sum_{m=0}^{N-1}\sum_{m’=0}^{N-1} f(x_m)g(x_{m’})N\delta_{x_{m’},x_l-x_m} \\
&=\sum_{m=0}^{N-1}f(x_m)g(x_l-x_m)\Delta x \\
&=(f\ast g)(x_l)
\end{align}
\)

よって畳み込みの定理通りに戻るため、離散フーリエ変換を介して畳み込みの計算ができます。

※注意
\(\Delta x, \Delta \tilde{k}\)が掛けられていない方の離散フーリエ変換と畳み込みでは元に戻りません。畳み込みの計算時に係数倍違って出てきます。

数値計算で使われている離散フーリエ変換アルゴリズム (例: intel®MKL)


数値計算では、上記の離散フーリエ変換の形、指数の肩に\(\tilde{k},x\)を含む物はあまり使われません。
そして、数値計算上での離散フーリエ変換の定義方法は複数存在します。
マニュアルを見てください。

ここでは、intel®マス・カーネル・ライブラリ(MKL)に使われている離散フーリエ変換アルゴリズムについて言及します。

MKLの変換はどうやら
\(
\begin{align}
(-1)^n f(\tilde{k}_n)&=\sum_{m=0}^{N-1}f(x_m)e^{-i2\pi nm/N} \\
f(x_m)&=\sum_{m=0}^{N-1}(-1)^n f(\tilde{k}_n)e^{i2\pi nm/N}
\end{align}
\)

という定義に従い変換しているようです。規格化因子\(\frac{1}{N}\)はユーザー任せです。
ここで、
\(f(\tilde{k}_n)\)は本来のフーリエ変換後の関数、
\(f(x)\)は\(x\)空間での関数
を表します。

”しているようです”と書いたのは、マニュアルと異なっているからです。マニュアルでは\((-1)^n\)など存在しませんが、実際試してみると、こういう変換をしているからです。

MKLではどういう\(\tilde{k}\)を想定しているかを考えたところ、以下の状況のようです。

\(x_m\) \(\tilde{k}_n\)
添え字 \(m=0,1,2,\cdots,N-1\) \(n=0,1,2,\cdots,N-1\)
区間 \(\displaystyle [-a,a]\) \(\displaystyle \left[0,N\Delta \tilde{k}\right]\)
刻み幅 \(\displaystyle \Delta x=\frac{2a}{N}\) \(\displaystyle \Delta \tilde{k}=\frac{1}{2a}\)
分点の位置 \(
\begin{eqnarray}
x_m=\left\{
\begin{aligned}
& m\Delta x ~~~~(m=0,\cdots,\frac{N}{2}-1) \\
& m\Delta x-\frac{1}{\Delta \tilde{k}} (m=\frac{N}{2},\cdots,N-1)
\end{aligned}
\right.
\end{eqnarray}
\)
\(\displaystyle \tilde{k}_n=n\Delta \tilde{k} ~~~~(n=0,\cdots,N-1)\)

ただし、\(x_m, \tilde{k}_n\)は以下の周期境界条件を見たします。
\(
\begin{align}
& x_m+N\Delta x=x_m\\
& \tilde{k}_n+N\Delta \tilde{k}=\tilde{k}_n
\end{align}
\)

区間[-a,a]で定義された\(x_m\)の順番がおかしく感じますが、これは周期境界条件を用いて、
区間[0,2a]で定義された\(x_m=m\Delta x ~~~~(m=0,\cdots,N-1)\)
と見ても良いです。
なので、実用上どうやって関数\(f(x)\)を入れればよいか?に困ることは無いでしょう。
この上記離散フーリエ変換を考えると、上で提示したMKLの離散フーリエ変換が導くことが出来、プログラム上でも一致します。

MKLの離散フーリエ変換の導出)

離散フーリエ変換時の注意点

上に記した、余分な係数\((-1)^n\)のためにMKLを利用して求められた離散フーリエ変換の値は、本来のフーリエ変換後の関数と異なります。
具体的に例を示しましょう。

\(e^{-\alpha x^2}\)のフーリエ変換は
\(
\displaystyle \int_{-\infty}^{\infty} e^{-\alpha x^2}e^{-i2\pi\tilde{k} x}dx=\sqrt{\frac{\pi}{\alpha}}e^{-\pi^2\tilde{k}^2/\alpha}
\)

です。離散フーリエ変換を定義通り計算したものと、MKLを用いて離散フーリエ変換すると以下のようになります。
左から\(f(x), \mathcal{F}[f(x)](\tilde{k}),\mathcal{F}^{-1}[\mathcal{F}[f(x)]](x)\)になっています。

ここで、因子\((-1)^n\)のために交互に符号が変わっていることが分かるでしょう。
後でやりますが、このせいで畳み込みの計算をする際には気を付けなければなりません。

計算コードはこちら。

MKLのルーチンは優秀です。誤差が余りたまらないようです。
定義通り、愚直に作ったプログラムでは計算速度が遅いだけでなく、計算精度も悪くなります。
下の画像は、離散フーリエ変換→逆変換を複数回繰り返した時、元の関数とどれだけずれるかを表しています。

定義通りの離散フーリエ変換を\(10^n\)回繰り返すと\(n\)に対し線形にあがっていきます。この誤差は丸め誤差に起因しているのだと思います。しかし、MKLのルーチンには起こっていません。どういうアルゴリズムか知りませんが、優秀なのでしょう。

スポンサーリンク


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

MKLの離散フーリエ変換ルーチンを使って畳み込みを計算します。
例として、以下の関数を考えましょう。
\(
\begin{align}
h(x)&=\int \tau^5 e^{-\tau^2}\cdot e^{-4(x-\tau)^2}d\tau \\
&=\frac{1}{3125}\sqrt{\frac{\pi}{5}}x(1024x^4+1600x^2+375)e^{-\frac{4}{5}x^2}
\end{align}
\)

ここで、
\(
\begin{align}
f(\tau)=\tau^5 e^{-\tau^2}
g(\tau)=e^{-4\tau^2}
\end{align}
\)
です。

このページで紹介した定義の畳み込みでは証明した通り、
\(
h(x_m)=\mathcal{F}^{-1}\left[\mathcal{F}[f(x)](\tilde{k})\cdot\mathcal{F}[g(x)](\tilde{k})\right](x_m)
\)

で計算できます。
しかし、MKLの離散フーリエ変換(\(\mbox{DFT}_{MKL}\))の場合は簡単にはできません。計算は、
\(
\displaystyle (f\ast g)(x_m)=\frac{\Delta x}{N}\mathcal{F}_{MKL}^{-1}[\mathcal{F}_{MKL}[f(x)]\cdot\mathcal{F}_{MKL}[g(x)] (-1)^n](x_m)
\)

としなければなりません。\((-1)^n\)は順方向フーリエ変換が2回あり、MKLの逆方向離散フーリエ変換に必要なこの係数が落ちてしまうため入れなければなりません。また、MKLのルーチンでは係数\(\Delta x, \Delta\tilde{k}\)を含んでいないので計算しておく必要があります。

intel®MKLを用いる際、
順方向離散フーリエ変換にはDftiComputeForward、
逆方向離散フーリエ変換にはDftiComputeBackward
を使うと、擬コードはこのように書けるかと思います。

\(
\begin{align}
& Status = DftiComputeForward(hand,f) \\
& Status = DftiComputeForward(hand,g) \\
& h(\tilde{k}_n)=z(\tilde{k}_n)\times g(\tilde{k}_n)\times (-1)^n ~~~~~(\mbox{for all } k_n) \\
& Status = DftiComputebackward(hand,h) \\
& h(x_m)=h(x_m)\times \frac{\Delta x}{N} ~~~~~(\mbox{for all } x_m)\\
\end{align}
\)

\(\displaystyle \Delta x \Delta\tilde{k}=\frac{1}{N}\)の関係を用いて最後にまとめて定数倍しています。

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

という結果を得ます。

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

ニュートン写像に現れる綺麗な画像

ニュートン写像はニュートン=ラフソン法によって、ある初期値から適当な値へ収束するまでの回数でその関数を特徴づけるものです。
フラクタルと似ている考えであることを明記しておきます。

ニュートン写像


\(
\displaystyle N_f(z)=z-\frac{f(z)}{f'(z)}
\)

の繰り返しによって適当な初期値からスタートし、収束させていきます。このページでは画像が主です。説明の方はリーマンのゼータ関数の複素力学系をご覧ください。

基本はこれだけです。簡単な複素関数でも相当複雑になります。
早速適当な関数を用意して、ニュートン写像を見ていきましょう。


スポンサーリンク


まずは基本的な多項式を見ます。
\(
f(z)=z^2-1
\)

\(
f(z)=z^3-1
\)

\(
f(z)=z^4-1
\)

\(
\displaystyle f(z)=z^5-1, z^6-1, z^7-1, z^10-1
\)

少し複雑な関数の場合
\(
\displaystyle f(z)=\ln z, e^z-1
\)

\(
\displaystyle f(z)=e^{1/z}-0.01
\)

\(
\displaystyle f(z)=ze^{z}\ln z-0.01
\)

\(
\displaystyle f(z)=z^5e^z-0.1
\)

\(
\displaystyle f(z)=z^4e^z\ln z-0.1
\)

\(
\displaystyle f(z)=sin(z)/z
\)

\(
\displaystyle f(z)=z^2 e^{-z^2}-1
\)

\(
\displaystyle f(z)=ln(z) e^{-z^2}-1
\)

最後にリーマンのゼータ関数

複素関数\(f(z)\)が多くの根を持つとき、非常に複雑で綺麗な画像が得られるようです。
恐らく遠く離れたところからたまたま根のあたりに到達することが多くなるからでしょう。

スポンサーリンク

fortran90によるプログラム


fortran90によるプログラムを置いておきます。入力するべき場所は!=====の間と、関数fと導関数dfです。

program main
  implicit none
  integer,parameter::kmax=49
  double precision,parameter::eps=1d-8
  integer::i,j,k,Nx,Ny
  double precision::xa,xb,ya,yb,x,y,hx,hy
  complex(kind(0d0))::f,df,z,z1
  external::f,df

  !==================
  xa=-5d0; xb=5d0
  Nx=201

  ya=-5d0; yb=5d0
  Ny=201
  !==================
   
  hx=(xb-xa)/dble(Nx-1)
  hy=(yb-ya)/dble(Ny-1)
  do i=0,Nx-1
     x=xa+i*hx
     do j=0,Ny-1
        y=ya+j*hy
        z=dcmplx(x,y)
        do k=0,kmax
           z1=z-f(z)/df(z)
           if(abs((z-z1)/z).le.eps)exit
           z=z1
        enddo
        write(10,'(2e18.8e3,i5)')x,y,k
     enddo
     write(10,*)
     if(mod(i,100).eq.0)write(6,'(A,i6,A,i6)')"  ",i," / ",Nx
  enddo
 
  stop
end program main

function f(z)
  implicit none
  complex(kind(0d0))::f,z
 
  f=z**5*exp(z)-0.1d0
 
  return
end function f

function df(z)
  implicit none
  complex(kind(0d0))::df,z
 
  df=5*z**4*exp(z)+z**5*exp(z)
 
  return
end function df