第1種変形ベッセル関数

第1種変形ベッセル関数をFortranで書くとどのようなプログラムになるのでしょうか。
Fortran初心者です。よろしくお願いします。

No.5 回答者： stomachman 回答日時： 2013/09/20 23:01

　デキアイのライブラリをテストしてみると、用途に照らして精度が不足だったり余りに遅かったりということも、案外あるものです。

また、実行環境によっては、ライブラリを突っ込むスペースすらもない、ということだってあり、自作するのは必ずしもおかしなことではありません。

> Fortran初心者です

　いやいや、Fortranばかりでなく、数値計算の超初心者ですよね。

(0) 何のために第1種変形ベッセル関数を使うのか。ホントに必要なのか。

　「何のために計算するんだかきちんと考え直してみたら、そもそもそんな計算は要らなかった」だとか「もっと簡単な関数で充分代用（近似）できる」ということがしばしばあります。また時には、単に式を整理したらベッセル関数が消えちゃった、なんてマヌケな話もあります。

　ホントに要るんだという場合でも「高性能のハードウェア上でごくたまーに使う」というのと、「200円のマイコン上で毎秒何千回も」と言うのでは、考え方が全く違ってきます。

(1) 第1種変形ベッセル関数を、その都度計算する必要があるのか。

　たとえば、特定のαとxにおけるIα(x)を1回だけ使う、というのだったら、そもそもプログラムを書く必要がないでしょ？いや、そこまで極端でなくても、大抵の応用では、あらかじめ数表を作っておけば済みます。次に

(2) 第1種変形ベッセル関数をどんな範囲のαとxで使い、どれだけの精度が必要なのか。

　（「Iα(x)を計算して印刷しなさい」なんて練習問題は論外としますと、）その結果を使って何か意味のある計算をやる訳ですから、α、xの範囲と、関数の値に求められる精度が分かる筈です。数表にする場合にも、この情報は必須ですね。もちろん、αとxの範囲が狭く、さして精度が要らないのなら、小さな表を作っておけば済みます。
　αとxの範囲と必要な精度を見極めて、ようやく

(3) どんな計算式を使って第1種変形ベッセル関数を計算するか。その計算式を如何にして計算するか。

の検討ができる。
　たとえば「第1種変形ベッセル関数とは言っても、αとxがこの範囲で必要な精度がこれだけであるなら、双3次関数による近似だけで充分だ」なんてことが判明し、ほんの数行のコードで済んでしまう、なんてことも、しばしば（かなりしばしば）あります。
　そうも行かないという時には、大抵は繰り返し(loop)の計算になります。すると、単に計算式どおりに計算していく、というのでは良い性能（スピードと精度）が得られません。「計算式を如何にして計算するか」が問題です。たとえば、近似値を効率よく改良できる漸化式を考案したり、数値計算に伴う誤差を抑え込む工夫をしたり、無駄な計算をしないための打ち切りの判定条件を案出したりして、計算式からアルゴリズムを創り出す。（実は、ここが数値計算プログラミングの腕の見せ所であり、同じ精度を出すのにヘタと上手とでは計算時間が10000倍違う、なんてこともあります。）
　この段階では、どんなやり方がいいか、実際にテストプログラムを書いたり、表計算ソフトを使ったりして、結果を見て検討することになるでしょう。

(4) 実装方法を決め、コーディングする。

　時には、プログラミング言語に適した実装方法、つまり「書き方」を工夫しなくてはならないことがあります。Fortranは、版によってはイマドキのプログラミング言語なら当たり前に持っている機能がなかったりします。(なにしろ元は1955年製ですから。）

　ここまできてようやく「Fortranで書くとどのようなプログラムになる」かが決まる。だから、ご質問を見ただけでシロートだなとバレちゃうんです。

(5) テストベッド（FUNCTIONを呼び出すメインプログラム）を製作し、テストする。

　少なくとも2回、有効桁数(precision)を変えて計算してみて、両者の結果を比較します。必要な精度までどちらのテスト結果も一致していればいいんですが、さもなければ、有効桁数の多い方の計算ですら充分な精度が出ていないおそれがあります。もちろん、いくつかの(α,x)について、「正解」を用意する必要もありますね。

　というわけで、ご自身の質問の意味がお分かりになったところで、(0)から順に始めましょう。

この回答へのお礼

ご指摘とおりです。数値計算についても初心者です。今一度、自分の中で整理したいと考えております。

お礼日時：2013/09/23 19:08

No.4 ベストアンサー 回答者： spring135 回答日時： 2013/09/20 22:10

＃3です。

訂正です。


DIMENSION FAC(100)

DO 10 K=1,100,1

I=1

DO 20 J=1,K,1

I=I*J

20 CONTINUE

FAC(K)=I

10　CONTINUE

S=0

L=10

DO 100 I=0,L,1

DS=(X/2)**(2*I)/FAC(I)/FAC(N+I)

S=S+DS

100 CONTINUE

EIN=S*(X/2)**N

RETURN

END

Lを100より大きくするときはFACをそれに合わせて作ってください。

L=170ぐらいまではこのままでできます。

この回答へのお礼

ありがとうございました。早速試させていただきます・

お礼日時：2013/09/25 18:36

No.3 回答者： spring135 回答日時： 2013/09/20 22:03

整数(n)次の第１種変形ベッセル関数EI(n,x)は

EI(n,x)=(x/2)^n*[lim(j→∞)Σ(i=0,j)(x/2)^(2i)/i!/(n+i)!

で表されます。

FORTRANは最近まったく使わないので間違っているかもしれませんが一応書いてみます。

無限級数ですが収束するということは有限項でも間に合うということでL項まで計算してみて

Lを増減してあたりをみます。

階乗を最初に作っておきます。

DIMENSION FAC(100)

DO 10 K=1,100,1

I=1

DO 20 J=1,K,1

I=I*J

20 CONTINUE

FAC(K)=I

10　CONTINUE

S=0

L=10

DO 100 I=0,L,1

DS=(X/2)**(2*I)/I!/(N+I)!

S=S+DS

100 CONTINUE

EIN=S*(X/2)**N

RETURN

END

あちこち訂正しながらやってみてください。

変形でないベッセルを散々扱かったので10項ぐらいで３ケタは合うでしょう。

100項やれば５ケタぐらい合います。

No.2 回答者： Tacosan 回答日時： 2013/09/20 10:58

「ライブラリ等は使わないで組めれば」と考える理由が分からない. ライブラリ関数を使えるなら, その方がいいに決まってる.

わざわざ時間と手間をかけてバグの入る可能性を増やすことなど愚の骨頂.

No.1 回答者： spring135 回答日時： 2013/09/20 02:20

ライブラリーにあると思いますが、級数展開で数値会を求めたいということですか。

この回答への補足

ありがとうございます。ライブラリ等は使わないで組めればと考えております。よろしくお願いします

補足日時：2013/09/20 07:03