Romberg积分

在数值微分的最后我们提到了 Richardson 外展。 但由于数值微分本身并不稳定，Richardson 外展并没有发挥很大作用，我们也没有提供对应的例程。但数值积分不一样，它并不需要除以本身就很微小的采样步长 \(h\)， 所以数值稳定，Richardson 外展能在这个问题上发挥极大的作用。Romberg 积分就是这样的数值积分方法， 它在本身精度并不高的梯形法则上通过简单的加权求和显著提高了算法的计算精度。

本文我们从数学原理出发，伯努利数、伯努利多项式、欧拉-麦克劳林求和公式，获得关于采样步长 \(h\) 的级数形式，进而可以应用 Richardson 外展得到 Romberg 积分方法。 具体实现参考例程。

1. 伯努利数

首先我们来看一个经典的算术问题，对自然数的 \(k\) 次幂求和:

$$ S_{k}(n)=1^{k}+2^{k}+3^{k}+\cdots +n^{k} $$

当 \(k = 1,2,3\) 时分别对应着等差数列求和 \(S_1(n) = n(n+1)/2\)、平方和公式 \(S_2(n) = n(n+1)(2n+1)/6\)、立方和公式 \(S_3(n) = \left[n(n+1)/2\right]^2\)。 当 \(k = 0\) 时，有 \(S_0(n) = n\)。当幂次 \(k\) 不断增大时，想要找到一个简单的求和公式并不容易。Johann Faulhaber 曾手算到了 17 次幂。 后来伯努利(Jacob Benoulli)在观察低次幂公式的时候，敏锐地捕捉到了一组系数，能够写出任意幂次的求和公式:

$$ \begin{equation}\label{f1} S_{k}(n)=\frac{1}{k+1}\sum _{j=0}^{k}(-1)^j{k+1 \choose j}B_{j}n^{k+1-j} \end{equation} $$

上式中系数 \(B_j\) 就是伯努利数，给定首项 \(B_0 = 1\)，其余伯努利数都可以通过如下的递推公式求出:

$$ \begin{equation}\label{f2} \sum_{k=0}^n {n+1 \choose k} B_k = 0 \Longrightarrow B_n = -\frac{1}{n+1}\sum_{k=0}^{n+1} {n+1 \choose k} B_k \end{equation} $$

下表是前 12 的伯努利数，从中我们可以发现一个规律，除了第一个奇数项 \(B_1 = -\frac{1}{2}\) 之外，其余的奇数项比如 \(B_3, B_5, \cdots\) 都是 0。

有些材料可能觉得式(\(\ref{f1}\))中的 \((-1)^j\) 太麻烦了，他们把 \(B_1 = -\frac{1}{2}\) 的伯努利数列称为第一类伯努利数列，用符号 \(B_{n}^{-}\) 或 \(B_{n}\) 表示。 然后定义了第二类伯努利数列，用符号 \(B_{n}^{+}\) 标记，并规定 \(B_1^+ = +\frac{1}{2}\)。这样式(\(\ref{f1}\))就可以改写成如下的形式:

$$ \begin{equation}\label{f3} S_{k}(n)=\frac{1}{k+1}\sum _{j=0}^{k}{k+1 \choose j}B_{j}^+ n^{k+1-j} \end{equation} $$

2. 伯努利数的生成函数

根据伯努利数的递推式(\(\ref{f2}\))逐个向下计算 \(B_n\) 是很繁琐的。在各种契机下，欧拉构造出了伯努利数的生成函数。我们知道函数 \(e^t\) 的泰勒级数为:

$$ \begin{equation}\label{f4} e^t = 1 + t + \frac{1}{2!} t^2 + \cdots + \frac{1}{n!}t^n + \cdots = \sum_{n=0}^{\infty} \frac{t^n}{n!} \end{equation} $$

欧拉以 \(\{B_n\}\) 为系数，构造了一个无穷级数:

$$ \begin{equation}\label{f5} S(t) = B_0 + B_1 t + B_2 \frac{1}{2!} t^2 + \cdots + B_n \frac{1}{n!}t^n + \cdots = \sum_{n=0}^{\infty} B_n \frac{t^n}{n!} \end{equation} $$

