#include<iostream.h>
#include<iomanip.h>
#include<math.h>
double a,b,erel;
//函数声明
double f(double);
double simpson(double &h);
double tixing(double &h);
//主函数
void main()
{
double h,s,t1,t2,s1,s2,c1,c2,temp;int i,m;
cout<<"输入左区间:";
cin>>a;
cout<<"输入右区间:";
cin>>b;
cout<<"输入误差上限:";
cin>>erel;
h=b-a;
s1=h*(f(a)+4*f((a+b)/2)+f(b))/6;
temp=erel;
for(m=0;temp<1;m++) temp*=10;
//用simpson公式计算结果
s=simpson(h);
s2=s1/2+h*s/6;
for(i=0;fabs(s1-s2)>erel;i++)
{ h=h/2;
s1=s2;
s=simpson(h);
s2=0.5*s1+h*s/6;
}
cout<<endl<<"自适应simpson公式最后的结果为:"<<setprecision(m+1)<<s2<<endl<<"迭代的次数为:"<<i<<endl;
//用梯形公式计算结果
h=b-a;
t1=0.5*h*(f(a)+f(b));
s=tixing(h);
t2=0.5*t1+h/2*s;
for(i=0;fabs(t1-t2)>erel;i++)
{ h=h/2;
t1=t2;
s=tixing(h);
t2=0.5*t1+h/2*s;
}
cout<<endl<<"自适应梯形公式最后的结果为:"<<s2<<endl<<"迭代的次数为:"<<i<<endl;
//Romberg算法结果如下
cout<<"\nRomberg算法如下:";
h=b-a;
t1=0.5*h*(f(a)+f(b));
cout<<endl<<setw(20)<<setiosflags(ios::left)<<setprecision(m+1)<<t1;
i=1;
bu2:s=tixing(h);t2=0.5*t1+h/2*s;
cout<<endl<<setw(20)<<setprecision(m+1)<<t2;
s2=t2+(t2-t1)/3;
cout<<setw(20)<<setprecision(m+1)<<s2;
if(i==1)
goto bu5;
else goto bu6;
bu5:i=i+1;h=h/2;t1=t2;s1=s2;goto bu2;
bu6:c2=s2+(s2-s1)/15;
cout<<setw(20)<<setprecision(m+1)<<c2;
if(i==2)
{ c1=c2;
goto bu5;
}
else goto bu8;
bu8:if(fabs(c2-c1)<erel) goto bu9;
else
{ c1=c2;
goto bu5;
}
bu9:cout<<"\nRomberg算法最后结果为:"<<setprecision(m+1)<<c2<<endl;
}
//函数定义
double f(double x)
{ return(4/(x*x+1));
}
//步长为h时候的simpson求和
double simpson(double &h)
{ int k;double s,n;
n=(b-a)/h;
for(k=0,s=0;k<n;k++)
{ s+=2*f(a+(k+0.25)*h)-f(a+(k+0.5)*h)+2*f(a+(k+0.75)*h);
}
return s;
}
//步长为h时候的梯形求和
double tixing(double &h)
{ int k;double s,n;
n=(b-a)/h;
for(k=0,s=0;k<n;k++)
{ s+=f(a+(k+0.5)*h);
}
return s;
}
