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

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

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
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<n;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<m;i++)
{x=a+i*h2;
if(x==0) y=2;
else y=2*(sin(x)/x);
s=s+y;}
for(i=0;i<m;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;}}。