如果知道误差怎样随步长变化，就能把两个近似组合起来，消去首个误差项。Richardson 外推做这件事；Romberg 把同一操作反复应用于复合梯形公式。
本章以 lecture note 第 4.5–4.6、5.10–5.11 节为主要依据。微分时目标值可以是 f ′ ( x ) f'(x) f ′ ( x ) I I I R k , j R_{k,j} R k , j k k k j j j [1] [1] S. Rojas, “Lecture Notes on Computational Mathematics,” 2025. Course lecture note distributed with MTH2051; local source course-lecture-notes.pdf., [2] [2] M. U. School of Mathematics, “MTH2051: Introduction to Computational Mathematics: Course Lecture Notes and Study Workbooks,” 2026. Monash University course materials and study workbooks, Semester 2, 2026., [3] [3] R. L. Burden and J. D. Faires, Numerical Analysis , 9th ed. Brooks/Cole, Cengage Learning, 2011.
Richardson：消去已知幂次
定理 一般外推公式
假设 p > 0 p>0 p > 0 q > p q>p q > p
Q ( h ) = L + c h p + O ( h q ) , Q(h)=L+ch^p+O(h^q), Q ( h ) = L + c h p + O ( h q ) , 其中 c c c h h h
R ( h ) = 2 p Q ( h / 2 ) − Q ( h ) 2 p − 1 . R(h)=\frac{2^pQ(h/2)-Q(h)}{2^p-1}. R ( h ) = 2 p − 1 2 p Q ( h /2 ) − Q ( h ) . 则 R ( h ) = L + O ( h q ) R(h)=L+O(h^q) R ( h ) = L + O ( h q )
证明
在 h / 2 h/2 h /2 c h p / 2 p ch^p/2^p c h p / 2 p 2 p 2^p 2 p 2 p − 1 2^p-1 2 p − 1 O ( h q ) O(h^q) O ( h q )
这不是“自动提升两阶”的定理。单侧差分一般含相邻幂次，q = p + 1 q=p+1 q = p + 1 q = p + 2 q=p+2 q = p + 2
例如 f ∈ C 5 f\in C^5 f ∈ C 5
D c e n ( x ; h ) = f ′ ( x ) + h 2 6 f ′ ′ ′ ( x ) + O ( h 4 ) . D_{\mathrm{cen}}(x;h)=f'(x)+\frac{h^2}{6}f'''(x)+O(h^4). D cen ( x ; h ) = f ′ ( x ) + 6 h 2 f ′′′ ( x ) + O ( h 4 ) . 于是
4 D c e n ( x ; h / 2 ) − D c e n ( x ; h ) 3 = f ′ ( x ) + O ( h 4 ) . \frac{4D_{\mathrm{cen}}(x;h/2)-D_{\mathrm{cen}}(x;h)}3=f'(x)+O(h^4). 3 4 D cen ( x ; h /2 ) − D cen ( x ; h ) = f ′ ( x ) + O ( h 4 ) . 单侧前向差分在 f ∈ C 3 f\in C^3 f ∈ C 3 f ′ ( x ) + h f ′ ′ ( x ) / 2 + O ( h 2 ) f'(x)+hf''(x)/2+O(h^2) f ′ ( x ) + h f ′′ ( x ) /2 + O ( h 2 ) 2 D f w d ( x ; h / 2 ) − D f w d ( x ; h ) 2D_{\mathrm{fwd}}(x;h/2)-D_{\mathrm{fwd}}(x;h) 2 D fwd ( x ; h /2 ) − D fwd ( x ; h )
改写外推公式，可看见精细近似的修正量：
R ( h ) = Q ( h / 2 ) + Q ( h / 2 ) − Q ( h ) 2 p − 1 . R(h)=Q(h/2)+\frac{Q(h/2)-Q(h)}{2^p-1}. R ( h ) = Q ( h /2 ) + 2 p − 1 Q ( h /2 ) − Q ( h ) . 后一个量估计的是 L − Q ( h / 2 ) L-Q(h/2) L − Q ( h /2 )
微分的递归外推
讲义第 4.6 节将微分外推记为 R k ( h ) R_k(h) R k ( h ) R k , j R_{k,j} R k , j R 0 ( h ) = D ( h ) R_0(h)=D(h) R 0 ( h ) = D ( h )
R k ( h ) = 2 p + 2 ( k − 1 ) R k − 1 ( h / 2 ) − R k − 1 ( h ) 2 p + 2 ( k − 1 ) − 1 . R_k(h)=\frac{2^{p+2(k-1)}R_{k-1}(h/2)-R_{k-1}(h)}{2^{p+2(k-1)}-1}. R k ( h ) = 2 p + 2 ( k − 1 ) − 1 2 p + 2 ( k − 1 ) R k − 1 ( h /2 ) − R k − 1 ( h ) . 命题 有限偶次展开下的递归精度
设固定整数 K ≥ 1 K\geq1 K ≥ 1 p > 0 p>0 p > 0
D ( h ) = f ′ ( x ) + ∑ r = 0 K − 1 c r h p + 2 r + O ( h p + 2 K ) . D(h)=f'(x)+\sum_{r=0}^{K-1}c_rh^{p+2r}+O(h^{p+2K}). D ( h ) = f ′ ( x ) + r = 0 ∑ K − 1 c r h p + 2 r + O ( h p + 2 K ) . 系数均不随 h h h 0 ≤ k ≤ K 0\leq k\leq K 0 ≤ k ≤ K R k ( h ) = f ′ ( x ) + O ( h p + 2 k ) R_k(h)=f'(x)+O(h^{p+2k}) R k ( h ) = f ′ ( x ) + O ( h p + 2 k )
证明
第 k k k h p + 2 r h^{p+2r} h p + 2 r ( 2 p + 2 ( k − 1 ) − 2 p + 2 r ) / ( 2 p + 2 r ( 2 p + 2 ( k − 1 ) − 1 ) ) (2^{p+2(k-1)}-2^{p+2r})/(2^{p+2r}(2^{p+2(k-1)}-1)) ( 2 p + 2 ( k − 1 ) − 2 p + 2 r ) / ( 2 p + 2 r ( 2 p + 2 ( k − 1 ) − 1 )) r = k − 1 r=k-1 r = k − 1 k k k
中心一阶差分取 p = 2 p=2 p = 2
梯形公式的偶次幂展开为何成立
不能把一个形式上的无限级数当作普遍定理。下面只证明有限项展开，足够支持任意固定列的 Romberg 分析。
定理 有限阶 Euler–Maclaurin 展开
取整数 m ≥ 1 m\geq1 m ≥ 1 f ∈ C 2 m + 2 [ a , b ] f\in C^{2m+2}[a,b] f ∈ C 2 m + 2 [ a , b ] h = ( b − a ) / n h=(b-a)/n h = ( b − a ) / n
T ( h ) = I + ∑ r = 1 m B 2 r ( 2 r ) ! [ f ( 2 r − 1 ) ( b ) − f ( 2 r − 1 ) ( a ) ] h 2 r + O ( h 2 m + 2 ) . T(h)=I+\sum_{r=1}^{m}
\frac{B_{2r}}{(2r)!}
\bigl[f^{(2r-1)}(b)-f^{(2r-1)}(a)\bigr]h^{2r}
+O(h^{2m+2}). T ( h ) = I + r = 1 ∑ m ( 2 r )! B 2 r [ f ( 2 r − 1 ) ( b ) − f ( 2 r − 1 ) ( a ) ] h 2 r + O ( h 2 m + 2 ) . 这里 B 2 r B_{2r} B 2 r B 2 = 1 / 6 B_2=1/6 B 2 = 1/6 B 4 = − 1 / 30 B_4=-1/30 B 4 = − 1/30
证明
定义 Bernoulli 多项式 B 0 ( t ) = 1 B_0(t)=1 B 0 ( t ) = 1 r ≥ 1 r\geq1 r ≥ 1
B r ′ ( t ) = r B r − 1 ( t ) , ∫ 0 1 B r ( t ) d t = 0. B_r'(t)=rB_{r-1}(t),\qquad \int_0^1 B_r(t)\,dt=0. B r ′ ( t ) = r B r − 1 ( t ) , ∫ 0 1 B r ( t ) d t = 0. 积分与归一化唯一确定它们。前几个是
B 1 ( t ) = t − 1 2 , B 2 ( t ) = t 2 − t + 1 6 , B_1(t)=t-\frac12,\quad B_2(t)=t^2-t+\frac16, B 1 ( t ) = t − 2 1 , B 2 ( t ) = t 2 − t + 6 1 , B 3 ( t ) = t 3 − 3 2 t 2 + 1 2 t , B 4 ( t ) = t 4 − 2 t 3 + t 2 − 1 30 . B_3(t)=t^3-\frac32t^2+\frac12t,\quad
B_4(t)=t^4-2t^3+t^2-\frac1{30}. B 3 ( t ) = t 3 − 2 3 t 2 + 2 1 t , B 4 ( t ) = t 4 − 2 t 3 + t 2 − 30 1 . 递推和零均值给出 B r ( 1 ) = B r ( 0 ) B_r(1)=B_r(0) B r ( 1 ) = B r ( 0 ) r ≥ 2 r\geq2 r ≥ 2 B r ( 1 − t ) = ( − 1 ) r B r ( t ) B_r(1-t)=(-1)^rB_r(t) B r ( 1 − t ) = ( − 1 ) r B r ( t ) r ≥ 3 r\geq3 r ≥ 3 B r = B r ( 0 ) B_r=B_r(0) B r = B r ( 0 )
将每个多项式按网格周期延拓，记为 B ~ r ( x ) = B r ( { ( x − a ) / h } ) \widetilde B_r(x)=B_r(\{(x-a)/h\}) B r ( x ) = B r ({( x − a ) / h })
T ( h ) − I = h ∫ a b B ~ 1 ( x ) f ′ ( x ) d x . T(h)-I=h\int_a^b\widetilde B_1(x)f'(x)\,dx. T ( h ) − I = h ∫ a b B 1 ( x ) f ′ ( x ) d x . 这一式也可直接验证：每块的边界贡献为两端函数值平均，积分项为块积分除以 h h h r ≥ 2 r\geq2 r ≥ 2 2 m 2m 2 m
T ( h ) − I = ∑ r = 1 m B 2 r h 2 r ( 2 r ) ! [ f ( 2 r − 1 ) ( b ) − f ( 2 r − 1 ) ( a ) ] − h 2 m ( 2 m ) ! ∫ a b B ~ 2 m ( x ) f ( 2 m ) ( x ) d x . T(h)-I=
\sum_{r=1}^{m}\frac{B_{2r}h^{2r}}{(2r)!}
\bigl[f^{(2r-1)}(b)-f^{(2r-1)}(a)\bigr]
-\frac{h^{2m}}{(2m)!}
\int_a^b\widetilde B_{2m}(x)f^{(2m)}(x)\,dx. T ( h ) − I = r = 1 ∑ m ( 2 r )! B 2 r h 2 r [ f ( 2 r − 1 ) ( b ) − f ( 2 r − 1 ) ( a ) ] − ( 2 m )! h 2 m ∫ a b B 2 m ( x ) f ( 2 m ) ( x ) d x . 再分部积分两次，最后的余项成为一个 h 2 m + 2 h^{2m+2} h 2 m + 2 h 2 m + 2 h^{2m+2} h 2 m + 2 O ( h 2 m + 2 ) O(h^{2m+2}) O ( h 2 m + 2 )
特别地，在 f ∈ C 6 f\in C^6 f ∈ C 6
T ( h ) = I + h 2 12 [ f ′ ( b ) − f ′ ( a ) ] − h 4 720 [ f ′ ′ ′ ( b ) − f ′ ′ ′ ( a ) ] + O ( h 6 ) . T(h)=I+\frac{h^2}{12}[f'(b)-f'(a)]
-\frac{h^4}{720}[f'''(b)-f'''(a)]+O(h^6). T ( h ) = I + 12 h 2 [ f ′ ( b ) − f ′ ( a )] − 720 h 4 [ f ′′′ ( b ) − f ′′′ ( a )] + O ( h 6 ) . 端点导数匹配可能让某些系数为零；导数奇异则可能让整个展开不适用。
Romberg 表：行负责采样，列负责外推
讲义允许任意初始均匀网格，记其步长为 h h h h = b − a h=b-a h = b − a
h k = b − a 2 k , R k , 0 = T ( h k ) , h_k=\frac{b-a}{2^k},\qquad R_{k,0}=T(h_k), h k = 2 k b − a , R k , 0 = T ( h k ) , R k , j = R k , j − 1 + R k , j − 1 − R k − 1 , j − 1 4 j − 1 , 1 ≤ j ≤ k . R_{k,j}=R_{k,j-1}
+\frac{R_{k,j-1}-R_{k-1,j-1}}{4^j-1},
\qquad 1\leq j\leq k. R k , j = R k , j − 1 + 4 j − 1 R k , j − 1 − R k − 1 , j − 1 , 1 ≤ j ≤ k . 定理 固定列的精度
对固定 j ≥ 0 j\geq0 j ≥ 0 f ∈ C 2 j + 2 [ a , b ] f\in C^{2j+2}[a,b] f ∈ C 2 j + 2 [ a , b ] k → ∞ k\to\infty k → ∞
R k , j = I + O ( h k 2 j + 2 ) . R_{k,j}=I+O(h_k^{2j+2}). R k , j = I + O ( h k 2 j + 2 ) . 证明
j = 0 j=0 j = 0 j ≥ 1 j\geq1 j ≥ 1 m = j m=j m = j h 2 , … , h 2 j h^2,\ldots,h^{2j} h 2 , … , h 2 j O ( h 2 j + 2 ) O(h^{2j+2}) O ( h 2 j + 2 ) h 2 r h^{2r} h 2 r 4 r 4^r 4 r s s s
4 s − 4 r 4 s − 1 \frac{4^s-4^r}{4^s-1} 4 s − 1 4 s − 4 r 在 r = s r=s r = s s = 1 , … , j s=1,\ldots,j s = 1 , … , j j j j
第一外推列正是细网格 Simpson 公式，由权重恒等式 可知，不只是两种方法精度恰好相同。
复用已有采样
网格减半时，旧节点保留，只计算新中点。对 k ≥ 1 k\geq1 k ≥ 1
R k , 0 = 1 2 R k − 1 , 0 + h k ∑ i = 1 2 k − 1 f ( a + ( 2 i − 1 ) h k ) . R_{k,0}=\frac12R_{k-1,0}
+h_k\sum_{i=1}^{2^{k-1}}f\bigl(a+(2i-1)h_k\bigr). R k , 0 = 2 1 R k − 1 , 0 + h k i = 1 ∑ 2 k − 1 f ( a + ( 2 i − 1 ) h k ) . 证明就是把细网格梯形和拆成旧节点与新节点两部分。完成第 k k k 2 k + 1 2^k+1 2 k + 1
算法 1 Romberg 表
Require: 函数 f f f [ a , b ] [a,b] [ a , b ] K K K
1: R 0 , 0 ← ( b − a ) ( f ( a ) + f ( b ) ) / 2 R_{0,0}\gets (b-a)(f(a)+f(b))/2 R 0 , 0 ← ( b − a ) ( f ( a ) + f ( b )) /2
2: for k ← 1 k\gets1 k ← 1 to K K K do
3: h ← ( b − a ) / 2 k h\gets(b-a)/2^k h ← ( b − a ) / 2 k
4: R k , 0 ← R k − 1 , 0 / 2 + h ∑ i = 1 2 k − 1 f ( a + ( 2 i − 1 ) h ) R_{k,0}\gets R_{k-1,0}/2+h\sum_{i=1}^{2^{k-1}}f(a+(2i-1)h) R k , 0 ← R k − 1 , 0 /2 + h ∑ i = 1 2 k − 1 f ( a + ( 2 i − 1 ) h )
5: for j ← 1 j\gets1 j ← 1 to k k k do
6: R k , j ← R k , j − 1 + ( R k , j − 1 − R k − 1 , j − 1 ) / ( 4 j − 1 ) R_{k,j}\gets R_{k,j-1}+(R_{k,j-1}-R_{k-1,j-1})/(4^j-1) R k , j ← R k , j − 1 + ( R k , j − 1 − R k − 1 , j − 1 ) / ( 4 j − 1 )
7: end for
8: end for
9: return R R R
交互实验与停止判断
这个交互图需要启用 JavaScript。
在指数例子中，查看第一列到对角线如何逼近 e − 1 e-1 e − 1 k k k 2 k 2^k 2 k j j j j j j
常用停止检查是相邻对角值的差是否小于给定尺度，但这仍是实践中的判断。还应设置最大行数，检查非有限值，并留意浮点精度下的停滞。外推使用负系数，可能放大输入噪声；即使数学上的截断项继续减小，计算误差也可能不再改善。
例如 ∫ 0 1 e x d x = e − 1 \int_0^1 e^x\,dx=e-1 ∫ 0 1 e x d x = e − 1 1.859140914 1.859140914 1.859140914 1.718861152 1.718861152 1.718861152 1.718282688 1.718282688 1.718282688 1.718281829 1.718281829 1.718281829 自适应 Simpson 。
参考文献
[1] S. Rojas, “Lecture Notes on Computational Mathematics,” 2025. Course lecture note distributed with MTH2051; local source course-lecture-notes.pdf. ↩ [2] M. U. School of Mathematics, “MTH2051: Introduction to Computational Mathematics: Course Lecture Notes and Study Workbooks,” 2026. Monash University course materials and study workbooks, Semester 2, 2026. ↩ [3] R. L. Burden and J. D. Faires, Numerical Analysis , 9th ed. Brooks/Cole, Cengage Learning, 2011. ↩
评论