学校の課題でベッセル(Bessel）関数 0次をＣ言語でプログラミングして、ｘ→∞ にし...

学校の課題でベッセル(Bessel）関数 0次をＣ言語でプログラミングして、ｘ→∞　にしたときのV0(x)の値を求めろと課題を出されました。

#include <stdio.h>

#include <math.h>



main()

{

double k=1,a2,a3,f,b;

int kurikaeshi=10; //kurikaeshiはk→∞の範囲

double k1[90],c1[90],s,e,a1[90],b1[90],l1[90],f1[90],x,end[90],j;

int i,a,a4;





k1[0] = 1;

a1[0] = 1;

c1[0] = 1;

for(x=2.00;x<kurikaeshi;x+0.01){



//階乗する分母の式----------------------------------------------------------



for(i=0;i<kurikaeshi;i++){

if(i == 0) k1[i] = 1;

else{

k=k*i;

k1[i] = k * k;}

// printf("k! * k! = %e \n",k1[i]);





//----分子の計算------------------------------------------------------------

if(i == 0) a1[i] = 1,f1[i] = 1;

else{

a = pow(-1,i);

a1[i] = a;

f1[i] = a1[i] / k1[i];

}

// printf("k=%d (-1)^k / k!*k! = %e \n",i,f1[i]);



//----(x/2)^2k--------------------------------------------------------------

if(x == 0) a3 = 1.0;

else{ if(i == 0 ) a3 = 1.0;

else{

a2 = x / 2.0;

a4 = 2 * i;

a3 = pow(a2,a4);

}}

c1[i] = a3;

// printf("(x/2)^2k = %e \n",c1[i]);





//----最後手前の計算------------------------------------------------------------

if(i == 0) l1[i] = 1;

else{

b = f1[i] * c1[i];

l1[i] = b + l1[i-1];

}

printf("x=%e k=%d %e \n\n",x,i,l1[i]);



j = l1[i];

while(fabs(j) < fabs(j) * 1e-25) {

end[i] = j;

printf("☆ x=0 y=%e \n",end[i]);

}

}





}

return 0;

}



と記述してみたのですが上手く計算されず途中でエラーとなってしまいます。

どうか解決する方法を教えてはもらえませんか？

お願いします。



参考ＵＲＬはhttp://www8.plala.or.jp/randomism/besselj/　です。

0次なのでn=0 で計算しています。

//_ベッセル関数：



//_試作品です：（級数計算を改善）

/*******************************

J0(x)は自前の関数Bess0で求めるようにしました。

テストメインでは、

① J0(x)の零点付近の様子を確かめ

② 適切な区間を指定して二分探索法で求めます

最後に、printfでの出力がついています。

*******************************/

// 二分探索法について、３行の簡潔な手順説明がここにあります↓

// http://epp.eps.nagoya-u.ac.jp/~sirono/suti02/lect1203/lect1203.html

// 他の話題もありますが、大きく二分探索法と出ていますので見つけ易いです。



#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;

}

void_main(){

____int_____iter,_loop,_maxlp=1000;

____double__x,_J0,_eps=1.e-12;

____double__xa,_xb,_Ja,_Jb;



____printf("1.0<x<3.0_でのJ0(x)の分布\n");

____printf("_x_________N___J0\n");

____for(x=1.;x<=3.;x+=0.5)_{

________J0_=_Bess0(_x,_eps,_maxlp,_&loop_);

________printf("%f_%3d_%9f\n",x,loop,J0);

____}

____printf("\n");

____printf("1.0<x<3.0_でのJ0(x)=0の根（J0(x)の零点）\n");

____xa_=_1.0;

____xb_=_3.0;

____Ja_=_Bess0(_xa,_eps,_maxlp,_&loop_);

____Jb_=_Bess0(_xb,_eps,_maxlp,_&loop_);

____iter_=_1;

____do_{____//_Binary_search_for_J0(x)_=_0.

________x__=_0.5*(xa+xb);

________J0_=_Bess0(_x,_eps,_maxlp,_&loop_);

________printf("trial_%2d:_at_x=%f,_J0=%9.1e\n",iter,x,J0);

________if_(_fabs(J0)_<_eps_)_break;

________if_(_J0*Ja_>=_0_)_xa_=_x;

________else______________xb_=_x;

____}_while_(_++iter_<_maxlp_);

____printf("\n");

____if(_iter_<_maxlp_)_{

________printf("Converged_at_%dth_trial_of_half_way_search.\n",iter);

________printf("J0(x)_=_%15.12f,_at_x_=_%15.12f\n",J0,x);

____}

____return;

}



実行結果：



1.0<x<3.0_でのJ0(x)の分布

_x_________N___J0

1.000000___8__0.765198

1.500000___9__0.511828

2.000000__10__0.223891 <=== ここまでプラス

2.500000__11_-0.048384 <=== ここからマイナス

3.000000__12_-0.260052



２分探索法での経過： <== 逐次近似の経過は普通出力しない、出すのは勝手

1.0<x<3.0_でのJ0(x)=0の根（J0(x)の零点）

trial__1:_at_x=2.000000,_J0=_2.2e-001

trial__2:_at_x=2.500000,_J0=-4.8e-002

trial__3:_at_x=2.250000,_J0=_8.3e-002

trial__4:_at_x=2.375000,_J0=_1.6e-002

trial__5:_at_x=2.437500,_J0=-1.7e-002

trial__6:_at_x=2.406250,_J0=-7.4e-004

trial__7:_at_x=2.390625,_J0=_7.4e-003

trial__8:_at_x=2.398438,_J0=_3.3e-003

trial__9:_at_x=2.402344,_J0=_1.3e-003

trial_10:_at_x=2.404297,_J0=_2.7e-004

trial_11:_at_x=2.405273,_J0=-2.3e-004

trial_12:_at_x=2.404785,_J0=_2.1e-005

trial_13:_at_x=2.405029,_J0=-1.1e-004

trial_14:_at_x=2.404907,_J0=-4.2e-005

trial_15:_at_x=2.404846,_J0=-1.1e-005

trial_16:_at_x=2.404816,_J0=_5.1e-006

trial_17:_at_x=2.404831,_J0=-2.8e-006

trial_18:_at_x=2.404823,_J0=_1.2e-006

trial_19:_at_x=2.404827,_J0=-8.1e-007

trial_20:_at_x=2.404825,_J0=_1.8e-007

trial_21:_at_x=2.404826,_J0=-3.1e-007

trial_22:_at_x=2.404826,_J0=-6.7e-008

trial_23:_at_x=2.404825,_J0=_5.6e-008

trial_24:_at_x=2.404826,_J0=-5.5e-009

trial_25:_at_x=2.404826,_J0=_2.5e-008

trial_26:_at_x=2.404826,_J0=_1.0e-008

trial_27:_at_x=2.404826,_J0=_2.3e-009

trial_28:_at_x=2.404826,_J0=-1.6e-009

trial_29:_at_x=2.404826,_J0=_3.5e-010

trial_30:_at_x=2.404826,_J0=-6.2e-010

trial_31:_at_x=2.404826,_J0=-1.3e-010

trial_32:_at_x=2.404826,_J0=_1.1e-010

trial_33:_at_x=2.404826,_J0=-1.4e-011

trial_34:_at_x=2.404826,_J0=_4.7e-011

trial_35:_at_x=2.404826,_J0=_1.7e-011

trial_36:_at_x=2.404826,_J0=_1.6e-012

trial_37:_at_x=2.404826,_J0=-5.9e-012

trial_38:_at_x=2.404826,_J0=-2.2e-012

trial_39:_at_x=2.404826,_J0=-2.8e-013



Converged_at_39th_trial_of_half_way_search.

J0(x)_=_-0.000000000000,_at_x_=__2.404825557696 <=== J0(x)はほぼゼロ