mnewt1d.m

nonlinear eq routines in matlab

function [x, ithist, iflag] = mnewt1d( f, x, m, tolf, tolx, maxit )
%
%  function [x, ithist, iflag] = mnewt1d( f, x, m, tolf, tolx, maxit )
%
%  mnewt1d attempts to compute a root of f using the modified Newton method.
%
%  Input parameters:
%    f       name of a matlab function that evaluates 
%            f and its derivative.
%    x       initial iterate
%    m       estimate of the multiplicity of the root.
%            m = 1 : Newton's method
%            m = 0 : adaptive modified Newton's method
%    tolf    stopping tolerance (optional. Default tolf = 1.e-7)
%            Newton's method stops if  |f(x)| < tolf
%    tolx    stopping tolerance (optional. Default tolx = 1.e-7)
%            Newton's method stops if  |s| < tolx, 
%            where s = -f(x)/f'(x)  is the Newton step.
%    maxit   maximum number of iterations (optional. Default maxit = 100)
%
%
%  Output parameters:
%    x       approximation of the solution. 
%    ithist  array with the iteration history
%            The i-th row of ithist contains  [it, x, fx, s. m]
%    ifag    return flag
%            iflag =  0  |f(x)| <= tolf 
%            iflag =  1  iteration terminated because maximum number of 
%                        iterations was reached. |f(x)| > tolf 
%
%  Matthias Heinkenschloss
%  Department of Computational and Applied Mathematics
%  Rice University
%  Jan 17, 2002
%

% set tolerances if necessary
if( nargin <= 4 ) tolf = 1.e-7; tolx = 1.e-7; maxit = 100; end
if( nargin <= 5 ) tolx = 1.e-7; maxit = 100; end
if( nargin <= 6 ) maxit = 100; end


   
it       = 0;
iflag    = 0;
[fx,fpx] = feval(f, x);
s        = 2*tolx;
m_in     = m;
if( m_in <= 0 ) m = 1; end

while( it < maxit & abs(s) > tolx & abs(fx) > tolf )
   
   if( m_in <= 0 & it >= 2 ) 
       m = (ithist(it,2)-ithist(it-1,2))...
           /(2*ithist(it,2)-x-ithist(it-1,2));
       m = max([1,m]);
   end
   s  = - fx/fpx;
   
   ithist(it+1,:) = [it, x, fx, s, m];
   
   x  = x+m*s;
   it = it+1;
   [fx,fpx] = feval(f, x);
end


% check why the Newton's method truncated and set iflag
if( abs(fx) > tolf )
    % Newton's method truncated because the maximum number of iterations
    % was reached
    iflag = 1;
    return
else
    % Newton's method truncated because abs(fx) <= tolf
    % print info for last iteration
    ithist(it+1,:) = [it, x, fx, 0, 0];
end