友意白雑記帳

だいたい自分用の覚え書き

複素数積分の数値計算用コード(.f90)

ずいぶん昔にお遊びで組んだコードとノートを発掘した。せっかくなので公開。

ソースコード↓(Fortran 90)

********************************

program MAIN !PROGRAM to integrate the complex function on the closed-circle contour.
implicit none
integer :: Nct = 100 != # of points of the closed-circle contour.
complex(kind(0d0)) :: ai = (0.d0, 1.d0)
complex(kind(0d0)) :: zp = (3.d0, -2.d0) != pole of the function.
double precision :: pi = 3.14159265359d0

integer :: i
double precision :: r, dth, theta
complex(kind(0d0)) :: W, z1, z0
complex(kind(0d0)), dimension(:), allocatable :: CS !contour coordinates.
complex(kind(0d0)), external :: fc !function of interest.

10034 allocate(CS(0:Nct)) ; CS = (0.d0, 0.d0)

dth = 2.d0*pi/dble(Nct)
r = 9.d0
CS(0) = cmplx(r)
do i = 1, Nct-1
theta = dble(i)*dth
CS(i) = r*cos(theta) + ai*r*sin(theta) !"x+iy" of the i-th contour point.
end do
CS(Nct) = CS(0) !to close the contour.
CS = CS + zp !origin ==> pole for the contour.
!! do i = 0, Nct
!! write(6,*) CS(i)
!! end do


W = (0.d0, 0.d0) !RESULT = integral value.
do i = 1, Nct
z0 = CS(i-1)
z1 = CS(i)
!! W = W + fc(z1,zp) * (z1-z0) !==> WRONG. Imag(W) gets (crazily) a finite value.
W = W + 0.5d0*(fc(z1,zp) + fc(z0,zp)) * (z1-z0)
end do

W = W*(-ai)*0.5d0/pi
write(6,"(' Real(W) & Imag(W) /(2*pi*i) =',2(3X,G19.11))") real(W), aimag(W)

20034 deallocate(CS)
end program MAIN


!--- Complex function to be integrated.
pure function fc(z,p) result (ff)
implicit none
complex(kind(0d0)), intent(IN) :: z, p
complex(kind(0d0)) :: ff
ff = z-p
!! ff = ff*ff
ff = 1.d0 / ff
!! ff = ff**3 ; ff = 1.d0 / ff
end function fc

********************************

 

そしてこちらがノート。テキトーにも程がある、、、。

f:id:Tomoishiro:20210824235546j:plain