龙贝格求积公式也称为逐次分半加速法。它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上，构造出一种加速计算积分的方法。 作为一种外推算法，它在不增加计算量的前提下提高了误差的精度。
在等距基点的情况下，用计算机计算积分值通常都采用把区间逐次分半的方法进行。这样，前一次分割得到的函数值在分半以后仍可被利用,且易于编程。
一、变步长的梯形法
求
<1> 首先在上用梯形公式，得
<2>将二等分，步长，用复化梯形公式，得
<3>将四等分，步长 ，用复化梯形公式，得
<n>将将2n 等分，步长 ，用复化梯形公式，得
问题：
每一步是否能用到上一步得计算结果？区间划分到什么程度，才满足精度要求？
分析：
<1>、n等分，步长，
<2>、2n 等分，步长，记的中点(图中红线)为。
则积分节点为： （分为原节点 和 新增节点）
(注意：此是计算时的)
算法总结：
计算由递推公式判断 （为误差限），成立则返回并输出。
变步长的梯形法代码（C++)：
//变步长梯形积分double T2n(double(*f)(double), double a, double b, double error) {int n = 1; double h = (b - a) / n;double Tn = (f(a) + f(b)) * h / 2.0; //计算 n=1 时的积分，粗略值double T2n;//步长变为一半后的积分值bool key = false;do { //至少循环一次double Hn = 0.0;for (int k = 0; k < n; k++) {double x = a + (k + 0.5) * h;Hn += h * f(x);}T2n = (Tn + Hn) / 2.0;//判断误差if (fabs(Tn - T2n) < error) {key = true;}else {Tn = T2n;n *= 2;h /= 2;}} while (!key);return T2n;}
二、龙贝格算法
准确值数值计算近似值 为误差1、记为.
(序列)
还存在误差
2、求公式的余项，可得：
（柯特斯序列）
还存在误差
3、求柯特斯公式的余项，可得：
（龙贝格序列）
加权平均公式：
三、代码实现
例： ,程序结果为：
具体代码为(没有安装armadillo矩阵计算库的点这里)：
//#include<iostream>//#include<armadillo>//#include<iomanip>#include"romberg.h"using namespace std;using namespace arma;int main(){const double eps = 1e-6;int size = 50;Romberg.zeros(size, 4);Coefficient << 4 / 3.0 << 1 / 3.0 << endr<< 16 / 15.0 << 1 / 15.0 << endr<< 64 / 63.0 << 1 / 63.0 << endr;double a =1, b=2;cout << "请输入积分下限, a = ";cin >> a;cout << "请输入积分上限, b = ";cin >> b;cout << endl;double result = getSn(a, b, eps);cout << "积分结果为：" << setprecision(8) << result << endl;}
头文件为：
#pragma once#include<iostream>#include<math.h>#include<armadillo>#include<iomanip>#include<vector>using namespace std;using namespace arma;mat Romberg; //用于存放T、S、 C、Rmat Coefficient; //存放 计算系数// 记录坐标值vector<double> x_vec;vector<double> y_vec;// 积分函数double fun(double x) {double y = log(1 + pow(x, 2)) / (1 + pow(x, 2));return y;}//存储对应坐标void GreatXY(double a, double b, int n) {x_vec.clear();y_vec.clear();double h = (b - a) / n;for (int k = 0; k <= n; k++) {double x = a + h * k;double y = fun(x);x_vec.push_back(x);y_vec.push_back(y);}}/*变步梯形积分a:积分下限b:积分上限n:节点数*/double getTn(double a, double b, int n){GreatXY(a, b, n);double h = (b - a) / n;double Tn = h / 2.0 * (fun(a) + fun(b));for (int k = 1; k < n; k++) {double x = a + k * h;Tn += h * fun(x);}return Tn;}void FillMatrix(double &a, double &b) {//填充Tnfor (int k = 0; k <= 4; k++) { //先计算k =0,1,2,3,4的情况，即将S1、S2计算出来，在|S1-S2|不符合误差情况时在更新矩阵。Romberg(k, 0) = getTn(a, b, pow(2, k));}//填充Sn(i=1)、Cn(i=2)、Rn(i=3) //编程思想的行列，从零开始for (int i = 1; i <= 3; i++) { //列for (int j = i; j <= 4; j++) { //行Romberg(j, i) = Romberg(j, i - 1)*Coefficient(i - 1, 0) - Romberg(j - 1, i - 1)*Coefficient(i - 1, 1);}}}double getSn(double a, double b,double eps){FillMatrix(a, b);//判断误差if (fabs(Romberg(4, 3) - Romberg(3, 3)) < eps) {return Romberg(4, 3);}else {for (int k = 5;; k++) {Romberg(k, 0) = getTn(a, b, pow(2, k));//TRomberg(k, 1) = 4 / 3.0 *Romberg(k, 0) - 1 / 3.0*Romberg(k - 1, 0); //SRomberg(k, 2) = 16 / 15.0*Romberg(k, 1) - 1 / 15.0*Romberg(k - 1, 1);//CRomberg(k, 3) = 64 / 63.0 *Romberg(k, 2) - 1 / 63.0*Romberg(k - 1, 2);//R//再次比较if (fabs(Romberg(k, 3) - Romberg(k-1, 3)) < eps) {return Romberg(k, 3);break;}if (k >= 50) {cout << "循环次数太多，请更换程序算法";return -1;break;}}}}
