わたしは複素変数に対する(第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さんがこの質問をご覧になられた場合、
私のやり方のどこがどう悪いのかご教授くださると有り難いです。
長い質問ですが、是非お願い致します。
A 回答 (2件)
- 最新から表示
- 回答順に表示
No.2
- 回答日時:
絶対値の小さい引数に対しては通常のべき級数展開、絶対値の大きい引数に対しては漸近展開を用いれば高速に計算できるでしょう。
「ニューメリカルレシピ イン シー」という本やその他にベッセル関数を計算するプログラムはあると思います。
No.1
- 回答日時:
高速と言えるかどうかは分かりませんが、前の回答にも書いたように、ベッセル関数の積分表示を実部と虚部に分けて数値積分してしまえばよいのではないでしょうか。
Jn(z)
=((z/2)^n/√πΓ(n+1/2))×∫[0~π]exp(iz cosθ)sin^(2n)θdθ
例えば1次のベッセル関数実部を数値積分の台形公式で計算するRの関数の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
で計算したものと一致していますから、たぶん正しいでしょう。(上のサイトでどのように計算しているかは知りません)
ご解答ありがとうございます。
積分表示は私も試したことがあります。
確かに正確な値が計算できます。
しかし、・・・、積分するのに時間がかかりますね。
ベッセル関数の値が大量に必要なため、
これでは全く使い物にならないのです。
お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!
似たような質問が見つかりました
- C言語・C++・C# C言語 3 2022/10/04 15:07
- その他(プログラミング・Web制作) FORTRAN77の配列(除算) 2 2023/02/01 14:34
- その他(ソフトウェア) F-BASICで計算中の実行が中途で勝手に止まり、大変困っています。 2 2023/03/02 16:15
- その他(コンピューター・テクノロジー) 50台の織機から回転数を取得・集計しモニターに表示したい 2 2022/11/05 15:48
- Visual Basic(VBA) ファイル全てを .xlsm に変更したところ、プログラムが途中で落ちてしまっています 17 2022/12/07 12:03
- C言語・C++・C# numpyスライス機能を使った数値計算 2 2023/05/08 16:01
- C言語・C++・C# このプログラミングの問題を教えてほしいです。 キーボードからデータ数nとn個のデータを入力し、平均値 3 2022/12/19 22:51
- Excel(エクセル) Excel(エクセル)でフィルター抽出後、非表示の行を計算しないで、合計を算出する方法 【内容】 添 4 2023/01/30 17:17
- C言語・C++・C# このプログラミング誰か教えてくれませんか 1 2022/06/02 15:27
- その他(プログラミング・Web制作) VBA 1 2023/01/19 16:19
このQ&Aを見た人はこんなQ&Aも見ています
おすすめ情報
このQ&Aを見た人がよく見るQ&A
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
eのマイナス無限大乗
-
「割る」と「割りかえす」の違い
-
30パーセントオフで371円だった...
-
楕円の円周の長さの計算の仕方...
-
10進法で時間の計算で30分が0.5...
-
公共工事の現場管理費率(%)...
-
分数の計算で分子が0になったら...
-
映画を1.3倍速で見た時の時間計...
-
冪乗の計算について教えてください
-
面積から辺の長さを出す計算式
-
プール計算って何ですか?
-
中学生の数学を習う順番に並べ...
-
半径の計算方法を教えてください。
-
積分のエクセル計算式を教えて...
-
ベランダプールの重量
-
一個当たり15秒の製品を1時間で...
-
2割負担の計算。
-
(x^2-x+1)^10の展開式におけるx...
-
パーセントの算出のしかたを教...
-
無理関数と二次関数の交点
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
eのマイナス無限大乗
-
公共工事の現場管理費率(%)...
-
30パーセントオフで371円だった...
-
分数の計算で分子が0になったら...
-
「割る」と「割りかえす」の違い
-
10進法で時間の計算で30分が0.5...
-
楕円の円周の長さの計算の仕方...
-
面積から辺の長さを出す計算式
-
プール計算って何ですか?
-
積分のエクセル計算式を教えて...
-
1/6n(n+1)(2n+1)+1/2n(n+1) の...
-
(2√2+1)(√2-2)の計算の仕方教...
-
中学生の数学を習う順番に並べ...
-
袋のサイズから容量を計算する方法
-
映画を1.3倍速で見た時の時間計...
-
一個当たり15秒の製品を1時間で...
-
2割負担の計算。
-
2の12乗、32乗・・・とい...
-
エクセルで日数を年数に置き換...
-
3・2+6・3+9・4+.....+3n(n...
おすすめ情報