C++实现龙贝格求积分算法

1. 算法原理简介

T(0)=(b-a)(f(a)+f(b))/2；
H(0)=(b-a)·f( (a+b)/2 )；
T(i)=(T(i-1)+H(i-1))/2。
T(i) 的语义是将积分区间做 2i 等分时复化梯形求积分公式算出的近似值。
H(i) 的语义是将积分区间做 2i 等分时，将每个小区间的长度乘该小区间中点处函数值的乘积进行累加求和的结果。

生成 i>=j 的积分近似值表 R(i, j）具体计算如下：.

• R(i,0)=T(i)=(T(i-1)+H(i-1))/2，i>=0，T(i)为区间逐次减半递推梯形求积分公式算出的结果；
• R(i,1)=S(i)，i>=1，S(i)为由T(i)、T(i-1)联合递推得出的区间逐次减半递推辛普森求积分公式算出的结果；
• R(i,2) = B(i)，i>=2，B(i) 为由S(i)、S(i-1)联合递推得出布尔求积分公式算出的结果；
  归纳总结有如下递推关系式：

由上面递推关系式计算积分近似值 R(i, j）构成三角表，循环结束条件 |R(i, i)-R(i-1, i-1)|< eps成立时，则以 R(i, i) 作为最终的积分近似值，否则继续循环计算。

2. 应用实例

用龙贝格积分法计算积分：

3. 程序代码

#include <iostream>
#include <cmath>
using namespace std;
#define N 20
#define MAX 10     //数组存的最大行数
#define a 0    //积分上限
#define b 1      //积分下限
#define eps 1e-6 //精度

double f(double x)//所求积分公式
{	
    return 4 / (1 + x * x);
}

double ComputeT(double aa, double bb, long int n)//复化梯形公式
{	
    int i;	
    double sum = 0, h = (bb - aa) / n;	
    for (i = 1; i < n; i++)	
    {		
        sum += f(aa + i * h);	
    }
    sum += (f(aa) + f(bb)) / 2;
    return h * sum;
}

double f2(double x)
{	
    return x * x;
}

int main()
{	
    cout << "被积函数 f = 4 / (1 + x * x)" << endl;	
    cout << "积分区间为（0，1）" << endl;	
    cout << "精度为1e-6" << endl;	
    int i;	
    long int n = N, m = 0;	
    double T[MAX + 1][2];	T[0][1] = ComputeT(a, b, n);	n *= 2;	
    for (m = 1; m < MAX; m++)	
    {		
        for (i = 0; i < m; i++)		
        {			
            T[i][0] = T[i][1];		
        }		
        T[0][1] = ComputeT(a, b, n);		
        n *= 2;		
        for (i = 1; i <= m; i++) //T的m(h)		
        {			
            T[i][1] = T[i - 1][1] + (T[i - 1][1] - T[i - 1][0]) / (pow(2, 2 * m) - 1);		
        }		
        if ((T[m - 1][1] < T[m][1] + eps) && (T[m - 1][1] > T[m][1] - eps))		
        {			
            cout << "积分值为：" << T[m][1] << endl;			
            return 0;		
        }	
    }
    cout << "此题无解！" << endl;	
    return 0;
} 

4.运行结果