「プログラミングと数値計算」カテゴリーアーカイブ

ガウス=クロンロッド求積法

2016年7月22日 sikino コメントする

ガウス=クロンロッド求積法と呼ばれる数値積分法です。

ガウス=クロンロッド求積法の説明
fortran90によるプログラム
参考文献

説明

ガウス=ルジャンドル求積法の補助的な存在、というとらえ方でいいと思います。
通常はガウス=ルジャンドル求積法の誤差推定のために使われるらしいです。

積分
$
\displaystyle \int_{-1}^{1} f(x) dx
$
を、
$
\displaystyle \int_{-1}^{1} f(x) dx \approx \sum_{i=1}^{2N+1} \omega_i f(x_i)+\mathcal{O}(x^{3N+1})
$
として近似します。ここで、$\omega_i$は点$x_i$での重みを意味します。
$\omega_i, x_i$は任意に決められず、$N$をいくつにするかによって強制的に決定されます。

ガウス=クロンロッド求積法は次数($2N+1$)を用いたとき、被積分関数が
$x^{3N+1}$以下の多項式であれば厳密な積分値を返しますCHAPTER 5 Numerical Quadrature。
Nは次数。分点数と一致。

重要な事ですが、被積分関数が多項式ではない場合、あまりよくない積分結果になってしまいます。
例えば、積分
$
\displaystyle \int_0^b \sqrt{x} dx
$
を考えると、被積分関数が多項式で書けないので、精度はガクッと落ちます。
この理由は、被積分関数の微分が$x=0$で発散するためであり、高次の積分法になるほど顕著に表れます。

また、ルンゲ現象に代表される、高次で発生する様々な問題を排除するため、通常は積分区間を複数に分けて、低次の積分法を組み合わせて計算を行います。

[adsense1]

fortran90によるプログラム

以下のプログラムは15次ガウス=クロンロッド求積法のプログラムです。
これは、積分
$
\displaystyle \int_{-1}^{1}e^{-x} dx = 2\sinh(1) \approx 2.350402387287602913764764
$
を計算するプログラムとなっています。

program main
implicit none
integer::i,N
double precision,allocatable::x(:),w(:)
double precision::f,s
external::f

N=15
allocate(x(1:N),w(1:N)); x=0d0; w=0d0
call GaussKronrod15(x,w)
s=0d0
do i=1,N
s=s+w(i)*f(x(i))
enddo
write(6,*)s

deallocate(x,w)

stop
end program main

function f(x)
implicit none
double precision,intent(in)::x
double precision::f

f=exp(-x)

return
end function f

subroutine GaussKronrod15(x,w)
!Gauss-Kronrod Quadrature Nodes and Weights
!http://www.advanpix.com/2011/11/07/gauss-kronrod-quadrature-nodes-weights/
implicit none
double precision,intent(out)::x(1:15),w(1:15)

integer::i
integer,parameter::N=15

x=0d0; w=0d0

x( 8) = 0d0
x( 9) = 2.077849550078984676006894037732449d-1
x(10) = 4.058451513773971669066064120769615d-1
x(11) = 5.860872354676911302941448382587296d-1
x(12) = 7.415311855993944398638647732807884d-1
x(13) = 8.648644233597690727897127886409262d-1
x(14) = 9.491079123427585245261896840478513d-1
x(15) = 9.914553711208126392068546975263285d-1
do i=1,7
x(i)=-x(N-i+1)
enddo

w( 8) = 2.094821410847278280129991748917143d-1
w( 9) = 2.044329400752988924141619992346491d-1
w(10) = 1.903505780647854099132564024210137d-1
w(11) = 1.690047266392679028265834265985503d-1
w(12) = 1.406532597155259187451895905102379d-1
w(13) = 1.047900103222501838398763225415180d-1
w(14) = 6.309209262997855329070066318920429d-2
w(15) = 2.293532201052922496373200805896959d-2
do i=1,7
w(i)=w(N-i+1)
enddo

return
end subroutine GaussKronrod15

実行結果は
2.35040238728760
となり、厳密値
2.350402387287602913
と全桁一致します。

区間変更

定義そのままでは$[-1\sim 1]$なので、区間を変更しましょう。
区間変更は$y=\frac{b-a}{2}x+\frac{b+a}{2}$と変数変換をすれば、

$
\begin{align}
\int_a^b f(x) dx &= \int_{-1}^{1} f(y) \frac{b-a}{2}dy \\
&\approx \sum_{i=1}^N \left(\frac{b-a}{2}\omega_i\right) f(\frac{b-a}{2}x+\frac{b+a}{2}) \\
&= \sum_{i=1}^N \omega_i’ f(x_i’)
\end{align}
$
と掛けます。これをプログラムします。
ここでは、積分
$
\displaystyle \int_{1}^{2.4}e^{-x} dx = e^{-1}-e^{-2.4} \approx 0.27716148788202981
$
を計算します。

program main
implicit none
integer::i,N
double precision,allocatable::x(:),w(:)
double precision::f,s
external::f

N=15
allocate(x(1:N),w(1:N)); x=0d0; w=0d0
call GaussKronrod15ab(1d0,2.4d0,x,w)
s=0d0
do i=1,N
s=s+w(i)*f(x(i))
enddo
write(6,*)s

deallocate(x,w)

stop
end program main

function f(x)
implicit none
double precision,intent(in)::x
double precision::f

f=exp(-x)

return
end function f

subroutine GaussKronrod15ab(a,b,x,w)
!Gauss-Kronrod Quadrature Nodes and Weights
!http://www.advanpix.com/2011/11/07/gauss-kronrod-quadrature-nodes-weights/
implicit none
double precision,intent(in)::a,b
double precision,intent(out)::x(1:15),w(1:15)

integer::i
integer,parameter::N=15

x=0d0; w=0d0

x( 8) = 0d0
x( 9) = 2.077849550078984676006894037732449d-1
x(10) = 4.058451513773971669066064120769615d-1
x(11) = 5.860872354676911302941448382587296d-1
x(12) = 7.415311855993944398638647732807884d-1
x(13) = 8.648644233597690727897127886409262d-1
x(14) = 9.491079123427585245261896840478513d-1
x(15) = 9.914553711208126392068546975263285d-1
do i=1,7
x(i)=-x(N-i+1)
enddo

w( 8) = 2.094821410847278280129991748917143d-1
w( 9) = 2.044329400752988924141619992346491d-1
w(10) = 1.903505780647854099132564024210137d-1
w(11) = 1.690047266392679028265834265985503d-1
w(12) = 1.406532597155259187451895905102379d-1
w(13) = 1.047900103222501838398763225415180d-1
w(14) = 6.309209262997855329070066318920429d-2
w(15) = 2.293532201052922496373200805896959d-2
do i=1,7
w(i)=w(N-i+1)
enddo

x=0.5d0*((b-a)*x+(a+b))
w=0.5d0*(b-a)*w

return
end subroutine GaussKronrod15ab

となります。

実行結果は
0.277161487882030
となり、厳密値
0.27716148788202981
と全桁一致します。数値計算なので、最後のところは四捨五入されています。

61次ガウス=クロンロッド求積法

61次ガウス=クロンロッド求積法のコードも置いておきます。
これほどの高次は余り使わず、低次のものを組み合わせたほうが良いことを再び書いておきます。

積分
$
\displaystyle \int_{-3}^{20}e^{-x} dx = e^{3}-e^{-20} \approx 20.0855369211265141
$
を計算します。

program main
implicit none
integer::i,N
double precision,allocatable::x(:),w(:)
double precision::f,s
external::f

N=61
allocate(x(1:N),w(1:N)); x=0d0; w=0d0
call GaussKronrod61ab(-3d0,20d0,x,w)
s=0d0
do i=1,N
s=s+w(i)*f(x(i))
enddo
write(6,*)s

deallocate(x,w)

stop
end program main

function f(x)
implicit none
double precision,intent(in)::x
double precision::f

f=exp(-x)

return
end function f

subroutine GaussKronrod61ab(a,b,x,w)
!Gauss-Kronrod Quadrature Nodes and Weights
!http://www.advanpix.com/2011/11/07/gauss-kronrod-quadrature-nodes-weights/
implicit none
double precision,intent(in)::a,b
double precision,intent(out)::x(1:61),w(1:61)

integer::i
integer,parameter::N=61

x=0d0; w=0d0

x(31:61) = &
(/0.000000000000000000000000000000000d0 &
, 5.147184255531769583302521316672257d-2 &
, 1.028069379667370301470967513180006d-1 &
, 1.538699136085835469637946727432559d-1 &
, 2.045251166823098914389576710020247d-1 &
, 2.546369261678898464398051298178051d-1 &
, 3.040732022736250773726771071992566d-1 &
, 3.527047255308781134710372070893739d-1 &
, 4.004012548303943925354762115426606d-1 &
, 4.470337695380891767806099003228540d-1 &
, 4.924804678617785749936930612077088d-1 &
, 5.366241481420198992641697933110728d-1 &
, 5.793452358263616917560249321725405d-1 &
, 6.205261829892428611404775564311893d-1 &
, 6.600610641266269613700536681492708d-1 &
, 6.978504947933157969322923880266401d-1 &
, 7.337900624532268047261711313695276d-1 &
, 7.677774321048261949179773409745031d-1 &
, 7.997278358218390830136689423226832d-1 &
, 8.295657623827683974428981197325019d-1 &
, 8.572052335460610989586585106589439d-1 &
, 8.825605357920526815431164625302256d-1 &
, 9.055733076999077985465225589259583d-1 &
, 9.262000474292743258793242770804740d-1 &
, 9.443744447485599794158313240374391d-1 &
, 9.600218649683075122168710255817977d-1 &
, 9.731163225011262683746938684237069d-1 &
, 9.836681232797472099700325816056628d-1 &
, 9.916309968704045948586283661094857d-1 &
, 9.968934840746495402716300509186953d-1 &
, 9.994844100504906375713258957058108d-1 /)

do i=1,30
x(i)=-x(N-i+1)
enddo
w(31:61)= &
(/5.149472942945156755834043364709931d-2 &
, 5.142612853745902593386287921578126d-2 &
, 5.122154784925877217065628260494421d-2 &
, 5.088179589874960649229747304980469d-2 &
, 5.040592140278234684089308565358503d-2 &
, 4.979568342707420635781156937994233d-2 &
, 4.905543455502977888752816536723817d-2 &
, 4.818586175708712914077949229830459d-2 &
, 4.718554656929915394526147818109949d-2 &
, 4.605923827100698811627173555937358d-2 &
, 4.481480013316266319235555161672324d-2 &
, 4.345253970135606931683172811707326d-2 &
, 4.196981021516424614714754128596976d-2 &
, 4.037453895153595911199527975246811d-2 &
, 3.867894562472759295034865153228105d-2 &
, 3.688236465182122922391106561713597d-2 &
, 3.497933802806002413749967073146788d-2 &
, 3.298144705748372603181419101685393d-2 &
, 3.090725756238776247288425294309227d-2 &
, 2.875404876504129284397878535433421d-2 &
, 2.650995488233310161060170933507541d-2 &
, 2.419116207808060136568637072523203d-2 &
, 2.182803582160919229716748573833899d-2 &
, 1.941414119394238117340895105012846d-2 &
, 1.692088918905327262757228942032209d-2 &
, 1.436972950704580481245143244358001d-2 &
, 1.182301525349634174223289885325059d-2 &
, 9.273279659517763428441146892024360d-3 &
, 6.630703915931292173319826369750168d-3 &
, 3.890461127099884051267201844515503d-3 &
, 1.389013698677007624551591226759700d-3/)
do i=1,30
w(i)=w(N-i+1)
enddo

x=0.5d0*((b-a)*x+(a+b))
w=0.5d0*(b-a)*w

return
end subroutine GaussKronrod61ab

実行結果
20.0855369211265
厳密値
20.0855369211265141

[adsense2]

fortran90で、任意の次数の求積点を知りたければ、f90のプログラムが
KRONROD Gauss-Kronrod Quadrature Rules -Source Codes in
Fortran90
にて公開されています。

参考文献

Gauss-Kronrod Quadrature Nodes and Weights

プログラミングと数値計算

unformattedによるファイル入出力(fortran90)

2016年7月10日 sikino コメントする

fortran90でファイルの入出力に関する設定です。

fortranでデータの書き出しには、
書式付きと書式なしがあります。

書式付きはそこでプログラムがすべて終わる時
書式なしはそのデータを使って別のプログラムを動かす時
と使い分けられると思います。

[adsense1]

書式付き

利点
・テキストデータなので、結果を人間がわかりやすい形で出力出来る。
・ミスを発見しやすい
欠点
・データの書き出しに時間がかかる
・データ量が多くなる(書式なしの約3倍)

書式付きは、例えば

open(10,file="test.d")
write(10,'(A,i5,e15.5e2)')"aaa",j,alpha

見たいに書式を指定して書ださせることです。

書式なし

利点
・書き出し時間を短くできる
・データ量が少なくなる(書式付の約1/3倍)
・簡潔にコードが書ける
欠点
・バイナリファイルなので出力結果を直接見ることが出来ない。

書式なしのプログラムでは、

open(10,file="test.bin",form="unformatted")
write(10)"aaa",j,alpha

の様に書きます。

実際のプログラムでは、

program main
implicit none
integer,parameter::N1=4
integer,parameter::N2=7

integer::i,j
double precision::A(1:N1),B(1:N2)
complex(kind(0d0))::CC(1:N1,1:N2)

do i=1,N1
A(i)=100d0+dble(i)
enddo
do j=1,N2
B(j)=dble(j)
enddo
do i=1,N1
do j=1,N2
CC(i,j)=A(i)+B(j)
enddo
enddo

!q open(21,file="unformatted.bin",form="unformatted")
!q write(21)N1,N2
!q do i=1,N1
!q do j=1,N2
!q write(21)A(i),B(j),CC(i,j)
!q enddo
!q enddo
!q close(21)

