复制内容到剪贴板
代码:
#include<iostream.h>
#include<math.h>
float f(float x)
{
return 1/(1+x*x);
}
float Romberg(float a,float b,float (*f)(float),float epsilon)
{
int n=1,k;
float h=b-a,x,temp;
float T1,T2,S1,S2,C1,C2,R1,R2;
T1=(b-a)/2*((*f)(a)+(*f)(b));
while(1)
{
temp=0;
for(k=0;k<=n-1;k++)
{
x=a+k*h+h/2;
temp+=(*f)(x);
}
T2=(T1+temp*h)/2;
if(fabs(T2-T1)<epsilon)
return T2;
S2=T2+(T2-T1)/3.0;
if(n==1)
{
T1=T2;
S1=S2;
h/=2;
n*=2;
continue;
}
C2=S2+(S2-S1)/15;
if(n==2)
{
C1=C2;
T1=T2;
S1=S2;
h/=2;
n*=2;
continue;
}
R2=C2+(C2-C1)/63;
if(n==4)
{
R1=R2;
C1=C2;
T1=T2;
S1=S2;
h/=2;
n*=2;
continue;
}
if(fabs(R2-R1)<epsilon)
return R2;
R1=R2;
C1=C2;
T1=T2;
S1=S2;
h/=2;
n*=2;
}
}
void main()
{
float epsilon=5e-6;
cout<<"R="<<Romberg(0,1,f,epsilon)<<endl;
}
