4.4龙贝格（Romberg）求积公式 +本研授课

龙贝格（Romberg）求积公式 深度讲解

龙贝格求积公式是数值积分中高效高精度的自适应加速算法，核心是在复合梯形公式递推二分的基础上，通过理查森（Richardson）外推技巧，用低精度的梯形值线性组合出高精度的积分结果，完美解决了复合梯形公式收敛慢、高精度需求下计算量过大的问题，是工程中高精度数值积分的首选算法之一。

---

一、算法背景：复合梯形公式的痛点与递推化

1. 复合梯形公式的收敛痛点

复合梯形公式具备稳定、易实现的优点，但误差阶仅为\(O(h^2)\)，收敛速度极慢。例如例4.5中的正弦积分\(\int_0^1 \frac{\sin x}{x}dx\)，要达到7位有效数字，需要将区间二分10次，计算1025个节点的函数值，计算成本极高。

要提升效率，核心思路是：不额外增加大量函数值计算，仅通过已有结果的线性组合消除低阶误差项，直接提升精度阶数。而实现这个思路的前提，是先把复合梯形公式改造成可递推的形式，避免步长二分时重复计算函数值。

2. 梯形法的递推化（步长二分递推公式）

区间二分的核心特点

将积分区间\([a,b]\)先做\(n\)等分，步长\(h=\frac{b-a}{n}\)，节点\(x_k = a+kh\)，得到复合梯形值\(T_n\)。
若将步长减半（区间二分）得到\(2n\)等分，此时原有的\(n+1\)个节点全部保留，仅在每个原小子区间\([x_k,x_{k+1}]\)中新增一个中点\(x_{k+1/2} = \frac{x_k+x_{k+1}}{2}\)，无需重复计算原有节点的函数值，仅需计算\(n\)个新增中点的函数值。

递推公式推导

原\(n\)等分的复合梯形公式：

\[T_n = \frac{h}{2}\left[ f(a) + 2\sum_{k=1}^{n-1}f(x_k) + f(b) \right] \]

二分后\(2n\)等分的复合梯形公式，可拆分为原有节点和新增中点两部分，化简后得到梯形法的步长二分递推公式：

\[T_{2n} = \frac{1}{2}T_n + \frac{h}{2}\sum_{k=0}^{n-1}f(x_{k+1/2}) \tag{4.17} \]

递推公式的核心优势

• 避免重复计算：二分后仅需计算新增中点的函数值，原有节点结果全部复用，计算量减半；
• 实现简单：可通过循环不断二分区间，自动递推计算梯形值，无需提前确定等分数，便于实现精度自适应控制。

3. 递推公式实例（例4.5）

以正弦积分\(I=\int_0^1 \frac{\sin x}{x}dx\)为例，补充定义\(f(0)=1\)，递推计算过程如下：

1. 初始步长\(h=1\)（\(n=1\)），\(T_1 = \frac{1}{2}[f(0)+f(1)] = 0.9207355\)；
2. 第一次二分（\(n=2\)），新增中点\(x=1/2\)，递推得\(T_2 = \frac{1}{2}T_1 + \frac{1}{2}f(1/2) = 0.9397933\)；
3. 第二次二分（\(n=4\)），新增中点\(x=1/4,3/4\)，递推得\(T_4 = \frac{1}{2}T_2 + \frac{1}{4}[f(1/4)+f(3/4)] = 0.9445135\)；
4. 持续二分的结果如下表：

二分次数\(k\)	等分数\(n=2^k\)	梯形值\(T_n\)
1	2	0.9397933
2	4	0.9445135
3	8	0.9456909
4	16	0.9459850
5	32	0.9460586
10	1024	0.9460830

可以看到，要达到7位有效数字，复合梯形公式需要二分10次、计算1025个函数值，收敛速度极慢，因此必须引入外推加速技术。

---

二、核心技术：理查森（Richardson）外推算法

外推算法的核心思想是：利用不同步长下的近似值，通过线性组合消去误差展开式中的低阶项，直接得到更高阶精度的近似值，无需额外计算函数值，仅通过代数运算就能提升精度。

1. 梯形值的误差展开定理

