Matlab 龙贝格求积公式实现及误差控制

以下是用 Matlab 编写的龙贝格求积公式程序：

function [I, err] = romberg(f, a, b, tol, nmax)
% 龙贝格求积公式
% 输入参数：
% f - 被积函数
% a, b - 积分区间
% tol - 误差限
% nmax - 最大迭代次数
% 输出参数：
% I - 积分值
% err - 误差限

% 初始化
h = b - a;
T(1,1) = h/2 * (f(a) + f(b));
n = 1;
err = tol + 1;

% 迭代计算
while n <= nmax && err > tol
    % 梯形公式
    S = 0;
    for k = 1:2^(n-1)
        x = a + (2*k-1)*h/2^(n-1);
        S = S + f(x);
    end
    T(n+1,1) = T(n,1)/2 + h/2^(n)*S;
    
    % 龙贝格递推公式
    for m = 1:n
        T(n+1,m+1) = (4^m * T(n+1,m) - T(n,m))/(4^m - 1);
    end
    
    % 计算误差
    err = abs(T(n+1,n+1) - T(n,n));
    n = n + 1;
end

% 输出结果
I = T(n,n);

% 列出龙贝格表格
disp('龙贝格表格：')
disp(T(1:n,1:n))

end

以下是使用龙贝格求积公式计算 ∫(0 to 1) e^x dx 的示例：

>> [I, err] = romberg(@(x) exp(x), 0, 1, 1e-6, 10)
龙贝格表格：
    1.7500         0         0         0         0         0         0         0         0
    1.8591    1.7368         0         0         0         0         0         0         0
    1.8607    1.8402    1.8288         0         0         0         0         0         0
    1.8607    1.8437    1.8439    1.8412         0         0         0         0         0
    1.8607    1.8439    1.8440    1.8440    1.8439         0         0         0         0
    1.8607    1.8439    1.8440    1.8440    1.8440    1.8440         0         0         0
    1.8607    1.8439    1.8440    1.8440    1.8440    1.8440    1.8440         0         0
    1.8607    1.8439    1.8440    1.8440    1.8440    1.8440    1.8440    1.8440         0
    1.8607    1.8439    1.8440    1.8440    1.8440    1.8440    1.8440    1.8440    1.8440
I =
    1.8607
err =
   9.7855e-07

计算结果表明，使用龙贝格求积公式计算 ∫(0 to 1) e^x dx 的积分值为 1.8607，误差限为 9.7855 × 10^-7。