function [] = chap_6_sec_6_prob_1()
%
% 
% Written by:
% -- 
% John L. Weatherwax                2007-04-23
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all;  clc; 

% THE EXACT SOLUTIONS: 
% 
% -for the phi iterations:
% 
Tm = (1-sqrt( 1 - 4*2*0.0315 ))/4; % <- x_{n+1}=phi(x_n) converges to this root
Tp = (1+sqrt( 1 - 4*2*0.0315 ))/4; 

%
% -for the (X,Y) iterations:
% 
XSQRm = (0.41-sqrt( 0.41^2 - 4*0.0049 ))/2; 
XSQRp = (0.41+sqrt( 0.41^2 - 4*0.0049 ))/2; 
Xmm   = -sqrt(XSQRm); Ymm = -0.07/Xmm; 
Xmp   = +sqrt(XSQRm); Ymp = -0.07/Xmp; 
Xpm   = -sqrt(XSQRp); Ypm = -0.07/Xpm; 
Xpp   = +sqrt(XSQRp); Ypp = -0.07/Xpp; 

% plot all true roots: 
% 
fh=figure; hold on; grid on; 
plot( Xmm, Ymm, 'kx', 'MarkerSize', 10 );
plot( Xmp, Ymp, 'kx', 'MarkerSize', 10 ); 
plot( Xpm, Ypm, 'kx', 'MarkerSize', 10 ); 
plot( Xpp, Ypp, 'kx', 'MarkerSize', 10 ); 
title( 'locus of all possible zeros of f' ); 

% perform all iterations (and tabulate errors):
% 

N_INFTY = 15; 

% assign the initial conditions for our two by two example 
% and our phi function iterations 
% 
x0 = 0.1; y0 = -0.6; t0 = 0; 

% plot the iterates:
% 
fi=figure; plot( x0, y0, 'b.', 'MarkerSize', 15 ); hold on; grid on; 
%figure(fi); plot( x0, y0, 'gs', 'MarkerSize', 15 ); hold on; grid on; 
title( 'the iterates' ); 

% Iterate each scheme:
%
xn_K=x0; yn_K=y0; 
xn_F=x0; yn_F=y0; 
tn  =t0; 
errs=zeros(N_INFTY,3); 
%fprintf('');
fprintf('ii=%3d: K=(%12.6f,%12.6f,%16.6g) F=(%12.6f,%12.6f,%16.6g) T=(%12.6f,%16.6g)\n',0,...
        xn_K,yn_K,norm( [xn_K-Xmp,yn_K-Ymp] ),xn_F,yn_F,norm( [xn_F-Xmp,yn_F-Ymp] ),tn,norm( [t0-Tm] ));
for ii=1:N_INFTY,

  % iterate our modified Newton-Raphson method:
  xn_Kp1 = xn_K + (1/7)*(   xn_K^2 +   yn_K^2 + 12*xn_K*yn_K + 0.43 ); 
  yn_Kp1 = yn_K + (1/7)*( 6*xn_K^2 + 6*yn_K^2 +  2*xn_K*yn_K - 2.32 );
  %[ xn_Kp1, yn_Kp1 ] 
  %errs(ii,1) = norm( [ xn_K-xn_Kp1, yn_K-yn_Kp1 ] ); 
  errs(ii,1) = norm( [ xn_K-Xmp, yn_K-Ymp ] );
  xn_K=xn_Kp1; yn_K=yn_Kp1; 
  figure(fi); plot( xn_K, yn_K, 'x', 'MarkerSize', 10, 'MarkerEdgeColor', [ 0 1 0 ], 'MarkerFaceColor', [ 0 1 0 ] ); 

  % iterate our fixed point method:
  xn_Fp1 = ( xn_F*yn_F + xn_F + 0.07 ); 
  yn_Fp1 = ( xn_F^2 + yn_F^2 + yn_F - 0.41 );
  %[ xn_Fp1, yn_Fp1 ] 
  %errs(ii,2) = norm( [ xn_F-xn_Fp1, yn_F-yn_Fp1 ] ); 
  errs(ii,2) = norm( [ xn_F-Xmp, yn_F-Ymp ] ); 
  xn_F=xn_Fp1; yn_F=yn_Fp1; 
  figure(fi); plot( xn_F, yn_F, 'd', 'MarkerSize', 10, 'MarkerEdgeColor', [ 0 0 1 ], 'MarkerFaceColor', [ 0 0 1 ] ); drawnow; 
  
  % iterate our phi function:
  tnp1 = 2*tn^2 + 0.0315;
  tn=tnp1; 
  errs(ii,3) = norm( [ Tm - tnp1 ] ); 

  fprintf('ii=%3d: K=(%12.6f,%12.6f,%16.6g) F=(%12.6f,%12.6f,%16.6g) T=(%12.6f,%16.6g)\n',ii,...
          xn_K,yn_K,errs(ii,1),xn_F,yn_F,errs(ii,2),tn,errs(ii,3));  

  %pause; 
end

% plot the error observing the various convergence: 
%
figure; semilogy( errs+eps, 'o' ); grid on; 
legend( gca, 'Mod Newton', 'Functional', 'Phi' ); 
%saveas( gcf, '../WriteUp/Pictures/chap_6_sec_6_prob_1_conv.eps' );