2007年8月16日木曜日

ベッセル関数

いまやってる計算でベッセル関数を計算しないといけない。
いろいろあさってると、Numerical recipes in Fortran77のサンプルコードが見れるサイトがあった。
この本、洋書しかないし、買うと一万ちかいのでこういうサイトは大変ありがたい。
ざっと読んだところ、ベッセル関数の引数の大きさによって近似多項式を変えている模様。
cosとかsinとかで展開する方法を使ってるプログラムは結構見かけるが、どうやらこの三角関数を使った展開はある程度引数が大きくないと不正確ならしいので、場合わけしたほうがよいみたい。
載ってたサンプルコードをちょいちょいいじって、0次のベッセル関数をプロットするようなのを作ってみたのでメモっとく。また使うかも知れんし。それにしても なんちゅーかひねりのないプログラムやなぁ。計算できりゃえーやって感じでいつも作ってるから、美しさからは程遠いものしか作れない。fortranはグラフィックの関数を持ってないのでそのまんまグラフ書くとかは無理なのもいけてないよな。

implicit none
REAL*8 bessj0,x
integer i,N
c Returns the Bessel function J0(x) for any real x.
REAL*8 ax,xx,z
DOUBLE PRECISION p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,
* r5,r6,s1,s2,s3,s4,s5,s6,y !We’ll accumulate polynomials in double precision.
SAVE p1,p2,p3,p4,p5,q1,q2,q3,q4,q5,r1,r2,r3,r4,r5,r6,
* s1,s2,s3,s4,s5,s6
DATA p1,p2,p3,p4,p5/1.d0,-.1098628627d-2,.2734510407d-4,
* -.2073370639d-5,.2093887211d-6/, q1,q2,q3,q4,q5/-.1562499995d-1,
* .1430488765d-3,-.6911147651d-5,.7621095161d-6,-.934945152d-7/
DATA r1,r2,r3,r4,r5,r6/57568490574.d0,-13362590354.d0,
& 651619640.7d0,
* -11214424.18d0,77392.33017d0,-184.9052456d0/,
* s1,s2,s3,s4,s5,s6/57568490411.d0,1029532985.d0,
* 9494680.718d0,59272.64853d0,267.8532712d0,1.d0/
N=10
do i=1,100
x=0.0d0+1.0d0/dble(N)*dble(i)
if(dabs(x).lt.8.)then !Direct rational function fit.
y=x**2
bessj0=(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6)))))
* /(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))
else !Fitting function (6.5.9).
ax=dabs(x)
z=8./ax
y=z**2
xx=ax-.785398164
bessj0=sqrt(.636619772/ax)*(cos(xx)*(p1+y*(p2+y*(p3+y*(p4+y
* *p5))))-z*sin(xx)*(q1+y*(q2+y*(q3+y*(q4+y*q5)))))
endif
write(*,*) x,bessj0
end do

END

0 コメント: