複素数変数のベッセル関数

わたしは複素変数に対する(第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さんがこの質問をご覧になられた場合、
私のやり方のどこがどう悪いのかご教授くださると有り難いです。


長い質問ですが、是非お願い致します。

No.2 回答者： grothendieck 回答日時： 2008/11/20 21:59

絶対値の小さい引数に対しては通常のべき級数展開、絶対値の大きい引数に対しては漸近展開を用いれば高速に計算できるでしょう。

「ニューメリカルレシピ　イン　シー」という本やその他にベッセル関数を計算するプログラムはあると思います。

No.1 回答者： grothendieck 回答日時： 2008/11/18 00:33

高速と言えるかどうかは分かりませんが、前の回答にも書いたように、ベッセル関数の積分表示を実部と虚部に分けて数値積分してしまえばよいのではないでしょうか。

Jn(z)
=((z/2)^n/√πΓ(n+1/2))×∫[0～π]exp(iz cosθ)sin^(2n)θdθ
例えば１次のベッセル関数実部を数値積分の台形公式で計算するＲの関数の1例は

complexbessel <- function(x,y)
{
h <- 0.00001*3.1415926358
s <- 0
a <- 0
for(i in 1:100000){
s<- s+x*(exp(-y*cos(a))*cos(x*cos(a))*(sin(a))^2
+ exp(-y*cos(a+h))*cos(x*cos(a+h))*(sin(a+h))^2 )/2
s<- s-y*(exp(-y*cos(a))*sin(x*cos(a))*(sin(a))^2
+ exp(-y*cos(a+h))*sin(x*cos(a+h))*(sin(a+h))^2 )/2

a <- h*i }
s*0.00001
}
これでJ1(1.7+0.5i)の実部は0.6301477となりました。
http://keisan.casio.jp/has10/SpecExec.cgi
で計算したものと一致していますから、たぶん正しいでしょう。（上のサイトでどのように計算しているかは知りません）

この回答へのお礼

ご解答ありがとうございます。
積分表示は私も試したことがあります。
確かに正確な値が計算できます。

しかし、・・・、積分するのに時間がかかりますね。
ベッセル関数の値が大量に必要なため、
これでは全く使い物にならないのです。

お礼日時：2008/11/18 17:56