複素数変数のベッセル関数
わたしは複素変数に対する(第1種)ベッセル関数を計算せねばなりません。
そこでnetlibのhttp://www.netlib.org/amos/から、
D.E. Amos作のFORTRANプログラム一式を、
cbesj.f plus dependencies
をクッリクしてダウンロードしました
(cbesj.fがBessel関数を計算するプログラム、その他はその付属品)。
そして、cbesj.fの冒頭に書いてある使用方法に従い、
わたしは次のようなメインプログラムを書きました。
----------------------------------------------------------
PROGRAM MAIN
IMPLICIT NONE
COMPLEX CY(1), Z
REAL FNU
INTEGER IERR, KODE, N, NZ
Z = (1.2, 0.5)
FNU = 1.0
KODE = 1
N = 1
CALL CBESJ(Z, FNU, KODE, N, CY, NZ, IERR)
WRITE(*,*) IERR
WRITE(*,*) CY(1)
end
----------------------------------------------------------
(このメインプログラムは、次数1、変数(実部1.2,虚部0.5)の(第1種)ベッセル関数の値を計算するためのものです。)
makefileを作成しコンパイルを実行すると、コンパイルは成功しa.outができました。
(※makefileを作らずとも、次のようにしてもよい。
g77 main.f cbesj.f cbinu.f i1mach.f r1mach.f casyi.f cbuni.f cmlri.f cseri.f cuoik.f cwrsk.f xerror.f cuni1.f cuni2.f gamln.f cuchk.f cunhj.f cunik.f cbknu.f crati.f cairy.f ckscl.f cshch.f cacai.f cs1s2.f)
しかし、a.outを実行するとエラーコードIERR=4が返ってきます。
このエラーコードの意味は、cbesj.fの冒頭の説明によると
IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTATION BECAUSE OF COMPLETE LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION
だそうです。でも、何が悪いのか全くわかりません。いろいろメインプログラムを変えてみても、
やはりエラーコードIERR=4が返ってくるだけで、全く受け付けてくれません。
以上に関して、私のやり方のどこがどう悪いのかご教授ください。
また、別のライブラリを用いる方法でもよいので、C言語やFORTRANで(第1種)ベッセル関数
を高速計算する方法をご存知の方は、その方法を詳しくご教授ください。
なお、http://oshiete1.goo.ne.jp/qa2300727.html
において回答者:grothendieckさんが、
やはりAmosのプログラムを用いて複素変数ベッセル関数を計算されているようです。
しかし、そこには詳細なやり方は説明されていないので、計算素人の私には使い方がわかりません。
もしもgrothendieckさんがこの質問をご覧になられた場合、
私のやり方のどこがどう悪いのかご教授くださると有り難いです。
長い質問ですが、是非お願い致します。
お礼
返信遅れまして申し訳ありません。 つまり、ベッセル関数とは、軸対称を用いる場合に使われる変数rを用いた関数で、“p=sin(kx)→p=Jn(Kr)”と変換したものということですね。 なんとかこの関数に対してのイメージが固まった気がします。 ご回答、ありがとうございました。