数值分析复习：Richardson外推和Romberg算法

本篇文章适合个人复习翻阅，不建议新手入门使用
本专栏：数值分析复习 的前置知识主要有：数学分析、高等代数、泛函分析

本节继续考虑数值积分问题

Richardson外推

命题：复合梯形公式的另一形式
设 $f\in C^{\infty }[a,b]$，记 $I={\int }_a^{b}f(x)\mathrm{d}x$ ，将复合梯形公式记为
$$T(h)=\frac{h}{2}\sum _{i=0}^{n-1}[f(x_i)+f(x_{i+1})]$$

则 $$T(h)=I+{\alpha }_1{h}^2+{\alpha }_2{h}^4+\cdots +{\alpha }_l{h}^{2l}+\cdots$$

其中 ${\alpha }_l(l=1,2,\dots )$ 与 $h$ 无关

证明
设 $x_{i+\frac{1}{2}}=\frac{x_i+x_{i+1}}{2},i=0,1,\dots ,n-1$

考虑 $f(x)$ 在 $x=x_{i+\frac{1}{2}}$ 处的Taylor展开公式
$$f(x)=f(x_{i+\frac{1}{2}})+f^{\prime }(x_{i+\frac{1}{2}})(x-x_{i+\frac{1}{2}})+\frac{f^{\prime \prime }(x_{i+\frac{1}{2}})}{2!}(x-x_{i+\frac{1}{2}}{)}^2+\cdots$$

若对上述 Taylor 公式代入 $x=x_i,x=x_{i+1}$，则得
$$f(x_{i+1})=f(x_{i+\frac{1}{2}})+f^{\prime }(x_{i+\frac{1}{2}})\frac{h}{2}+\frac{f^{\prime \prime }(x_{i+\frac{1}{2}})}{2!}(\frac{h}{2}{)}^2+\cdots$$$$f(x_i)=f(x_{i+\frac{1}{2}})+f^{\prime }(x_{i+\frac{1}{2}})(-\frac{h}{2})+\frac{f^{\prime \prime }(x_{i+\frac{1}{2}})}{2!}(-\frac{h}{2}{)}^2+\cdots$$

两式加和，得到
$$\frac{f(x_i)+f(x_{i+1})}{2}=f(x_{i+\frac{1}{2}})+\frac{h^2}{8}{f}^{\prime \prime }(x_{i+\frac{1}{2}})+\cdots$$

等式两端求和，乘以 $h$ 得到
$$T(h)=h\sum _{i=0}^{n-1}f(x_{i+\frac{1}{2}})+\frac{h^3}{8}\sum _{i=0}^{n-1}{f}^{\prime \prime }(x_{i+\frac{1}{2}})+\cdots \text{\hspace{1em}(1)}$$

另一方面，对Taylor公式从 $x_i$ 到 $x_{i+1}$ 进行积分，得到
$${\int }_{x_i}^{x_{i+1}}f(x)\mathrm{d}x=h\cdot f(x_{i+\frac{1}{2}})+\frac{f^{\prime }(x_{i+\frac{1}{2}})}{2}[(\frac{h}{2}{)}^2-(-\frac{h}{2}{)}^2]+\frac{f^{\prime \prime }(x_{i+\frac{1}{2}})}{6}[(\frac{h}{2}{)}^3-(-\frac{h}{2}{)}^3]+\cdots$$

等式两端求和得

$$I=\sum _{i=0}^{n-1}{\int }_{x_i}^{x_{i+1}}f(x)\mathrm{d}x=h\sum _{i=0}^{n-1}f(x_{i+\frac{1}{2}})+\frac{h^3}{24}\sum _{i=0}^{n-1}{f}^{\prime \prime }(x_{i+\frac{1}{2}})+\cdots \text{\hspace{1em}(2)}$$

结合（1）（2）式，可得
$$T(h)=I+\frac{h^3}{12}\sum _{i=0}^{n-1}{f}^{\prime \prime }(x_{i+\frac{1}{2}})+\cdots \text{\hspace{1em}(3)}$$

类似（2）式的推导，可得
$${\int }_a^b{f}^{\prime \prime }(x)\mathrm{d}x=h\sum _{i=0}^{n-1}{f}^{\prime \prime }(x_{i+\frac{1}{2}})+\frac{h^3}{24}\sum _{i=0}^{n-1}{f}^{(4)}(x_{i+\frac{1}{2}})+\cdots$$

结合 ${\int }_a^b{f}^{\prime \prime }(x)\mathrm{d}x=f^{\prime }(b)-f^{\prime }(a)$，可将（3）式化为
$$T(h)=I+{\alpha }_1{h}^2+h^5{c}_4\sum _{i=0}^{n-1}{f}^{(4)}(x_{i+\frac{1}{2}})+\cdots$$

重复上述操作，考虑 ${\int }_a^b{f}^{(4)}(x)\mathrm{d}x$，消去 $h^5$ 的项，得到 $h^4$ 的项，继续重复操作，可得
$$T(h)=I+{\alpha }_1{h}^2+{\alpha }_2{h}^4+\cdots +{\alpha }_l{h}^{2l}+\cdots$$

