学步园

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

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