复化梯形公式，复化抛物线公式和Romberg求积法的算法程序：
以下程序均定义误差限为1*10^-5;
1)复化梯形公式：
#include
#include
#define e 1e-5
#define a 0 //积分下限a
#define b 1 //积分上限b
#define f(x) (4/(1+(x*x))) //被积函数f(x)
int main()
{
int i,n;
double h,t0,t,g;
n=1; //赋初值
h=(double)(b-a)/2;
t=h*(f(a)+f(b));
do
{
t0=t;
g=0;
for (i=1;i<=n;i++)
g+=f((a+(2*i-1)*h));
t=(t0/2)+(h*g); //复化梯形公式
n*=2;
h/=2;
}
while (fabs(t-t0)>e); //自定义误差限e
printf("%.8lf",t); //输出积分的近似值
return 0;
}
2)复化抛物线公式：
#include
#include
#define e 1e-5
#define a 0 //积分下限a
#define b 1 //积分上限b
#define f(x) (4/(1+(x*x))) //被积函数f(x)
int main()
{
int i,n;
double f1,f2,f3,h,s0,s;
f1=f(a)+f(b); //赋初值
f2=f(((double)(b+a)/2));
f3=0;
s=((double)(b-a)/6)*(f1+4*f2);
n=2;
h=(double)(b-a)/4;
do //复化抛物线算法
{
f2+=f3;
s0=s;
f3=0;
for (i=1;i<=n;i++)
f3+=f((a+(2*i-1)*h));
s=(h/3)*(f1+2*f2+4*f3);
n*=2;
h/=2;
}
while (fabs(s-s0)>e); //自定义误差限
printf("%.8lf",s);
return 0;
}
3)Romberg求积法：
#include
#include
#define e 1e-5
#define a 0 //积分下限a
#define b 1 //积分上限b
#define f(x) (4/(1+(x*x))) //被积函数f(x)
double t[100][100];
int main()
{
int n,k,i,m;
double h,g,p;
h=(double)(b-a)/2;
t[0][0]=h*(f(a)+f(b));
k=1;
n=1;
do //Romberg算法
{
g=0;
for (i=1;i<=n;i++)
g+=f((a+((2*i-1)*h)));
t[k][0]=(t[k-1][0]/2)+(h*g);
for (m=1;m<=k;m++)
{
p=pow(4,(double)(m));
t[k-m][m]=(p*t[k-m+1][m-1]-t[k-m][m-1])/(p-1);
}
m-=1;
h/=2;
n*=2;
k+=1;
}
while (fabs(t[0][m]-t[0][m-1])>e); //自定义误差限e
printf("%.8lf",t[0][m]);
return 0;
}
给定精度，定义误差限为1*10^-5，分别求出步长的先验估计值：用复化梯形公式计算，要求h<0. 007746。用复化抛物线公式计算，要求h<0.131607。而实际上，在运用后验估计的程序中，以相同的精度，得到的复化梯形步长为h= 0.003906，得到的复化抛物线步长为h=0.0625，它们大致上分别为先验估计值的一半，符合要求。
在实践中比较上体三种算法的计算量，当取相同精度1*10^-5时，复化梯形调用函数f257次，孵化抛物线公式调用函数f17次，Romberg算法调用函数 f17次，从计算量上看，后两者较小。