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

複素関数のメモ

複素関数論のメモです。

  1. 微分
  2. 積分
  3. 参考文献

微分


コーシー・リーマンの関係式(まとめ)


複素関数がコーシー・リーマンの関係式を満たすと、その複素関数は通常の変数通りに微分することが出来ます。複素数\(z=x+iy\)を引数に持つ複素関数\(f(z)=f(x+iy)\)が、実数値関数\(u(x,y), v(x,y)\)を用いて

と書けているとき、コーシー・リーマンの関係式

を満たすならば、\(f(z)\)は正則関数と呼ばれます。ここで、
\(f_x(x,y)=\frac{\partial f}{\partial x}, f_y(x,y)=\frac{\partial f}{\partial y}\)を意味します。式(2)を満たす場合、複素関数\(f(z)\)の複素数\(z\)による微分は

と書くことができ、式(3), (4)は同じ値となります。

コーシー・リーマンの関係式(導出)


複素関数の引数の数は1つであり、それは複素数です。すなわち1変数関数です。故に、複素関数の微分を定義しようとすれば、
\begin{align}
\frac{df(z)}{dz}
\end{align}
という量は1つに定まるはずです。しかし実数の空間で生活している我々には、1つの複素数を実部と虚部に分けて考えてしまう癖のようなものがあります。そのせいで、複素関数は平面に記述される関数であり、実部方向と虚部方向をそれぞれ独立に考えることが出来る2変数関数だと勘違いしてしまいますが、そうではありません。

すなわち、2変数の実数値関数と1変数の複素関数は全く異なる関数であり、同じ様に扱ってはいけません。しかし、我々は1次元上に複素数を表現する方法を持っておらず、本来、1変数の複素関数を2変数の実数値関数のように表現しなければイメージすることが出来ません。この場合、複素関数の表現に何らかの制約がかかることがイメージできるでしょう。その制約が正則である、という条件であり、コーシー・リーマンの関係式なのです。

では実際に本来1変数関数の複素関数を実軸と虚軸という、あたかも2変数関数のように拡張して考えて、微分を考えてみましょう。

前提として、複素関数の微分
\begin{align}
\frac{df(z)}{dz}
\end{align}
は一意に定まるとします。この下で2変数関数として\(f(z)=f(x,y)\)を拡張した場合、微分したい点への近づき方が二通りあることに気が付きます。それは、\(x\)方向と\(y\)方向ということです。2変数関数であれば、偏微分の事です。どんな方向から近づいても微分値は一意に定まると考えているので、両結果は等しくなければなりません。
微分は実数値関数と同様に差分の極限として定義して、
\begin{align}
\frac{df(z)}{dz}\equiv \lim_{\Delta z\to 0}\frac{f(z+\Delta z)-f(z)}{\Delta z}
\end{align}
と書けるとします。ここで、\(\Delta z = \Delta x+i\Delta y\)を意味しており、\(\Delta x, \Delta y\)は実数です。

点\(z=(x, y)\)に実数軸に沿って近づく場合、\(\Delta y=0\)を考えればよいので、微分は差分の極限として
\begin{align}
\frac{df(z)}{dz}=\lim_{\Delta x\to 0}\frac{f(x+\Delta x,y)-f(x,y)}{\Delta x}
\end{align}
と書かれます。ここで\(h\)は実数です。また、虚数軸に沿って近づく場合\(\Delta x=0\)を考えればよいので、
\begin{align}
\frac{df(z)}{dz}=\lim_{\Delta y\to 0}\frac{f(x,y+\Delta y)-f(x,y)}{i\Delta y}
\end{align}
と書かれるはずです。複素関数の微分は一意に決まることから、二つの近づき方を考えたとしても一致していなければなりません。複素関数が\(f(z)=u(x,y)+iv(x,y)\)と実部と虚部があらわに分離されて書かれていたとするならば、
\begin{align}
&\lim_{\Delta x\to 0}\frac{f(x+\Delta x,y)-f(x,y)}{\Delta x}\\
&=\lim_{\Delta x\to 0}\frac{[u(x+\Delta,y)+iv(x+\Delta,y)]-[u(x,y)+iv(x,y)]}{\Delta x} \\
&=\lim_{\Delta x\to 0}\frac{[u(x+\Delta,y)-u(x,y)]+i[v(x+\Delta,y)-v(x,y)]}{\Delta x} \\
&=\frac{\partial u}{\partial x}+i\frac{\partial v}{\partial x} \\
\end{align}
と書けます。また、虚軸方向は
\begin{align}
&\lim_{\Delta y\to 0}\frac{f(x,y+\Delta y)-f(x,y)}{i\Delta y}\\
&=\lim_{\Delta y\to 0}\frac{[u(x,y+\Delta y)+iv(x,y+\Delta y)]-[u(x,y)+iv(x,y)]}{i\Delta y} \\
&=\lim_{\Delta y\to 0}\frac{-i[u(x,y+\Delta y)-u(x,y)]+[v(x,y+\Delta y)-v(x,y)]}{\Delta y} \\
&=-i\frac{\partial u}{\partial y}+\frac{\partial v}{\partial y}
\end{align}
と書くことが出来ます。
両者はもともと同じものであるとしてきましたから、実部と虚部で比較します。
すると、