然后他把 \(S(t)\) 和展开后的 \(e^t - 1\) 相乘得到:

$$ \begin{equation}\label{f6} S(t)(e^t - 1) = \left(B_0 + B_1 t + B_2 \frac{1}{2!} t^2 + \cdots \right) \cdot \left(t + \frac{t^2}{2!} + \frac{t^3}{3!} + \cdots \right) \end{equation} $$

将上式展开会发现，除了首项 \(t\) 的系数是 \(B_0 = 1\) 之外，其余各项的系数都是 \(\sum_{k=0}^n {n+1 \choose k} B_k = 0\)，这正是伯努利数递推式。于是有:

$$ \begin{equation}\label{f7} S(t) \cdot (e^t - 1) = t \Longleftrightarrow \frac{t}{e^t - 1} = \sum_{n=0}^{\infty} B_n \frac{t^n}{n!} \end{equation} $$

上式正是伯努利数的生成函数。说它是生成函数，并不是说通过它能直接算出 \(B_n\)，而是说它把一个无穷的离散的数列打包进一个连续的函数。 这样面对很多离散问题，我们就可以换个角度从连续函数来分析它们。

3. 伯努利多项式

据说欧拉为了研究离散求和 \(\sum f(k)\) 与连续求积 \(\int f(x)dx\) 之间的差距，他定义了一个多项式 \(B_n(x)\) 满足导数关系:

$$ \begin{equation}\label{f8} \frac{d}{dx} B_n(x) = n B_{n-1}(x) \end{equation} $$

并且满足积分关系:

$$ \begin{equation}\label{f9} \int_0^1 B_n(x) dx = 0 \end{equation} $$

定义首项 \(B_0(x) = 1\) 欧拉写出了几个多项式:

$$ \begin{aligned} B_0(x) & = 1 \\ B_1(x) & = x - \frac{1}{2} \\ B_2(x) & = x^2 - x + \frac{1}{6} \\ B_3(x) & = x^3 - \frac{3}{2}x^2 + \frac{1}{2}x \\ \end{aligned} $$

这些多项式的常数项 \(B_n(0)\) 恰好就是伯努利数 \(B_n\)，欧拉将这些多项式称为伯努利多项式，整理成通式有:

$$ \begin{equation}\label{f10} B_n(x) = \sum_{k=0}^n {n \choose k} B_k x^{n-k} \end{equation} $$

于是乎，欧拉有定义了一个级数 \(S(x,t) = \sum_{n=0}^{\infty} B_n(x) \frac{t^n}{n!}\)。将上式(\(\ref{f10}\))代入并展开， 利用柯西乘积公式 $$ \left(\sum _{k=0}^{\infty }a_{k}t^{k}\right)\left(\sum _{m=0}^{\infty }b_{m}t^{m}\right)=\sum _{n=0}^{\infty }\left(\sum _{k=0}^{n}a_{k}b_{n-k}\right)t^{n} $$ 有:

$$ \begin{aligned} S(x,t) & = \sum_{n=0}^{\infty}\left[\sum_{k=0}^{n}{n \choose k}B_{k}x^{n-k}\right]\frac{t^n}{n!} \\ & = \sum_{n=0}^{\infty}\left[\sum_{k=0}^{n}\frac{n!}{k!(n-k)!}B_{k}x^{n-k}\right]\frac{t^{n}}{n!} = \sum_{n=0}^{\infty}\left[\sum_{k=0}^{n}\left(\frac{B_{k}}{k!}\right)\left(\frac{x^{n-k}}{(n-k)!}\right)\right]t^{n} \\ & = \underbrace{\left(\sum _{k=0}^{\infty }B_{k}\frac{t^{k}}{k!}\right)}_{\frac{t}{e^t-1}} \underbrace{\left(\sum _{m=0}^{\infty }\frac{(xt)^{m}}{m!}\right)}_{e^{xt}} \end{aligned} $$

这正是伯努利多项式的生成函数:

$$ \begin{equation}\label{f11} \frac{t e^{xt}}{e^t - 1} = \sum_{n=0}^{\infty}B_n(x) \frac{t^n}{n!} \end{equation} $$

