本篇文章适合个人复习翻阅，不建议新手入门使用
本专栏：数值分析复习 的前置知识主要有：数学分析、高等代数、泛函分析
本节继续考虑数值积分问题
Richardson外推
命题：复合梯形公式的另一形式
设 f ∈ C ∞ [ a , b ] f\in C^{\infty}[a,b] f∈C∞[a,b]，记 I = ∫ a b f ( x ) d x I=\int_a^bf(x)\mathrm{d}x I=∫abf(x)dx ，将复合梯形公式记为
T ( h ) = h 2 ∑ i = 0 n − 1 [ f ( x i ) + f ( x i + 1 ) ] T(h)=\frac{h}{2}\sum\limits_{i=0}^{n-1}[f(x_i)+f(x_{i+1})] T(h)=2hi=0∑n−1[f(xi)+f(xi+1)]
则 T ( h ) = I + α 1 h 2 + α 2 h 4 + ⋯ + α l h 2 l + ⋯ T(h)=I+\alpha_1h^2+\alpha_2h^4+\cdots+\alpha_lh^{2l}+\cdots T(h)=I+α1h2+α2h4+⋯+αlh2l+⋯
其中 α l ( l = 1 , 2 , … ) \alpha_l(l=1,2,\dots) αl(l=1,2,…) 与 h h h 无关
证明
设 x i + 1 2 = x i + x i + 1 2 , i = 0 , 1 , … , n − 1 x_{i+\frac{1}{2}}=\frac{x_i+x_{i+1}}{2},i=0,1,\dots,n-1 xi+21=2xi+xi+1,i=0,1,…,n−1
考虑 f ( x ) f(x) f(x) 在 x = x i + 1 2 x=x_{i+\frac{1}{2}} x=xi+21 处的Taylor展开公式
f ( x ) = f ( x i + 1 2 ) + f ′ ( x i + 1 2 ) ( x − x i + 1 2 ) + f ′ ′ ( x i + 1 2 ) 2 ! ( x − x i + 1 2 ) 2 + ⋯ f(x)=f(x_{i+\frac{1}{2}})+f'(x_{i+\frac{1}{2}})(x-x_{i+\frac{1}{2}})+\frac{f''(x_{i+\frac{1}{2}})}{2!}(x-x_{i+\frac{1}{2}})^2+\cdots f(x)=f(xi+21)+f′(xi+21)(x−xi+21)+2!f′′(xi+21)(x−xi+21)2+⋯
若对上述 Taylor 公式代入 x = x i , x = x i + 1 x=x_{i},x=x_{i+1} x=xi,x=xi+1，则得
f ( x i + 1 ) = f ( x i + 1 2 ) + f ′ ( x i + 1 2 ) h 2 + f ′ ′ ( x i + 1 2 ) 2 ! ( h 2 ) 2 + ⋯ f(x_{i+1})=f(x_{i+\frac{1}{2}})+f'(x_{i+\frac{1}{2}})\frac{h}{2}+\frac{f''(x_{i+\frac{1}{2}})}{2!}(\frac{h}{2})^2+\cdots f(xi+1)=f(xi+21)+f′(xi+21)2h+2!f′′(xi+21)(2h)2+⋯ f ( x i ) = f ( x i + 1 2 ) + f ′ ( x i + 1 2 ) ( − h 2 ) + f ′ ′ ( x i + 1 2 ) 2 ! ( − h 2 ) 2 + ⋯ f(x_i)=f(x_{i+\frac{1}{2}})+f'(x_{i+\frac{1}{2}})(-\frac{h}{2})+\frac{f''(x_{i+\frac{1}{2}})}{2!}(-\frac{h}{2})^2+\cdots f(xi)=f(xi+21)+f′(xi+21)(−2h)+2!f′′(xi+21)(−2h)2+⋯
两式加和，得到
f ( x i ) + f ( x i + 1 ) 2 = f ( x i + 1 2 ) + h 2 8 f ′ ′ ( x i + 1 2 ) + ⋯ \frac{f(x_i)+f(x_{i+1})}{2}=f(x_{i+\frac{1}{2}})+\frac{h^2}{8}f''(x_{i+\frac{1}{2}})+\cdots 2f(xi)+f(xi+1)=f(xi+21)+8h2f′′(xi+21)+⋯
等式两端求和，乘以 h h h 得到
T ( h ) = h ∑ i = 0 n − 1 f ( x i + 1 2 ) + h 3 8 ∑ i = 0 n − 1 f ′ ′ ( x i + 1 2 ) + ⋯ (1) T(h)=h\sum\limits_{i=0}^{n-1}f(x_{i+\frac{1}{2}})+\frac{h^3}{8}\sum\limits_{i=0}^{n-1}f''(x_{i+\frac{1}{2}})+\cdots\tag 1 T(h)=hi=0∑n−1f(xi+21)+8h3i=0∑n−1f′′(xi+21)+⋯(1)
另一方面，对Taylor公式从 x i x_i xi 到 x i + 1 x_{i+1} xi+1 进行积分，得到
∫ x i x i + 1 f ( x ) d x = h ⋅ f ( x i + 1 2 ) + f ′ ( x i + 1 2 ) 2 [ ( h 2 ) 2 − ( − h 2 ) 2 ] + f ′ ′ ( x i + 1 2 ) 6 [ ( h 2 ) 3 − ( − h 2 ) 3 ] + ⋯ \int_{x_i}^{x_{i+1}}f(x)\mathrm{d}x=h\cdot f(x_{i+\frac{1}{2}})+\frac{f'(x_{i+\frac{1}{2}})}{2}[(\frac{h}{2})^2-(-\frac{h}{2})^2]+\frac{f''(x_{i+\frac{1}{2}})}{6}[(\frac{h}{2})^3-(-\frac{h}{2})^3]+\cdots ∫xixi+1f(x)dx=h⋅f(xi+21)+2f′(xi+21)[(2h)2−(−2h)2]+6f′′(xi+21)[(2h)3−(−2h)3]+⋯
等式两端求和得
I = ∑ i = 0 n − 1 ∫ x i x i + 1 f ( x ) d x = h ∑ i = 0 n − 1 f ( x i + 1 2 ) + h 3 24 ∑ i = 0 n − 1 f ′ ′ ( x i + 1 2 ) + ⋯ (2) I=\sum\limits_{i=0}^{n-1}\int_{x_i}^{x_{i+1}}f(x)\mathrm{d}x=h\sum\limits_{i=0}^{n-1}f(x_{i+\frac{1}{2}}) +\frac{h^3}{24}\sum\limits_{i=0}^{n-1}f''(x_{i+\frac{1}{2}}) +\cdots\tag 2 I=i=0∑n−1∫xixi+1f(x)dx=hi=0∑n−1f(xi+21)+24h3i=0∑n−1f′′(xi+21)+⋯(2)
结合（1）（2）式，可得
T ( h ) = I + h 3 12 ∑ i = 0 n − 1 f ′ ′ ( x i + 1 2 ) + ⋯ (3) T(h)=I+\frac{h^3}{12}\sum\limits_{i=0}^{n-1}f''(x_{i+\frac{1}{2}})+\cdots\tag 3 T(h)=I+12h3i=0∑n−1f′′(xi+21)+⋯(3)
类似（2）式的推导，可得
∫ a b f ′ ′ ( x ) d x = h ∑ i = 0 n − 1 f ′ ′ ( x i + 1 2 ) + h 3 24 ∑ i = 0 n − 1 f ( 4 ) ( x i + 1 2 ) + ⋯ \int_a^bf''(x)\mathrm{d}x=h\sum\limits_{i=0}^{n-1}f''(x_{i+\frac{1}{2}})+\frac{h^3}{24}\sum\limits_{i=0}^{n-1}f^{(4)}(x_{i+\frac{1}{2}})+\cdots ∫abf′′(x)dx=hi=0∑n−1f′′(xi+21)+24h3i=0∑n−1f(4)(xi+21)+⋯
结合 ∫ a b f ′ ′ ( x ) d x = f ′ ( b ) − f ′ ( a ) \int_a^bf''(x)\mathrm{d}x=f'(b)-f'(a) ∫abf′′(x)dx=f′(b)−f′(a)，可将（3）式化为
T ( h ) = I + α 1 h 2 + h 5 c 4 ∑ i = 0 n − 1 f ( 4 ) ( x i + 1 2 ) + ⋯ T(h)=I+\alpha_1h^2+h^5c_4\sum\limits_{i=0}^{n-1}f^{(4)}(x_{i+\frac{1}{2}})+\cdots T(h)=I+α1h2+h5c4i=0∑n−1f(4)(xi+21)+⋯
重复上述操作，考虑 ∫ a b f ( 4 ) ( x ) d x \int_a^bf^{(4)}(x)\mathrm{d}x ∫abf(4)(x)dx，消去 h 5 h^5 h5 的项，得到 h 4 h^4 h4 的项，继续重复操作，可得
T ( h ) = I + α 1 h 2 + α 2 h 4 + ⋯ + α l h 2 l + ⋯ T(h)=I+\alpha_1h^2+\alpha_2h^4+\cdots+\alpha_lh^{2l}+\cdots T(h)=I+α1h2+α2h4+⋯+αlh2l+⋯
定义：Richardson外推
从低阶精度格式的截断误差的渐近展开式出发，做简单线性计算从而得到高阶精度格式的方法称为Richardson外推
例：
考虑复合梯形公式 T ( h ) T(h) T(h) 满足的式子
T ( h ) = I + α 1 h 2 + α 2 h 4 + ⋯ + α l h 2 l + ⋯ T(h)=I+\alpha_1h^2+\alpha_2h^4+\cdots+\alpha_lh^{2l}+\cdots T(h)=I+α1h2+α2h4+⋯+αlh2l+⋯
此时截断误差量级为 O ( h 2 ) O(h^{2}) O(h2)
取步长为 h 2 \frac{h}{2} 2h，则有
T ( h 2 ) = I + α 1 h 2 4 + α 2 h 4 16 + ⋯ + α l h 2 l 2 2 l + ⋯ T(\frac{h}{2})=I+\alpha_1\frac{h^2}{4}+\alpha_2\frac{h^4}{16}+\cdots+\alpha_l\frac{h^{2l}}{2^{2l}}+\cdots T(2h)=I+α14h2+α216h4+⋯+αl22lh2l+⋯
结合这两个式子，消去 h 2 h^{2} h2项，得
4 T ( h 2 ) − T ( h ) 3 = I − 1 4 α 2 h 4 + ⋯ + α l 3 ( 1 2 2 l − 1 ) h 2 l + ⋯ \frac{4T(\frac{h}{2})-T(h)}{3}=I-\frac{1}{4}\alpha_2h^4+\cdots+\frac{\alpha_l}{3}(\frac{1}{2^{2l}}-1)h^{2l}+\cdots 34T(2h)−T(h)=I−41α2h4+⋯+3αl(22l1−1)h2l+⋯
记 T 1 ( h ) = 4 T ( h 2 ) − T ( h ) 3 T_1(h)=\frac{4T(\frac{h}{2})-T(h)}{3} T1(h)=34T(2h)−T(h)，且
T 1 ( h ) = I + β 2 h 4 + β 3 h 6 + ⋯ + β l h 2 l + ⋯ T_1(h)=I+\beta_2h^4+\beta_3h^6+\cdots+\beta^lh^{2l}+\cdots T1(h)=I+β2h4+β3h6+⋯+βlh2l+⋯
若用 T 1 ( h ) T_1(h) T1(h) 估计 I I I ，则截断误差量级提高到 O ( h 4 ) O(h^{4}) O(h4)
类似地，可继续做……
注：只要截断误差可表示为 h h h 的幂级数，均可使用 Richardson外推提高精度
Romberg（龙贝格）算法
在上述对复合梯形公式的截断误差进行Richardson外推的过程中，记复合梯形公式 T 0 ( h ) = T ( h ) = h 2 ∑ i = 0 n − 1 [ f ( x i ) + f ( x i + 1 ) ] T_0(h)=T(h)=\frac{h}{2}\sum\limits_{i=0}^{n-1}[f(x_i)+f(x_{i+1})] T0(h)=T(h)=2hi=0∑n−1[f(xi)+f(xi+1)]
加速一次（即进行一次Richardson外推）后的估计式记为
T 1 ( h ) = 4 T ( h 2 ) − T ( h ) 3 T_1(h)=\frac{4T(\frac{h}{2})-T(h)}{3} T1(h)=34T(2h)−T(h)
记加速 n n n 次的估计式为 T n ( h ) T_n(h) Tn(h)，则有递推式
T n ( h ) = 4 n 4 n − 1 T n − 1 ( h 2 ) − 1 4 n − 1 T n − 1 ( h ) T_n(h)=\frac{4^n}{4^n-1}T_{n-1}(\frac{h}{2})-\frac{1}{4^n-1}T_{n-1}(h) Tn(h)=4n−14nTn−1(2h)−4n−11Tn−1(h)
若记 T m ( k ) = T m ( h 2 k ) , k = 0 , 1 , 2 , … T_m^{(k)}=T_m(\frac{h}{2^k}),k=0,1,2,\dots Tm(k)=Tm(2kh),k=0,1,2,…，则有递推式
T n ( k ) = 4 n 4 n − 1 T n − 1 ( k + 1 ) − 1 4 n − 1 T n − 1 ( k ) T_n^{(k)}=\frac{4^n}{4^n-1}T_{n-1}^{(k+1)}-\frac{1}{4^n-1}T_{n-1}^{(k)} Tn(k)=4n−14nTn−1(k+1)−4n−11Tn−1(k)
定理：
设被积函数 f ( x ) f(x) f(x) 充分光滑
- lim k → ∞ T m ( k ) = I \lim\limits_{k\to\infty}T_m^{(k)}=I k→∞limTm(k)=I
- lim m → ∞ T m ( k ) = I \lim\limits_{m\to\infty}T_m^{(k)}=I m→∞limTm(k)=I
注：证明略去，第一个结论说明当节点数目无穷多时， T m ( k ) T_m^{(k)} Tm(k) 收敛于准确的积分值；第二个结论说明随着Richardson外推的进行， T m ( k ) T_m^{(k)} Tm(k) 也收敛于准确的积分值
上述递推式和收敛定理给出了如下的Romberg算法
定义：Romberg算法
对预先给定的精度 ε \varepsilon ε，求 I = ∫ a b f ( x ) d x I=\int_a^bf(x)\mathrm{d}x I=∫abf(x)dx 的近似值，算法如下：
初始取 k = 0 , m = 0 , h = b − a k=0,m=0,h=b-a k=0,m=0,h=b−a
- 代入梯形公式，求 T 0 ( k ) ( k = 0 , 1 , 2 , … ) T_0^{(k)}(k=0,1,2,\dots) T0(k)(k=0,1,2,…)
- 加速一次，由递推公式求 T 1 ( k ) T_1^{(k)} T1(k)
- 直至 ∣ T k ( 0 ) − T k − 1 ( 0 ) ∣ < ε |T_k^{(0)}-T_{k-1}^{(0)}|<\varepsilon ∣Tk(0)−Tk−1(0)∣<ε，则取 T k ( 0 ) ≈ I T_{k}^{(0)}\approx I Tk(0)≈I
注：具体求解顺序如下表
参考书籍：《数值分析》李庆扬 王能超 易大义 编