编程实现数值积分

 

2010-11-05 09:52:43|  分类： 数值计算|字号 订阅

在区间 不大时 , 用梯形公式、辛卜生公式计算定积分是简单实用的 , 但当区间 较大时 , 用梯形公式、辛卜生公式计算定积分达不到精确度要求 . 为了提高计算的精确度，我们将 [a,b] 区间n等分，在每个小区间上应用梯形公式、辛卜生公式计算定积分，然后将其结果相加，这样就得到了复化梯形公式和复化辛卜生公式。 1. 复化梯形公式 将积分区间 等分 , 设 , 则节点为  对每个小区间上应用梯形公式 , 然后将其结果相加，则得 (3.14) 称 (3.14) 式为复化梯形公式 . 当 在 [a,b] 上有连续的二阶导数时，则复化梯形公式 (3.14) 的余项推导如下： 因为 所以在区间 [a,b] 上公式 (3.14) 的误差为 又因为 在区间 [a,b] 上连续，由连续函数的性质知，在区间 [a,b] 上存在一点 ， 于是 （ 3.15 ）

 

 

复化梯形公式，复化抛物线公式和Romberg求积法的算法程序：

以下程序均定义误差限为1*10^-5;

1)复化梯形公式：

#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)))           //被积函数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 <stdio.h>
#include <math.h>
#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 <stdio.h>
#include <math.h>
#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次，从计算量上看，后两者较小。