定义：Richardson外推
从低阶精度格式的截断误差的渐近展开式出发，做简单线性计算从而得到高阶精度格式的方法称为Richardson外推

例：
考虑复合梯形公式 $T(h)$ 满足的式子
$$T(h)=I+{\alpha }_1{h}^2+{\alpha }_2{h}^4+\cdots +{\alpha }_l{h}^{2l}+\cdots$$

此时截断误差量级为 $O(h^2)$

取步长为 $\frac{h}{2}$，则有
$$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$$

结合这两个式子，消去$h^2$项，得
$$\frac{4T(\frac{h}{2})-T(h)}{3}=I-\frac{1}{4}{\alpha }_2{h}^4+\cdots +\frac{{\alpha }_l}{3}(\frac{1}{2^{2l}}-1)h^{2l}+\cdots$$
记 $T_1(h)=\frac{4T(\frac{h}{2})-T(h)}{3}$，且
$$T_1(h)=I+{\beta }_2{h}^4+{\beta }_3{h}^6+\cdots +{\beta }^l{h}^{2l}+\cdots$$
若用 $T_1(h)$ 估计 $I$ ，则截断误差量级提高到 $O(h^4)$
类似地，可继续做……

注：只要截断误差可表示为 $h$ 的幂级数，均可使用 Richardson外推提高精度

Romberg（龙贝格）算法

Romberg算法用以解决问题：对预先给定的精度 $\epsilon$，求 $I={\int }_a^{b}f(x)\mathrm{d}x$ 的近似值

记梯形公式为
$$T(1,1)=\frac{(b-a)}{2}[f(a)+f(b)]$$

设对区间 $[a,b]$ 作 $m$ 次等分得到的复合梯形公式为 $T(m+1,1)$
设对 $T(s,t)$ 进行一次 Richardson外推得到的积分为 $T(s+1,t+1)$

由Richardson外推过程，可归纳地证明
$$T(s+1,t+1)=\frac{4^{n}T(s+1,t)-T(s,t)}{4^n-1}$$

故可设计算法依次计算如下表格

由 Euler-Maclaurin 公式，可保证当 $t\to \infty$ 时，$T(t,t)$ 收敛于积分的真实值 $I$

随着节点加密（即 $s\to \infty$），$T(s,t)$ 也收敛于 $I$

故对预先设置的误差限 $\epsilon$，用
max ⁡ { ∣ T ( s + 1 , s + 1 ) − T ( s , s ) ∣ , ∣ T ( s + 1 , s + 1 ) − T ( s + 1 , s ) ∣ } < ε \max\{|T(s+1,s+1)-T(s,s)|,|T(s+1,s+1)-T(s+1,s)|\}<\varepsilon max{∣T(s+1,s+1)−T(s,s)∣,∣T(s+1,s+1)−T(s+1,s)∣}<ε

代替
$$|T(s,s)-I|<\epsilon$$

此时 $T(s+1,s+1)$ 即为 $I$ 的近似值

matlab代码如下，其中节点从1开始计数 $x_1,x_2,\dots ,x_{n+1}$

%程序：龙贝格算法求数值积分
%输入：被积函数f,积分上限b,下限a,要求的误差限tol
%输出：T数表，误差估计err
%     用到的梯形公式函数值个数cnt（即节点个数）
%     数值积分结果ans
function [q,cnt] = romberg(f,a,b,tol)

% 计算初始梯形公式T(1,1)
n=1;        %节点个数 $n+1$
h = b-a;    %初始步长
fv = [f(a);f(b)];     %初始的梯形公式所需的函数值
T(1,1) = .5*h*(fv(1)+fv(2));
cnt = 2;         % 返回函数值的个数

%计算梯形公式
err = 100*abs(tol); k=0;    % 初始化 err 使其大于误差，启动循环
while err > tol &k < 10
    k = k+1; n=2*n; h=h/2;  % 节点数翻倍，步长减半
    fvnew = zeros(n+1,1);   % 初始化，预备储存位置
    fvnew(1:2:n+1) = fv(1:n/2+1);%定义边界点x1,x(n+1)处的函数值
    for i =2:2:n           %定义其余节点处的函数值
        fvnew(i)= f(a+(i-1)*h);
    end
    cnt = cnt +(n/2);       %计算函数值个数，即节点个数
    fv = fvnew;
    %计算新的梯形公式
    trap = .5*(fv(1)+fv(n+1));
    for i= 2:n            % 估计积分
        trap = trap + fv(i);
    end
    T(k+1,1)=h*trap;

    % 作Richardson外推
    for j=2:k+1
        T(k+1,j) = ((4^(j-1))*T(k+1,j-1)-T(k,j-1))/(4^(j-1)-1);
    end
    q = T(k+1,k+1);

    % 估计误差
    err =max([abs(q-T(k,k));abs(q-T(k+1,k))]);
end
disp(T(1:k+1,1:k+1));
disp(['err = ' num2str(err, '%g')]);
disp(['cnt = ' num2str(cnt)]);

参考书籍：《数值分析》李庆扬 王能超 易大义 编