第二类整数阶贝塞尔函数(诺伊曼函数)

#include "stdio.h" #include "math.h" /************************************************************ 第一类贝塞尔函数:贝塞尔函数 *************************************************************/ double Bessel(int n,double x) { int i,m; double t,y,z,p,q,s,b0,b1; static double a[6]={ 57568490574.0,-13362590354.0, 651619640.7,-11214424.18,77392.33017, -184.9052456}; static double b[6]={ 57568490411.0,1029532985.0, 9494680.718,59272.64853,267.8532712,1.0}; static double c[6]={ 72362614232.0,-7895059235.0, 242396853.1,-2972611.439,15704.4826, -30.16036606}; static double d[6]={ 144725228443.0,2300535178.0, 18583304.74,99447.43394,376.9991397,1.0}; static double e[5]={ 1.0,-0.1098628627e-02, 0.2734510407e-04,-0.2073370639e-05, 0.2093887211e-06}; static double f[5]={ -0.1562499995e-01, 0.1430488765e-03,-0.6911147651e-05, 0.7621095161e-06,-0.934935152e-07}; static double g[5]={ 1.0,0.183105e-02, -0.3516396496e-04,0.2457520174e-05, -0.240337019e-06}; static double h[5]={ 0.4687499995e-01, -0.2002690873e-03,0.8449199096e-05, -0.88228987e-06,0.105787412e-06}; t=fabs(x); if (n<0) n=-n; if (n!=1) { if (t<8.0) { y=t*t; p=a[5]; q=b[5]; for (i=4; i>=0; i--) { p=p*y+a[i]; q=q*y+b[i]; } p=p/q; } else { z=8.0/t; y=z*z; p=e[4]; q=f[4]; for (i=3; i>=0; i--) { p=p*y+e[i]; q=q*y+f[i]; } s=t-0.785398164; p=p*cos(s)-z*q*sin(s); p=p*sqrt(0.636619772/t); } } if (n==0) return(p); b0=p; if (t<8.0) { y=t*t; p=c[5]; q=d[5]; for (i=4; i>=0; i--) { p=p*y+c[i]; q=q*y+d[i];} p=x*p/q; } else { z=8.0/t; y=z*z; p=g[4]; q=h[4]; for (i=3; i>=0; i--) { p=p*y+g[i]; q=q*y+h[i]; } s=t-2.356194491; p=p*cos(s)-z*q*sin(s); p=p*x*sqrt(0.636619772/t)/t; } if (n==1) return(p); b1=p; if (x==0.0) return(0.0); s=2.0/t; if (t>1.0*n) { if (x<0.0) b1=-b1; for (i=1; i<=n-1; i++) { p=s*i*b1-b0; b0=b1; b1=p;} } else { m=(n+(int)sqrt(40.0*n))/2; m=2*m; p=0.0; q=0.0; b0=1.0; b1=0.0; for (i=m-1; i>=0; i--) { t=s*(i+1)*b0-b1; b1=b0; b0=t; if (fabs(b0)>1.0e+10) { b0=b0*1.0e-10; b1=b1*1.0e-10; p=p*1.0e-10; q=q*1.0e-10; } if ((i+2)%2==0) q=q+b0; if ((i+1)==n) p=b1; } q=2.0*q-b0; p=p/q; } if ((x<0.0)&&(n%2==1)) p=-p; return(p); } /************************************************************* 第二类贝塞尔函数:诺伊曼函数 *************************************************************/ double Neumann(int n,double x) { int i; double y,z,p,q,s,b0,b1; static double a[6]={ -2.957821389e+9,7.062834065e+9, -5.123598036e+8,1.087988129e+7,-8.632792757e+4, 2.284622733e+2}; static double b[6]={ 4.0076544269e+10,7.452499648e+8, 7.189466438e+6,4.74472647e+4,2.261030244e+2,1.0}; static double c[6]={ -4.900604943e+12,1.27527439e+12, -5.153438139e+10,7.349264551e+8,-4.237922726e+6, 8.511937935e+3}; static double d[7]={ 2.49958057e+13,4.244419664e+11, 3.733650367e+9,2.245904002e+7,1.02042605e+5, 3.549632885e+2,1.0}; static double e[5]={ 1.0,-0.1098628627e-02, 0.2734510407e-04,-0.2073370639e-05, 0.2093887211e-06}; static double f[5]={ -0.1562499995e-01, 0.1430488765e-03,-0.6911147651e-05, 0.7621095161e-06,-0.934935152e-07}; static double g[5]={ 1.0,0.183105e-02, -0.3516396496e-04,0.2457520174e-05, -0.240337019e-06}; static double h[5]={ 0.4687499995e-01, -0.2002690873e-03,0.8449199096e-05, -0.88228987e-06,0.105787412e-06}; if (n<0) n=-n; if (x<0.0) x=-x; if (x==0.0) return(-1.0e+70); if (n!=1) { if (x<8.0) { y=x*x; p=a[5]; q=b[5]; for (i=4; i>=0; i--) { p=p*y+a[i]; q=q*y+b[i];} p=p/q+0.636619772*Bessel(0,x)*log(x); } else { z=8.0/x; y=z*z; p=e[4]; q=f[4]; for (i=3; i>=0; i--) { p=p*y+e[i]; q=q*y+f[i];} s=x-0.785398164; p=p*sin(s)+z*q*cos(s); p=p*sqrt(0.636619772/x); } } if (n==0) return(p); b0=p; if (x<8.0) { y=x*x; p=c[5]; q=d[6]; for (i=4; i>=0; i--) { p=p*y+c[i]; q=q*y+d[i+1];} q=q*y+d[0]; p=x*p/q+0.636619772*(Bessel(1,x)*log(x)-1.0/x);; } else { z=8.0/x; y=z*z; p=g[4]; q=h[4]; for (i=3; i>=0; i--) { p=p*y+g[i]; q=q*y+h[i];} s=x-2.356194491; p=p*sin(s)+z*q*cos(s); p=p*sqrt(0.636619772/x); } if (n==1) return(p); b1=p; s=2.0/x; for (i=1; i<=n-1; i++) { p=s*i*b1-b0; b0=b1; b1=p; } return(p); } main() { int n; double x; double y0,y1,y2; FILE * fp = fopen("D://data.txt","w"); fprintf(fp,"x Neumann(0) Neumann(1) Neumann(2)/n"); n = 2 ; for (x=0.5;x<20;x+=0.02) { y0 = Neumann(0,x); y1 = Neumann(1,x); y2 = Neumann(2,x); fprintf(fp,"%6.3f %6.3f %6.3f %6.3f/n",x,y0,y1,y2); } fprintf(fp,"/n"); fclose(fp); }  

 

对得到的TXT数据作图如下所示：