という関係式が導けるという訳です。これがコーシー・リーマンの関係式と呼ばれるもので、1変数関数の複素関数を実部、虚部という2変数関数に無理矢理拡張した場合に現れる条件です。
複素関数の微分が定義できるとき、複素関数の実部虚部は、この関係式を満たしていなければなりません。ある範囲においてすべての点で複素関数の微分が定義できることを、正則であると言ったり、滑らかである、解析的であると言います。
”複素関数の滑らか”は通常の”滑らか”と意味が異なることに注意しましょう。通常は不連続点が無いことを指しますが、複素関数の場合はそれ以上の制約(コーシー・リーマンの関係式)を含んでいます。

また、ここでは述べませんが、複素関数は1度微分できると何回でも微分することが出来ます。
さらに微分が定義できる正則関数であれば定義域を超えた領域へ関数が拡張できる場合があります。これが解析接続と呼ばれるテクニックです。コーシー・リーマンの関係式を制約として今まで捉えてきたのですが、この制約があるからこそまだ見ぬ定義域外をこの関係式を手掛かりに関数を拡張できるのです。

例として、究極的に滑らかそうに見えるけど微分できない関数を上げましょう。
\begin{align}
f(z)=z=x+iy
\end{align}
はコーシー・リーマンの関係式を満たすため、複素関数として微分が定義でき、正則です。しかし、その複素共役を取った関数
\begin{align}
f(z)=z^*=x-iy
\end{align}
はコーシー・リーマンの関係式を満たさないため、複素関数として微分が定義できず、正則ではありません。近づき方によって微分の値が変わってしまうのです。

通常の微分が可能という意味を明確にするため、例を挙げます。\(f(z)\)が\(z\)の多項式で書かれている場合、微分は

として求める事が出来ます。通常と同じですね。これは、複素数の成分を持たない場合には実軸上で考えたものと同等になるので、実数軸上の微分を自然に拡張したようになっています。

積分


特異点周りの積分


端的にまとめます。
点\(a\)を囲うように複素平面上の閉じた積分経路の方向を\(C\)で指定するとします。\(z\)の多項式で表現される関数\(f(z)=(z-a)^n\)の積分は、

となります。点\(a\)の周りを反時計回りに1周回る場合、\(C\)はパラメータ\(t\)を用いて\(z(t)\)と書くと、
\begin{align}
z(t)=a+re^{it},~(0\le t \le 2\pi)
\end{align}
と書けます。\(r\)は積分経路の半径を意味します。式(6)の計算をすれば、

となります。\(n=-1\)の時、右辺の積分は\(2\pi\)に等しい事は明らかです。それ以外の\(n\ne -1\)ならば、

となります。よって、

が導かれます。

留数定理


一つ前の、ある点周りの積分をもう少し一般的に考えましょう。\(z=a\)を中心とした閉曲線に沿った複素関数\(f(z)\)の複素積分\(\oint_C f(z)dz\)は

と書けます。ここで\(\text{Res}[f,a]\)は留数と呼ばれる量です。これから留数について述べていきます。

\(\text{Res}[f,a]\)は、\(f(z)\)が点\(a\)周りで

とローラン展開出来るとします。この時、\(a_{-1}\)を留数と呼びます。すなわち、

を意味します。

では、留数を実際に求めていきましょう。\(f(z)\)が点\(a\)周りでローラン展開可能な場合、\(n\)の下限が重要になります。下限を極の位数\(k\)と表現します。\(k\)は、

