romberg.m

龙贝格计算程序

function R = romberg(f, a, b, n) 
format long 
% ROMBERG -- Compute Romberg table integral approximation. 
% 
% SYNOPSIS: 
% R = romberg(f, a, b, n) 
% 
% DESCRIPTION: 
% Computes the complete Romberg table approximation to the integral 
% 
% / b 
% I = | f(x) dx 
% / a 
% 
% PARAMETERS: 
% f - The integrand. Assumed to be a function callable as 
% y = f(x) 
% with `x' in [a, b]. 
% a - Left integration interval endpoint. 
% b - Right integration interval endpoint. 
% n - Maximum level in Romberg table. 
% 
% RETURNS: 
% R - Romberg table. Represented as an (n+1)-by-(n+1) lower 
% triangular matrix of integral approximations. 
% 
% SEE ALSO: 
% TRAPZ, QUAD, QUADL. 

% NOTE: all indices adjusted for MATLAB's one-based indexing scheme. 

% Pre-allocate the Romberg table. Avoids subsequent re-allocation which 
% is often very costly. 
R = zeros([n + 1, n + 1]); 

% Initial approximation. Single interval trapezoidal rule. 
R(0+1, 0+1) = (b - a) / 2 * (feval(f, a) + feval(f, b)); 

% First column of Romberg table. Increasingly accurate trapezoidal 
% approximations. 
for i = 1 : n, 
    h = (b - a) / 2^i; 
    s = 0; 
    for k = 1 : 2^(i-1), 
        s = s + feval(f, a + (2*k - 1)*h); 
    end 
    R(i+1, 0+1) = R(i-1+1, 0+1)/2 + h*s; 
end 

% Richardson extrapolation gives remainder of Romberg table. 
% 
% Note: The table is computed by columns rather than the more 
% traditional row version. The reason is that this prevents frequent 
% (and needless) re-computation of the `fac' quantity. 
% 
% Moreover, MATLAB matrices internally use ``column major'' ordering so 
% this version is less harmful to computer memory cache systems. This 
% reason is an implementational detail, though, and less important in 
% introductory courses such as MA2501. 
for j = 1 : n, 
    fac = 1 / (4^j - 1); 
    for m = j : n, 
        R(m+1, j+1) = R(m+1, j-1+1) + fac*(R(m+1, j-1+1) - R(m-1+1, j-1+1)); 
    end 
end