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

ディガンマ関数\(\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つに分けて、

  1. \(10 \lt {\rm Re}(z)\),      漸近展開
  2. \(1 \lt {\rm Re}(z) \lt 10\),    漸近展開+漸化式
  3. \(-9 \lt {\rm Re}(z) \lt 1\),    漸近展開+漸化式+相反公式
  4. \({\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
スポンサーリンク

参考文献

[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


「ディガンマ関数の数値計算」への4件のフィードバック

  1. ガンマ関数やポリガンマ関数は助変数に対する評価が多いので、メモ化をしてしまうのもひとつだと思います。

    1. 「助変数に対する評価」や「メモ化」とは何を指しているんでしょうか?具体的なアイデアをお教えください。

  2. 例えば、超幾何関数やベッセル関数は展開式に助変数のΓ関数を含みます。
    引数をキーとして関数値をペアにしてmapを作ってしまえば、同じ引数に対する計算は表引きだけで済む様になります。また差が整数の引数の場合には漸化式の結果が直接に関数値になります。特に精度を重視するときは効果大です。

  3. 正直なところ、何の話をしているのか全く分かりません。ここはディガンマ関数についてのページです。
    ディガンマ関数に限った場合、上記プログラムのどの文章辺りで、もしくは数式のどこで「メモ化」なる工夫が使えるのでしょうか?お教え願えればと思います。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です