を満たすときの値として定義できます。

もし、関数\(f(z)\)が\(z=a\)で\(k\)位の極であれば、\(z=a\)周りで

と展開できます。
留数は\(a_{-1}\)なので、\(a_{-1}\)について式変形を行っていけば、

と計算できます。

参考文献


原惟行、松永秀章著、『複素解析入門』(第3刷2010, 初版2007)

Hyper-dual numbersによる二階偏微分の計算

Hyper-dual numbersと呼ばれる、実数を拡張した考えを導入すると導関数が計算できます。

あらかじめ、引数がHyper-dual numbersである時の数々の関数の定義を実装しておけば、その関数の組み合わせで作られる関数の一階導関数、二階導関数のほぼ厳密な答えを得ることが出来ます。
この考えに従う関数の微分方法は、Forward型の自動微分と呼ばれます。

  1. Dual number
    1. Dual numberと導関数
  2. Hyper-dual numbers
    1. Hyper-dual numbersの演算
    2. Hyper-dual numbersのプログラム
    3. ヘッセ行列
  3. 参考文献

Dual number


Dual numberという考えがあります[2,3,5]。

これは実数を拡張する、という考えで複素数に似た考え方です。

良く知られている実数の拡張方法の一つは、複素数です。
通常の実数に\(i^2=-1\)を満たす\(i\)という数を付加するのが複素空間です。

拡張の方法は何も\(i\)だけではありません。
例えば、\(\epsilon^2=0, (\epsilon\ne 0)\)を満たす\(\epsilon\)という数を付加することもできます。

この\(\epsilon\)は複素数ではありません
複素平面上でこの性質を満たす数は無いことからも、この\(\epsilon\)を追加するということは新しい方向への数の拡張です。
実際、複素平面上で\(\epsilon^2=0, (\epsilon\ne 0)\)を満たす数があるのかを探しても
\(
\begin{align}
(a+ib)^2=a^2-b^2+i2ab=0
\end{align}
\)

であるので、これを満たすのは\(a=b=0\)しか存在せず、複素平面上の数ではないことが分かります。

そんな新しい方向\(\epsilon\)の平面で定義された数をDual number(二重数、または双対数)と呼びます。

二重数\(a\)は実数部と非実数部から構成されており、
\(
a=a_0+a_1\epsilon
\)

と書くことが出来ます。ここで、\(a_0, a_1\)は実数であり、 \(\epsilon\)は、虚数単位\(i\)に倣って二重数単位とでも名付けておきましょう。

Dual numberと導関数


二重数は関数の導関数と大きな関係がある事を示しましょう。
テーラー展開を行います。
関数\(f\)の\(x\)周りの展開は刻み幅\(\Delta\)を用いて
\(\displaystyle
f(x+\Delta)=f(x)+f'(x) \Delta+\frac{1}{2!}f^{\prime\prime}(x)\Delta^2+\cdots
\)

と記述することが出来ます。もし仮に\(\Delta\)が純非実数だとしましょう。すなわち、\(\Delta=h\epsilon\)とします。ここで、\(h\)は実数です。テーラー展開の式に代入すれば、
\(
\begin{eqnarray}
f(x+h\epsilon)&=& f(x)+f'(x) h\epsilon+\frac{1}{2!}f^{\prime\prime}(x){h}^2\epsilon^2+\cdots \\
&=& f(x)+f'(x) h\epsilon
\end{eqnarray}
\)

となるわけです。注記しますが、上の式は\(h\)の1次で打ち切っているのではなく、厳密にイコールが成り立っています。
この式が言っているのは、関数\(f(a), (a\text{は二重数})\)を計算し、その非実数部を\(h\)で割ると関数の導関数になっているということです。

二重数の非実数部を取り出す関数\(\text{Dual}\)を定義します。
すると、導関数は
\(\displaystyle
f'(x)=\frac{1}{h}\text{Dual}(f(x+h\epsilon))
\)

として得られます。
言葉で書けば、二重数空間に拡張した\(f(x+h\epsilon)\)を計算すると、その非実数部に導関数が現れる、ということです。

二重数のプログラムはJeffrey Fikeさんによる[4]にありますので、そちらをご参考ください。

Hyper-dual numbers


