Matlab 龙贝格求积公式：误差控制与表格生成

以下是使用 Matlab 编写的龙贝格求积公式和表格：

function [I, err] = romberg(f, a, b, nmax, tol)
% f: 被积函数
% a, b: 积分区间
% nmax: 最大迭代次数
% tol: 误差容限

% 初始化数组
R = zeros(nmax, nmax);
h = b - a;
R(1,1) = h/2 * (f(a) + f(b));

% 计算龙贝格递推公式
for i = 2:nmax
    h = h/2;
    sumf = 0;
    for k = 1:(2^(i-2))
        x = a + (2*k-1)*h;
        sumf = sumf + f(x);
    end
    R(i,1) = 1/2*R(i-1,1) + h*sumf;
    for j = 2:i
        R(i,j) = (4^(j-1)*R(i,j-1) - R(i-1,j-1))/(4^(j-1) - 1);
    end
    err = abs(R(i,i) - R(i-1,i-1));
    if err < tol
        break;
    end
end

% 返回积分值和误差
I = R(i,i);
err = err;

% 输出龙贝格表格
for k = 1:i
    fprintf('R(%d,1) = %f\n', k, R(k,1));
    for j = 2:k
        fprintf('R(%d,%d) = %f\n', k, j, R(k,j));
    end
    fprintf('\n');
end

end

例如，对于被积函数 f(x) = exp(-x^2) 在区间 [0,1] 上求积分，可以使用以下代码调用 romberg 函数：

f = @(x) exp(-x.^2);
a = 0;
b = 1;
nmax = 10;
tol = 1e-6;

[I, err] = romberg(f, a, b, nmax, tol);

输出结果为：

R(1,1) = 0.683939
R(2,1) = 0.746824
R(2,2) = 0.755183
R(3,1) = 0.750064
R(3,2) = 0.753750
R(3,3) = 0.754710

...

积分值为 I = 0.754710，误差为 err = 9.5245e-07。