4. 欧拉-麦克劳林求和公式

在伯努利多项式的基础上，人们证明了 Euler-Maclaurin 求和公式。该公式揭示了离散求和与连续求积之间的关系。 记 \(a,b\) 为两个整数，那么在区间 \([a,b]\) 对一个足够光滑的函数 \(f(x)\) 积分 \(\int_a^b f(x) dx\) 与累加运算 \(\sum_{i=a}^b f(x)\) 有如下的关系:

$$ \begin{equation}\label{f12} \sum_{i=a}^{b}f(i) = \int_{a}^{b}f(x) dx + \frac{f(a)+f(b)}{2} + \sum_{k=1}^{m}\frac{B_{2k}}{(2k)!}\left(f^{(2k-1)}(b)-f^{(2k-1)}(a)\right) + R_{m} \end{equation} $$

上式中 \(B_{2k}\) 正是我们在前面讨论的伯努利数。因为在伯努利数列中除了 \(B_1\) 之外的奇数项都为 0，所以上式中只保留了其中的偶数项。 \(R_m\) 是一个余项。根据伯努利多项式 \(B_n(x)\) 的导数关系式(\(\ref{f8}\))，我们引入 \(B_1(x) = x - \frac{1}{2}\)，对 \(f(x)\) 进行在单位区间 \([0, 1]\) 上分部积分有:

$$ \int_0^1 f(x) dx = \int_0^1 f(x)B_1'(x) dx = \left . f(x)B_1(x)\right |_0^1 - \int_0^1 f'(x)B_1(x) dx $$

于是有

$$ \begin{equation}\label{f13} \int_0^1 f(x) dx = \frac{f(1) + f(0)}{2} - \int_0^1 f'(x)B_1(x) dx \end{equation} $$

我们对上式中的第二项再做一次分部积分有:

$$ \int_0^1 f'(x)B_1(x) dx = \frac{1}{2}\int_0^1 f'(x)B_2'(x) dx = \frac{1}{2}\left(\left . f'(x)B_2(x) \right |_0^1 - \int_0^1 f''(x)B_2(x) dx\right) $$

因为伯努利多项式的积分性质式(\(\ref{f9}\))，有 \(B_n(1) = B_n(0) = B_n\)，代入上式有

$$ \begin{equation}\label{f14} \int_0^1 f'(x)B_1(x) dx =\frac{B_2}{2}\left(f'(1) - f'(0)\right) - \frac{1}{2}\int_0^1 f''(x)B_2(x) dx \end{equation} $$

将式(\(\ref{f14}\))代入式(\(\ref{f13}\))，有

$$ \begin{equation}\label{f15} \int_0^1 f(x) dx = \frac{f(1) + f(0)}{2} - \frac{B_2}{2}\left(f'(1) - f'(0)\right) + \frac{1}{2}\int_0^1 f''(x)B_2(x) dx \end{equation} $$

不断对最后一项重复分部积分的操作，有

$$ \begin{equation}\label{f16} \int_0^1 f(x) dx = \frac{f(1) + f(0)}{2} - \sum_{k=1}^m\frac{B_{2k}}{(2k)!}\left(f^{(2k-1)}(1) - f^{(2k-1)}(0)\right) + R_m \end{equation} $$

将上式整理一下就得到了单位区间 \([0,1]\) 上的 Euler-Maclaurin 求和公式:

$$ \begin{equation}\label{f17} f(1) + f(0) = \int_0^1 f(x) dx + \frac{f(1) + f(0)}{2} + \sum_{k=1}^m\frac{B_{2k}}{(2k)!}\left(f^{(2k-1)}(1) - f^{(2k-1)}(0)\right) + R_m \end{equation} $$

