2012年5月4日金曜日

ベッセル関数 (J0(x) と J1(x)) の数表 (例えば x は 0~10 の範囲で 0.1 刻み) を...

ベッセル関数 (J0(x) と J1(x)) の数表 (例えば x は 0~10 の範囲で 0.1 刻み) を出力する C のプログラムを作れ。という課題なのですができません教えて下さい







> itagi_juniorさん



知恵袋過去を探したら

detail.chiebukuro.yahoo.co.jp/qa/question_detail/q1112026261

がありました。





あとGSLと言うライブラリにもベッセル関数含まれてます。

それのソース見れば、答になってるはず。



http://ftp.jaist.ac.jp/pub/GNU/gsl/gsl-1.15.tar.gz



展開してそれらしい名前のファイルが以下

gsl-1.15/specfunc/bessel_J0.c

gsl-1.15/specfunc/bessel_J1.c



答えあわせも使えると思います。



よろしくお願い致します。







****追記1

あと、何重にも同じ投稿されているので、他消していただけないでしょうか。

よろしくお願い致します。





*** 補足分



私も良く見てませんが、テイラー展開(かな?)した項を計算していって、項の値がある一定以下になるまで計算してるってパターンじゃないのでしょうか?



後で過去のプログラムとGS見ときます。何かわかったら追記します。



よろしくお願い致します。





*** 追記2

上のプログラム後半削って、最初の方にあるループの範囲を変更したらJ0だと思う。

数値未確認、GSLと比較したら検算になると思います。

後は以下を頑張ってください。



1)検算



2)検算合ってたら、下のプログラムを改造してJ1を計算する



よろしくお願い致します。



--- em163.c

#include <stdio.h>

#include <math.h>



double Bess0( // J0(x) を計算し、戻り値で返す

double x, //[i] 独立変数

double eps, //[i] 打ち切り誤差

int max, //[i] 級数項数の上限

int* lp) //[o] 計算した項数

{

int k;

double fct, dJ, J0;



k = 0; J0 = 1;

dJ = 1; fct = -pow(x/2,2.);

do {

if ( ++k > max ) break;

dJ *= fct/(k*k);

J0 += dJ;

} while ( fabs(dJ) > eps );

*lp = k;

return J0;

}

int main(){

int loop, maxlp=1000;

double x, J0, eps=1.e-12;



printf("1.0<x<10.0 でのJ0(x)の分布\n");

printf(" x N J0\n");

for(x=1.;x<=10.;x+=0.1) {

J0 = Bess0( x, eps, maxlp, &loop );

printf("%f %3d %9f\n",x,loop,J0);

}

printf("\n");

return 0;

}

---

0 件のコメント:

コメントを投稿