function [sol, it_hist, ierr] = nsol(x,f,tol,parms)% Newton solver, locally convergent % solver for f(x) = 0%% Hybrid of Newton, Shamanskii, Chord%% C. T. Kelley, November 26, 1993%% This code comes with no guarantee or warranty of any kind.%% function [sol, it_hist, ierr] = nsol(x,f,tol,parms)%% inputs:%        initial iterate = x%		 function = f%        tol = [atol, rtol] relative/absolute%			error tolerances%	 parms = [maxit, isham, rsham]%			maxit = maxmium number of iterations%				default = 40%		isham, rsham: The Jacobian matrix is% 		computed and factored after isham%               updates of x or whenever the ratio%		of successive infinity norms of the%               nonlinear residual exceeds rsham.%			isham = 1, rsham = 0 is Newton's method,%			isham = -1, rsham = 1 is the chord method,% 			isham = m, rsham = 1 is the Shamanskii method%      	       defaults = [40, 1000, .5]%% output:%	sol = solution%	it_hist = infinity norms of nonlinear residuals%			for the iteration%	ierr = 0 upon successful termination%	ierr = 1 if either after maxit iterations%             the termination criterion is not satsified%             or the ratio of successive nonlinear residuals%             exceeds 1. In this latter case, the iteration%				is terminted.%%% internal parameter:%       debug = turns on/off iteration statistics display as%               the iteration progresses%% Requires: diffjac.m, dirder.m%% Here is an example. The example computes pi as a root of sin(x)% with Newton's method and plots the iteration history.%%%  x=3; tol=[1.d-6, 1.d-6]; params=[40, 1, 0];%  [result, errs, it_hist] = nsol(x, 'sin', tol, params);%  result%  semilogy(errs)%%% set the debug parameter, 1 turns display on, otherwise off%debug=1;%% initialize it_hist, ierr, and set the iteration parameters%ierr = 0;maxit=40;isham=1000;rsham=.5;if nargin == 4maxit=parms(1); isham=parms(2); rsham=parms(3);endrtol=tol(2); atol=tol(1);it_hist=[];n = length(x);fnrm=1;itc=0;%% evaluate f at the initial iterate% compute the stop tolerance%f0= feval(f,x);fnrm=norm(f0,inf);it_hist=[it_hist,fnrm];fnrmo=1;itsham=isham;stop_tol=atol+rtol*fnrm;%% main iteration loop%while(fnrm > stop_tol & itc < maxit)%% keep track of the ratio (rat = fnrm/frnmo)% of successive residual norms and % the iteration counter (itc)%    rat=fnrm/fnrmo;    outstat(itc+1, :) = [itc fnrm rat];    fnrmo=fnrm;     itc=itc+1;%% evaluate and factor the Jacobian% on the first iteration, every isham iterates, or% if the ratio of successive residual norm is too large%    if(itc == 1 | rat > rsham | itsham == 0)        itsham=isham;        [l, u] = diffjac(x,f,f0);    end    itsham=itsham-1;%% compute the step%    tmp = -l\f0;    step = u\tmp;    xold=x;    x = x + step;    f0= feval(f,x);    fnrm=norm(f0,inf);    it_hist=[it_hist,fnrm];    rat=fnrm/fnrmo;    if debug==1        disp([itc fnrm rat])    end    outstat(itc+1, :)=[itc fnrm rat];%% if residual norms increase, terminate, set error flag%    if rat >= 1        ierr=1;        sol=xold;        disp('increase in residual')        disp(outstat)        return;    end% end whileendsol=x;if debug==1    disp(outstat)end%% on failure, set the error flag%if fnrm > stop_tol    ierr = 1;end