理查森外推法——龙贝格积分

理查森外推法

        根据积分估计值本身，可以构造出方法来改进数值积分的结果。这些方法一般被称为理查森外推法。

以复合梯形法则为例进行推导

复合梯形法则的估计值和误差可以表示为：

其中I为积分精确值，I(h)为积分近似值，E为截断误差。用两个步长计算结果为：

其中E的表达式为：

假设二阶导为与步长无关的常数：

代入可得：

当h2 = h1/2时:

龙贝格积分公式

Matlab实现

function [q, ea, iter] = romberg(func, a, b, es, maxit, varargin)
% 应用理查森外推法计算积分
%input:
% func；被积函数
% a, b：积分上下限
% es：期望误差
% maxit：最大迭代步数
%output:
% q：积分结果
% ea：实际误差
% iter：实际迭代步数
if nargin<4||isempty(es), es = 0.000001; end
if nargin<5||isempty(maxit), maxit = 50; end
n = 1;
I(1, 1) = trap(func, a, b, n, varargin{:});
iter = 0;
while iter<maxit
    iter = iter+1;
    n = 2^iter;
    I(iter+1, 1) = trap(func, a, b, n, varargin{:});
    for k = 2:iter+1
        j = 2+iter-k;
        I(j, k) = (4^(k-1)*I(j+1,k-1) - I(j,k-1)) / (4^(k-1)-1);
    end
    ea = abs((I(1,iter+1)-I(2,iter)) / I(1,iter+1)) * 100;
    if ea<=es, break; end
end
q = I(1, iter+1);
end