Romberg求积法是一种数值分析中的求积算法, 它利用了外推的思想方法, 以提高代数精度. 这种方法通过不断缩小积分区间, 将区间划分为更小的子区间, 并计算每个子区间的积分值, 从而逼近更为准确的求积值. 这个过程是通过计算一系列的近似值, 并逐步改进这些近似值, 以得到更为精确的结果. 在实现过程中, 通常采用逐次分半加速法, 将积分区间逐次分割为更小的子区间. 前一次分割得到的函数值在分半之后仍可被利用, 且易于编程实现.

算法说明

首先我们有复化梯形公式的余项:
$$I-T(h)=\sum _{k=1}^m\frac{B_{2k}}{(2k)!}[f^{(2k-1)}(a)-f^{(2k-1)}(b)]h^{2k}+r_{m+1}$$

其中, $B_{2k}$为Bernoulli常数,
$$r_{m+1}:=-\frac{B_{2m+2}}{\left(2m+2\right)!}(b-a)f^{\left(2m+2\right)}\left(\eta \right)h^{2m+2},\eta \in (a,b)$$

此外, 我们还有Richardson外推法如下:
假设$F(h)$逼近$F(0)$的余项为
$$F(h)-F(0)=a_1{h}^{p_1}+a_2{h}^{p_2}+\cdots$$
其中$p_1<p_2<\cdots$, $a_k$是与$h$无关的非零常数, 则由
$$\{\begin{array}{l}{F}_0(h)=F(h)\\ F_{k+1}(h)=\frac{q^{p_k}{F}_k(h/q)-F_k(h)}{q^{p_k}-1},k=0,1,2,\cdots \end{array}$$
定义的序列{Fk(h)}\{F_k(h)\}{Fk​(h)}有
$$F_n(h)-F(0)=a_{n+1}^{(n)}{h}^{p_{n+1}}+a_{n+2}^{(n)}{h}^{p_{n+2}}+a_{n+3}^{(n)}{h}^{p_{n+3}}+\cdots$$
其中, $a_{n+k}^{(n)}(k=1,2,\cdots )$与$h$无关, $q>1$.

利用复化梯形余项公式和Richardson外推法可得Romberg求积法:
$$\{\begin{array}{l}{T}_1^{(0)}=\frac{b-a}{2}[f(a)+f(b)]\\ T_1^{(i)}=\frac{1}{2}{T}_1^{(i-1)}+\frac{b-a}{2^i}\sum _{j=1}^{2^{i-1}}f\left[a+\frac{2j-1}{2^i}(b-a)\right],i=1,2,\cdots \\ T_{m+1}^{(k-1)}=\frac{4^m{T}_m^{(k)}-T_m^{(k-1)}}{4^m-1},k=1,2,\cdots ,i,m=1,2,\cdots \end{array}$$

对给定的精度标准$\epsilon$, 我们可由
$$|T_m^{(0)}-T_{m-1}^{(0)}|<\epsilon$$
作为迭代终止的标准.

算法实现

首先包含头文件如下:

#include <functional>
#include <math.h>
using namespace std;
typedef function<double(double)> func;

接着定义一个“T”表类:

class Romberg_Table
{
    double *arr;
    unsigned n, k;

public:
    Romberg_Table(const func &f, const double &a, const double &b, const double &fa, const double &fb, unsigned n0) : n(n0), k(0)
    {
        if (!n0)
            throw "n不能为零！";
        double *p(arr = new double[n0]), h(b - a);
        *p = h * (fa + fb) / 2;
        unsigned i(n0);
        while (--i)
        {
            double pre_h(h), x(a + (h /= 2)), fx(f(x)), &I(*p);
            if (isnan(fx))
            {
                delete[] arr;
                throw "区间内存在奇点！";
            }
            while ((x += pre_h) < b)
                if (isnan(fx += f(x)))
                {
                    delete[] arr;
                    throw "区间内存在奇点！";
                }
            *++p = I / 2 + fx * h;
        }
    }
    ~Romberg_Table()
    {
        delete[] arr;
    }
    // 外推, 注意n不能等于1
    void extrapolation() noexcept
    {
        unsigned k4(1 << (++k << 1)), k41(k4 - 1), m(--n); // k4表示4^k
        double *p(arr);
        *p = (k4 * *(p + 1) - *p) / k41;
        while (--m)
        {
            ++p;
            *p = (k4 * *(p + 1) - *p) / k41;
        }
    }
    // 输出当前值
    inline double &get_value() noexcept
    {
        return arr[n - 1];
    }
    // 判断是否达到精度要求
    inline bool judge_accuracy(const double &e) noexcept
    {
        if (n == 1)
            return false;
        return abs(arr[n - 2] - arr[n - 1]) >= e;
    }
};

接着实现Romberg求积法:

/*
 * Romberg求积法(端点可能为奇点)
 * f :被积函数
 * a :积分下限
 * b :积分上限
 * fa:积分下限处的函数值
 * fb:积分上限处的函数值
 * e :精度
 * n :复化次数
 */
double Romberg(func f, const double &a, const double &b, const double &fa, const double &fb, const double &e = 1e-8, const unsigned &n = 16)
{
    if (a >= b)
        throw "区间上限小于区间下限！";
    Romberg_Table table(f, a, b, fa, fb, n);
    while (table.judge_accuracy(e))
        table.extrapolation();
    return table.get_value();
}
/*
 * Romberg求积法
 * f :被积函数
 * a :积分下限
 * b :积分上限
 * e :精度
 * n :复化次数
 */
double Romberg(func f, const double &a, const double &b, const double &e = 1e-8, const unsigned &n = 16)
{
    double fa(f(a));
    if (isnan(fa))
        throw "积分下限为奇点！";
    double fb(f(b));
    if (isnan(fb))
        throw "积分上限为奇点！";
    return Romberg(f, a, b, fa, fb, e, n);
}

其中, 对于积分上下限为奇点的情形提供了手动输入积分上下限函数值的接口.

测试案例

计算积分$I={\int }_0^1\sqrt{x}\ln x\mathrm{d}x$.

#include <stdio.h>
// Romberg实现代码略
int main()
{
    auto f = [](double x) -> double
    {
        return sqrt(x) * log(x);
    };
    try
    {
        printf("%f", Romberg(f, 0., 1., 0., 0., 1e-6, 15));
    }
    catch (const char *s)
    {
        printf(s);
    }
    return 0;
}

最终得到结果为-0.444444, 计算过程如下: