rombf.m

来自「matlab的数学物理方程数值算法源程序。这是"Numerical Method」· M 代码 · 共 40 行

function R = rombf(a,b,N,func,param)
%  Function to compute integrals by Romberg algorithm
%  R = rombf(a,b,N,func,param)
%  Inputs
%    a,b    Lower and upper bound of the integral
%    N      Romberg table is N by N
%    func   Name of integrand function in a string such as
%           func='errintg'.  The calling sequence is func(x,param)
%    param  Set of parameters to be passed to function
%  Output 
%     R     Romberg table; Entry R(N,N) is best estimate of
%           the value of the integral

%* Compute the first term R(1,1)
h = b - a;        % This is the coarsest panel size
np = 1;           % Current number of panels
R(1,1) = h/2 * (feval(func,a,param) + feval(func,b,param));

%* Loop over the desired number of rows, i = 2,...,N
for i=2:N

  %* Compute the summation in the recursive trapezoidal rule
  h = h/2;          % Use panels half the previous size
  np = 2*np;        % Use twice as many panels
  sumT = 0;
  for k=1:2:np-1    % This for loop goes k=1,3,5,...,np-1
    sumT = sumT + feval(func, a + k*h, param);
  end

  %* Compute Romberg table entries R(i,1), R(i,2), ..., R(i,i)
  R(i,1) = 1/2 * R(i-1,1) + h * sumT;   
  m = 1;
  for j=2:i
    m = 4*m;
    R(i,j) = R(i,j-1) + (R(i,j-1) - R(i-1,j-1))/(m-1);
  end
end
return;