复化梯形公式和复化辛卜生公式

实验三：分别用复化梯形公式和复化辛卜生公式计算f(x)=sin(x)/x的积分，并与准确值比较判断精度。

#include

#include

#include

double Trapezoid(float,float,float,int);

double Simpson(float,float,float,int);

double Integral(float,float);

int SigDigT(double,double,int);

void Euler(float a[50],float b);

float f(float);

float Bisection(float,float);

main()

{char q;

double d1,d2;

float a,b,h1,h2;

int i,m,n;

printf("\t\t\t实验三：数值积分\n");

printf("实验题目：分别用复化梯形公式和复化辛卜生公式计算f(x)=sin(x)/x的积分，并与准确值比较判断精度。\n\n");

printf("请输入定积分的上下限：\n");

printf("定积分下限：a=");

scanf("%f",&a);

printf("定积分上限：b=");

scanf("%f",&b);

if(a==b)

{printf("积分上下限相等，定积分的值为0。\n");

goto sign1;}

else

printf("请输入复化梯形公式划分的份数：");

scanf("%d",&n);

printf("请输入复化辛卜生公式划分的份数：");

scanf("%d",&m);

h1=(b-a)/n;

h2=(b-a)/m;

printf("复化梯形公式的步长为：h=%f\n",h1);

printf("复化辛朴生公式的步长为：h=%f\n",h2);

printf("复化梯形公式的结果是：\nT=%f\n",Trapezoid(a,b,h1,n));

printf("复化辛卜生公式的结果是：\nS=%f\n",Simpson(a,b,h2,m));

printf("定义计算该公式的结果是：\nI=%f\n",Integral(a,b));

d1=fabs(Integral(a,b)-Trapezoid(a,b,h1,n));

d2=fabs(Integral(a,b)-Simpson(a,b,h2,m));

printf("复化梯形公式的误差是：\nd=%f\n",d1);

printf("复化辛卜生公式的误差是：\nd=%f\n",d2);

if(d1>d2) printf("复化梯形公式的精度高于复化辛卜生公式的精度\n");

else if(d1==d2) printf("复化辛卜生公式和复化梯形公式的精度相等\n");

else printf("复化辛卜生公式的精度高于复化梯形公式的精度\n");

for(i=0;i<2;i++)

{if(i==0) printf("划分%d份的复化梯形公式的有效位数是小数点后的第%d位\n",n,SigDigT(d1,d2,i));

else printf("划分%d份的复化辛卜生公式的有效位数是小数点后的第%d位\n",m,SigDigT(d1,d2,i));}

sign1:

printf("\n是否继续?(y/n)\n");

printf("注：选y重新输入，选n返回主界面\n");

getchar();

q=getchar();

if(q=='y')

system("cls");

else if(q=='n')

{system("cls");

main();}

else {printf("您输入的字符不可识别，请重新输入!\n");

goto sign1;}}

double Trapezoid(float a,float b,float h1,int n)

{double x,y,t=0;

int i;

for(i=1;i

{x=a+i*h1;

if(x==0) y=2;

else y=2*(sin(x)/x);

t=t+y;}

if(a==0) t=t+1+sin(b)/b;

else if(b==0) t=t+1+sin(a)/a;

else t=t+sin(a)/a+sin(b)/b;

t=(h1/2)*t;

return t;}

double Simpson(float a,float b,float h2,int m)

{double x,y,s=0;

int i;

for(i=1;i

{x=a+i*h2;

if(x==0) y=2;

else y=2*(sin(x)/x);

s=s+y;}

for(i=0;i

{x=a+((1+2*i)*h2)/2;

if(x==0) y=4;

else y=4*(sin(x)/x);

s=s+y;}

if(a==0) s=s+1+sin(b)/b;

else if(b==0) s=s+1+sin(a)/a; else s=s+sin(a)/a+sin(b)/b;

s=(h2/6)*s;

return s;}

double Integral(float a,float b) {double h,x,y,I=0;

int n=500000;

h=(b-a)/500000;

for(n=0;n<500000;n++)

{x=a+n*h;

if(x==0) y=1;

else y=sin(x)/x;

I=I+y*h;}

return I;}

SigDigT(double d1,double d2,int i) {int j=0,k=0;

double r=1,t;

if(i==0)

{do{r=r*(0.1);

t=d1/r;

j++;}while(t<0.1);

return j;}

else

{do{r=r*(0.1);

t=d2/r;

j++;}while(t<0.1);

return j;}}