function chandranewt
close all
clear all
N     = 100;
c     = .9;
Tr    = .000001;
Ta    = Tr;
tol   = [ Ta, Tr ];
parmsc = [ 40, -1, 1 ]; 
parmsn = [ 40,  1, 0 ];
x = zeros(N,1);
for i = 1:N
    x(i) = (i - .5)/N; 
end

[ hc, it_histc, ierrc ] = nsol(x, @Chandra, tol, parmsc);
[ hn, it_histn, ierrn ] = nsol(x, @Chandra, tol, parmsn);
semilogy((1:length(it_histc)),it_histc,'r*',(1:length(it_histn)),it_histn,'b-')
figure
plot(x, hc, x, hn)
%----------------------------------------------------------------------
function h = Chandra(x)
N = length(x'); c = .9;
h=zeros(N,1);
for i = 1:N
    mu(i) = (i - .5)/N;
end
for i = 1:N
    sum = 0;
    for j = 1:N
        sum = sum + mu(i)*x(j)/(mu(i) + mu(j));
    end
    den = 1 - c*sum/(2*N);
    h(i) = x(i) - 1/den;
end   
%----------------------------------------------------------------------
function [ sol, it_hist, ierr ] = nsol ( x, f, tol, parms )

% Newton solver, locally convergent 
% solver for f(x) = 0
%
% Hybrid of Newton, Shamanskii, Chord
%
% 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
%----------------------------------------------------------------------
function [ l, u ] = diffjac ( x, f, f0 )

%% DIFFJAC computes and factors the jacobian of a function.
%
%  Discussion:
%
%    A forward difference estimate is used.
%
%    The routine DIRDER is needed to compute the directional derivative
%    that defines each column of the jacobian estimate.
%
%  Parameters:
%
%    Input, real X(N), the point at which the jacobian is estimated.
%
%    Input, function F ( X ), the function whose jacobian is to be estimated.
%
%    Input, real F0(N), the value of F(X).
%
%    Output, real L(N,N), U(N,N), the LU factorization of the 
%    (approximate) jacobian.
%
  n = length ( x );

  for j = 1 : n
    zz = zeros ( n, 1 );
    zz(j) = 1;
    jac(:,j) = dirder ( x, zz, f, f0 );
  end

  [ l, u ] = lu ( jac );
%----------------------------------------------------------------------
function z = dirder ( x, w, f, f0 )

%% DIRDER estimates the directional derivative of a function F.
%
%  Parameters:
%
%    Input, real X(N), the point at which the derivative is to be estimated.
%
%    Input, real W(N), the direction.
%
%    Input, real F(X), the function whose derivative is to be estimated.
%
%    Input, real F0(N), the value F(X).
%
%    Output, real Z(N), the estimate for the directional derivative of F.
%

%
% Hardwired difference increment.
%
  epsnew = 1.0E-07;

  n = length ( x );
%
%  Scale the step.
%
  if ( norm ( w ) == 0.0 )
    z = zeros ( n, 1 );
    return
  end

  epsnew = epsnew / norm ( w );

  if ( 0.0 < norm ( x ) )
    epsnew = epsnew * norm ( x );
  end
%
%  DEL and F1 could share the same space if storage
%  is more important than clarity.
%
  del = x + epsnew * w;
  f1 = feval ( f, del );
  z = ( f1 - f0 ) / epsnew;

%----------------------------------------------------------------------
%----------------------------------------------------------------------