open(21,file="unformatted.bin",form="unformatted")
write(21)N1,N2
write(21)A,B,CC
close(21)

stop
end program main

といった感じです。

書式なしの読み込み

書式なしのファイルを読み込むためにはまったく同じ形で読み込ませればいいです。
例えば、上の形で書いた場合、読み込みは

program main
implicit none
integer::i,j,N1,N2

double precision,allocatable::A(:),B(:)
complex(kind(0d0)),allocatable::CC(:,:)

open(21,file="unformatted.bin",form="unformatted")
read(21)N1,N2
allocate(A(1:N1),B(1:N2),CC(1:N1,1:N2))
!q do i=1,N1
!q do j=1,N2
!q read(21) A(i),B(j),CC(i,j)
!q enddo
!q enddo
read(21) A,B,CC
close(21)

do i=1,N1
do j=1,N2
write(6,'(4f10.5)')A(i),B(j),CC(i,j)
enddo
enddo

stop
end program main

とすればいいです。
”!q”のコメントアウトは２つのプログラムに共通しています。
なのでどういう書き方でも、同じ入出力の形を指定してあげれば良いです。

[adsense2]

プログラミングと数値計算, 数学

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

2016年7月4日 sikino 4件のコメント

ディガンマ関数$\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

プログラミングと数値計算, 数学

ベルヌーイ数の数値計算

2016年7月3日 sikino コメントする

これがベルヌーイ数を求める最速のアルゴリズムだと思います。
直接代入。これに限ります。

n=0~100までサポートしています。
実用的にはこの範囲で足りるんじゃないかなと思います。

値はwolfram alphaより。
Table[bernoulli b(i), {i,0,100}] with 17 digits

program main
implicit none
integer::i
double precision::B
external::B

do i=0,30
write(6,'(i5,f30.17)')i,B(i)
enddo

stop
end program main

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
double precision::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

B(256)まで増やしたのはこちら。

一応逐次的に求める事が出来るアルゴリズムがありまして実装しましたが、全然良くないです[1]。
4倍精度演算を使って$B_{30}$を求めようとするとたったの8桁しか一致しません。
倍精度演算のみだと$B_{14}$を求めようとするとたったの6桁しか一致しません。

▼クリックでこの場に展開

参考文献

[1]Algorithmic description -Bernoulli number, wikipedia

[adsense2]

プログラミングと数値計算, 数学

リーマンのゼータ関数の数値計算コード(複素平面)

2016年6月2日 sikino 8件のコメント

ここではゼータ関数を数値計算で得るためのアルゴリズムとfortran90によるプログラムを載せています。

ゼータ関数の計算方法

ゼータ関数は通常
・$\rm{Re}(s)\gt 1$に対しては直接計算を行い、$\rm{Re}(s)\le 1$の領域に対しては相反公式
$
\displaystyle \zeta(s)=2^{s}\pi^{s-1}\sin\left(\frac{\pi s}{2}\right)\Gamma(1-s)\zeta(1-s)
$
を使って求める事が推奨されます。
ガンマ関数を使うことを避けたい場合、全複素平面($s\ne 1$)で有効なゼータ関数の表示
$
\displaystyle
\zeta(s)=\frac{1}{1-2^{1-s}}
\sum_{n=0}^{\infty} \frac{1}{2^{n+1}}
\sum_{k=0}^{n} {n \choose k} \frac{(-1)^k}{(k+1)^s}
$
を使うのが良いでしょう。

ゼータ関数の級数表示

ゼータ関数は
$
\displaystyle
\zeta(s)=\sum_{n=1}^{\infty} \frac{1}{n^s} ,~~~(1\lt s)
$
と定義されます。が、これはlogの収束であるため非常に収束が遅いです。
Euler-Knopp変換
$\displaystyle
\sum_{n=0}^\infty a_n = \sum_{n=0}^{\infty} \frac{1}{(1+q)^{n+1}}
\sum_{j=0}^{n}{n \choose j} q^{n-j}a_j
$
(q=1のときEuler変換)によって早く収束する級数に変換されます[4]。
より詳しい証明は[5]。

Riemann Zeta Function — from Wolfram MathWorld, Eq.(21)
によると、全領域(s=1以外)で収束するリーマンのゼータ関数の級数表示は、
$
\displaystyle
\zeta(s)=\frac{1}{1-2^{1-s}}
\sum_{n=0}^{\infty} \frac{1}{2^{n+1}}
\sum_{k=0}^{n} {n \choose k} \frac{(-1)^k}{(k+1)^s}
$
と書けます。もう一つ形式があるのですが、そちらは収束が遅そうな形式だったので用いていません。

数値計算的には2の乗数とかはあまり計算したくありません。そこで、漸化式を考えます。
$
\begin{align}
\zeta(s) &=\sum_{n=0}^{\infty} b_n \\
b_n &= \sum_{k=0}^{n}a_k^n \\
a_k^n &= \frac{1}{1-2^{1-s}}\frac{1}{2^{n+1}}{n \choose k}\frac{(-1)^k}{(k+1)^s}
\end{align}
$
と$b_n, a_k^n$を定義します。$a_k^n\;\; (n\ge k)$を計算するために漸化式を考えると

$
\begin{align}
a_k^n&=\frac{1}{2}\left(\frac{n}{n-k}\right)a_{k}^{n-1} \\
a_k^n&=\frac{k-n-1}{k}\left(\frac{k}{k+1}\right)^s a_{k-1}^n \\
a_0^0&=\frac{1}{2}\frac{1}{1-2^{1-s}}
\end{align}
$

と導く事が出来ます。
計算手順としては、

としていけばよいでしょう。
計算量を考えて赤の手順を少なくするようにしています。

fortran90によるプログラム

本当のゼータ関数の値との相対誤差$\delta$は約
$\delta\lt 3\times 10^{-14} (Re(s)\lt 1)$,
$\delta\lt 3\times 10^{-15} (1 \lt Re(s)) $
負の領域の方が精度が悪いのは、相反公式から来ています。
相反公式のみを4倍精度演算で行えば
$\delta\lt 7\times 10^{-15} (Re(s) \lt 1) $
にまで精度が上がります。目立って精度はあがりません。

2017/01/12に大幅な更新をしました。
以前のものと比べ変更点は、
・4倍精度を使用していない
・計算速度は以前の54倍
になっているので、以前のものを使っている方は更新することを推奨します。
ガンマ関数に関しては複素ガンマ関数の数値計算をご参照ください。

function zeta(s)
implicit none
complex(kind(0d0)),intent(in)::s
complex(kind(0d0))::zeta
!
!Checked : ifort (success : version 16.0.2 20160204)
! : gfortran(success : version 4.8.4)
! : gfortran(success : version 4.4.7)
!
!http://slpr.sakura.ne.jp/qp/
! Author : sikino
! date: 2016/06/02, (yyyy/mm/dd)
! : 2016/06/03
! : 2016/06/04
! : 2016/06/05
! : 2016/12/25
! : 2017/01/12
! : 2017/01/14
!
integer::n,k
integer,parameter::nmax=1000
double precision::nq,nq1,ln1,ln2
double precision,parameter::eps=1d-15
double precision,parameter::pi=3.14159265358979324d0
complex(kind(0d0))::ds,is,d,s0,tc,tt,qzeta
complex(kind(0d0))::a(0:Nmax),b(0:nmax),cgamma
external::cgamma
!complex*32::zgamma
!external::zgamma

ds=dble(s); is=imag(s)
if(ds.eq.1d0.and.is.eq.0d0)then
zeta=1d+300
elseif(ds.eq.0d0.and.is.eq.0d0)then
zeta=-0.5d0
else
if(real(s).lt.1d0)then
s0=1d0-s
else
s0=s
endif

tt=1d0/(1d0-2d0**(1d0-s0))
qzeta=0.5d0*tt

a=cmplx(0d0,0d0,kind=8); b=cmplx(0d0,0d0,kind=8)
b(0)=0.25d0*tt
do n=1,nmax
nq=n
nq1=nq+1d0

a(0)=b(0)
do k=1,n-1
a(k)=b(k)*nq/(n-k)
enddo
!----
ln2=log(dble(n)/dble(n+1))
if(n.eq.1)then
d=s0*ln2
else
d=d*ln2/ln1
endif
a(n)=-a(n-1)*exp(d)/nq
ln1=ln2
!----

tc=cmplx(0d0,0d0,kind=8)
do k=0,n
tc=tc+a(k)
enddo

qzeta=qzeta+tc

if(abs(qzeta).ge.1d0)then
if(abs(tc/qzeta).le.eps)exit
else
if(abs(tc).le.eps)exit
endif

b=0.5d0*a
enddo
if(n.eq.nmax+1)write(6,*)" ***worning, didn't converge at zeta"

if(real(s).lt.1d0)then
zeta=(2d0**s)*pi**(-s0)*sin(0.5d0*pi*s)*qzeta*cgamma(s0)
!zeta=cmplx((2q0**s)*acos(-1q0)**(-s0)*sin(0.5q0*acos(-1q0)*s)*qzeta*zgamma(s0*1q0),kind=8)
else
zeta=qzeta
endif
end if

return
end function zeta

