romberg.c

自己编写的龙贝格算法

#include <stdio.h> 
#include <conio.h> 
#include <math.h> 
float f(float x) 
{ 
return (x*x*x+sin(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;h/=2;n*=2; 
} 
} 
main() 
{ 
float epsilon=1e-5; 
printf("R=%f",Romberg(0.3,0.8,f,epsilon));/*在Romberg(0,1,f,epsilon))中a=0,b=1,f为所调用子函数，epsilon为误差精度*/ 
getch(); 
}