#include <stdio.h>
#include <conio.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;
S1=S2;
T1=T2;
h/=2;
n*=2;
}
}
void main()
{
float epsilon=5e-6;
printf("R=%f",Romberg(0,1,f,epsilon));
getch();
}
用龙贝格法计算积分
最新推荐文章于 2023-01-08 11:11:59 发布