c语言编程-微积分程序

辛普森法

#include

#include

#defineF(X)(4.0/(1+X*X))

staticfloatEPS=1.0E-14;

staticintCOUNT=1;

staticdoublea=0.0,b=1.0;

staticdoubleM[32],T[32],S[32],C[32],R[32],E[32];

doubleGETM(intK)

{unsignedlongj,n=1;

doublex=0,y=0,step=0;

for(j=0;j<K;j++)n*=2;

step=(b-a)/n;

x=a+step/2;

for(j=0;j<n;j++){y+=F(x);x+=step;}

return(y*step);

}

Simpson(floatEPS)

{intk=0;

T[0]=(b-a)*(F(a)+F(b))/2.0;

S[0]=T[0];

E[0]=1.0E10;

for(k=0;k<20;k++)

{M[k]=GETM(k);

T[k+1]=(M[k]+T[k])/2.0;

S[k+1]=(4.0*T[k+1]-T[k])/3.0;

COUNT++;

E[k+1]=fabs(S[k+1]-S[k]);

if(k<3)continue;

if(E[k+1]<EPS)break;

}

return;

}

ShawResult()

{intk;

system(cls);

printf(\nKM[K]T[K]S[K]);

printf(E[K]);

printf(\n-------------------------------------------);

printf(-----------------------------------);

for(k=0;k<COUNT;k++)

printf(\n%2d%20.15lf%20.15lf%20.15lf%12.4e,k,M[k],T[k],S[k],E[k]);

printf(\n--------------------------------------------);

printf(-----------------------------------);

getch();

system(cls);

return;

}

SaveResult()

{intk;

FILE*fp;

fp=fopen(Simpson.htm,w);

if(!fp)return;

fprintf(fp,);

fprintf(fp,);

fprintf(fp,);

fprintf(fp,\nKM[K]T[K]S[K]);

fprintf(fp,E[K]);

fprintf(fp,\n----------------------------------------------);

fprintf(fp,----------------------------------);

for(k=0;k<COUNT;k++)

fprintf(fp,\n%2d%20.15f%20.15f%20.15f%12.4le,k,M[k],T[k],S[k],E[k]);

fprintf(fp,\n----------------------------------------------);

fprintf(fp,----------------------------------);

fprintf(fp,);

fclose(fp);

return;

}

main()

{Simpson(EPS);

ShawResult();

SaveResult();

}