对于一般的整数 \(a, b \in \mathbb{Z}\) 且 \(a < b\)，我们可以将区间 \([a,b]\) 拆分为 \(b-a\) 个长度为 1 的子区间 \([a, a+1], [a+1, a+2], \dots, [b-1, b]\)。 定积分 \(\int_a^b f(x) dx\) 就可以写成这些子区间定积分的和。然后我们分别对每个子区间重复上述单位区间上的分部积分，经过一系列的消减操作之后，就可以得到式(\(\ref{f12}\))。 更一般情况下，我们面对的是一个任意的区间 \([a,b]\)，即 \(a,b\) 不会是整数。此时，我们可以引入一个步长 \(h\)，把整个区间 \([a,b]\) 等分成 \(n\) 份。 步长为 \(h = \frac{b-a}{n}\)。此时我们就会得到一个带步长的 Euler-Maclaurin 公式: $$ \begin{equation}\label{f18} \sum _{i=0}^{n}f(a+ih)=\frac{1}{h}\int _{a}^{b}f(x)\,dx+\frac{f(a)+f(b)}{2}+\sum _{k=1}^{m}\frac{B_{2k}}{(2k)!}h^{2k-1}\left(f^{(2k-1)}(b)-f^{(2k-1)}(a)\right)+R_{m} \end{equation} $$ 这个公式正是本文要找的，将Richardson外展应用于数值积分的理论依据。 基于此对梯形法则进行Richardson外展的方法就是 Romberg 积分方法。

5. Romberg积分

在数值积分中，为了计算 \(\int_a^b f(x) dx\)，我们把区间 \([a,b]\) 进行 \(n\) 等分，步长为 \(h = \frac{b-a}{n}\)，节点为 \(x_j = a + jh\)。 对式(\(\ref{f18}\))做如下变形:

$$ \begin{equation}\label{f19} \underbrace{h\left[\frac{1}{2}f(a)+\sum_{j=1}^{n-1}f(x_{j})+\frac{1}{2}f(b)\right]}_{T(h)} = \int_{a}^{b}f(x) dx + \underbrace{\sum_{k=1}^{m}\frac{B_{2k}}{(2k)!}h^{2k}\left[f^{(2k-1)}(b)-f^{(2k-1)}(a)\right]}_{B(h)} + O(h^{2m+2}) \end{equation} $$

上式左侧就是组合分段的梯形公式，将之记为 \(T(h)\)。右侧的第二项是一个关于 \(h\) 的级数，给定区间 \([a,b]\) 其中 \(\frac{B_{2k}}{(2k)!} [f^{(2k-1)}(b) - f^{(2k-1)}(a)]\) 是一个与步长 \(h\) 无关的常数，我们可以用符号 \(\beta_k\) 来替代它，可以得到一个很漂亮的形式:

$$ \begin{equation}\label{f20} T(h) = \int_a^b f(x)dx + \beta_1 h^2 + \beta_2 h^4 + \cdots + \beta_m h^{2m} + O(h^{2m+2}) \end{equation} $$

写出了级数，我们就可以套用Richardson外展的套路，先计算 \(N(h)\)，再减半步长计算 \(N(h/2)\)，最后加权消元。

$$ \begin{equation}\label{f21} \begin{cases} T(h) & = \int_a^b f(x)dx + \beta_1 h^2 + \beta_2 h^4 + \cdots \\ T(\frac{h}{2}) & = \int_a^b f(x)dx + \frac{1}{4}\beta_1 h^2 + \frac{1}{16}\beta_2 h^4 + \cdots \end{cases} \Longrightarrow \frac{4T(\frac{h}{2}) - T(h)}{3} = \int_a^b f(x)dx - \frac{1}{4}\beta_2 h^4 - \frac{5}{16}\beta_3h^6 + \cdots \end{equation} $$

如上式所示，通过简单的加权求和，我们可以直接消除 \(h^2\) 项，得到 Simpson 序列。如果持续相同的操作，还可以依次消除 \(h^4, h^6, \cdots\)。 整个消除过程可以用一张表来规范化的表示出来，就得到了 Romberg 积分。

我们用 \(R(k,j)\) 来表示表中第 \(k\) 行 \(j\) 列的元素。那么第一列直接按照复合梯形法则依次按照减半步长 \(h\) 计算。后续的外推列可以直接按照如下递推公式得出:

$$ R(k,j) = \frac{4^j \cdot R(k,j-1) - R(k-1, j-1)}{4^j-1} $$

具体实现参考 例程。

6. 完