Fortranでベッセル関数のべき級数展開が正しく計算できない!!!
タイトルの通りです。ベッセル関数のべき級数展開は画像に添付した通りです。
この計算をFortranで行っているのですが、式中のzが30以上になると、正しく
計算できません。zが30以下のときは、ベッセル関数の値と完全に一致しています。
どうも、(z/2)^2nの項が大きくなりすぎることが原因です。
これの解決法をどなたかご存知ですか?
もしご存知でしたらご教授願います。
ちなみに∞はn=100までで計算をしています。
ほんと死にそうです。。。助けて下さい。
No.8
- 回答日時:
ん~, 少なくとも Horner の方法では桁落ちしますね. まあ, 交項級数なのでその可能性は想定できたわけですが.... 他の方法でも, 部分和の大きさに対して最終的な結果は非常に小さいので, やっぱり落ちてると思います.
ベッセル関数を計算するルーチンには公開されているものもありますが, 大きな z に対しては漸近展開に切り替えてますね. Numerical Recipes だと |z| > 8 で切り替えてます.
No.5
- 回答日時:
>(1)1/N!の計算結果の呼び出し(別のループで計算済み)
>(2)(A/2)^2Nの計算
これではだめです。
(A/2)^2N/N!(N=1,2,...)
を計算して足し合わせるようにしてください。大きな数と小さな数を作らないようにするのがコツです。
漸近式の計算結果と比較してみましたか
No.3
- 回答日時:
まず「正しく計算できない」とだけ述べるのではなく, 「期待する結果」と「実際に得られた結果」を書くようにしてください. 例えば「z の値がこれこれのときに本当はこうならなきゃいけないんだけど得られた結果はこうである」と書いてください.
また, 数値計算を行うにあたっては「どのようなプログラムで計算したのか」は重要な情報となります. ですので, できる限りあなたが用いたプログラムを出すようにしてください.
ということで本題だけど, 交項級数なので #2 でも言われるようにどうしても精度は悪くなりがちです. できれば努力と根性で交項級数じゃないようにした方がいい. あるいは, n の上限を 100 にしているのだから単なる多項式とみなせるので, Horner の方法なんかを使ってもいいかもしれない.
この回答への補足
ご指摘有難うございます。
プログラムですが、簡単に書きますと、
*******************************
DO 1 Z=1,40,1 (ベッセル級数のzの値)
SUM=0
DO 2 N=0,NN,1 (無限級数の重ね合わせ NN:無限級数の上限(100))
(1)1/N!の計算結果の呼び出し(別のループで計算済み)
(2)(A/2)^2Nの計算
(3)(-1)^Nの計算
DSUM=(1)*(2)*(3)
SUM=SUM+DSUM
2 CONTINUE
結果出力(Z=○のときにJ0(A)は▲)
1 CONTINUE
********************************
そして、結果ですが、
A__J0(A)の真値__級数による解
1__0.76519____0.76519
・
32__0.13807____0.1380465
40__7.37E-03___-0.398868
42__-1.15E-01___3.65E+00
44__8.63E-02___ -25.0222
46__3.94E-02___ 294.5768
・
50__5.58E-02___ -1.03E+04
・
60__-9.15E-02___2.40E+08
のようになります。
1~31程度 ほぼ一致
32~40程度 僅かにずれ
それ以降 爆発的に誤差が大きくなる
これで再度検討していただけたらありがたいです。
No.1
- 回答日時:
>どうも、(z/2)^2nの項が大きくなりすぎることが原因です。
数値が大きくなりすぎることに対処するには、Σの中身の対数をとるなどして、数値を小さくして計算すると良いと思います。
L(n)=log|(-1)^n/n!^2 (z/2)^(2n)} とします。(logは常用対数とします。)
すると、L(n)は次のように表されます。
L(n)= log{(z/2)^(2n)/n!^2} =2n log(z/2)-2 log(n!)
これで計算すれば、数値が扱える範囲を超えることは防げると思います。
あとは、J_0(n)はつぎのように足し併せて計算できます。
J_0(n) =Σ[n=0→∞] (-1)^n 10^L(n)
この回答への補足
直ちに回答して頂き有難うございます。
さて、対数の件ですが、既に試したのですがやはりzが同じぐらいの値のときに
計算結果がおかしくなってしまいました。
底を大きくしたりも試したのですが、うまくいきません。
他に考えられることってなんかありますか?
現在は(z/2)^2nをどうにか簡単にできないかと考えてます。
(1/n!)^2を(1/n!)(1/n!)と分離して(z/2)^n(z/2)^nと分離したものと
打ち消しても同様の結果でした。orz...
お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!
似たような質問が見つかりました
- Excel(エクセル) エクセルのSUM関数について 4 2023/04/18 10:37
- 数学 「f(z)=1/(z^2-1)に関して ローラン展開を使う場合、マクローリン展開を使う場合、テイラー 3 2022/08/27 19:56
- モニター・ディスプレイ LOW/HIGH-PASS FILTER passive 8 2022/04/22 11:17
- Excel(エクセル) エクセル 関数について質問です。 2 2022/10/03 11:14
- 工学 以前、線形代数からフーリエ級数展開を導く上で 式v=(v, e1)e1+(v, e2)e2+…+(v 6 2022/06/29 17:24
- 数学 『数は実在するのか』 6 2023/06/04 15:15
- Excel(エクセル) エクセル 自動計算 1 2023/01/30 13:28
- その他(ソフトウェア) F-BASICで計算中の実行が中途で勝手に止まり、大変困っています。 2 2023/03/02 16:15
- Excel(エクセル) エクセルの関数で、間違っている物はないのでしょうか 4 2022/04/14 00:44
- 数学 数2 分数の計算について 写真の問題はいちいち展開したあとに計算しないといけないのでしょうか 例えば 6 2023/04/29 16:37
おすすめ情報
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
らせんRの計算の仕方
-
「再帰的」の意味を教えてください
-
数学の問題です。
-
中3の有効数字の範囲の問題で √...
-
10^0.2 = 1.58489319246111の計...
-
いつもありがとうございます、 ...
-
信頼区間90%は何σ?
-
電力ケーブルのインピーダンス...
-
1512の1/5乗
-
倍率とデシベルの計算式
-
ガウスの消去法
-
2の128乗の計算方法
-
円周率はなぜ3.14で計算するよ...
-
算数の時計算で「1時と2時の間...
-
3割の計算
-
至急【高校数学A n進法の四則...
-
時間を100進法であらわしたい。
-
10,000百万円っていくらですか?
-
「強度」は高い?強い?
-
穴が開く? 空く? 明く?
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
3割の計算
-
信頼区間90%は何σ?
-
減少率計算式教えて下さい
-
倍率とデシベルの計算式
-
「再帰的」の意味を教えてください
-
らせんRの計算の仕方
-
電力ケーブルのインピーダンス...
-
2の128乗の計算方法
-
常用対数の求め方 log10の2は約...
-
1.01の12乗の計算
-
電卓の機能の名称が分からない...
-
10^0.2 = 1.58489319246111の計...
-
実生活で役に立つ数学ってあり...
-
指数計算を教えてください
-
2次関数って何の仕事で必要な...
-
計算の方法を教えてください。 ...
-
円周率の計算式って何ですか?
-
100の101乗と101の100乗ではど...
-
1512の1/5乗
-
三角形の面積・・・ヘロンと座...
おすすめ情報