龙贝格积分

龙贝格积分（Romberg integration）把逐次减半步长的复合梯形公式组织成一个三角表，再用理查森外推逐层消去主导误差项。它利用光滑函数的梯形误差通常按偶次步长展开这一结构，而不是仅把更多采样值作平均。

以 $I={\int }_0^1{x}^4\mathrm{d}x=1/5$ 为教学算例。梯形网格从一个区间细化为两个、四个区间，近似值依次为 $1/2,9/32,113/512$。仅靠细化还未精确；对这三项作两层外推，会得到恰好 $1/5$。

梯形网格复用旧节点

令 $h_k=(b-a)/2^k$，以 $T_k$ 记步长 $h_k$ 的复合梯形近似。第一项为 $T_0=(b-a)(f(a)+f(b))/2$。细化时旧节点权重减半，只需补上原各小区间的中点，因此 $$T_k=\frac{1}{2}{T}_{k-1}+h_k\sum _{i=1}^{2^{k-1}}f(a+(2i-1)h_k).$$ 主例有 $T_0=1/2$，然后 $$T_1=\frac{1}{4}+\frac{1}{2}{\left(\frac{1}{2}\right)}^4=\frac{9}{32},$$ $$T_2=\frac{1}{2}\frac{9}{32}+\frac{1}{4}[{\left(\frac{1}{4}\right)}^4+{\left(\frac{3}{4}\right)}^4]=\frac{113}{512}.$$ 每次只计算新出现的奇数网格点，不必重新评价全部旧节点。

消去偶次误差的递推证明

在足够光滑的条件下，梯形误差有展开 $$T(h)=I+c_2{h}^2+c_4{h}^4+\cdots .$$ 例如要按此展开论证到指定有限阶，需相应高阶导数在整个闭区间存在并连续，不能只检查节点处函数值有限。比较步长 $h$ 与 $2h$，两项的主导误差分别为 $c_2{h}^2$ 和 $4c_2{h}^2$，所以 $$\frac{4T(h)-T(2h)}{3}=I+O(h^4).$$ 系数四来自二次误差在步长加倍后变为四倍；分母三则使常数项仍保持为 $I$。

定义 $R_{k,0}=T_k$。第 $j$ 层消去 $h^{2j}$ 项，递推为 $$R_{k,j}=\frac{4^j{R}_{k,j-1}-R_{k-1,j-1}}{4^j-1},1\le j\le k.$$ 其正确性可作归纳：若上一层在步长 $h_k$ 下为 $I+C{h}_k^{2j}+O(h_k^{2j+2})$，同层上一行的步长为 $2h_k$，主项系数变为 $4^{j}C$。按上式相减便消去该项，常数项又除以 $4^j-1$ 保持不变，余项为 $O(h_k^{2j+2})$。因此每加一列，在光滑展开有效的前提下提高两阶，而非由某一次结果恰好很准推断。

主例的完整三角表

第一列是梯形值，其余列依次计算为 $$R_{1,1}=\frac{4(9/32)-1/2}{3}=\frac{5}{24},R_{2,1}=\frac{4(113/512)-9/32}{3}=\frac{77}{384},$$ $$R_{2,2}=\frac{16(77/384)-5/24}{15}=\frac{1}{5}.$$

每个新条目使用左侧本行与上一行的同列条目；第二次外推的系数是十六。

对这个四次多项式，梯形误差实际上恰为 $T(h)-I=h^2/3-h^4/30$。第一列之后只有这两个误差幂，两层外推把它们都消掉，因此精确。这一特殊精确性不表示一般函数也会在四个小区间后得到精确积分。

变式。 将被积函数改成 $x^2$，则 $T_0=1/2,T_1=3/8$，一次外推就得到 $(4\cdot 3/8-1/2)/3=1/3$。因为此时梯形误差只有 $h^2$ 项，一层已经足够。

固定层数的Python实现

函数返回全部三角表，输入层数 $L\ge 1$ 表示行编号零到 $L-1$，最细网格有 $2^{L-1}$ 个小区间。它不进行自适应停止，也不改写用户数据；每个采样值必须有限。区间端点允许相同，此时直接返回零积分表，不调用被积函数。

Python 3

from math import fsum, isfinite

def romberg_table(function, a, b, levels):
    a, b = float(a), float(b)
    if not isfinite(a) or not isfinite(b) or b < a:
        raise ValueError("finite ordered endpoints required")
    if type(levels) is not int or levels < 1:
        raise ValueError("levels must be a positive integer")
    if a == b:
        return [[0.0]*(k+1) for k in range(levels)]
    def sample(x):
        value = float(function(x))
        if not isfinite(value):
            raise ValueError("nonfinite integrand sample")
        return value
    table = [[(b-a)*(sample(a)+sample(b))/2]]
    for k in range(1, levels):
        h = (b-a) / 2**k
        new_sum = fsum(sample(a+(2*i-1)*h) for i in range(1, 2**(k-1)+1))
        row = [table[k-1][0]/2 + h*new_sum]
        for j in range(1, k+1):
            factor = 4**j
            row.append((factor*row[j-1]-table[k-1][j-1])/(factor-1))
        table.append(row)
    return table

table = romberg_table(lambda x: x**4, 0, 1, 3)
assert table[0][0] == 1/2
assert table[1][0] == 9/32
assert table[2][0] == 113/512
assert abs(table[2][2] - 1/5) < 1e-14
assert abs(romberg_table(lambda x: x*x, 0, 1, 2)[1][1] - 1/3) < 1e-14
assert romberg_table(lambda x: 1/0, 2, 2, 2) == [[0.0], [0.0, 0.0]]
cusp = romberg_table(lambda x: abs(x-0.5), 0, 1, 2)
assert cusp[1][0] == 1/4 and abs(cusp[1][1] - 1/6) < 1e-14
print(table)

到第 $L-1$ 行共需 $2^{L-1}+1$ 次函数评价，旧节点不重算；外推表有 $O(L^2)$ 个条目，因此时间为 $O(2^L+L^2)$，保存完整表占 $O(L^2)$ 空间。若只保留前后两行可减少表存储，但本实现保留全部步骤便于复核。

光滑性失效与停止误差

取 $f(x)=|x-1/2|$，真积分为 $1/4$。两小区间的梯形公式已恰好精确，但粗网格 $T_0=1/2$，所以外推值却为 $(4\cdot 1/4-1/2)/3=1/6$，反而变差。原因是中点处有尖角，跨越不同网格的误差不满足所用的固定偶次展开。

相邻对角条目的差常用作经验误差指标，却不是无条件的严格误差界。端点奇异、强振荡、舍入误差或函数评价误差都可能破坏外推改进；必要时应先分割奇点、作变量替换，或选择适合该结构的求积公式。