関数論　学術図書出版　遠木幸成・阪井章先生(複素積分の数値計算法)

$z=u(x,y)+iv(x,y)$とするとき、$Cauchy-Riemann$の公式は
$$u_{x}(x,y)=v_{y}(x,y), u_{y}(x,y)=-v_{x}(x,y)$$
となることです。下記のコードは$x^{2}+2xyi-y^2$について計算したものです。

```r

#cauchy-Riemanの定理

f=function(x,y){
return(x^2+2*x*y*1i-y^2)  
}

h=0.0001

x=2;y=3

#xの微分
fx=(f(x+h,y)-f(x,y))/h

#yの微分
fy=(f(x,y+h)-f(x,y))/h

#一致するはず(=fx)
-1i*fy

```

複素領域での微分の応用として、ヤコビアンを使用したニュートン法です。

```r

#複素領域のニュートン法

f=function(x){
f1=3*x[1]+2*x[2]-1+3*1i
f2=5*x[1]+3*x[2]-4+5*1i
return(c(f1,f2))  
}

X=rep(1+1i,2)

ite=10^3
eta=0.01
h=0.001

for(l in 1:ite){

dF=array(0,dim=c(length(X),length(X))) 
for(j in 1:length(X)){
Xsub=X
Xsub[j]=Xsub[j]+h+h*1i
dF[,j]=(f(Xsub)-f(X))/(h+h*1i)
}  
 
X=X-eta*solve(dF)%*%f(X)
print(sum(f(X)))   
  
}

```

$Cauchy$の積分定理 : 関数$f(z)$は領域$D$において正則、単一閉曲線$C$は周も内部も$D$に属するとする。そのとき、
$$\int_{C}f(z)dz=0$$
を使って、$$\int_{C}\frac{e^{iz}}{z}dz$$を計算していきます。ここでは積分領域は$R1:\epsilon\rightarrow R$、$R2:半径Rの0\rightarrow 2\piの円弧$、$R3:-R\rightarrow -\epsilon$、$R4:半径\epsilonの2\pi\rightarrow 0の円弧$上で積分します。それぞれの積分領域は$\epsilon \rightarrow 0、R \rightarrow \infty$とするとそれぞれ$R1+R3$上の積分では$\pi i$、$R2$、$R4$では$0$となります。ここでは数値積分を行うため、
$$\int_{C}f(z)dz \approx \sum_{k:曲線上の座標の細分} f(z_{k})(z_{k+1}-z_{k})$$
を使って計算しました。こういった方法で計算できるとき$Stieltjes$積分可能といいます。




```r

#Cauchyの積分定理
#複素積分をRで計算する　p.46

ep=10^(-3)
R=1000

f=function(Z){
return(exp(Z*1i)/Z)  
}


#半径epの積分
n=10^5
sita=-seq(-pi,0,1/n)

z=ep*exp(sita*1i)
z2=z[2:length(z)]
z1=z[1:(length(z)-1)]
f1=sum(f(z2)*(z2-z1))

#半径Rの積分
n=10^5
sita=seq(0,pi,1/n)

z=R*exp(sita*1i)
z2=z[2:length(z)]
z1=z[1:(length(z)-1)]
f2=sum(f(z2)*(z2-z1))

#-R~-epの積分
n=10^6
z=-ep-(-ep-(-R))*seq(1,n,1)/n
z2=z[2:length(z)]
z1=z[1:(length(z)-1)]
f3=-sum(f(z2)*(z2-z1))

#R~epの積分
n=10^6
z=ep+(R-ep)*seq(1,n,1)/n
z2=z[2:length(z)]
z1=z[1:(length(z)-1)]
f4=sum(f(z2)*(z2-z1))

#cauchyの積分定理を確認
f1+f2+f3+f4

```

$Cauchy$の積分公式 : 単一閉曲線$C$の内部および周上で$f(z)$が正則ならば$C$の内部の1点$a$(極 : 非正則点)に対して

$$f(a)=\frac{1}{2\pi i} \int_{C}\frac{f(z)}{z-a}dz$$

を数値積分によって計算します。$f(z)=e^{az}(aは任意の数)$とします。

$$\int_{|z|=1}\frac{e^{az}}{z}dz=2\pi i(\frac{1}{2\pi i}\int_{|z|=1}\frac{e^{az}}{z-0}dz)=2\pi ie^{0}=2\pi i$$

となることを同じく数値積分により確認してみました。



```r

#Cauchyの積分公式
#原点で極がある場合
a=2
f=function(z){
return(exp(a*z)/z)
}


#|z|=1、すなわち単位円盤上の閉曲線で積分する

#半径Rの積分
n=10^5
sita=seq(0,2*pi,1/n)
R=1
z=R*exp(sita*1i)
z2=z[2:length(z)]
z1=z[1:(length(z)-1)]
#積分値
fval=sum(f(z2)*(z2-z1))



```

また、次の定理により導関数もすぐに求まります：
単一閉曲線上$C$の内部および周上で$f(z)$が正則ならば$C$の内部の1点$a$に対して
$$f^{(n)}(a)=\frac{n!}{2\pi i} \int_{C}\frac{f(z)}{(z-a)^{n+1}}dz$$
で求まる。ただし、$f^{(n)}(z)$は$n$階の導関数。<font color="red">この方法でも実関数の高次微分が計算可能。</font>

下記プログラムでは$$\frac{(n-1)!}{2\pi i} \int_{|z|=1}\frac{e^{bz}}{(z-0)^n}dz$$
を数値積分により計算しています。ちなみに

$$\int_{|z|=1}\frac{e^{bz}}{z^n}dz=\frac{2\pi ib^{n-1}}{(n-1)!}$$

です。Laurent展開の係数はこの導関数により求まるので、数値計算でも計算可能だと思われます(無限級数なのでどこまで展開するかによるが)。

```r

#n階微分の導関数

a=2
N=1 #微分回数
A=0 #極
f=function(z){
value=exp(a*z)  
return(value/((z-A)^N))  
}

#|z|=1、すなわち単位円盤上の閉曲線で積分する

#半径Rの積分
n=10^5
sita=seq(0,2*pi,1/n)
R=1
z=R*exp(sita*1i)
z2=z[2:length(z)]
z1=z[1:(length(z)-1)]
fval=sum(f(z2)*(z2-z1))

#解析解(fval)
if(N>1){
2*pi*1i*(a^(N-1))/(prod(c(1:(N-1))))
}else{
2*pi*1i*(a^(N-1))
}

#実際の導関数値は
if(N>1){
fval*(prod(c(1:(N-1)))/(2*pi*1i))
}else{
fval*(1/(2*pi*1i))  
}

```

$Taylor$展開の係数は$f(z)=\sum_{n=0}^{\infty}a_n(z-a)^n$に対し、

$$a_n=\frac{1}{2\pi i}\int_{C}\frac{f(z)}{(z-a)^{n+1}}dz$$

として求められる。$e^z$の係数は$a_n=\frac{1}{n!}$であることを数値計算により求めた。

```r

#Taylor展開の係数

coefficients=c()

for(k in 1:10){

#n階微分の導関数  
  
N=k+1 #微分回数
A=0 #極
f=function(z){
value=exp(z)  
return(value/((z-A)^N))  
}

#正則領域は|z|<inf. |z|=1、すなわち単位円盤上の閉曲線で積分する

#半径Rの積分
n=10^5
sita=seq(0,2*pi,1/n)
R=1
z=R*exp(sita*1i)
z2=z[2:length(z)]
z1=z[1:(length(z)-1)]
fval=sum(f(z2)*(z2-z1))

df=fval*(1/(2*pi*1i))  
coefficients=c(coefficients,df)

}

#exp(z)のTaylor展開の係数が順に計算される(1/n!)
coefficients


```

次は留数を数値計算で求めてみる。留数とは、

$f(z)$が$0<|z-a|<r$で正則であるとき、$0<p<r$に対して、

$$\frac{1}{2\pi i}\int_{|z-a|=p}f(z)dz$$

を$z=a$における$f(z)$の留数という。数値計算では、

$$\int_{|z-a|=p}f(z)dz \approx \sum_{k:曲線上の座標の細分} f(z_{k})(z_{k+1}-z_{k})$$

$z_{k}=a+pe^{i\theta_k},  \space \space \theta_kは0\sim 2\piの間に値をとる。$閉曲線上に細分を作って計算する。積分領域の半径 $p$ を決めないと計算できないが、収束半径を超えないよう、十分小さくとればいいと思う。

```r

#留数定理について

f=function(z){
return((z^2+z+1)/((z-1)*z^2))  
}

#収束半径の中に積分領域の半径が収まっていればいいから十分小さくとる(0.5)
#気になる場合は収束半径を計算してもいい

#半径Rの積分
#極がz0=1のとき
z0=1
n=10^5
sita=seq(0,2*pi,1/n)
R=0.5
z=z0+R*exp(sita*1i)
z2=z[2:length(z)]
z1=z[1:(length(z)-1)]
fval=sum(f(z2)*(z2-z1))
#z=1での留数の値
fval/(2*pi*1i)

```