#include "stdio.h" #include "math.h" double first_modified_Bessel(int n,double x) { int i,m; double t,y,p,b0,b1,q; static double a[7]={ 1.0,3.5156229,3.0899424,1.2067492, 0.2659732,0.0360768,0.0045813}; static double b[7]={ 0.5,0.87890594,0.51498869, 0.15084934,0.02658773,0.00301532,0.00032411}; static double c[9]={ 0.39894228,0.01328592,0.00225319, -0.00157565,0.00916281,-0.02057706, 0.02635537,-0.01647633,0.00392377}; static double d[9]={ 0.39894228,-0.03988024,-0.00362018, 0.00163801,-0.01031555,0.02282967, -0.02895312,0.01787654,-0.00420059}; if (n<0) n=-n; t=fabs(x); if (n!=1) { if (t<3.75) { y=(x/3.75)*(x/3.75); p=a[6]; for (i=5; i>=0; i--) p=p*y+a[i]; } else { y=3.75/t; p=c[8]; for (i=7; i>=0; i--) p=p*y+c[i]; p=p*exp(t)/sqrt(t); } } if (n==0) return(p); q=p; if (t<3.75) { y=(x/3.75)*(x/3.75); p=b[6]; for (i=5; i>=0; i--) p=p*y+b[i]; p=p*t; } else { y=3.75/t; p=d[8]; for (i=7; i>=0; i--) p=p*y+d[i]; p=p*exp(t)/sqrt(t); } if (x<0.0) p=-p; if (n==1) return(p); if (x==0.0) return(0.0); y=2.0/t; t=0.0; b1=1.0; b0=0.0; m=n+(int)sqrt(40.0*n); m=2*m; for (i=m; i>0; i--) { p=b0+i*y*b1; b0=b1; b1=p; if (fabs(b1)>1.0e+10) { t=t*1.0e-10; b0=b0*1.0e-10; b1=b1*1.0e-10; } if (i==n) t=b0; } p=t*q/b1; if ((x<0.0)&&(n%2==1)) p=-p; return(p); } double second_modified_Bessel(int n,double x) { int i; double y,p,b0,b1; static double a[7]={ -0.57721566,0.4227842,0.23069756, 0.0348859,0.00262698,0.0001075,0.0000074}; static double b[7]={ 1.0,0.15443144,-0.67278579, -0.18156897,-0.01919402,-0.00110404,-0.00004686}; static double c[7]={ 1.25331414,-0.07832358,0.02189568, -0.01062446,0.00587872,-0.0025154,0.00053208}; static double d[7]={ 1.25331414,0.23498619,-0.0365562, 0.01504268,-0.00780353,0.00325614,-0.00068245}; if (n<0) n=-n; if (x<0.0) x=-x; if (x==0.0) return(1.0e+70); if (n!=1) { if (x<=2.0) { y=x*x/4.0; p=a[6]; for (i=5; i>=0; i--) p=p*y+a[i]; p=p-first_modified_Bessel(0,x)*log(x/2.0); } else { y=2.0/x; p=c[6]; for (i=5; i>=0; i--) p=p*y+c[i]; p=p*exp(-x)/sqrt(x); } } if (n==0) return(p); b0=p; if (x<=2.0) { y=x*x/4.0; p=b[6]; for (i=5; i>=0; i--) p=p*y+b[i]; p=p/x+first_modified_Bessel(1,x)*log(x/2.0); } else { y=2.0/x; p=d[6]; for (i=5; i>=0; i--) p=p*y+d[i]; p=p*exp(-x)/sqrt(x); } if (n==1) return(p); b1=p; y=2.0/x; for (i=1; i<n; i++) { p=b0+i*y*b1; b0=b1; b1=p; } return(p); } main() { int n; double x; double y0,y1; FILE * fp = fopen("D://data.txt","w"); fprintf(fp,"x K0(x) K1(x)/n"); for (x=0.5;x<6;x+=0.01) { y0 = second_modified_Bessel(0,x); y1 = second_modified_Bessel(1,x); fprintf(fp,"%6.3f %6.3f %6.3f/n",x,y0,y1); } fprintf(fp,"/n"); fclose(fp); } 对得到的TXT数据作图如下所示: