clear;
clc;
echo on;

% Tutorial 8

% In this tutorial we show the different converging speed
% between modified Newton's method and ordinary Newton's method.

% Let's look at the equation (x-1)*(x-sqrt(2))^3=0. Obviously, this
% equation has one repeated root 'sqrt(2)', with multiplicity of 3. Thus,
% if we use a ordinary Newton's method to find root, it will converge
% linearly. It is possible to devise an accelerated scheme, however.  If
% your root has a multiplicity of m (m is 3, above, for example), a
% modified rule with:
%
%        x_i+1=x_i-m*(f(x)/f'(x)) 
%
% again converges quadratically.  Effectively we are amplifying the
% Newton's method correction to account for the multiplicity.

% We demonstrate this strategy:
ft8 = @(x) (x-1).*(x-sqrt(2)).^3;

ft8prime = @(x) -(-4*x+sqrt(2)+3).*(sqrt(2)-x).^2;
% We set the interval:
a=.75;
b=1.75;
pause
% First, let's plot the function:
y=[a:(b-a)/50:b];
figure(1)
plot(y,ft8(y));
xlabel('y')
ylabel('f(y)')
grid
pause
hold on

% Let's first look at the ordinary Newton's method, adopting codes from
% example 23...
tol=1e-12; % Set the tolerance... 
xlast=a;
x(1)=b;
i=1;
figure(1)
plot(y,ft8(y))
xlabel('y')
ylabel('f(y)')
grid on
hold on
plot(x(i),ft8(x(i)),'o')
pause(1)
while min(abs(ft8(x(i))),abs(x(i)-xlast))>tol,
  x(i+1)=x(i)-ft8(x(i))/ft8prime(x(i)); % Note, here we use ordinary Newton's method.
  % Thus there is no coefficient before f(x(i))/fprime(x(i).
  xlast=x(i);
  i=i+1;
  figure(1)
  plot(x(i),ft8(x(i)),'o')
  pause(1)
  echo off
end
echo on
hold off
pause 
% How many iteration we used:  
i
pause

% Then, using the similar strategy, let's look at the modified Newton's
% method, by introducing a coefficient before f(x(i))/fprime(x(i) in the
% root finding formula. Note here the coefficient m=3, because we have the
% root sqrt(2) with multiplicity of 3. 
tol=1e-12;
xlast=a; 
xnew(1)=b;
j=1;
figure(1)
plot(y,ft8(y))
xlabel('y')
ylabel('f(y)')
grid on
hold on
plot(xnew(j),ft8(xnew(j)),'*')
pause(1)
while min(abs(ft8(xnew(j))),abs(xnew(j)-xlast))>tol,
  xnew(j+1)=xnew(j)-3*ft8(xnew(j))/ft8prime(xnew(j)); % Note here m = 3.
  xlast=xnew(j);
  j=j+1;
  figure(1)
  plot(xnew(j),ft8(xnew(j)),'*')
  pause(1)
  echo off
end
echo on
hold off
pause
% How many iteration we used this time:  
j
% See it is much smaller! 
pause

% Now let's look at the error as a function of iteration of this two
% method.
answer = sqrt(2); % The actual root sqrt(2)
[m,nN]=size(x);
[m,nModifiedN]=size(xnew);
figure(2)
plot([1:nN],log10(abs(x-answer)),...
     [1:nModifiedN],log10(abs(xnew-answer)))
xlabel('iteration')
ylabel('error')
grid on
legend('ordinary Newton', 'modified Newton','Location', 'NorthEast') ;
% We can see that the the error of modified Newton's method converge much
% faster than ordinary Newton's method in this case. This is the consequence
% of a quadratic convergence.
pause

% To estimate the rate of convergence r, We make the log-log plot of
% error(i+1) vs error(i).  The rate r is the slope of this plot as the error
% gets small (e.g., ignore the first point for the modified Newton's method).
figure(3)
loglog(abs(xnew(1:nModifiedN-1)-answer),abs(xnew(2:nModifiedN)-answer),...
     abs(x(1:nN-1)-answer),abs(x(2:nN)-answer))
xlabel('error(i)')
ylabel('error(i+1)')
grid on
% We can add in the asymptotic behavior for this function. 
% The constant in the quadratic convergence is
% C = f(m+1)th/f(m)th / (m(m+1)).  In this case, m=3 and we have the ratio:
C = 1/3/(2^.5-1)
% so let's add it in!
hold on
ei=abs(xnew(1:nModifiedN-1)-answer);
plot(ei,ei.^2*C,'r--')
hold off
legend('modified Newton', 'ordinary Newton','asymptotic formula','Location','NorthWest') ;
% The slope of the ordinary Newton's method error result is about half of
% the quadratic modified Newton's method result. 

echo off;