ディガンマ関数の数値計算

ディガンマ関数\(\psi(z)\)は

\(
\displaystyle \psi(z)=\frac{d}{dz} \ln\Gamma(z)=\frac{\Gamma'(z)}{\Gamma(z)}
\)
として定義されます。

数値計算でこれを求めます。
用いる式は、
ディガンマ関数の漸近展開式[1],
\(
\displaystyle \psi(z)\sim \ln z -\frac{1}{2z}-\sum_{n=1}^{\infty} \frac{B_{2n}}{2nz^{2n}},\;\; z\to\infty (|{\rm arg} z| \lt \pi )
\)
ここで\(B_{2n}\)はベルヌーイ数[2]
と
漸化式[1]
\(
\displaystyle \psi(z+1)=\psi(z)+\frac{1}{z}
\)
相反公式[1]
\(
\psi(1-z)=\psi(z)+\pi \cot(\pi z)
\)
を利用します。

具体的には領域を4つに分けて、

\(10 \lt {\rm Re}(z)\), 漸近展開
\(1 \lt {\rm Re}(z) \lt 10\), 漸近展開+漸化式
\(-9 \lt {\rm Re}(z) \lt 1\), 漸近展開+漸化式+相反公式
\({\rm Re}(z) \lt -9\), 漸近展開+相反公式

を利用しています。

program main
implicit none
integer::i,j
double precision::x,y
complex(kind(0d0))::z,digamma
external::digamma

do i=0,200
x=-10.05d0+i*0.1d0
do j=0,200
y=-10.05d0+j*0.1d0
z=dcmplx(x,y)
write(10,'(4e20.10e2)')x,y,dble(digamma(z)),dimag(digamma(z))
enddo
write(10,*)
enddo

stop
end program main

!--------------------------

function digamma(z)
! digamma function for complex plane
! sikinote, http://slpr.sakura.ne.jp/qp/
! author : sikino
! date : 2016/07/04 (yyyy/mm/dd)
! date : 2016/07/07 (yyyy/mm/dd)
implicit none
complex(kind(0d0)),intent(in)::z
complex(kind(0d0))::digamma

integer::i,j,n,m,as,key
integer,parameter::nmax=100
double precision::eps=1d-13
double precision::pi=4d0*atan(1d0)
complex(kind(0d0))::c,w,t,ds,s
double precision::B
external::B

as=10

if(dble(z).ge.dble(as))then
w=z
key=0
elseif(dble(z).ge.1d0)then
w=z+dble(as)
key=1
elseif(dble(z).ge.-dble(as-1))then
w=1d0-z+dble(as)
key=2
else
w=1d0-z
key=3
endif

! Asymptotic expansion
s=log(w)-0.5d0/w
do n=1,nmax
m=2*n
ds=B(m)/(m*w**m)
s=s-ds
if(abs(ds)/abs(s).le.eps)exit
enddo
if(n.eq.nmax+1)then
write(6,*)"Worning! did not converge at digamma"
endif

if(key.eq.1.or.key.eq.2)then
!Recurrence relation
do i=1,as
s=s+1d0/(1d0-w)
w=w-1d0
enddo
endif
if(key.eq.2.or.key.eq.3)then
! Reflection Formula
s=s-pi/tan(pi*z)
endif

digamma=s

return
end function digamma
!--------------------------
function B(n)
! 1st bernoulli number B(n), 0 <= n <= 100
! from wolfram alpha
! " Table[bernoulli b(i), {i,0,100}] with 17 digits "
implicit none
integer,intent(in)::n
double precision::B,A(0:100)

if(n.lt.0.or.n.gt.100)then
write(6,*)"Not expected n at bernoulli number function. "
write(6,*)"program stop"
stop
endif

A(0:100)=(/1d0 &
, -0.50000000000000000d0, 0.16666666666666667d0, 0d0, -0.033333333333333333d0, 0d0 &
, 0.023809523809523810d0, 0d0, -0.033333333333333333d0, 0d0, 0.075757575757575758d0, 0d0 &
, -0.25311355311355311d0, 0d0, 1.1666666666666667d0, 0d0, -7.0921568627450980d0, 0d0 &
, 54.971177944862155d0, 0d0, -529.12424242424242d0, 0d0, 6192.1231884057971d0, 0d0 &
, -86580.253113553114d0, 0d0, 1.4255171666666667d+6, 0d0, -2.7298231067816092d+7, 0d0 &
, 6.0158087390064237d+8, 0d0, -1.5116315767092157d+10, 0d0, 4.2961464306116667d+11, 0d0 &
, -1.3711655205088333d+13, 0d0, 4.8833231897359317d+14, 0d0, -1.9296579341940068d+16, 0d0 &
, 8.4169304757368262d+17, 0d0, -4.0338071854059455d+19, 0d0, 2.1150748638081992d+21, 0d0 &
, -1.2086626522296526d+23, 0d0, 7.5008667460769644d+24, 0d0, -5.0387781014810689d+26, 0d0 &
, 3.6528776484818123d+28, 0d0, -2.8498769302450882d+30, 0d0, 2.3865427499683628d+32, 0d0 &
, -2.1399949257225334d+34, 0d0, 2.0500975723478098d+36, 0d0, -2.0938005911346378d+38, 0d0 &
, 2.2752696488463516d+40, 0d0, -2.6257710286239576d+42, 0d0, 3.2125082102718033d+44, 0d0 &
, -4.1598278166794711d+46, 0d0, 5.6920695482035280d+48, 0d0, -8.2183629419784576d+50, 0d0 &
, 1.2502904327166993d+53, 0d0, -2.0015583233248370d+55, 0d0, 3.3674982915364374d+57, 0d0 &
, -5.9470970503135448d+59, 0d0, 1.1011910323627978d+62, 0d0, -2.1355259545253501d+64, 0d0 &
, 4.3328896986641192d+66, 0d0, -9.1885528241669328d+68, 0d0, 2.0346896776329074d+71, 0d0 &
, -4.7003833958035731d+73, 0d0, 1.1318043445484249d+76, 0d0, -2.8382249570693707d+78/)

B=A(n)

return
end function B

[adsense2]

参考文献

[1]Abramowitz, M. and I.A. Stegun, {\it Handbook of Mathematical Functions}, Dover Publications, Ninth Printing, 1970., P.258-259, Eq.6.3.18(漸近展開式), Eq.6.3.5(漸化式), Eq.6.3.7(相反公式)
[2]Bernoulli Number -wolfram mathworld