#include <stdio.h>
#include <conio.h>
#include <math.h>
float f(float x)
{
return exp(x)*cos(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=5e-6;
printf("R=%f",Romberg(0,1,f,epsilon));/*在Romberg(0,1,f,epsilon))中a=0,b=1,f為所調用子函數，epsilon為誤差精度*/
getch();
}
將a=0和b=3.14159,esplion=0.0001帶入的話答案應該是:-12.070289;此答案也可以用手工驗證的
f(x)=exp(x)cosx 的積分結果為:-0.5*e^(pie)-0.5=-12.07031,可以說結果還是符合要求的哦.
用復合梯形公式為：
#include <stdio.h>
#include <math.h>
#define epsilon 1e-6
int n;
void main()
{
int i;
float s;
float f(float);
float AutoTrap(float(*)(float),float,float);
s=AutoTrap(f,0.0,1.0); /*在AutoTrap(f,0.0,1.0)中: f為所調用子函數，a=0.0,b=1.0*/
printf("T(%d)=%f\n",n,s);
getch();
}
float AutoTrap(float(*f)(float),float a,float b)
{
int i;
float x,s,h=b-a;
float t1,t2=h/2.0*(f(a)+f(b));
n=1;
do
{
s=0.0;
t1=t2;
for(i=0;i<=n-1;i++)
{
x=a+i*h+h/2;
s+=f(x);
}
t2=(t1+s*h)/2.0;
n*=2;
h/=2.0;
}
while(fabs(t2-t1)>epsilon);
return t2;
}
float f(float x)
{
return exp(x)*cos(x);
}
復合Simpson公式算法為：
#include <stdio.h>
#include <conio.h>
#include <math.h>
void main()
{
int i,n=2;
float s;
float f(float);
float Simpson(float(*)(float),float,float,int);
for(i=0;i<=2;i++)
{
s=Simpson(f,0,1,n); /*在Simpson(f,0,1,n)中:f為所調用子函數, a=0,b=1，n為次數*/
printf("s(%d)=%f\n",n,s);
n*=2;
}
getch();
}
float Simpson(float(*f)(float),float a,float b,int n)
{
int k;
float s,s1,s2=0.0;
float h=(b-a)/n;
s1=f(a+h/2);
for(k=1;k<=n-1;k++)
{
s1+=f(a+k*h+h/2);
s2+=f(a+k*h);
}
s=h/6*(f(a)+4*s1+2*s2+f(b));
return s;
}
float f(float x)
{
return exp(x)*cos(x);
}
以上程序都是本人從書上摘抄且全部運行通過！！！