Fortranでベッセル関数のべき級数展開が正しく計算できない！！！

Fortranでベッセル関数のべき級数展開が正しく計算できない！！！


タイトルの通りです。ベッセル関数のべき級数展開は画像に添付した通りです。

この計算をFortranで行っているのですが、式中のzが30以上になると、正しく

計算できません。ｚが30以下のときは、ベッセル関数の値と完全に一致しています。

どうも、(z/2)^2nの項が大きくなりすぎることが原因です。

これの解決法をどなたかご存知ですか？

もしご存知でしたらご教授願います。

ちなみに∞はn=100までで計算をしています。

ほんと死にそうです。。。助けて下さい。

No.1 回答者： Mr_Holland 回答日時： 2010/05/17 16:04

＞どうも、(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)

この回答への補足

直ちに回答して頂き有難うございます。

さて、対数の件ですが、既に試したのですがやはりｚが同じぐらいの値のときに

計算結果がおかしくなってしまいました。

底を大きくしたりも試したのですが、うまくいきません。

他に考えられることってなんかありますか？

現在は（ｚ/2）^2nをどうにか簡単にできないかと考えてます。

（１/n!）^2を（１/n!）（１/n!）と分離して（ｚ/2）^n（ｚ/2）^nと分離したものと

打ち消しても同様の結果でした。orz...

補足日時：2010/05/17 16:30

No.2 回答者： -ok 回答日時： 2010/05/17 17:28

よくわかりませんが、ひょっとして桁落ちしていませんか？

正の項と負の項があるので、気になります。

また、収束性が悪くなっているということはありませんか？

この回答への補足

回答ありがとうございます。

桁落ちはしていることはしてるのですが、それが原因なのかもわからない状況です。

補足日時：2010/05/18 10:48

No.3 回答者： Tacosan 回答日時： 2010/05/17 18:17

まず「正しく計算できない」とだけ述べるのではなく, 「期待する結果」と「実際に得られた結果」を書くようにしてください. 例えば「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程度　僅かにずれ
それ以降　　爆発的に誤差が大きくなる


これで再度検討していただけたらありがたいです。

補足日時：2010/05/17 19:27

No.4 回答者： -ok 回答日時： 2010/05/17 20:13

＃２です。

桁落ちが原因であるかどうかをみるために、各Nでの項の値を調べることをお勧めします。

No.5 回答者： spring135 回答日時： 2010/05/17 21:02

＞(1)1/N!の計算結果の呼び出し（別のループで計算済み）

＞(2)(A/2)^2Nの計算

これではだめです。
(A/2)^2N/N!(N=1,2,...)
を計算して足し合わせるようにしてください。大きな数と小さな数を作らないようにするのがコツです。

漸近式の計算結果と比較してみましたか

No.6 回答者： spring135 回答日時： 2010/05/17 21:20

NO.5です。

もっと丁寧に書きます。
Sn=Π(k=1→n)(z/2k)^2
J0=Σ(n=1→∞)(-1)^n*Sn
とします。

この回答への補足

回答有難うございます。

さっそく試してみましたが、ｚが３８ぐらいまで計算できるようになりました。

目から鱗の方法でした。

しかし、まだまだ大きなｚについて出来るようになりたいと考えております。

これはもう諦めた方が賢明ですかね？

漸近展開については比較はしておりません。

また試してみたいと思います。

補足日時：2010/05/18 11:13

No.7 ベストアンサー 回答者： Tacosan 回答日時： 2010/05/18 00:20

「これでいけば簡単」という方針は #3 に書いておいたんだけどなぁ....

よく読めば「プログラムを要求している」だけじゃないと気付く, はず.

この回答への補足

コメント足らずで申し訳ございません。

Horner の方法などはいろいろ不都合がありまして元々考えてはいませんでした。

べき級数展開をするしかない状況で計算が出来ないので非常に困っております。

補足日時：2010/05/18 10:51

この回答へのお礼

honerをためしましたが、結局桁落ちしているみたいです。

お礼日時：2010/05/18 12:52

No.8 回答者： Tacosan 回答日時： 2010/05/19 01:02

ん～, 少なくとも Horner の方法では桁落ちしますね. まあ, 交項級数なのでその可能性は想定できたわけですが.... 他の方法でも, 部分和の大きさに対して最終的な結果は非常に小さいので, やっぱり落ちてると思います.

ベッセル関数を計算するルーチンには公開されているものもありますが, 大きな z に対しては漸近展開に切り替えてますね. Numerical Recipes だと |z| > 8 で切り替えてます.