さて、ここまでで二重数の概念を簡単に説明し、新しい実数の拡張を行いました。
Dual Numberのままでは高階導関数は得られません。なぜなら、導関数の二次以降は二重数の性質\(\epsilon^2=0\)によってゼロになるからです。

そこで、若干工夫します。
Dual Numberを二種類用意することで二階微分を得ることが出来ます[2,3]。
すなわち、
\(\begin{gather}
a=a_0+a_1\epsilon_1+a_2\epsilon_2+a_3\epsilon_1\epsilon_2\\
\epsilon_1^2=\epsilon_2^2=(\epsilon_1\epsilon_2)^2=0,~\epsilon_1\ne 0,~~\epsilon_2\ne 0,~~\epsilon_1\epsilon_2 \ne 0
\end{gather}
\)

と実数を拡張します。この様に拡張した数\(a\)をHyper-dual numbersと呼びます[2,3]。

Hyper-dual numbersの演算


Hyper-dual numbersである
\(
\begin{align}
a&=a_0+a_1\epsilon_1+a_2\epsilon_2+a_3\epsilon_1\epsilon_2\\
b&=b_0+b_1\epsilon_1+b_2\epsilon_2+b_3\epsilon_1\epsilon_2
\end{align}
\)

を用意します。和、積、商はそれぞれ
\(
\begin{align}
a+b&=(a_0+b_0)+(a_1+b_1)\epsilon_1+(a_2+b_2)\epsilon_2+(a_3+b_3)\epsilon_1\epsilon_2\\
ab&=a_0b_0+(a_0b_1+a_1b_0)\epsilon_1+(a_0b_2+a_2b_0)\epsilon_2+(a_0b_3+a_1b_2+a_2b_1+a_3b_0)\epsilon_1\epsilon_2\\
\frac{1}{a}&=\frac{1}{a_0}-\frac{a_1}{a_0^2}\epsilon_1-\frac{a_2}{a_0^2}\epsilon_2+\left(-\frac{a_3}{a_0^2}+\frac{2a_1a_2}{a_0}\right)\epsilon_1\epsilon_2
\end{align}
\)
と定義されます[2]。任意の関数は
\(
f(x)=f(x_0)+x_1f'(x_0)\epsilon_1+x_2f'(x_0)\epsilon_2
+\left(x_3f'(x_0)+x_1x_2f^{\prime\prime}(x_0)\right)\epsilon_1\epsilon_2
\)

として計算することが出来ます。なので、この結果から\(\epsilon_1\epsilon_2\)の係数として関数の二階微分が得られます。

Hyper-dual numbersのプログラム


実際にプログラムを組みましょう。
Hyper-dual numbersの型を持つ変数はFortran90では定義できません。
なので、構造体を利用して自分で型と、演算を定義します。

Fortran90のプログラムは以下の通りになるかと思います。
基本的な四則演算、基本的な初等関数のHyper-dual numbersの演算をモジュールとして書いています。

上のモジュールと下のメインプログラムを一緒にコンパイルすることにより、関数の二階微分が得られます。

例として
2変数関数
\(\displaystyle
f(x,y)=\frac{\text{ln}(xy^2)e^x}{\sqrt{\sin^3{x}+\sin^3{y}}}
\)

の偏微分
\(\begin{align}
\frac{\partial f}{\partial x},~~\frac{\partial f}{\partial y},~~\frac{\partial^2 f}{\partial x\partial y}
\end{align}
\)

を得ることを考えます。
\(
f(x+1\epsilon_1,y+1\epsilon_2)
\)

を計算すると
\(
f(x+1\epsilon_1,y+1\epsilon_2)=a_0+a_1\epsilon_1+a_2\epsilon_2+a_3\epsilon_1\epsilon_2
\)

のように実数係数\(a_0, a_1, a_2, a_3\)が得られます。
すると、
\(
\begin{align}
\frac{\partial f}{\partial x}=a_1 \\
\frac{\partial f}{\partial y}=a_2 \\
\frac{\partial^2 f}{\partial x\partial y} =a_3
\end{align}
\)

として偏微分が得られます。
ちなみに、二階微分が欲しい場合は
\(
f(x+1\epsilon_1+1\epsilon_2, y)
\)

を考えると
\(
\begin{align}
\frac{\partial f}{\partial x}=a_1 \\
\frac{\partial f}{\partial x}=a_2 \\
\frac{\partial^2 f}{\partial x^2} =a_3
\end{align}
\)

として得られます。\(a_1,~a_2\)はどちらを採用しても構いません。

プログラムでは変数の型Hyperdualを持つ入力変数を\(\text{xH,yH}\),出力を\(\text{rH}\)と置いています。

program main
  use Hyperdualmod
  implicit none
  type(Hyperdual)::xH,yH,rH
 
  xH%x0 = 0.3d0 ! real part
  xH%x1 = 1d0   ! unreal part \epsilon_1
  xH%x2 = 0d0   ! unreal part \epsilon_2
  xH%x3 = 0d0   ! unreal part \epsilon_1\epsilon_2

  yH%x0 = 0.4d0 ! real part
  yH%x1 = 0d0   ! unreal part \epsilon_1
  yH%x2 = 1d0   ! unreal part \epsilon_2
  yH%x3 = 0d0   ! unreal part \epsilon_1\epsilon_2

  write(6,'(4f23.16)')xH%x0,xH%x1,xH%x2,xH%x3
  write(6,'(4f23.16)')yH%x0,yH%x1,yH%x2,yH%x3  
 
  rH = log(xH*yH**2)*exp(xH)/sqrt(sin(xH)**3+cos(yH)**3)
  !rH = asin(2d0*xH)*acos(yH)/atan(xH*yH)
  !rH = xH**yH
 
  write(6,'(4f23.16)')rH%x0,rH%x1,rH%x2,rH%x3  

  stop
end program main

ヘッセ行列


二階偏微分の計算が出来たので、ヘッセ行列が簡単に計算できます。
ルーチンを作れば、以下の通りになるかと思います。
下のプログラムは3変数関数
\(\displaystyle
f(x,y,z)=\exp(xy)\tan(z)
\)

の\(x=-2,~y=3,~z=1\)におけるヘッセ行列を計算します。

program main
  use Hyperdualmod
  implicit none
  integer::N
  double precision::x,y,z,f
  double precision,allocatable::nabla(:),Hesse(:,:)
  double precision,allocatable::w(:)
  external::func

  N=3

  allocate(nabla(1:N),Hesse(1:N,1:N))
  nabla=0d0
  Hesse=0d0
  allocate(w(1:N))
  w=0d0

  x = -2d0
  y =  3d0
  z =  1d0
  w(1) = x
  w(2) = y
  w(3) = z
  call Hessian(N,w,func,nabla,Hesse)

  f=exp(x*y)*tan(z)
  write(6,'(1e24.16)')nabla(1)
  write(6,'(1e24.16)')y*f
  write(6,'(1e24.16)')nabla(2)
  write(6,'(1e24.16)')x*f
  write(6,'(1e24.16)')nabla(3)
  write(6,'(1e24.16)')exp(x*y)/(cos(z)**2)
  write(6,'(3e24.16)')Hesse(1,1:3)
  write(6,'(3e24.16)')f*y**2, f*(1+x*y), f*y/(cos(z)*sin(z))
  write(6,'(3e24.16)')Hesse(2,1:3)
  write(6,'(3e24.16)')f*(1+x*y), x**2*f, f*x/(cos(z)*sin(z))
  write(6,'(3e24.16)')Hesse(3,1:3)
  write(6,'(3e24.16)')f*y/(cos(z)*sin(z)), f*x/(cos(z)*sin(z)), f*2d0/(cos(z)**2)

  stop
end program main

subroutine func(N,x,f)
  use Hyperdualmod
  implicit none
  integer::N
  type(Hyperdual),intent(in)::x(1:N)
  type(Hyperdual),intent(out)::f
 
  f = exp(x(1)*x(2))*tan(x(3))
 
  return
end subroutine func

subroutine Hessian(N,x,func,nabla,Hesse)
  use Hyperdualmod
  implicit none
  integer,intent(in)::N
  double precision,intent(in)::x(1:N)
  double precision,intent(out)::nabla(1:N)
  double precision,intent(out)::Hesse(1:N,1:N)
  external::func

  integer::i,j,k
  type(Hyperdual),allocatable::xH(:)
  type(Hyperdual)::rH
 
  allocate(xH(1:N))
  do i=1,N
     xH(i)%x0 = 0d0
     xH(i)%x1 = 0d0
     xH(i)%x2 = 0d0
     xH(i)%x3 = 0d0
  enddo
 
  do i=1,N
     xH(i)%x0 = x(i)
     xH(i)%x3 = 0d0
  enddo

  do i=1,N
     do k=1,N
        xH(k)%x1 = 0d0
        xH(k)%x2 = 0d0
     enddo
     xH(i)%x1 = 1d0
     xH(i)%x2 = 1d0
     
     call func(N,xH,rH)
     nabla(i) = rH%x1
     Hesse(i,i) = rH%x3
  enddo

  do i=1,N
     do j=i+1,N
        do k=1,N
           xH(k)%x1 = 0d0
           xH(k)%x2 = 0d0
        enddo
        xH(i)%x1 = 1d0
        xH(i)%x2 = 0d0
        xH(j)%x1 = 0d0
        xH(j)%x2 = 1d0
        call func(N,xH,rH)
        Hesse(i,j) = rH%x3
        Hesse(j,i) = Hesse(i,j)
     enddo
  enddo  
 
  return
end subroutine Hessian

Hyperdual.f90にモジュールを入れ、メインプログラムをmain.f90に入れたとすると、以下の結果を得ます。

> gfortran Hyperdual.f90 main.f90  
> ./a.out
  0.1158128336233602E-01
  0.1158128336233602E-01
 -0.7720855574890680E-02
 -0.7720855574890680E-02
  0.8491012233306163E-02
  0.8491012233306162E-02
  0.3474385008700806E-01 -0.1930213893722670E-01  0.2547303669991849E-01
  0.3474385008700806E-01 -0.1930213893722670E-01  0.2547303669991849E-01
 -0.1930213893722670E-01  0.1544171114978136E-01 -0.1698202446661233E-01
 -0.1930213893722670E-01  0.1544171114978136E-01 -0.1698202446661233E-01
  0.2547303669991849E-01 -0.1698202446661233E-01  0.2644793608458059E-01
  0.2547303669991849E-01 -0.1698202446661233E-01  0.2644793608458058E-01

実行結果の奇数行目はHyper-dual Numberによる計算結果、偶数行目は解析解を表します。
また、6行目までは一階微分、7行目以降はヘッセ行列を表します。

参考文献


[1]関数
\(\displaystyle
f(x,y)=\frac{\text{ln}(xy^2)e^x}{\sqrt{\sin^3{x}+\sin^3{y}}}
\)

の偏微分\(\frac{\partial^2 f}{\partial x\partial y}\)
の計算
https://www.wolframalpha.com/input/?i=D%5Bln(x*y%5E2)*e%5E(x)%2Fsqrt(sin%5E3(x)%2Bcos%5E3(y)),+y,x%5D

[2]Jeffrey A. Fike and Juan J. Alonso,~”The Development of Hyper-Dual Numbers for Exact Second-Derivative Calculations”, 49th AIAA Aerospace Sciences Meeting(2011)
http://adl.stanford.edu/hyperdual/fike_aiaa-2011-886_slides.pdf,
J. A. Fike and J. J. Alonso. The Development of Hyper-Dual Numbers for Exact Second Derivative Calculations. AIAA paper 2011-886, 49th AIAA Aerospace Sciences Meeting, January 4-7, 2011.http://adl.stanford.edu/hyperdual/Fike_AIAA-2011-886.pdf

[3]JeffreyA.Fike,~”Derivative Calculations Using Hyper-Dual Numbers”, Sandia National Laboratories (2016)
https://www.osti.gov/servlets/purl/1368722

[4]Jeffrey Fike, Aerospace Design Lab, http://adl.stanford.edu/hyperdual/

[5]松本佳彦, 新しい数をつくる, (2018) http://ymatz.net/assets/docs/20180629-jtpa-slide-mod

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

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

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

まとめ

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

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

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

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

実際に、LAPACKのLU分解を使用したプログラムの例を載せておきます。

高次のニュートン法?


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コードは以下のもので、
https://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}\theta(x)e^{- ikx} dx=\int_{0}^{\infty}e^{-ikx}dx
\end{align}
\)

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

\(k\ne0\)の場合


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

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

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

\(k=0\)の場合


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

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

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

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

特異点の周り方の指定

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

と記述されます。

演算子としての表現


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

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

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

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

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

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

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

結論


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

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

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

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

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

図では

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

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

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

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

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

が成立します。

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


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

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

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

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

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

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

注意


これは私の解釈です。

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

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

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

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