复化梯形公式和复化辛卜生公式
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
实验三:分别用复化梯形公式和复化辛卜生公式计算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;}}