满意答案

辛普森法

#include

#include

#define F(X) (4.0/(1+X*X))

static float EPS = 1.0E-14;

static int COUNT=1;

static double a=0.0, b=1.0;

static double M[32],T[32],S[32],C[32],R[32],E[32];

double GETM(int K)

{ unsigned long j,n=1;

double x=0,y=0,step=0;

for(j=0;j

step = (b-a)/n;

x = a+step/2;

for(j=0;j

return (y*step);

}

Simpson(float EPS)

{ int k=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]

}

return;

}

ShawResult()

{ int k;

system("cls");

printf("\n K M[K] T[K] S[K]");

printf(" E[K]");

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

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

for(k=0;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()

{ int k;

FILE * fp;

fp=fopen("Simpson.htm","w");

if(!fp) return;

fprintf(fp,"

");

fprintf(fp,"

");

fprintf(fp,"

fprintf(fp,"\n K M[K] T[K] S[K]");

fprintf(fp," E[K]");

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

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

for(k=0;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();

}

10分享举报