外推的前提是近似值的误差可表示为步长\(h\)的幂级数，对于复合梯形公式，有如下定理：

定理4.4 若\(f(x) \in C^\infty[a,b]\)（无穷阶连续可导），则步长\(h\)对应的复合梯形值\(T(h)\)可展开为：

\[T(h) = I + \alpha_1 h^2 + \alpha_2 h^4 + \alpha_3 h^6 + \dots + \alpha_l h^{2l} + \dots \tag{4.18} \]

其中系数\(\alpha_1,\alpha_2,\dots\)与步长\(h\)无关，仅和被积函数\(f(x)\)有关。

这个定理是外推的核心基础：梯形值的误差仅包含\(h\)的偶次幂项，我们可以通过不同步长的\(T(h)\)和\(T(h/2)\)，逐步消去低阶误差项。

2. 逐次外推过程

第一次外推：梯形值→辛普森值（消去\(h^2\)项）

我们有步长\(h\)和\(h/2\)对应的梯形值展开式：

\[T(h) = I + \alpha_1 h^2 + \alpha_2 h^4 + \dots \]

\[T\left(\frac{h}{2}\right) = I + \alpha_1 \cdot \frac{h^2}{4} + \alpha_2 \cdot \frac{h^4}{16} + \dots \]

用\(4\times T(h/2) - T(h)\)消去\(h^2\)项，两边除以3，记结果为\(S(h)\)：

\[S(h) = \frac{4T(h/2) - T(h)}{3} = I + \beta_1 h^4 + \beta_2 h^6 + \dots \tag{4.20} \]

此时\(S(h)\)的误差阶从\(O(h^2)\)提升到\(O(h^4)\)，结果与复合辛普森公式完全一致，仅通过线性运算就实现了精度的跨越式提升。

第二次外推：辛普森值→柯特斯值（消去\(h^4\)项）

同理，用\(S(h)\)和\(S(h/2)\)继续外推，消去\(h^4\)项，得到柯特斯值\(C(h)\)：

\[C(h) = \frac{16S(h/2) - S(h)}{15} = I + \gamma_1 h^6 + \gamma_2 h^8 + \dots \tag{4.22} \]

误差阶提升至\(O(h^6)\)，与复合柯特斯公式结果一致。

第三次外推：柯特斯值→龙贝格值（消去\(h^6\)项）

继续外推得到龙贝格值\(R(h)\)：

\[R(h) = \frac{64C(h/2) - C(h)}{63} = I + \delta_1 h^8 + \dots \tag{4.23} \]

误差阶达到\(O(h^8)\)，精度已极高。

外推的通用规律

第\(m\)次外推的统一公式为：

\[T_m(h) = \frac{4^m T_{m-1}(h/2) - T_{m-1}(h)}{4^m - 1} \]

其中\(T_m(h)\)表示\(m\)次外推的结果，每外推一次，误差阶提升2个量级。

---

三、龙贝格（Romberg）求积算法

将上述外推过程标准化，就得到了完整的龙贝格算法，它是一个递推式的表格计算方法，便于编程实现。

1. 离散化递推公式

在实际计算中，用二分次数\(k\)标记步长\(h_k = \frac{b-a}{2^k}\)，记：

• \(T_0^{(k)}\)：二分\(k\)次后的梯形值（0次外推）；
• \(T_m^{(k)}\)：基于二分\(k\)次的结果，做\(m\)次外推得到的值。

则离散化的龙贝格递推公式为：

\[T_m^{(k)} = \frac{4^m T_{m-1}^{(k+1)} - T_{m-1}^{(k)}}{4^m - 1}, \quad k=1,2,\dots; \ m=1,2,\dots \tag{4.25} \]

2. 龙贝格T表与计算步骤

龙贝格算法通过构造T表实现，T表的每一行对应一次区间二分，每一列对应一次外推加速，对角线元素是精度最高的结果。T表结构如下：

