第一類修正貝塞爾函數

#include "stdio.h"<br />#include "math.h"<br />/*****************************************************************<br /> 第一類變型貝塞爾函數/<br /> 第一類修正貝塞爾函數/<br /> 第一類變態貝塞爾函數/<br /> 第一類虛宗量貝塞爾函數/<br /> 雙曲型貝塞爾函數<br />******************************************************************/<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 />main()<br />{<br />int n;<br />double x;<br />double y0,y1,y2,y3,y4;<br />FILE * fp = fopen("D://data.txt","w");<br />fprintf(fp,"x I0(x) I1(x) I2(x) I3(x) I4(x)/n");<br />for (x=0.0;x<6;x+=0.01)<br />{<br />y0 = first_modified_Bessel(0,x);<br />y1 = first_modified_Bessel(1,x);<br />y2 = first_modified_Bessel(2,x);<br />y3 = first_modified_Bessel(3,x);<br />y4 = first_modified_Bessel(4,x);<br />fprintf(fp,"%6.3f %6.3f %6.3f %6.3f %6.3f %6.3f/n",x,y0,y1,y2,y3,y4);<br />}<br />fprintf(fp,"/n");<br />fclose(fp);<br />}<br /> 

 

對得到的TXT資料作圖如下所示：