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 == 4 maxit=parms(1); isham=parms(2); rsham=parms(3); end rtol=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 while end sol=x; if debug==1 disp(outstat) end % % on failure, set the error flag % if fnrm > stop_tol ierr = 1; end