function cgamma(z)
! sikinote (http://slpr.sakura.ne.jp/qp)
! author : sikino
! date : 2017/01/09
! 2017/01/10

implicit none
complex(kind(0d0)),intent(in)::z
complex(kind(0d0))::cgamma

integer::i,n,M
double precision,parameter::pi=3.14159265358979324d0
double precision,parameter::e1=0.3678794411714423216d0
double precision,parameter::ln2pi2=0.91893853320467274d0
double precision,parameter::sq2pi=2.5066282746310005d0
double precision::a(1:30),AB(1:14),dz,iz
complex(kind(0d0))::q,s,w,r,q1,q2

dz=dble(z)
iz=imag(z)

if(dz.eq.nint(dz).and.iz.eq.0d0)then
if(dz.gt.0d0)then
s=dcmplx(1d0,0d0)
n=nint(dz)
do i=2,n
s=s*i
enddo
else
s=1d+300
endif
cgamma=s
else
q=z
if(dz.lt.0d0)q=1d0-z

if(abs(iz).lt.1.45d0)then
! For x=arb, |y|\lt 1.45
n=int(dble(q))
s=1d0
do i=1,n
s=s*(q-i)
enddo
w=q-dble(n)

a(1) = 1d0
a(2) = 0.57721566490153286d0
a(3) =-0.65587807152025388d0
a(4) =-0.42002635034095236d-1
a(5) = 0.16653861138229149d0
a(6) =-0.42197734555544337d-1
a(7) =-0.96219715278769736d-2
a(8) = 0.72189432466630995d-2
a(9) =-0.11651675918590651d-2
a(10)=-0.21524167411495097d-3
a(11)= 0.12805028238811619d-3
a(12)=-0.20134854780788239d-4
a(13)=-0.12504934821426707d-5
a(14)= 0.11330272319816959d-5
a(15)=-0.20563384169776071d-6
a(16)= 0.61160951044814158d-8
a(17)= 0.50020076444692229d-8
a(18)=-0.11812745704870201d-8
a(19)= 0.10434267116911005d-9
a(20)= 0.77822634399050713d-11
a(21)=-0.36968056186422057d-11
a(22)= 0.51003702874544760d-12
a(23)=-0.20583260535665068d-13
a(24)=-0.53481225394230180d-14
a(25)= 0.12267786282382608d-14
a(26)=-0.11812593016974588d-15
a(27)= 0.11866922547516003d-17
a(28)= 0.14123806553180318d-17
a(29)=-0.22987456844353702d-18
a(30)= 0.17144063219273374d-19

M=int(11.3*abs(w)+13)
if(M.gt.30)M=30
r=a(M)
do i=M-1,1,-1
r=r*w+a(i)
enddo

cgamma=s/(r*w)
else
! For x=arb, |y|\gt 1.45
s=1d0
if(abs(q).lt.9d0)then
do i=0,7
s=s*(q+i)
enddo
s=1d0/s
q=q+8d0
endif

AB(1) = 0.83333333333333333d-1
AB(2) =-0.27777777777777778d-2
AB(3) = 0.79365079365079365d-3
AB(4) =-0.59523809523809524d-3
AB(5) = 0.84175084175084175d-3
AB(6) =-0.19175269175269175d-2
AB(7) = 0.64102564102564103d-2
AB(8) =-0.29550653594771242d-1

q1=1d0/q
q2=q1*q1

r=AB(8)
do i=7,1,-1
r=r*q2+AB(i)
enddo
cgamma=s*exp((q-0.5d0)*log(q)-q+ln2pi2+r*q1)
endif
if(dz.lt.0d0)cgamma=pi/(cgamma*sin(pi*z))
endif

return
end function cgamma

ifort version 16.0.2 20160204, 12.1.0 20110811
gfortran version 4.8.4
で動くことを確かめています。

こちらに4倍精度用のゼータ関数のコードを置いておきます。
精度は何点かで試しましたが、30~32桁一致するくらいです。

▼ここクリックでこの場に4倍精度ゼータ関数のコードを展開

function zzeta(s)
implicit none
complex*32,intent(in)::s
complex*32::zzeta
!Checked : ifort (success : version 16.0.2 20160204)
! : gfortran(success : version 4.8.4)
! : gfortran(fail : version 4.4.7)
!
!http://slpr.sakura.ne.jp/qp/
! Author : sikino
! date: 2017/01/12

integer,parameter::nmax=2000
integer::n,k
real*16::nq,nq1,ln1,ln2
real*16,parameter::eps=1q-32
real*16,parameter::pi=3.14159265358979323846264338327950288q0
complex*32::ds,is,d,s0,tc,tt,qzeta
complex*32::a(0:Nmax),b(0:nmax)
complex*32::zgamma
external::zgamma

ds=dble(s); is=imag(s)
if(ds.eq.1q0.and.is.eq.0q0)then
zzeta=1q+3000
elseif(ds.eq.0q0.and.is.eq.0q0)then
zzeta=-0.5q0
else
if(real(s).lt.1q0)then
s0=1q0-s
else
s0=s
endif

tt=1q0/(1q0-2q0**(1q0-s0))
qzeta=0.5q0*tt

a=cmplx(0q0,0q0,kind=16); b=cmplx(0q0,0q0,kind=16)
b(0)=0.25q0*tt
do n=1,nmax
nq=n
nq1=nq+1q0

a(0)=b(0)
do k=1,n-1
a(k)=b(k)*nq/(n-k)
enddo
!----
ln2=log(nq/nq1)
if(n.eq.1)then
d=s0*ln2
else
d=(ln2/ln1)*d
endif
a(n)=-a(n-1)*exp(d)/nq
ln1=ln2
!----

tc=cmplx(0q0,0q0,kind=16)
do k=0,n
tc=tc+a(k)
enddo

qzeta=qzeta+tc
if(abs(qzeta).ge.1q0)then
if(abs(tc/qzeta).le.eps)exit
else
if(abs(tc).le.eps)exit
endif

b=0.5q0*a
enddo
if(n.eq.nmax+1)write(6,*)" ***worning, didn't converge at zzeta"

if(real(s).lt.1q0)then
zzeta=(2q0**s)*pi**(-s0)*sin(0.5q0*pi*s)*qzeta*zgamma(s0*1q0)
else
zzeta=qzeta
endif
end if

return
end function zzeta

function zgamma(z)
! sikinote (http://slpr.sakura.ne.jp/qp)
! author : sikino
! date : 2017/01/09
! 2017/01/10

implicit none
complex*32,intent(in)::z
complex*32::zgamma

integer::i,n,M
real*16,parameter:: pi=3.14159265358979323846264338327950288q0
real*16,parameter:: e1=0.36787944117144232159552377016146087q0
real*16,parameter::ln2pi2=0.91893853320467274178032973640561764q0
real*16,parameter:: sq2pi=2.50662827463100050241576528481104525q0
real*16::a(1:50),AB(1:30),dz,iz,aw
complex*32::q,s,w,r,q1,q2

dz=real(z); iz=imag(z)
if(dz.eq.nint(dz).and.iz.eq.0q0)then
if(dz.gt.0q0)then
s=dcmplx(1q0,0q0)
n=nint(dz)
do i=2,n
s=s*i
enddo
else
s=1q+3000
endif
zgamma=s
else
q=z
if(dz.lt.0q0)q=1q0-z

if(abs(imag(q)).lt.1.35q0)then
! For x=arb, |y|\lt 1.35
n=int(real(q))
s=1q0
do i=1,n
s=s*(q-i)
enddo
w=q-real(n)

a(1) = 1q0
a(2) = 0.577215664901532860606512090082402431q0
a(3) =-0.655878071520253881077019515145390481q0
a(4) =-0.420026350340952355290039348754298187q-1
a(5) = 0.166538611382291489501700795102105236q0
a(6) =-0.421977345555443367482083012891873913q-1
a(7) =-0.962197152787697356211492167234819898q-2
a(8) = 0.721894324666309954239501034044657271q-2
a(9) =-0.116516759185906511211397108401838867q-2
a(10)=-0.215241674114950972815729963053647806q-3
a(11)= 0.128050282388116186153198626328164323q-3
a(12)=-0.201348547807882386556893914210218184q-4
a(13)=-0.125049348214267065734535947383309224q-5
a(14)= 0.113302723198169588237412962033074494q-5
a(15)=-2.056338416977607103450154130020572837q-7
a(16)= 6.116095104481415817862498682855342867q-9
a(17)= 5.002007644469222930055665048059991303q-9
a(18)=-1.181274570487020144588126565436505578q-9
a(19)= 1.043426711691100510491540332312250191q-10
a(20)= 7.782263439905071254049937311360777226q-12
a(21)=-3.696805618642205708187815878085766237q-12
a(22)= 5.100370287454475979015481322863231803q-13
a(23)=-2.058326053566506783222429544855237420q-14
a(24)=-5.348122539423017982370017318727939949q-15
a(25)= 1.226778628238260790158893846622422428q-15
a(26)=-1.181259301697458769513764586842297831q-16
a(27)= 1.186692254751600332579777242928674071q-18
a(28)= 1.412380655318031781555803947566709037q-18
a(29)=-2.298745684435370206592478580633699260q-19
a(30)= 1.714406321927337433383963370267257067q-20
a(31)= 1.337351730493693114864781395122268023q-22
a(32)=-2.054233551766672789325025351355733797q-22
a(33)= 2.736030048607999844831509904330982015q-23
a(34)=-1.732356445910516639057428451564779799q-24
a(35)=-2.360619024499287287343450735427531031q-26
a(36)= 1.864982941717294430718413161878666894q-26
a(37)=-2.218095624207197204399716913626860384q-27
a(38)= 1.297781974947993668824414486330595180q-28
a(39)= 1.180697474966528406222745415509988346q-30
a(40)=-1.124584349277088090293654674261653627q-30
a(41)= 1.277085175140866203990206677750563204q-31
a(42)=-7.391451169615140823461289330001078362q-33
a(43)= 1.134750257554215760954165252401478595q-35
a(44)= 4.639134641058722029944804894433875948q-35
a(45)=-5.347336818439198875077418068570325956q-36
a(46)= 3.207995923613352622861238213439296964q-37
a(47)=-4.445829736550756882101576639628564322q-39
a(48)=-1.311174518881988712901096702293607192q-39
a(49)= 1.647033352543813886818056435068597372q-40
a(50)=-1.056233178503581218582902184722885523q-41

aw=abs(w)
if(aw.le.0.42q0)then
M=int(35q0*aw+20q0)+1
else
M=int(15.5q0*aw+28.2q0)+1
endif
if(M.gt.50)M=50
r=a(M)
do i=M-1,1,-1
r=r*w+a(i)
enddo

zgamma=s/(r*w)
else
! For x=arb, |y|\ge 1.35
s=1q0
if(abs(q).lt.18q0)then
do i=0,16
s=s*(q+i)
enddo
s=1q0/s
q=q+17q0
endif

AB(1) = 0.83333333333333333333333333333333333q-1
AB(2) =-0.27777777777777777777777777777777778q-2
AB(3) = 0.79365079365079365079365079365079365q-3
AB(4) =-0.59523809523809523809523809523809524q-3
AB(5) = 0.84175084175084175084175084175084175q-3
AB(6) =-0.19175269175269175269175269175269175q-2
AB(7) = 0.64102564102564102564102564102564103q-2
AB(8) =-0.29550653594771241830065359477124183q-1
AB(9) = 0.17964437236883057316493849001588940q0
AB(10)=-1.3924322169059011164274322169059011q0
AB(11)= 13.402864044168391994478951000690131q0
AB(12)=-156.84828462600201730636513245208897q0
AB(13)= 2193.1033333333333333333333333333333q0
AB(14)=-36108.771253724989357173265219242231q0
AB(15)= 691472.268851313067108395250775673468q0
AB(16)=-1.52382215394074161922833649588867805q7

q1=1q0/q; q2=q1*q1

r=AB(16)
do i=15,1,-1
r=r*q2+AB(i)
enddo

zgamma=s*exp((q-0.5q0)*log(q)-q+ln2pi2+r*q1)
endif
if(dz.lt.0q0)zgamma=pi/(zgamma*sin(pi*z))
end if

return
end function zgamma

ゼータ関数のデータ

計算結果のデータも置いておきます。
(11MB)riemanndata.tar.gz
(9MB)riemanndata.zip
それぞれ、中に
“riemannzeta_high_resolution.txt”
“riemannzeta_low_resolution.txt”
が入っています。
highの方は1列目から4列目にかけて
${\rm Re}(s):[-20\sim 20],{\rm Im}(s):[-50\sim 50],{\rm Re}(\zeta(s)),{\rm Im}(\zeta(s))$、0.1ずつ、小数点以下13桁で出力したもので、
lowの方は1列目から4列目にかけて
${\rm Re}(s):[-20\sim 20],{\rm Im}(s):[-50\sim 50],{\rm Re}(\zeta(s)),{\rm Im}(\zeta(s))$、0.2ずつ、小数点以下5桁で出力したものです。

x,y,z軸を${\rm Re}(s),{\rm Im}(s),{\rm Re}(\zeta(s))$にとりグラフにすると、

というグラフが得られます。

~~${\rm Re}(s)$が負で${\rm Im}(s)$が小さいところで発散が起こるようですので、注意をして使ってください。~~
何故倍精度でおかしい問題が起こっていたのか？今回の場合はオーバーフローでもなく、アンダーフローでもありません。
原因は桁落ちです。和と差によって激しく桁落ちしてしまうのです。

また、倍精度よりも計算が早いです。なぜならば目標精度までスムーズに収束してくれるからです。

相反公式を用いたので上記問題は発生しません。

[adsense1]

ゼータ関数の実部、虚部、絶対値のグラフ

ゼータ関数の実部、虚部、絶対値のグラフは以下の通りになります。

実部${\rm Re}(\zeta(s))$

カラーマップ(横軸:${\rm Re}(s)$,縦軸$:{\rm Im}(s)$,カラー:${\rm Re}(\zeta(s))$)

虚数${\rm Im}(\zeta(s))$

カラーマップ(横軸:${\rm Re}(s)$,縦軸$:{\rm Im}(s)$,カラー:${\rm Im}(\zeta(s))$)

絶対値$|\zeta(s)|$

カラーマップ(横軸:${\rm Im}(s)$,縦軸$:{\rm Re}(s)$,カラー:$|\zeta(s)|$)

対数プロット(横軸:${\rm Re}(s)$,縦軸$:{\rm Im}(s)$,カラー:$(\log{10}(|\zeta(s)|$)

（おまけ）ゼータ関数のニュートン写像(詳しくはリーマンのゼータ関数の複素力学系へ)

非自明なゼロ点

上記の計算方法(4倍精度のプログラム)で、リーマンの非自明なゼロ点を考えます。
「非自明なゼロ点」はゼータ関数の実数軸上には無いゼロ点のことです。
リーマン予想とは
ゼータ関数$\zeta(s)$の「非自明なゼロ点」が$s=\frac{1}{2}+ix$($x$は実数)上に全てある
というものです。$\zeta(s)=0$が満たされるときはゼータ関数の実部も虚部もゼロになるので、簡単に言えばゼータ関数の絶対値がゼロになる$|\zeta(s)|=0$と言うことです。
実際に$|\zeta(s)|$を可視化してみますとこうなります。

自明なゼロ点

僕は他のページでは自明な点に関する図を見たことがありません。非常に小さいからでしょう。
小さいものを見るには対数スケールで見るのが良いです。実際に書かせてみます。

確かに小さいですね。対数にしないと明らかに見えません。

絶対値$|\zeta(s)|$を上の図のように表示させるgnuplotスクリプトはこちら。

set xr[-15:10]
set yr[-40:40]
set zr[0:2]
set view 69, 73, 2, 0.3
set view equal xy
set ticslevel 0

set palette defined (-9 "purple",-6 "blue",-3 "turquoise", 0 "#f5f5f5", 3 "gold", 6 "orange", 9 "red")
set colorbox horiz user origin .1,.1 size .8,.04
set pm3d at b

splot "fort.11" u 1:2:(sqrt($3**2+$4**2)) every :3 w l lw 0.1

[adsense2]

クリティカルライン上での計算

クリティカルラインとは
$
\displaystyle s=\frac{1}{2}+ix
$
の線です。リーマン予想は、非自明なゼロ点が全てこの上にあるか？という問いかけであり、この計算と繋がっています。

倍精度型のアルゴリズム(zeta)では、$x\sim 400$程度まで計算することが出来ます。
これは、

相反公式の
$
\displaystyle \sin(\frac{\pi}{2}s)
$
の値がオーバーフローを起こすためです。

倍精度演算では約$10^{306}$まで扱えるので、
$
s=\sin^{-1}(10^{306})\frac{2}{\pi}\sim 1+i449
$
4倍精度では約$10^{4932}$まで扱えるので、
$
s=\sin^{-1}(10^{4932})\frac{2}{\pi}\sim 1+i7230
$
となります。
しかし、実際に計算してみると倍精度では確かに$x=450$程度まで計算可能ですが、4倍精度では$x=1600$程度までしか正しい値を返してくれません。どこかで何かが起こっているはずですが、解析はしてません。

実際に試してみますと、
倍精度では

となります。

ゼータ関数の導関数とプログラム

ゼータ関数の導関数も求めましょう。
導関数はゼータ関数の式をそのまま微分すればいいです。なぜなら、正則関数として解析接続されているからです。
実際に微分しますと、
$
\displaystyle
\zeta'(s)=\frac{1}{(2^s-2)^2}\sum_{n=0}^{\infty} \frac{1}{2^{n+1}} \sum_{k=0}^{n} {n \choose k} \frac{(-1)^k}{(k+1)^s}
\left[ -2^s \{(2^s-2)\ln(k+1)+\ln4\} \right]
$
変形して以下のように
$
\begin{align}
\zeta'(s) &=\sum_{n=0}^{\infty} b_n \\
b_n &= \sum_{k=0}^{n}a_k^n \\
a_k^n &= (-1)^k {n \choose k} \left[ \frac{-2^s (k+1)^{-s}\{(2^s-2)\ln(k+1)+\ln4\}}{(2^s-2)^2} \right]
\end{align}
$
となります。同じように漸化式を求めてやれば、
$
\begin{align}
a_k^n&=\frac{1}{2}\left(\frac{n}{n-k}\right)a_{k}^{n-1} \\
a_k^n&=-\frac{k-n-1}{k}\left(\frac{k}{k+1}\right)^s\frac{(2^s-2)\ln(k+1)+\ln 4}{(2^s-2)\ln(k)+\ln 4} a_{k-1}^n \\
a_0^0&=\frac{1}{2}\frac{-2^s \ln 4}{(2^s-2)^2}
\end{align}
$

として求められます。
実際の数値計算では、2番目の漸化式は余り有能ではないかもしれません。
ゼータ関数を求めるときと同じようにやるならば、必要なのは$a_n^n$を$a_{n-1}^n$から求める作業です。漸化式を用いても直接求めても計算量はさほど変わらない気がします。

導関数の相反公式

相反公式ももちろん存在します[3]。
解析接続されているので、元の相反公式の左辺、右辺をそのまま微分することで求められます。
$
\begin{align}
-\zeta'(1-s)&=2(2\pi)^{-s}\left\{x\cos\left(\frac{\pi s}{2}\right)+y\sin\left(\frac{\pi s}{2}\right)\right\}\Gamma(s)\zeta(s) \\
&\hspace{3em} +2(2\pi)^{-s}\cos\left(\frac{\pi s}{2}\right)\left\{\Gamma(s)\zeta(s)\right\}’ \\
& x = -\ln(2\pi), \;\; y=-\frac{\pi}{2}
\end{align}
$

$
\Gamma'(s)=\Gamma(s)\psi(s)
$
ここで$\psi(z)$はディガンマ関数

確かめはwolfram alphaを用いて以下のコマンドで行いました。
ReplaceAll[D[Zeta[x], x], {x -> -6+2i}]

fortran90によるプログラム

fortran90によるプログラムはこんな感じになるかと思います。
相反公式で、ゼータ関数、倍精度ガンマ関数とディガンマ関数を用います。
ゼータ関数は上の、倍精度ガンマ関数とディガンマ関数はつけておきました。詳しくはディガンマ関数の数値計算をご覧ください。
ここに記述しているのは倍精度のみです。

▼ここクリックでゼータ関数の導関数のコード展開

function cdzeta(s)
implicit none
complex(kind(0d0)),intent(in)::s
complex(kind(0d0))::cdzeta
!
!Checked : ifort (success : version 16.0.2 20160204)
! : gfortran(success : version 4.8.4)
!
!http://slpr.sakura.ne.jp/qp/
! Author : sikino
! date: 2016/12/30
! : 2017/01/13
!
integer,parameter::nmax=100
integer::n,k
double precision::nq,nq1,x,y,g1,g2,ds,is
double precision,parameter::eps=1q-14
double precision,parameter::pi=3.14159265358979324d0
double precision,parameter::ln4=1.386294361119890619d0 ! ln(4)
complex(kind(0d0))::d,s0,tc,tt,ys,gs,zs,s1,s2,h1,h2
complex(kind(0d0))::a(0:Nmax),b(0:nmax)
complex(kind(0d0))::cgamma,czeta,digamma
external::cgamma,czeta,digamma

ds=dble(s); is=imag(s)
if(ds.eq.1d0.and.is.eq.0d0)then
cdzeta=1d+300
elseif(ds.eq.0d0.and.is.eq.0d0)then
cdzeta=-0.5d0*log(2d0*pi)
else
if(real(s).lt.1d0)then
s0=1d0-s
else
s0=s
endif

s1=2d0**s0
s2=s1-2d0
tt=-s1*ln4/(s2*s2)
cdzeta=0.5d0*tt

a=cmplx(0d0,0d0,kind=8); b=cmplx(0d0,0d0,kind=8)
b(0)=cmplx(0.25d0,0d0,kind=8)*tt
do n=1,nmax
nq=n
nq1=nq+1d0

a(0)=b(0)
do k=1,n-1
a(k)=b(k)*nq/(n-k)
enddo
!----
g2=log(nq/nq1)
if(n.eq.1)then
h1=ln4
d=s0*log(0.5d0)
else
d=d*g2/g1
endif
h2=s2*log(nq1)+ln4
a(n)=-a(n-1)*exp(d)*h2/(h1*nq)
g1=g2
h1=h2
!----

tc=cmplx(0d0,0d0,kind=8)
do k=0,n
tc=tc+a(k)
enddo
cdzeta=cdzeta+tc

if(abs(cdzeta).ge.1d0)then
if(abs(tc/cdzeta).le.eps)exit
else
if(abs(tc).le.eps)exit
endif

b=0.5d0*a
enddo
if(n.eq.nmax+1)write(6,*)" ***worning, didn't converge at cdzeta"

if(real(s).lt.1d0)then
x=-log(2d0*pi)
y=-pi*0.5d0
ys=-y*s0
zs=czeta(1d0-s)
gs=cgamma(s0)
cdzeta=-2d0*(2d0*pi)**(-s0)*(&
(x*cos(ys)+y*sin(ys))*gs*zs+cos(ys)*(digamma(1d0-s)*gs*zs+gs*cdzeta))
else
cdzeta=cdzeta
endif
end if

return
end function cdzeta

function cgamma(z)
! sikinote (http://slpr.sakura.ne.jp/qp)
! author : sikino
! date : 2017/01/09
! 2017/01/10

implicit none
complex(kind(0d0)),intent(in)::z
complex(kind(0d0))::cgamma

integer::i,n,M
double precision,parameter::pi=3.14159265358979324d0
double precision,parameter::e1=0.3678794411714423216d0
double precision,parameter::ln2pi2=0.91893853320467274d0
double precision,parameter::sq2pi=2.5066282746310005d0
double precision::a(1:30),AB(1:14),dz,iz
complex(kind(0d0))::q,s,w,r,q1,q2

dz=dble(z)
iz=imag(z)

if(dz.eq.nint(dz).and.iz.eq.0d0)then
if(dz.gt.0d0)then
s=dcmplx(1d0,0d0)
n=nint(dz)
do i=2,n
s=s*i
enddo
else
s=1d+300
endif
cgamma=s
else
q=z
if(dz.lt.0d0)q=1d0-z

if(abs(iz).lt.1.45d0)then
! For x=arb, |y|\lt 1.45
n=int(dble(q))
s=1d0
do i=1,n
s=s*(q-i)
enddo
w=q-dble(n)

a(1) = 1d0
a(2) = 0.57721566490153286d0
a(3) =-0.65587807152025388d0
a(4) =-0.42002635034095236d-1
a(5) = 0.16653861138229149d0
a(6) =-0.42197734555544337d-1
a(7) =-0.96219715278769736d-2
a(8) = 0.72189432466630995d-2
a(9) =-0.11651675918590651d-2
a(10)=-0.21524167411495097d-3
a(11)= 0.12805028238811619d-3
a(12)=-0.20134854780788239d-4
a(13)=-0.12504934821426707d-5
a(14)= 0.11330272319816959d-5
a(15)=-0.20563384169776071d-6
a(16)= 0.61160951044814158d-8
a(17)= 0.50020076444692229d-8
a(18)=-0.11812745704870201d-8
a(19)= 0.10434267116911005d-9
a(20)= 0.77822634399050713d-11
a(21)=-0.36968056186422057d-11
a(22)= 0.51003702874544760d-12
a(23)=-0.20583260535665068d-13
a(24)=-0.53481225394230180d-14
a(25)= 0.12267786282382608d-14
a(26)=-0.11812593016974588d-15
a(27)= 0.11866922547516003d-17
a(28)= 0.14123806553180318d-17
a(29)=-0.22987456844353702d-18
a(30)= 0.17144063219273374d-19

M=int(11.3*abs(w)+13)
if(M.gt.30)M=30
r=a(M)
do i=M-1,1,-1
r=r*w+a(i)
enddo

cgamma=s/(r*w)
else
! For x=arb, |y|\gt 1.45
s=1d0
if(abs(q).lt.9d0)then
do i=0,7
s=s*(q+i)
enddo
s=1d0/s
q=q+8d0
endif

AB(1) = 0.83333333333333333d-1
AB(2) =-0.27777777777777778d-2
AB(3) = 0.79365079365079365d-3
AB(4) =-0.59523809523809524d-3
AB(5) = 0.84175084175084175d-3
AB(6) =-0.19175269175269175d-2
AB(7) = 0.64102564102564103d-2
AB(8) =-0.29550653594771242d-1

q1=1d0/q
q2=q1*q1

r=AB(8)
do i=7,1,-1
r=r*q2+AB(i)
enddo
cgamma=s*exp((q-0.5d0)*log(q)-q+ln2pi2+r*q1)
endif
if(dz.lt.0d0)cgamma=pi/(cgamma*sin(pi*z))
endif

return
end function cgamma

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=50
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

注釈

他のページでの例えば
リーマンのゼータ関数で遊び倒そう (Ruby編)
等では、漸化式を使わないで直接計算しています。このため、余り項の数$n$の数を増やす事が出来ないこと、倍精度と思われるため、桁落ちの問題が回避できない問題が起こります。
実際、私も直接計算するコードをfortran90で書きましたが、項の数を増やせず収束しない事、また計算時間が2~3倍余計にかかります。

参考文献
Riemann Zeta Function — from Wolfram MathWorld, Eq.(21)
draw a horizontal gradient somewhere at the bottom of the graph. -ColorBox

[3]Tom M. Apóstol “Formulas for Higher Derivatives of the Riemann Zeta Function”, Math. of. Comp vol 44, 169, 1985, pp. 223-232.

[4]Madhava-Gregory-Leibniz級数の収束の加速(収束の加速1)

[5]リーマンゼータ関数の級数表示による解析接続インテジャーズ

プログラミングと数値計算, 数学

特殊関数

2016年5月29日 sikino 3件のコメント

以下に示すのは

・ラゲール陪多項式(Associated Laguerre Polynomials):$L_n^{(\alpha)}(x)$
・ルジャンドル陪多項式(Associated Legendre Polynomials):$P_l^m(x)$
・ヤコビ多項式(Jacobi Polynomials)$P_n^{(\alpha,\beta)}(x)$
・エルミート多項式(Hermite Polynomials):$H_n(x)$
とその1階微分、2階微分

のfortran90用のプログラムです。
複素平面で有効です。

このプログラムを使用して生じた問題に対する責任は一切取りません。
ミスがあるかもしれないことを念頭に置いて使ってください。
また、使えない場合がありましたらご連絡いただけると幸いです。

※2017/10/27
整数次数の、実数、複素引数エルミート多項式を追加しました。

特殊関数のモジュール”specialfunction”

▼ここクリックでモジュールが展開されます

module specialfunction
! sikinote, http://slpr.sakura.ne.jp/qp
! author : sikino
! date : 2016/05/29
! : 2017/10/27
!
! Reference
! Abramowitz & Stegun (1965), ISBN 978-0486612720
! (Numerical Method, 22.18, Page 789)
!
! -----------------------------------------
! laguerre(n,a,x) L_n^a(x)
! ...valid below combinations
! (n,a,x) = (I,I,D)
! (n,a,x) = (I,D,D)
! (n,a,x) = (I,C,C)
! 1st derivative laguerre1(n,a,x)
! 2nd derivative laguerre2(n,a,x)
! -----------------------------------------
! legendre(l,m,x) P_l^m(x)
! ...valid below combinations
! (l,m,x) = (I,I,D)
! (l,m,x) = (I,I,Z)
! 1st derivative legendre1(l,m,x)
! 2nd derivative legendre2(l,m,x)
! -----------------------------------------
! hermite(n,x) H_n(x)
! ...valid below combinations
! (n,x) = (I,D)
! (n,x) = (I,Z)
! 1st derivative hermite1(n,x)
! 2nd derivative hermite2(n,x)
! -----------------------------------------
! jacobi(n,alpha,beta,x) P_n^(alpha,beta)(x)
! ...valid below combinations
! (n,alpha,beta,x) = (I,D,D,C)
! (n,alpha,beta,x) = (I,C,C,C)
!
! 1st derivative jacobi1(n,alpha,beta,x)
! 2nd derivative jacobi2(n,alpha,beta,x)
! -----------------------------------------
!
! I -> integer
! D -> double precision
! C -> complex(kind(0d0))
!
implicit none
interface laguerre
module procedure &
ilaguerre, &
dlaguerre, &
claguerre
end interface laguerre
interface legendre
module procedure &
dalegendre, &
calegendre
end interface legendre
interface hermite
module procedure &
dhermite, &
chermite
end interface hermite
interface jacobi
module procedure &
djacobi, &
cjacobi
end interface jacobi

interface laguerre1
module procedure &
ilaguerre1, &
dlaguerre1, &
claguerre1
end interface laguerre1
interface legendre1
module procedure &
dalegendre1 , &
calegendre1
end interface legendre1
interface hermite1
module procedure &
dhermite1, &
chermite1
end interface hermite1
interface jacobi1
module procedure &
djacobi1, &
cjacobi1
end interface jacobi1

interface laguerre2
module procedure &
ilaguerre2, &
dlaguerre2, &
claguerre2
end interface laguerre2
interface legendre2
module procedure &
dalegendre2 , &
calegendre2
end interface legendre2
interface hermite2
module procedure &
dhermite2, &
chermite2
end interface hermite2
interface jacobi2
module procedure &
djacobi2, &
cjacobi2
end interface jacobi2

contains
function ilaguerre(n,a,x)
!
! a
! L (x)
! n
! a is integer
! not only 0 <= x < infty, but also negative value.
integer,intent(in)::n,a
double precision,intent(in)::x
double precision::ilaguerre,fx

integer::m
double precision::da,dn,dm,bm,cm,dnx

if(n.le.-1)then
ilaguerre=0d0
return
endif

dn=dble(n); da=dble(a)
!dnx=exp(dlgamma(da+dn+1d0)-dlgamma(dn+1d0)-dlgamma(da+1d0))
dnx=exp(log(gamma(da+dn+1d0))-log(gamma(dn+1d0))-log(gamma(da+1d0)))

fx=x
ilaguerre=1d0
do m=n,1,-1
dm=dble(m); cm=dm*(da+dm); bm=dble(n-m+1)
ilaguerre=1d0-ilaguerre*fx*bm/cm
enddo
ilaguerre=ilaguerre*dnx

return
end function ilaguerre
function dlaguerre(n,a,x)
!
! (a)
! L (x)
! n
!
! not only 0 <= x < infty, but also negative value.
integer,intent(in)::n
double precision,intent(in)::a
double precision,intent(in)::x
double precision::dlaguerre,fx

integer::m
double precision::dn,dm,bm,cm,dnx

if(n.le.-1)then
dlaguerre=0d0
return
endif

dn=dble(n)
!dnx=exp(dlgamma(a+dn+1d0)-dlgamma(dn+1d0)-dlgamma(a+1d0))
dnx=exp(log(gamma(a+dn+1d0))-log(gamma(dn+1d0))-log(gamma(a+1d0)))

fx=x
dlaguerre=1d0
do m=n,1,-1
dm=dble(m); cm=dm*(a+dm); bm=dble(n-m+1)
dlaguerre=1d0-dlaguerre*fx*bm/cm
enddo
dlaguerre=dlaguerre*dnx

return
end function dlaguerre
function claguerre(n,a,x)
!
! (a)
! L (x)
! n
!
! not only 0 <= x < infty, but also negative value.
integer,intent(in)::n
complex(kind(0d0)),intent(in)::a,x
complex(kind(0d0))::claguerre,fx,dnx,cm,dn

integer::m
double precision::dm,bm

if(n.le.-1)then
claguerre=dcmplx(0d0,0d0)
return
endif

dn=dcmplx(dble(n),0d0)
dnx=exp(clgamma(a+dn+1d0)-clgamma(dn+1d0)-clgamma(a+1d0))

fx=x
claguerre=dcmplx(1d0,0d0)
do m=n,1,-1
dm=dble(m); cm=dm*(a+dm); bm=dble(n-m+1)
claguerre=1d0-claguerre*fx*bm/cm
enddo
claguerre=claguerre*dnx

return
end function claguerre
!======================================
function dalegendre(l,m,x)
implicit none
integer,intent(in)::m,l
double precision,intent(in)::x
double precision::dalegendre
! m
! P (x)
! l
!
! not only -1 <= x <= 1, but also more larger analytical continuum
integer::am,q,i
double precision::fm,dq,dam2,fl
double precision::t1,t2

dalegendre=0d0
if(abs(m).gt.l)then
dalegendre=0d0
return
endif

am=abs(m)

fm=1d0
do i=1,2*am-1,2
fm=fm*dble(i)
enddo

dam2=dble(2*am-1)
if(l-am.eq.0)then
dalegendre=dble((-1)**am)*fm*(1d0-x*x)**(0.5d0*dble(am))
elseif(l-am.eq.1)then
t1=dble((-1)**am)*fm*(1d0-x*x)**(0.5d0*dble(am))
dalegendre=dble(2*am+1)*x*t1
else
t1=dble((-1)**am)*fm*(1d0-x*x)**(0.5d0*dble(am))
t2=dble(2*am+1)*x*t1
do q=2,l-am
dq=1d0/dble(q)
dalegendre=(dam2*dq+2d0)*x*t2-(dam2*dq+1d0)*t1
t1=t2
t2=dalegendre
enddo
endif
if(m.lt.0)then
fl=1d0
do i=l-am+1,l+am
fl=fl*dble(i)
enddo
dalegendre=dalegendre*dble((-1)**am)/fl
endif

return
end function dalegendre
function calegendre(l,m,x)
! m
! P (x)
! l
!
! |m| <= l
! not only -1 <= x <= 1, but also more larger analytical continuum
implicit none
integer,intent(in)::m,l
complex(kind(0d0)),intent(in)::x
complex(kind(0d0))::calegendre
integer::am,q,i
double precision::fm,dq,dam2,fl
complex(kind(0d0))::t1,t2

calegendre=dcmplx(0d0,0d0)
if(abs(m).gt.l)then
calegendre=dcmplx(0d0,0d0)
return
endif

am=abs(m)

fm=1d0
do i=1,2*am-1,2
fm=fm*dble(i)
enddo

dam2=dble(2*am-1)
if(l-am.eq.0)then
calegendre=dble((-1)**am)*fm*(1d0-x*x)**(0.5d0*dble(am))
elseif(l-am.eq.1)then
t1=dble((-1)**am)*fm*(1d0-x*x)**(0.5d0*dble(am))
calegendre=dble(2*am+1)*x*t1
else
t1=dble((-1)**am)*fm*(1d0-x*x)**(0.5d0*dble(am))
t2=dble(2*am+1)*x*t1
do q=2,l-am
dq=1d0/dble(q)
calegendre=(dam2*dq+2d0)*x*t2-(dam2*dq+1d0)*t1
t1=t2
t2=calegendre
enddo
endif
if(m.lt.0)then
fl=1d0
do i=l-am+1,l+am
fl=fl*dble(i)
enddo
calegendre=calegendre*dble((-1)**am)/fl
endif

return
end function calegendre
!===============================
function djacobi(n,alpha,beta,x)
implicit none
double precision,intent(in)::alpha,beta
integer,intent(in)::n
complex(kind(0d0)),intent(in)::x
complex(kind(0d0))::djacobi
!
! (alpha,beta)
! P (x)
! n
!
! Orthogonality
! 1
! / a b (a,b) (a,b) 2^(a+b+1) Gamma(n+a+1)Gamma(n+b+1)
! | (1-x) (1+x) P (x) P (x) dx = ---------- ------------------------- delta_nm
! / m n 2n+a+b+1 Gamma(n+a+b+1)Gamma(n+1)
! -1
!
! Defferential equation
!
! (1-x^2)y'' + (b-a-(a+b+2)x)y' + n(n+a+b+1)y = 0
!
integer::m
double precision::dnx,dn,bm,cm
complex(kind(0d0))::fx

djacobi=dcmplx(0d0,0d0)
dn=dble(n)
dnx=exp(log(gamma(alpha+dn+1d0))-log(gamma(dn+1d0))-log(gamma(alpha+1d0)))

fx=1d0-x
djacobi=dcmplx(1d0,0d0)
do m=n,1,-1
bm=dble(n-m+1)*(alpha+beta+dble(n+m))
cm=dble(2*m)*(alpha+dble(m))
djacobi=1d0-djacobi*fx*bm/cm
enddo
djacobi=dnx*djacobi

return
end function djacobi
!===============================
function cjacobi(n,alpha,beta,x)
implicit none
integer,intent(in)::n
complex(kind(0d0)),intent(in)::alpha,beta,x
complex(kind(0d0))::cjacobi
!
! (alpha,beta)
! P (x)
! n
!
! Orthogonality
! 1
! / a b (a,b) (a,b) 2^(a+b+1) Gamma(n+a+1)Gamma(n+b+1)
! | (1-x) (1+x) P (x) P (x) dx = ---------- ------------------------- delta_nm
! / m n 2n+a+b+1 Gamma(n+a+b+1)Gamma(n+1)
! -1
!
! Defferential equation
!
! (1-x^2)y'' + (b-a-(a+b+2)x)y' + n(n+a+b+1)y = 0
!
integer::m
complex(kind(0d0))::cnx,bm,cm,dn,fx

cjacobi=dcmplx(0d0,0d0)
dn=dcmplx(dble(n),0d0)
cnx=exp(clgamma(alpha+dn+1d0)-clgamma(dn+1d0)-clgamma(alpha+1d0))

fx=1d0-x
cjacobi=dcmplx(1d0,0d0)
do m=n,1,-1
bm=dble(n-m+1)*(alpha+beta+dble(n+m))
cm=dble(2*m)*(alpha+dble(m))
cjacobi=1d0-cjacobi*fx*bm/cm
enddo
cjacobi=cnx*cjacobi

return
end function cjacobi

!=======================================

function ilaguerre1(n,a,x)
implicit none
integer,intent(in)::n,a
double precision,intent(in)::x
double precision::ilaguerre1

ilaguerre1=-laguerre(n-1,a+1,x)

return
end function ilaguerre1
function dlaguerre1(n,a,x)
implicit none
integer,intent(in)::n
double precision,intent(in)::a,x
double precision::dlaguerre1

dlaguerre1=-laguerre(n-1,a+1d0,x)

return
end function dlaguerre1
function claguerre1(n,a,x)
implicit none
integer,intent(in)::n
complex(kind(0d0)),intent(in)::a,x
complex(kind(0d0))::claguerre1

claguerre1=-laguerre(n-1,a+1d0,x)

return
end function claguerre1
function ilaguerre2(n,a,x)
implicit none
integer,intent(in)::n,a
double precision,intent(in)::x
double precision::ilaguerre2

ilaguerre2=laguerre(n-2,a+2,x)

return
end function ilaguerre2
function dlaguerre2(n,a,x)
implicit none
integer,intent(in)::n
double precision,intent(in)::a,x
double precision::dlaguerre2

dlaguerre2=laguerre(n-2,a+2d0,x)

return
end function dlaguerre2
function claguerre2(n,a,x)
implicit none
integer,intent(in)::n
complex(kind(0d0)),intent(in)::a,x
complex(kind(0d0))::claguerre2

claguerre2=laguerre(n-2,a+2d0,x)

return
end function claguerre2

function dalegendre1(l,m,x)
implicit none
integer,intent(in)::m,l
double precision,intent(in)::x
double precision::dalegendre1

double precision::c0,c1,c2

dalegendre1=0d0
if(abs(m).gt.l)then
dalegendre1=0d0
return
endif

if(abs(m).eq.1)then
dalegendre1=(dble(l)*x*legendre(l,m,x)-dble(l+m)*legendre(l-1,m,x))/(x*x-1)
else
c0=1d0/dble(4*m+4)
c1=0.5d0*dble(l+m)*(1d0+dble(l)/dble(1-m*m))
c2=-dble((l+m)*(l-m+1)*(l+m-2)*(l+m-1))/dble(4*m-4)
dalegendre1=c0*legendre(l-1,m+2,x)+c1*legendre(l-1,m,x)+c2*legendre(l-1,m-2,x)
endif

return
end function dalegendre1
function calegendre1(l,m,x)
implicit none
integer,intent(in)::m,l
complex(kind(0d0)),intent(in)::x
complex(kind(0d0))::calegendre1

double precision::c0,c1,c2

calegendre1=dcmplx(0d0,0d0)
if(abs(m).gt.l)then
calegendre1=dcmplx(0d0,0d0)
return
endif

if(abs(m).eq.1)then
calegendre1=(dble(l)*x*legendre(l,m,x)-dble(l+m)*legendre(l-1,m,x))/(x*x-1)
else
c0=1d0/dble(4*m+4)
c1=0.5d0*dble(l+m)*(1d0+dble(l)/dble(1-m*m))
c2=-dble((l+m)*(l-m+1)*(l+m-2)*(l+m-1))/dble(4*m-4)
calegendre1=c0*legendre(l-1,m+2,x)+c1*legendre(l-1,m,x)+c2*legendre(l-1,m-2,x)
endif

return
end function calegendre1

function dalegendre2(l,m,x)
implicit none
integer,intent(in)::m,l
double precision,intent(in)::x
double precision::dalegendre2

double precision::c0,c1,c2

dalegendre2=0d0
if(abs(m).gt.l)then
dalegendre2=0d0
return
endif

if(abs(m).eq.1.or.abs(m).eq.3)then
dalegendre2=(2d0*x*legendre1(l,m,x)&
+dble((l+1)*(l+m))/dble(2*l+1)*legendre1(l-1,m,x)&
-dble(l*(l-m+1))/dble(2*l+1)*legendre1(l+1,m,x))/(1-x*x)
else
c0=1d0/dble(4*m+4)
c1=0.5d0*dble(l+m)*(1d0+dble(l)/dble(1-m*m))
c2=-dble((l+m)*(l-m+1)*(l+m-2)*(l+m-1))/dble(4*m-4)
dalegendre2=c0*legendre1(l-1,m+2,x)&
+c1*legendre1(l-1,m,x)+c2*legendre1(l-1,m-2,x)
endif

return
end function dalegendre2
function calegendre2(l,m,x)
implicit none
integer,intent(in)::m,l
complex(kind(0d0)),intent(in)::x
complex(kind(0d0))::calegendre2

double precision::c0,c1,c2

calegendre2=dcmplx(0d0,0d0)
if(abs(m).gt.l)then
calegendre2=dcmplx(0d0,0d0)
return
endif

if(abs(m).eq.1.or.abs(m).eq.3)then
calegendre2=(dble(l)*x*legendre(l,m,x)-dble(l+m)*legendre(l-1,m,x))/(x*x-1)
else
c0=1d0/dble(4*m+4)
c1=0.5d0*dble(l+m)*(1d0+dble(l)/dble(1-m*m))
c2=-dble((l+m)*(l-m+1)*(l+m-2)*(l+m-1))/dble(4*m-4)
calegendre2=c0*legendre(l-1,m+2,x)+c1*legendre(l-1,m,x)+c2*legendre(l-1,m-2,x)
endif

return
end function calegendre2

!=========================

function djacobi1(n,alpha,beta,x)
implicit none
double precision,intent(in)::alpha,beta
integer,intent(in)::n
complex(kind(0d0)),intent(in)::x
complex(kind(0d0))::djacobi1
double precision::c0

c0=log(gamma(alpha+beta+dble(n+2))-log(gamma(alpha+beta+dble(n+1)))-log(2d0))
djacobi1=c0*jacobi(n-1,alpha+1d0,beta+1d0,x)

return
end function djacobi1
function cjacobi1(n,alpha,beta,x)
implicit none
integer,intent(in)::n
complex(kind(0d0)),intent(in)::alpha,beta,x
complex(kind(0d0))::cjacobi1,c0

c0=clgamma(alpha+beta+dble(n+2))-clgamma(alpha+beta+dble(n+1))-log(2d0)
cjacobi1=c0*jacobi(n-1,alpha+1d0,beta+1d0,x)

return
end function cjacobi1

function djacobi2(n,alpha,beta,x)
implicit none
double precision,intent(in)::alpha,beta
integer,intent(in)::n
complex(kind(0d0)),intent(in)::x
complex(kind(0d0))::djacobi2
double precision::c0

c0=log(gamma(alpha+beta+dble(n+3)))-log(gamma(alpha+beta+dble(n+1)))-2d0*log(2d0)
djacobi2=c0*jacobi(n-2,alpha+2d0,beta+2d0,x)

return
end function djacobi2
function cjacobi2(n,alpha,beta,x)
implicit none
integer,intent(in)::n
complex(kind(0d0)),intent(in)::alpha,beta,x
complex(kind(0d0))::cjacobi2,c0

c0=clgamma(alpha+beta+dble(n+3))-clgamma(alpha+beta+dble(n+1))-2d0*log(2d0)
cjacobi2=c0*jacobi(n-2,alpha+2d0,beta+2d0,x)

return
end function cjacobi2

function dhermite(n,x)
implicit none
integer,intent(in)::n
double precision,intent(in)::x
double precision::dhermite

integer::m,k
double precision::dn,dm,bm,cm,dnx,fx

if(n.le.-1)then
dhermite=0d0
return
endif

fx=x*x
dhermite=1d0
if(mod(n,2).eq.0)then
k=n/2
dn=dble(k)
dnx=exp(log(gamma(2d0*dn+1d0))-log(gamma(dn+1d0)))
do m=k,1,-1
dm=dble(m); cm=dm*(2d0*dm-1d0); bm=2d0*dble(k-m+1)
dhermite=1d0-dhermite*fx*bm/cm
enddo
else
k=(n-1)/2
dn=dble(k)
dnx=2d0*x*exp(log(gamma(2d0*dn+2d0))-log(gamma(dn+1d0)))
do m=k,1,-1
dm=dble(m); cm=dm*(2d0*dm+1d0); bm=2d0*dble(k-m+1)
dhermite=1d0-dhermite*fx*bm/cm
enddo
endif
dhermite=dhermite*dnx
if(mod(k,2).eq.1)dhermite=-dhermite

return
end function dhermite

function dhermite1(n,x)
implicit none
integer,intent(in)::n
double precision,intent(in)::x
double precision::dhermite1

dhermite1=2*n*dhermite(n-1,x)

return
end function dhermite1

function dhermite2(n,x)
implicit none
integer,intent(in)::n
double precision,intent(in)::x
double precision::dhermite2

dhermite2=dble(4*n)*dble(n-1)*dhermite(n-2,x)

return
end function dhermite2

function chermite(n,z)
implicit none
integer,intent(in)::n
complex(kind(0d0)),intent(in)::z
complex(kind(0d0))::chermite

integer::m,k
double precision::dn,dm,bm,cm
complex(kind(0d0))::fx,dnx

if(n.le.-1)then
chermite=dcmplx(1d0,0d0)
return
endif

fx=z*z
chermite=dcmplx(1d0,0d0)
if(mod(n,2).eq.0)then
k=n/2
dn=dble(k)
dnx=exp(log(gamma(2d0*dn+1d0))-log(gamma(dn+1d0)))
do m=k,1,-1
dm=dble(m); cm=dm*(2d0*dm-1d0); bm=2d0*dble(k-m+1)
chermite=1d0-chermite*fx*bm/cm
enddo
else
k=(n-1)/2
dn=dble(k)
dnx=2d0*z*exp(log(gamma(2d0*dn+2d0))-log(gamma(dn+1d0)))
do m=k,1,-1
dm=dble(m); cm=dm*(2d0*dm+1d0); bm=2d0*dble(k-m+1)
chermite=1d0-chermite*fx*bm/cm
enddo
endif
chermite=chermite*dnx
if(mod(k,2).eq.1)chermite=-chermite

return
end function chermite

function chermite1(n,z)
implicit none
integer,intent(in)::n
complex(kind(0d0)),intent(in)::z
complex(kind(0d0))::chermite1

chermite1=2*n*chermite(n-1,z)

return
end function chermite1

function chermite2(n,z)
implicit none
integer,intent(in)::n
complex(kind(0d0)),intent(in)::z
complex(kind(0d0))::chermite2

chermite2=dble(4*n)*dble(n-1)*chermite(n-2,z)

return
end function chermite2

function clgamma(z)
implicit none
complex(kind(0d0)),intent(in)::z
complex(kind(0d0))::clgamma

clgamma=log(cgamma(z))

return
end function clgamma

function cgamma(z)
implicit none
complex(kind(0d0)),intent(in)::z
complex(kind(0d0))::cgamma
!
! sikinote (http://slpr.sakura.ne.jp/qp)
! author : sikino
! date : 2017/01/09
! 2017/01/10 optimize the number of a(M)
!
integer::i,n,M
double precision,parameter::pi=3.14159265358979324d0
double precision,parameter::e1=0.3678794411714423216d0
double precision,parameter::ln2pi2=0.91893853320467274d0
double precision,parameter::sq2pi=2.5066282746310005d0
double precision::a(1:30),AB(1:14),dz,iz
complex(kind(0d0))::q,s,w,r,q1,q2

dz=dble(z)
iz=imag(z)

if(dz.eq.nint(dz).and.iz.eq.0d0)then
if(dz.gt.0d0)then
s=dcmplx(1d0,0d0)
n=nint(dz)
do i=2,n-1
s=s*i
enddo
else
s=1d+300
endif
cgamma=s
else
q=z
if(dz.lt.0d0)q=1d0-z

if(abs(iz).lt.1.45d0)then
! For x=arb, |y|\lt 1.45
n=int(dble(q))
s=1d0
do i=1,n
s=s*(q-i)
enddo
w=q-dble(n)

a(1) = 1d0
a(2) = 0.57721566490153286d0
a(3) =-0.65587807152025388d0
a(4) =-0.42002635034095236d-1
a(5) = 0.16653861138229149d0
a(6) =-0.42197734555544337d-1
a(7) =-0.96219715278769736d-2
a(8) = 0.72189432466630995d-2
a(9) =-0.11651675918590651d-2
a(10)=-0.21524167411495097d-3
a(11)= 0.12805028238811619d-3
a(12)=-0.20134854780788239d-4
a(13)=-0.12504934821426707d-5
a(14)= 0.11330272319816959d-5
a(15)=-0.20563384169776071d-6
a(16)= 0.61160951044814158d-8
a(17)= 0.50020076444692229d-8
a(18)=-0.11812745704870201d-8
a(19)= 0.10434267116911005d-9
a(20)= 0.77822634399050713d-11
a(21)=-0.36968056186422057d-11
a(22)= 0.51003702874544760d-12
a(23)=-0.20583260535665068d-13
a(24)=-0.53481225394230180d-14
a(25)= 0.12267786282382608d-14
a(26)=-0.11812593016974588d-15
a(27)= 0.11866922547516003d-17
a(28)= 0.14123806553180318d-17
a(29)=-0.22987456844353702d-18
a(30)= 0.17144063219273374d-19

M=int(11.3*abs(w)+13)
if(M.gt.30)M=30
r=a(M)
do i=M-1,1,-1
r=r*w+a(i)
enddo

cgamma=s/(r*w)
else
! For x=arb, |y|\gt 1.45
s=1d0
if(abs(q).lt.9d0)then
do i=0,7
s=s*(q+i)
enddo
s=1d0/s
q=q+8d0
endif

AB(1) = 0.83333333333333333d-1
AB(2) =-0.27777777777777778d-2
AB(3) = 0.79365079365079365d-3
AB(4) =-0.59523809523809524d-3
AB(5) = 0.84175084175084175d-3
AB(6) =-0.19175269175269175d-2
AB(7) = 0.64102564102564103d-2
AB(8) =-0.29550653594771242d-1

q1=1d0/q
q2=q1*q1

r=AB(8)
do i=7,1,-1
r=r*q2+AB(i)
enddo
cgamma=s*exp((q-0.5d0)*log(q)-q+ln2pi2+r*q1)
endif
if(dz.lt.0d0)cgamma=pi/(cgamma*sin(pi*z))
endif

return
end function cgamma
end module specialfunction

プログラム例

呼び出す際は総称名
関数：Laguerre(n,a,x), legendre(l,m,x), jacobi(n,a,b,x)
1階微分：Laguerre1(n,a,x), legendre1(l,m,x), jacobi1(n,a,b,x)
2階微分：Laguerre2(n,a,x), legendre2(l,m,x), jacobi2(n,a,b,x)
を用いてください。

ラゲール陪多項式$L_n^{(a)}(x)$の引数の型(n,a,x)の組み合わせは、
(n,a,x) : 整数型、　整数型、　　倍精度型　:　出力は倍精度型
(n,a,x) : 整数型、倍精度型、　　倍精度型　:　出力は倍精度型
(n,a,x) : 整数型、倍精度型、複素倍精度型　:　出力は複素倍精度型

ルジャンドル陪多項式$P_l^m(x)$の引数の型(l,m,x)の組み合わせは、
(l,m,x) : 整数型、整数型、　　倍精度型　:　出力は倍精度型
(l,m,x) : 整数型、整数型、複素倍精度型　:　出力は複素倍精度型

ヤコビ多項式$P^{(a,b)}_n(x)$の引数の型(n,a,b,x)の組み合わせは、
(n,a,b,x) : 整数型、　　倍精度型、　　倍精度型、複素倍精度型　:　出力は複素倍精度型
(n,a,b,x) : 整数型、複素倍精度型、複素倍精度型、複素倍精度型　:　出力は複素倍精度型

というものが有効になっています。

実際に全部使っているプログラムはこちら。

program main
use specialfunction
implicit none

write(6,*)laguerre(4,2,4.78d0)
write(6,*)laguerre(3,2.5d0,4.78d0)
write(6,*)laguerre(8,dcmplx(1.6d0,3.01d0),dcmplx(5.64d0,49.1d0))

write(6,*)legendre(8,5,0.28d0)
write(6,*)legendre(6,2,dcmplx(1.8d0,2.7d0))

write(6,*)jacobi(3,2.15d0,1.85d0,dcmplx(4d0,8.1d0))
write(6,*)jacobi(3,dcmplx(1.2d0,0.5d0),dcmplx(5d0,0.3d0),dcmplx(4d0,8.1d0))

stop
end program main

引数は適当に入れています。
spfunc.f90にmoduleを入れ、
temp.f90に上の12行のコードを書きコンパイル→実行すると、

$ ifort spfunc.f90 temp.f90
$ ./a.out
3.2997056066666568
-8.4458666666671553E-002
( 454556598.95966828 , 376437598.22472727 )
9378.9803534037801
( -519656.04147187521 , 135464.19972749997 )
( -10936.767937499992 , -2092.2704999999974 )
( -18248.193121999979 , -10218.566324666655 )
$

と言う結果を得ます。
wolfram alphaと確かめてみましょう。計算で入れた値はそれぞれ、
$
\begin{align}
& L_4^2(4.78) \\
& L_3^{2.5}(4.78) \\
& L_8^{1.6+i3.01}(5.64+i49.1) \\
& P_8^{5}(0.28) \\
& P_6^{2}(1.8+i2.7) \\
& P_3^{(2.15,1.85)}(4+i8.1) \\
& P_3^{(1.2+i0.5,5+i0.3)}(4+i8.1) \\
\end{align}
$
です。

wolframの結果と比較したのがこちら。

下線部が上記プログラムとwolframとで一致した桁です。
まぁまぁいいんじゃないでしょうか。
$L_8^{1.6+i3.01}(5.64+i49.1)$の結果が著しく良くないのが気になります･･･使う際は安定な範囲を確かめてから使ってください。

以下、ここで示している特殊関数の定義と関係式です。

ラゲール陪多項式

定義域

$\;\;\;\;
0 \lt x \lt \infty
$
微分方程式

$\;\;\;\;
\displaystyle
\left[x\frac{d^2}{dx^2}+(\alpha+1-x)\frac{d}{dx}+n\right]L_m^{(\alpha)}(x)=0
$
具体例

$\;\;\;\;
\begin{align}
& L^{(0)}_0(x)=1 \\
& L^{(0)}_1(x)=-x+1 \\
& L^{(0)}_2(x)=\frac{1}{2}(x^2-4x+2) \\
& L^{(0)}_3(x)=\frac{1}{6}(-x^3+9x^2-18x+6) \\
& L^{(1)}_0(x)=1 \\
& L^{(1)}_1(x)=-x+2 \\
& L^{(1)}_2(x)=\frac{1}{2}\left(x^2-6x+6\right) \\
& L^{(1)}_3(x)=\frac{1}{6}\left(-x^3+12x^2-36x+24\right)
\end{align}
$
漸化式

$\;\;\;\;
\begin{align}
L_n^{(\alpha)}(x)&=\binom{n+\alpha}{n}a_0(x) \\
&\hspace{1em}a_{m-1}(x)=1-\frac{n-m+1}{m(\alpha+m)}\cdot x\cdot a_m(x) \\
&\hspace{5em}(m=n,n-1,n-2,\cdots,1,\;\;a_n(x)=1)
\end{align}
$
[1]
1階微分
$\;\;\;\;
\begin{align}
\frac{d}{dx}L_n^{(\alpha)}(x)=-L_{n-1}^{(\alpha+1)}(x)
\end{align}
$
2階微分
$\;\;\;\;
\begin{align}
\frac{d^2}{dx^2}L_n^{(\alpha)}(x)=L_{n-2}^{(\alpha+2)}(x)
\end{align}
$
直交性

$\;\;\;\;
\displaystyle
\int_{0}^{\infty} e^{-x}x^{\alpha}L^{(\alpha)}_n(x)L^{(\alpha)}_{n’}(x) dx =\frac{\Gamma(\alpha+n+1)}{n!}\delta_{n,n’},\;\;(\alpha>-1)
$
[2]

ルジャンドル陪多項式

定義域

$\;\;\;\;
1 \lt x \lt -1
$
微分方程式

$\;\;\;\;
\displaystyle (1-x^2)\frac{d^2}{dx^2}P_l^m(x)-2x\frac{d}{dx}P_l^m(x)+\left[l(l+1)-\frac{m^2}{1-x^2}P_l^m(x)\right]=0
$
具体例

$\;\;\;\;
\begin{align}
& P_0^0(x)=1 \\
& P_{1}^{-1}(x)=\frac{1}{2}(1-x^2)^{1/2} \\
& P_{1}^{0}(x)=x \\
& P_{1}^{1}(x)=-(1-x^2)^{1/2} \\
& P_{2}^{-2}(x)=\frac{1}{8}(1-x^2) \\
& P_{2}^{-1}(x)=\frac{1}{2}x(1-x^2)^{1/2} \\
& P_{2}^{0}(x)=\frac{1}{2}(3x^2-1) \\
& P_{2}^{1}(x)=-3(1-x^2)^{1/2} \\
& P_{2}^{2}(x)=3(1-x^2)
\end{align}
$
漸化式

用いている式
$\;\;\;\;
\begin{align}
& P_{|m|}^{|m|}(x)=(-1)^{|m|}(2|m|-1)!!(1-x^2)^{\frac{|m|}{2}} \\
& P_{l}^{-m}(x)=(-1)^m\frac{(l-m)!}{(l+m)!}P_l^m(x) \\
& P_{|m|}^{|m|+1}(x)=(2|m|+1)xP_{|m|}^{|m|}(x) \\
& P_{|m|}^{|m|+q}(x)=\left(\frac{2|m|-1}{q}+2\right)xP_{|m|}^{|m|+q-1}(x)-\left(\frac{2|m|-1}{q}+1\right)P_{|m|}^{|m|+q-2}(x)
\end{align}
$[3]
1階微分
$\;\;\;\;
\begin{align}
m=\pm 1&\; \\
&:\hspace{1em} \frac{d}{dx}P_l^m(x)=\frac{lx}{(x^2-1)}P_l^m(x)-\frac{l+m}{(x^2-1)} P_{l-1}^m(x) \\
m\ne\pm 1\; \\
&:\hspace{1em} \frac{d}{dx}P_l^m(x)=c_{2}P^{m+2}_{l-1}(x)+c_0 P^{m}_{l-1}(x)+c_{-2} P^{m-2}_{l-1}(x) \\
&\hspace{1em} c_{2}=\frac{1}{4(m+1)}\\
&\hspace{1em} c_{0}=\frac{l+m}{2}\left(1+\frac{l}{1-m^2}\right)\\
&\hspace{1em} c_{-2}=-\frac{(l+m)(l+m-1)(l+m-2)(l-m+1)}{4(m-1)}
\end{align}
$
※$m=\pm 1$の時、$x=\pm 1$で発散します。
2階微分
$\;\;\;\;
\begin{align}
m=\pm 1,\pm 3&\; \\
&:\hspace{1em} \frac{d^2}{dx^2}P_l^m(x)=\frac{1}{1-x^2}\left[2x\frac{dP_l^m(x)}{dx}+\frac{(l+1)(l+m)}{2l+1}\frac{dP_{l-1}^m(x)}{dx}-\frac{l(l-m+1)}{2l+1}\frac{dP_{l+1}^m(x)}{dx}\right] \\
m\ne\pm 1,\pm 3&\; \\
&:\hspace{1em} \frac{d^2}{dx^2}P_l^m(x)=c_{2}\frac{dP^{m+2}_{l-1}(x)}{dx}+c_0 \frac{dP^{m}_{l-1}(x)}{dx}+c_{-2} \frac{dP^{m-2}_{l-1}(x)}{dx}
\end{align}
$
※$m=\pm 1,\pm 3$の時、$x=\pm 1$で発散します。
直交性

$l$に対する直交性
$\;\;\;\;
\displaystyle
\int_{-1}^{1} P_{m}^{l}(x)P_{m}^{l’}(x) dx =\frac{2}{2l+1}\cdot \frac{(l+m)!}{(l-m)!}\delta_{l,l’}
$
[2]

ヤコビ多項式

定義域

$\;\;\;\;
-1 \lt x \lt 1
$
微分方程式

$\;\;\;\;
\displaystyle
\left[(1-x^2)\frac{d^2}{dx^2}+(\beta-\alpha-(\alpha+\beta+2)x)\frac{d}{dx}+n(n+\alpha+\beta+1)\right]P_n^{(\alpha,\beta)}(x)=0
$
漸化式

$\;\;\;\;
\begin{align}
P_n^{(\alpha,\beta)}(x)&=\binom{n+\alpha}{n}a_0(x) \\
&\hspace{1em}a_{m-1}(x)=1-\frac{(n-m+1)(\alpha+\beta+n+m)}{2m(\alpha+m)}(1-x)\cdot a_m(x)\\
&\hspace{5em}(m=n,n-1,n-2,\cdots,1,\;\;a_n(x)=1)
\end{align}
$
[1]
1階微分
$\;\;\;\;
\begin{align}
\frac{d}{dx}P_n^{(\alpha,\beta)}(x)=\frac{\Gamma(\alpha+\beta+n+2)}{2\Gamma(\alpha+\beta+n+1)}P_{n-1}^{(\alpha+1,\beta+1)}(x)
\end{align}
$
2階微分
$\;\;\;\;
\begin{align}
\frac{d^2}{dx^2}P_n^{(\alpha,\beta)}(x)=\frac{\Gamma(\alpha+\beta+n+3)}{2^2\Gamma(\alpha+\beta+n+1)}P_{n-2}^{(\alpha+2,\beta+2)}(x)
\end{align}
$
直交性

$\;\;\;\;
\begin{align}
\displaystyle
&\int_{-1}^{1} (1-x)^{\alpha}(1+x)^{\beta}P_{n}^{(\alpha,\beta)}(x)P_{n’}^{(\alpha,\beta)}(x)dx \\
&\hspace{4em}=\frac{2^{\alpha+\beta+1}}{2n+\alpha+\beta+1}
\frac{\Gamma(n+\alpha+1)\Gamma(n+\beta+1)}{n!\Gamma(n+\alpha+\beta+1)}\delta_{nm}\\
&(\alpha>-1,\ \beta>-1)
\end{align}
$[2]

参考文献

[1] Abramowitz & Stegun (1965), ISBN 978-0486612720 (Numerical Method, 22.18, Page 789)
[2] Abramowitz & Stegun (1965), ISBN 978-0486612720 (Numerical Method, 22.2, Page 774,775)
[3] Associated Legendre polynomials -wikipedia

テキストエディタ, プログラミングと数値計算

emacs上で再読み込み

2016年2月20日 sikino コメントする

emacs上でファイルを再読み込みしたいとします。
OSはLinux mint17.1です。

例えば、
emacsであるファイルを開いたまま、sedコマンドでファイルをいじると、そのいじったファイルは現在開いているemacsには反映されません。

これを変更してF5キーを押すだけでemacs上で再読み込みができるようにします。

ホームディレクトリに
.emacs
というファイルを作成し、その中に↓

(defun revert-buffer-no-confirm (&optional force-reverting)
"Interactive call to revert-buffer. Ignoring the auto-save
file and not requesting for confirmation. When the current buffer
is modified, the command refuses to revert it, unless you specify
the optional argument: force-reverting to true."
(interactive "P")
;;(message "force-reverting value is %s" force-reverting)
(if (or force-reverting (not (buffer-modified-p)))
(revert-buffer :ignore-auto :noconfirm)
(error "The buffer has been modified")))

;; reload buffer
(global-set-key (kbd "<f5>") 'revert-buffer-no-confirm)

と入力します。これを記述し、もう一度emacsを立ち上げ,sedコマンド等で外部から編集し、F5キーを押せば更新されます。

参考文献

Emacsで現在開いてるバッファを再読込する

プログラミングと数値計算

sedコマンドによるファイルの書き換え

2016年2月20日 sikino コメントする

シェルスクリプトによってファイルを書き換えます。

sedコマンドの概要
具体例一覧

sedコマンドの概要

対象とするファイル名”input.ini”の中身に

&input
a=7,
b=13,
c='./dat/',
&end

という文が書いてあるとします。
これのファイルの中身を端末からのコマンドで書き換えるためにはsedコマンドを利用します。

sedコマンドとは、置換用のコマンドです。
例えば、コマンドライン上で

sed -e "s/a=7/a=9/" ./input.ini

と入力すると”input.ini”の中身の”a=7″に一致する場所を”a=9″に置換して出力します。
ここで重要なのが、”a=7″に一致する場所だけなので、後ろにあるコンマは影響されません。
実際に実行してみると

$ sed -e "s/a=7/a=9/" ./input.ini
&input
a=9,
b=13,
c='./dat/',
&end

となってコマンドライン上に出力されます。
この時、書き換わった内容はコマンドラインのみに出力されるため、元のファイル自体に書き換えはされません。

元のファイルの書き換えを行いたければ、出力先を適当なファイル、例えば”tmp”というファイルにし、それを”input.ini”にリネームすればいいです。
すなわち、

sed -e "s/a=7/a=9/" ./input.ini > tmp
mv tmp input.ini

とすればいいのです。実際に実行して中身を見ますと

$ sed -e "s/a=7/a=9/" ./input.ini > tmp
$ mv tmp input.ini
$ cat input.ini
&input
a=9,
b=13,
c='./dat/',
&end

となります。

sedコマンドの概要はここで終了です。
以降は具体例とシェルスクリプトのコード(ファイル名”q.sh”),実行例を載せていきます。

特定の数値を変えたい場合

対象とするファイル”input.ini”

&input
a=7,
b=13,
c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
sed -e "s/a=7/a=9/" ./input.ini > tmp
mv tmp input.ini

実行例

$ sh q.sh
$ cat input.ini
&input
a=9,
b=13,
c='./dat/',
&end

特殊文字(/等)を含む文字列の変更

sedコマンドを使うために/を使います。ならば/を置換したい場合どうすればいいでしょう。
この時はsedコマンドで使っている区切りを表す/を別の特殊文字(例えば%)に置き換えればよいです。

対象とするファイル”input.ini”

&input
a=7,
b=13,
c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
sed -e "s%c='./dat/'%c='./data/'%" input.ini > tmp
mv tmp input.ini

実行例

$ sh q.sh
$ cat input.ini
&input
a=7,
b=13,
c='./data/',
&end

置換する数値の変数での指定

対象とするファイル”input.ini”

&input
a=7,
b=13,
c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
pr=5
sed -e "s/a=7/a=${pr}/" ./input.ini > tmp
mv tmp input.ini

実行例

$ sh q.sh
$ cat input.ini
&input
a=5,
b=13,
c='./dat/',
&ends

ある数値以降を全て変更したい

コンマも含めて消したいとします。
対象とするファイル”input.ini”

&input
a=7,
b=13,
c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
sed -e "s/c=.*/c=12345/" ./input.ini > tmp
mv tmp input.ini

実行例

$ sh q.sh
$ cat input.ini
&input
a=7,
b=13,
c=12345
&end

ファイルの空白行の削除

ファイルに含まれる空白行を全て削除します。
対象とするファイル”input.ini”

&input
a=7,

b=13,

c='./dat/',

&end

“q.sh”の中身

#!/bin/sh
sed -e '/^ *$/d' input.ini > tmp
mv tmp input.ini

実行例

$ sh q.sh
$ cat input.ini
&input
a=7,
b=13,
c='./dat/',
&end

ファイルの行頭に数値を入れる

対象とするファイル”input.ini”

&input
a=7,
b=13,
c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
i=1
while [ ${i} -le 1000 ]
do
sed -e "${i}s/^/${i} /" input.ini > temp
mv temp input.ini

i=`expr $i + 1`
done

実行例

$ sh q.sh
$ cat input.ini
1 &input
2 a=7,
3 b=13,
4 c='./dat/',
5 &end
6

“q.sh”中、1000をファイルに応じて変更してください。

※また、行番号を入れるだけの場合、catコマンドを利用して、

cat -b input.ini > tmp
mv tmp input.ini

とするのが簡単だと思います。
オプション”-b”は空行を含めず行番号を付けるという意味です。空行を含めたい場合は”-n”でできます。

ファイルの最後に改行を入れる

対象とするファイル”input.ini”

&input
a=7,
b=13,
c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
sed '$a\\n' input.ini > tmp
mv tmp input.ini

実行例

$ sh q.sh
$ cat input.ini
&input
a=7,
b=13,
c='./dat/',
&end

$

※どうしても2行分改行されてしまうようです。
sedコマンドで1行分の改行は、どうすればよいか分かりませんでした。

苦肉の策ですが、適当なファイル(例えばnn.txt)に改行だけを入力し、catコマンドでファイルを結合すると一応できます。
—nn.txt—

————
を作っておいて、
#!/bin/sh
cat input.ini nn.txt > tmp
mv tmp input.ini
とすれば”input.ini”の後ろに”nn.txt”が結合されるため、1行分の改行が可能です。

do文で何回も変更したい

対象とするファイル”input.ini”

&input
a=7,
b=13,
c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
for pr in 0 1 2
do
sed -e "s/a=.*/a=${pr},/" input.ini > tmp
mv tmp input.ini

cat input.ini
done

実行例

$ sh q.sh
&input
a=0,
b=13,
c='./dat/',
&end
&input
a=1,
b=13,
c='./dat/',
&end
&input
a=2,
b=13,
c='./dat/',
&end

全く違う変数を繰り返し代入したい

対象とするファイル”input.ini”

&input
a=7,
b=13,
c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
i=1
for j in 3d0 939d0
do
echo "count -> ${i}"
sed -e "s/a=.*/a=${j},/" ./input.ini > tmp
mv tmp input.ini

grep "a=" input.ini

i=`expr $i + 1`
done

実行例

$ sh q.sh
count -> 1
a=3d0,
count -> 2
a=939d0,
$

シェルスクリプトの中の”grep”は”a=”に一致する行を書き出させるものです。

整数の変数を規則的に何回も変更したい

対象とするファイル”input.ini”

&input
a=7,
b=13,
c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
i=1
while [ ${i} -le 3 ]
do
echo "count -> ${i}"
sed -e "s/a=.*/a=${i},/" ./input.ini > tmp
mv tmp input.ini

grep "a=" input.ini

i=`expr $i + 1`
done

実行例

$ sh q.sh
count -> 1
a=1,
count -> 2
a=2,
count -> 3
a=3,
$

小数点を含む変数を規則的に何回も変更したい

小数点を含む演算は”expr”コマンドではできません。
シェルスクリプトで数値計算を行う”bc”コマンドを用いましょう。
ここでは、変数$v$に、
$v=i\times0.7,\;\;(i=1,2,3)$
を代入する演算を考えましょう。

対象とするファイル”input.ini”

&input
a=3,
b=13,
c='./dat/',
&end

“q.sh”の中身

#!/bin/sh
i=1
while [ ${i} -le 3 ]
do
echo "count -> ${i}"

v=$(echo "scale=3; $i*0.70 " | bc | sed -e 's/^\./0./g')
echo "$v"

sed -e "s/a=.*/a=${v},/" ./input.ini > tmp
mv tmp input.ini

grep "a=" input.ini

i=`expr $i + 1`
done

実行例

$ sh q.sh
count -> 1
0.70
a=0.70,
count -> 2
1.40
a=1.40,
count -> 3
2.10
a=2.10,
$

あるディレクトリ以下にある全ファイル名を変数にしたい

状況

.
├── a.sh
├── dat
├── data
│   ├── asdfgh.d
│   ├── qwerty.d
│   └── zxcvbn.d
└── input

inputの中身

&input
filename="./data/xxxxx.d",
key = 0,
&end

a.shの中身

#!/bin/sh
i=1
for j in ./data/*
do
sed -e "s%filename=.*%filename="${j}",%" input > tmp
mv tmp ./dat/${j##*/}_ex
echo ${j##*/}
grep "filename=" ./dat/${j##*/}_ex
done

———————-

a.sh実行

$sh a.sh
asdfgh.d
filename="./data/asdfgh.d",
qwerty.d
filename="./data/qwerty.d",
zxcvbn.d
filename="./data/zxcvbn.d",

a.sh実行後の構成

.
├── a.sh
├── dat
│   ├── asdfgh.d_ex
│   ├── qwerty.d_ex
│   └── zxcvbn.d_ex
├── data
│   ├── asdfgh.d
│   ├── qwerty.d
│   └── zxcvbn.d
└── input

例えばdat/asdfgh.d_exの中身

&input
filename="./data/asdfgh.d",
key = 0,
&end

参考先

フィルタを使用した文字列操作 1 -UNIX & Linux コマンド・シェルスクリプトリファレンス
 sed の使い方 -サーバエンジニアの知恵袋
 制御構文 -シェルスクリプト入門
 How to pass results of bc to a variable – Ask Ubuntu
bcによる少数の演算 -SWEng TIPs

プログラミングと数値計算

openMPによる並列計算

2016年2月15日 sikino コメントする

まとめ

fortran90では、

program main
implicit none
integer::i,omp_get_thread_num

! number of threads for parallel computation
call omp_set_num_threads(2)

!$omp parallel do
do i=1,100
write(6,*) omp_get_thread_num()
enddo
!$omp end parallel do

stop
end program main

とすれば2個のスレッドで並列計算され、コンパイル時に
intel® fortranでは、

ifort -openmp main.f90

gfortranでは、

gfortran -fopenmp main.f90

とすれば良いです。

fortran90でのお話です。

openMPは複数のCPUを使って並列計算をするインターフェースのことです
C/C++,fortranで対応しています。
何も指定しない場合、計算は基本的に1つのコアのみを使います。
Linuxを使っている場合、恐らくデフォルトで使えるようになっているはずです。

並列計算を行うまでの大まかな手順は2つ。

プログラム内部に文章を書く
コンパイル時にオプションを追加
- intel®fortranコンパイラの場合
- gfortranコンパイラの場合

です。

※今使っているパソコンの並列計算できるCPUの数の確認は、

cat /proc/cpuinfo | grep processor

です。proc/cpuinfoの場所は使っているOSによって違う可能性があります。
”linux mint cpu 確認”
等とgoogle等で検索するのが手っ取り早いでしょう。

プログラム内部に文章を書く

プログラム内でopenMPを使用するdoループを挟むように、

!$omp parallel do
do i=1,100
call test(i)
enddo
!$omp end parallel do

と書きます。基本はこれだけです。
並列計算に使えるスレッドが2個である場合、並列計算というのは

スレッド0

do i=1,50
call test(i)
enddo

スレッド1

do i=51,100
call test(i)
enddo

と割り振る操作なのです。

上記のようにdoループを分割するため、注意しなければならない点として、
doループに用いる変数に対して別々に計算しても計算結果が変わらない事
が挙げられます。
上記例では、i=0から始めて、i=1,2,3,･･･と順次計算する必要は無く、i=5,10,1,3,･･･みたいな乱雑な順番でも計算結果が変わらない状況で使わなければなりません。
特に、i=2を計算したい時にi=1の計算結果を使う場合は並列化することは出来ません。

実際には同期したりとかで、出来なくはないですが面倒くさく、僕はやったことが無いので書きません。

もしも並列計算したいスレッド数を制限したい場合、

!$omp parallel do

を

!$omp parallel do num_threads(2)

のように指定するか(この例ではその部分の並列計算を2つのスレッドで行う)、
並列計算が始まる前のプログラムの初めの方に

call omp_set_num_threads(2)

と書きましょう。これを書くと全ての並列計算箇所が2個のスレッドで行われることになります。

また、関数omp_get_thread_num()はスレッドにつけられた番号を返す事が出来ます。引数は要りません。
なので、これを利用すると何個のスレッドが使われているのか直に確認できます。

補足ですが、MKL(マス・カーネル・ライブラリ)のルーチン内で使われる並列計算のスレッド数は

call mkl_set_num_threads(2)

と書くことで制御できます。

コンパイル時のオプション

intel® fortranでは、

ifort -openmp main.f90

gfortranでは、

gfortran -fopenmp main.f90

としてコンパイル、実行しましょう。

実例

例として

program main
implicit none
integer::i,omp_get_thread_num

! number of threads for parallel computation
call omp_set_num_threads(2)

!$omp parallel do
do i=1,4
write(6,*) omp_get_thread_num()
enddo
!$omp end parallel do

stop
end program main

を動かすと同出力されるか見てみましょう。
gfortranコンパイラで、使えるスレッドの数は2つです。
上記のコードを

gfortran -fopenmp main.f90

としてコンパイルし、実行すると

./a.out
1
0
0
1

という結果が得られます。
ここに出ている数字はスレッドの番号で、0,1と番号付けされたスレッドが確かに出力されている事が分かります。1,0,0,1の順番はどうでも良いです。実行するごとに変わるので、これは同時並行処理されているためでしょう。

call omp_set_num_threads(1)

として計算を行うと、

./a.out
0
0
0
0

と表示されることから、確かに使うスレッド数が制御されていることが分かります。

doループの並列処理を処理が終わった順に行う場合

doループ内で、あるときだけ計算が重くなる場合があります。
この場合、均等にdoループの変数を分けてしまうとCPU1つあたりの処理が大きく異なってしまいます。そこで、doループの変数を終わった順に処理していく指定方法があります。それがschedule構文です。以下のようにすると、i=1から計算が終わった順にiを増やしていきます。

program main
implicit none
integer::i,omp_get_thread_num,N

! number of threads for parallel computation
call omp_set_num_threads(2)

N=100
!$omp parallel do schedule(dynamic,1)
do i=1,N
write(6,*) omp_get_thread_num()
enddo
!$omp end parallel do

stop
end program main

とします。ただし、1回あたり計算時間がほとんどかからない場合、オーバーヘッド時間が無視できなくなるため、使うときは吟味してから使用するべきです。

また、並列化しても共通の変数を使うことができます。以下のようにすると、計算の進み具合がわかります。

program main
implicit none
integer::i,omp_get_thread_num,N,prog

! number of threads for parallel computation
call omp_set_num_threads(2)

N=100
prog=0
!$omp parallel do shared(prog) schedule(dynamic,1)
do i=1,N
write(6,*) omp_get_thread_num()
prog=prog+1
if(mod(prog,50).eq.0)write(6,'(A,i6,A,i6)')" ",prog," / ",N
enddo
!$omp end parallel do

stop
end program main

とします。

経験的に、openmpを使用する際はただ一つのサブルーチンを用いて、並列個所を簡単にするのが良いです。

!$omp parallel do
do i=1,N
a(i)=0
do j=i,N
a(i)=a(i)+1
enddo
a(i)=a(i)*i
enddo
!$omp end parallel do

と、!$~!$につらつらと書くのではなく、

!$omp parallel do
do i=1,N
call test(i,N,a(i))
enddo
!$omp end parallel do

として、サブルーチン１つにするだけの方がバグが出にくいです。

openMPの機能はこれだけではありません。
調べることをお勧めします。ここではこれ以上書きません。

参考文献

インテル® コンパイラーOpenMP*入門

プログラミングと数値計算

MATLABを用いた3次元等値面の描き方

2016年2月3日 whitearth コメントする

さて, 数値計算で得られた計算データを3次元で美しく描きたくなると思います.
なると思います.
MATLABで等値面を描くと, とても綺麗になるのでおすすめです.
ということで, 描き方について紹介したいと思います.
まぁ, 手始めに楕円体ですかね.

program main
implicit none

integer(8)::na,nb,nc,intangle
real(8)::x,y,z,V,theta
integer(8),parameter::no=50
real(8)::a,b,c
real(8)::ss(1:3)

Open(50,FILE='data.dat',STATUS='unknown')

a=0.6d0
b=1.2d0
c=0.6d0

do nc=-no,no
write(*,*)nc+no,"/",no*2
do na=-no,no
do nb=-no,no

ss(1)=dble(na)/dble(no)
ss(2)=dble(nb)/dble(no)
ss(3)=dble(nc)/dble(no)

V=(b**2d0 *c**2d0 *SS(1)**2d0)+(c**2d0 *a**2d0 *SS(2)**2d0)+(a**2d0 *b**2d0 *SS(3)**2d0) !楕円体の式

write(50,'(3f10.5,2f11.7)')ss(1),ss(2),ss(3),V
end do
end do
end do

close(50)

end do
end program main

今回は, 101x101x101で計算しますが, これらの範囲は自由に変更できます.
実際に使う際は, vの式を任意の式に変更するだけです.

x, y, zの回す順番だけは, このプログラムコードに合わせて下さい.

さて, data.datができました. 次に, MATLABで読み込めるようにファイルを変換します.

program main
implicit none

integer(8)::nx,ny,nz,xx,yy,zz
real(8),allocatable::x(:,:,:)
real(8),allocatable::y(:,:,:)
real(8),allocatable::z(:,:,:)
real(8),allocatable::v(:,:,:)

nx=50
ny=50
nz=50

open(62,file="data.dat",status="old")!読み込みファイル
open(64,file="x.txt",status="unknown")!出力ファイル
open(65,file="y.txt",status="unknown")!出力ファイル
open(66,file="z.txt",status="unknown")!出力ファイル
open(67,file="v.txt",status="unknown")!出力ファイル

allocate(x(-nx:nx,-ny:ny,-nz:nz))
allocate(y(-nx:nx,-ny:ny,-nz:nz))
allocate(z(-nx:nx,-ny:ny,-nz:nz))
allocate(v(-nx:nx,-ny:ny,-nz:nz))

do zz=-ny,ny
do yy=-nx,nx
do xx=-nz,nz
read(62,*)x(xx,yy,zz),y(xx,yy,zz),z(xx,yy,zz),v(xx,yy,zz)
end do
end do
end do

do xx=-nx,nx
write(64,*) x(xx,:,:)
write(65,*) y(xx,:,:)
write(66,*) z(xx,:,:)
write(67,*) v(xx,:,:)
end do

deallocate(x)
deallocate(y)
deallocate(z)
deallocate(v)

close(62)
close(64)
close(65)
close(66)
close(67)

end program

これで, x.txt, y.txt, z.txt, v.txtが出力されます.
nx, ny, nzの値は, 範囲設定を変えた場合は適宜それに合わせて書き換えて下さい.
さて, 準備ができたので, MATLABを起動しましょう.

x.txt, y.txt, z.txt, v.txtの階層に移動して,
次のコマンドを入力してこれらのファイルを読み込みます.
次のコマンドはスクリプトファイルにしちゃうと楽ですね.

xn=50*2+1 %xのデータ数
yn=50*2+1 %yのデータ数
zn=50*2+1 %zのデータ数
x=load('x.txt')
x2=reshape(x,xn,yn,zn)
y=load('y.txt')
y2=reshape(y,xn,yn,zn)
z=load('z.txt')
z2=reshape(z,xn,yn,zn)
v=load('v.txt')
v2=reshape(v,xn,yn,zn)

コードの内容は,
まず, 先ほどのデータを読み込みます.
次に, 読み込んだデータをちょっと整理してます.

さて, あとは描くだけです！！
次のコマンドを入力していきましょう.

E=0.5 %Eが描きたい等値面の値
p = patch(isosurface(x2,y2,z2,v2,E));
isonormals(x2,y2,z2,v2,p)
set(p,'FaceColor','green','EdgeColor','none'); %ここで色変更可
daspect([1,1,1])
view(30,30); axis([-1.5 1.5 -1.5 1.5 -1.5 1.5]) %axisで表示範囲を設定
camlight ('headlight') %光の当てる方向
camlight ('right') %光の当てる方向
lighting gouraud

これで, 綺麗な楕円体が描けると思います.