% demonstration of Newton's method on a system of equations
%
% Author: James C. Sutherland
% Date: August, 2009
%
% The equation set we are solving is:
%   x^2 + y^2 = 1
%   y = x^3 - x + 1
%
% we rewrite these as
%  /                        \
%  | f1 = x1^2 + x2^2 - 1   |
%  | f2 = x2 - x1^3 +x1 - 1 |
%  \                        /
%
% Note that we will need the Jacobian:
%
%   /               \
%   |  d f1   d f1  |
%   |  ----   ----  |   /                   \
%   |  d x1   d x2  |   |  2 x1       2 x2  |
%   |               | = |                   |
%   |  d f2   d f2  |   | -3 x1^2+1   1     | 
%   |  ----   ----  |   \                   /
%   |  d x1   d x2  |
%   \               /
%

clear; clc; close all;


atol = 1e-4;  % the absolute error tolerance - for convergence check

% build a function to calculate f1 and f2 as a function of x1 and x2
fval = @(x)( [x(1).^2 + x(2).^2 - 1;  x(2)-x(1)^3+x(1)-1] );

% build a function to calculate the jacobian
jacobian = @(x)( [2*(1),  2*x(2);  -3*x(1).^2+1,  1 ] );


x = [0;0];  % set the initial guess

res = atol+1;
iter=0;
fprintf('Iteration   x1      x2   Residual L2 norm\n');
fprintf('----------------------------------------------\n');

while res > atol
  
  J = jacobian(x);
  f = fval(x);
  
  dx = -J\f;
  
  x = x + dx;

  % calculate the L2 norm of the residual vector at the current solution
  res = norm( fval(x), 2 );
  iter=iter+1;
  
  fprintf('%5i    %-7.4f  %-7.4f %9.2e\n',iter,x(1),x(2),res);

end

