richextr.m

来自「用matlab写的一些数值算法」· M 代码 · 共 35 行

function  [y, info, R] = richextr(F,h0,q,p,tol)
% Richardson extrapolation of values given by 
%   F(h), h = h0, h0/q, h0/q^2, ...
% Assume that  F(h) = F(0) + a1*h^p(1) + ... + am*h^p(m) ,
% where  m = length(p) .
% Stop when two adjacent values in the same column differ 
% less than  tol .  
% info = [estimated error, row number, column number]

% Version 11.12.2003.  INCBOX

% Initialize the scheme
m = length(p);  R = zeros(m+2,m+1);
info = [inf  0  0];                          % best so far
h = h0;  R(1,1) = feval(F,h);
for  i = 1 : m+1
  h = h/q;   R(i+1,1) = feval(F,h);
  for  k = 1 : min(i,m)   
    d = R(i+1,k) - R(i,k);    
    if  abs(d) <= tol          % desired accuracy obtained
      y = R(i+1,k);  
      info = [abs(d)  i+1  k];
      if  nargout > 2            % return active part of R
        R = R(1:i+1,1:i);
      end
      return
    end
    if  abs(d) < info(1)              % update best so far
      info = [abs(d)  i+1  k];
    end
    R(i+1,k+1) = R(i+1,k) + d/(q^p(k) - 1);  % extrapolate
  end
end
% Required accuracy not obtained.  Return best estimate
y = R(info(2),info(3));