二分次数\(k\)	步长\(h_k\)	0次外推\(T_0^{(k)}\)（梯形值）	1次外推\(T_1^{(k)}\)（辛普森值）	2次外推\(T_2^{(k)}\)（柯特斯值）	3次外推\(T_3^{(k)}\)（龙贝格值）
0	\(b-a\)	\(T_0^{(0)}\)
1	\(\frac{b-a}{2}\)	\(T_0^{(1)}\)	\(T_1^{(0)}\)
2	\(\frac{b-a}{4}\)	\(T_0^{(2)}\)	\(T_1^{(1)}\)	\(T_2^{(0)}\)
3	\(\frac{b-a}{8}\)	\(T_0^{(3)}\)	\(T_1^{(2)}\)	\(T_2^{(1)}\)	\(T_3^{(0)}\)

具体计算步骤：

1. 初始计算：取步长\(h=b-a\)，计算初始梯形值\(T_0^{(0)} = \frac{h}{2}[f(a)+f(b)]\)，记二分次数\(k=0\)；
2. 区间二分与梯形值递推：将区间二分，用递推公式(4.17)计算新的梯形值\(T_0^{(k+1)}\)，填入T表第\(k+1\)行第1列；
3. 逐列外推加速：用递推公式(4.25)，从左到右计算第\(k+1\)行的所有外推值，填满T表的上三角部分；
4. 精度判断：检查对角线相邻两个值的绝对误差\(|T_m^{(0)} - T_{m-1}^{(0)}|\)，若小于预设精度\(\varepsilon\)，则终止计算，取\(T_m^{(0)}\)作为积分近似值；否则令\(k=k+1\)，返回步骤2继续二分。

3. 算法收敛性

• 若\(f(x)\)充分光滑，T表的每一列、对角线元素都收敛到积分真值\(I\)，且对角线元素的收敛速度远快于列元素；
• 即使\(f(x)\)光滑性较差，龙贝格算法依然收敛，仅收敛速度变慢，仍优于单纯的复合梯形公式。

---

四、算法实例（例4.6）

用龙贝格算法计算积分\(I=\int_0^1 x^{3/2}dx\)，被积函数\(f(x)=x^{3/2}\)在\([0,1]\)上仅一阶连续可导（光滑性较差），积分精确值为\(I=0.4\)。

龙贝格T表计算结果

二分次数\(k\)	\(T_0^{(k)}\)	\(T_1^{(k)}\)	\(T_2^{(k)}\)	\(T_3^{(k)}\)	\(T_4^{(k)}\)	\(T_5^{(k)}\)
0	0.500000
1	0.426777	0.402369
2	0.407018	0.400432	0.400303
3	0.401812	0.400077	0.400054	0.400050
4	0.400463	0.400014	0.400009	0.400009	0.400009
5	0.400118	0.400002	0.400002	0.400002	0.400002	0.400002

结果分析

即使被积函数光滑性较差，龙贝格算法仅通过5次二分（33个节点）就得到了6位有效数字的结果，收敛速度远快于复合梯形公式，充分体现了外推加速的优势。

---

五、龙贝格算法的核心优势与适用场景

核心优势

1. 精度高、收敛快：通过外推技术，误差阶从\(O(h^2)\)快速提升到\(O(h^8)\)甚至更高，用极少的函数值就能达到极高精度；
2. 计算效率高：二分过程复用所有已计算的函数值，无重复计算，外推仅为简单线性运算，计算成本极低；
3. 自适应精度控制：可通过对角线误差自动终止计算，无需提前确定区间等分数，便于编程实现自适应算法；
4. 稳定性好：所有外推系数均为正数，不会放大输入误差，计算过程绝对稳定；
5. 适用性广：无论被积函数是否充分光滑都能收敛，光滑性越好，收敛速度越快。

适用场景

• 工程中高精度数值积分计算，尤其是被积函数求值成本较高的场景（如实验测量数据、复杂函数计算）；
• 通用数值积分程序的核心算法，是绝大多数科学计算软件中数值积分的基础实现；
• 对积分精度有自适应要求的场景，可自动调整二分次数满足精度需求。

注意事项

• 龙贝格算法要求被积函数在积分区间上可积，且能方便地计算任意节点的函数值；
• 对于有奇点、振荡剧烈的被积函数，需要先做奇点处理或变量替换，再使用龙贝格算法。