%demo of newton's method.
%two parts - will do netwon's method for some starting value, with pause
%and graphics accomanying the first 10 interations.
%will also count up number iterations for convergence over a grid and plot
%convergence speed as a function of starting guess/approximation.

%clear fig and workspace
clf
clear
%grid step
grdstp=.001;
%resolution for determining zero 10e(resol)
resol=-10;
%region for convergence study (symmetric about zero)
mx=2;
mn=-mx;
%x vector of initial guesses/approximations
x=[mn:grdstp:mx]';
%num steps in x vector
nstp=length(x);
%fn to solve
y=x.^3-1;
plot(x,y)
grid
hold
pause

%max interations
cntmax=300;
%set number iterations to larger than max, to id init values that did not
%converge when done
cntmaxp1=cntmax+1;
numit=cntmaxp1(ones(nstp,1));
%find if one of initial trial solutions is good - if so - set iterations to
%zero, NaN out x and y. find indices that meet conditions, if there are any
%do the processing
sindx=find(abs(y)<10^resol);
nzero=length(sindx);
if ( nzero ~= 0 )
    soln=x(sindx);
    x(sindx)=NaN;
    y(sindx)=NaN;
    numit(sindx)=0;
end

%anlaytic expression for derivative
slp=3*x.^2;
%find places where newtons method will fail because derivative is zero and
%save them, remove from calc by setting x nad y to NaN.
zindx=find(abs(slp)<10^resol);
nzero=length(zindx);
if ( nzero ~= 0 )
    zslp=x(zindx);
    x(zindx)=NaN;
    y(zindx)=NaN;
end

%start building matrix with history of guesses/approximations - starts as
%uniform spacing on x axis, will replace each approximation with new
%approximation for next step.
xi(:,1)=x;
%yi does not need to be saved into matrix
yi=y;

%plot initial guess/approximation
xpt=1;
plot(x(xpt),y(xpt),'b+')
%interate looking for soln
for cnt=2:cntmax
%calc f(x) for current approximation
    yi=xi(:,cnt-1).^3-1;
%calc f'(x) at current approximation
    slp=3*xi(:,cnt-1).^2;
%test for slope=0 - remove these points from further procesing
    slp(slp<1e-5)=NaN;
%find new x such that y=slp*x
    adj=yi./slp;
%calc new approximantion by adjusting old approximation
    xi(:,cnt)=xi(:,cnt-1)-adj;
    
%plot up what is going on (old approximation, new approximation, line
%between them, line from new approx to f(x)
    plot(xi(xpt,cnt),0,'ro')
    plot(xi(xpt,cnt),xi(xpt,cnt)^3-1,'b+')
    plot([xi(xpt,cnt) xi(xpt,cnt-1)],[0 xi(xpt,cnt-1)^3-1],'g')
    plot([xi(xpt,cnt) xi(xpt,cnt)],[0 xi(xpt,cnt)^3-1],'m')
    
%now find approximations that have converged, record the number of
%iterations and set stuff to NaN to remove from processing loop
    tsty=abs(yi);
    indx=find(tsty<10^resol);
    nzero=length(indx);
    if ( nzero ~= 0 )
        numit(indx)=cnt;
        xi(indx,cnt)=NaN;
    end
%pause to see what is going on during first 10 interations
    if( cnt<10 )
        pause
    end
end
%done - max interations
%find points that did not converge, note them, set x and numit to NaN
ncindx=find(numit==cntmax+1);
nnc=length(ncindx);
if(nnc~=0)
    ncx=x(ncindx);
    x(ncindx)=NaN;
    numit(ncindx)=NaN;
end
%wait till ready to continue
pause
%clear plot
clf
%plot x (initial guess/approximation), and number interations for that
%guess/approximation to converge
%plot does not plot NaNs - so the bad points don't get plotted
plot(x,numit)
hold
%plot initial guesses/approximations that were the exact solution,
%had a zero derivative, never converged for unknown reason.
plot(zslp,zeros(length(zslp)),'ro')
plot(soln,zeros(length(soln)),'bo')
plot(ncx,zeros(length(ncx)),'r+')
%label and grid
xlabel('starting guess')
ylabel('number iterations to convergence')
grid