文章目录
- Richardson外推
- Romberg（龙贝格）算法
本篇文章适合个人复习翻阅，不建议新手入门使用
本专栏：数值分析复习 的前置知识主要有：数学分析、高等代数、泛函分析
本节继续考虑数值积分问题
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
注：具体求解顺序如下表
参考书籍：《数值分析》李庆扬 王能超 易大义 编
