復(fù)化梯形公式，復(fù)化拋物線(xiàn)公式和Romberg求積法的算法程序：

以下程序均定義誤差限為1*10^-5;

1)復(fù)化梯形公式：

#include <stdio.h>
#include <math.h>
#define e 1e-5
#define a 0                          //積分下限a
#define b 1                          //積分上限b
#define f(x) (4/(1+(x*x)))           //被積函數(shù)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);                //復(fù)化梯形公式
       n*=2;
       h/=2;
    }
    while (fabs(t-t0)>e);             //自定義誤差限e
    printf("%.8lf",t);                //輸出積分的近似值
  
return 0;
}

2)復(fù)化拋物線(xiàn)公式：

#include <stdio.h>
#include <math.h>
#define e 1e-5
#define a 0                          //積分下限a
#define b 1                          //積分上限b
#define f(x) (4/(1+(x*x)))           //被積函數(shù)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                                //復(fù)化拋物線(xiàn)算法
    {
                      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 <stdio.h>
#include <math.h>
#define e 1e-5
#define a 0                          //積分下限a
#define b 1                          //積分上限b
#define f(x) (4/(1+(x*x)))           //被積函數(shù)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，分別求出步長(zhǎng)的先驗(yàn)估計(jì)值：用復(fù)化梯形公式計(jì)算，要求h<0. 007746。用復(fù)化拋物線(xiàn)公式計(jì)算，要求h<0.131607。而實(shí)際上，在運(yùn)用后驗(yàn)估計(jì)的程序中，以相同的精度，得到的復(fù)化梯形步長(zhǎng)為h= 0.003906，得到的復(fù)化拋物線(xiàn)步長(zhǎng)為h=0.0625，它們大致上分別為先驗(yàn)估計(jì)值的一半，符合要求。

在實(shí)踐中比較上體三種算法的計(jì)算量，當(dāng)取相同精度1*10^-5時(shí)，復(fù)化梯形調(diào)用函數(shù)f257次，孵化拋物線(xiàn)公式調(diào)用函數(shù)f17次，Romberg算法調(diào)用函數(shù) f17次，從計(jì)算量上看，后兩者較小。