% f_Newton_Raphson_solution_
% of_a_2_Eqs_System_with_4_inputs.m
%                            07.10.2022
% Find the solutions of the intersection 
% of the functions (ff1) and  (ff2).
% See (Note 1)
% This statement is equivalent to the 
% solution of a two equations system.
% The script requires 4 inputs, in the 
% text marked with "Input"

clear; clc
% Maximum number of iterations
n = 10;
% Difference to stop the iteration
e = 0.0001; 

% 1.- Functions input ___________________
% To obtain the derivatives using diff,  | 
% the functions ff1 and ff2 needs to be  | 
% declared as symbolic                   |
syms ff1(x,y) ff2(x,y)                  %|
% Functions definition                   |
ff1(x,y) = x^2 + y^2 - 4;      % Input_1 |
ff2(x,y) = x*y - 1;            % Input_2 |
                                        %|
% Derivatives of ff1 and ff2 with respect| 
% to x and y are calculated using diff   |
          % In this, case the results are|
dff1x = diff(ff1,x);        % dff1x = 2*x|
dff1y = diff(ff1,y);        % dff1x = 2*y|
dff2x = diff(ff2,x);        % dff1x = y  |
dff2y = diff(ff2,y);        % dff1x = x  |
% _______________________________________|

% 2.-Input of initial assumptions________
% See note 1                             |
 x = 4;                       % Input_3  |
 y = 1;                       % Input_4  |
% _______________________________________|

% 3.- Solutions matrix___________________
% The matrix "Old" will contain initially|		
% the assumed values and then the values |		
% obtained by the iteration	             |
 Old = [x,y];                           %|
%________________________________________|
 
% 4.- Evaluation of functions ___________ 
% The functions ff1 and ff2 are evaluated|
% directly using the actual values of the|
% variables, since they are symbolic     |
                   % In this step:       |
f1 =ff1(x,y);                  %| f1 = 13|
f2 =ff2(x,y);                  %| f2 = 3 |
% _______________________________________|

% 5.- Start of iteration ________________              
% Using a for loop, to perform  iteration|
for i=1:n                               %|
% _______________________________________|
                                      
% 6.- Extraction of solution values _____
x = Old(1,1);  % First element  of row 1 |
y = Old(1,2);  % Second element of row 1 |
% _______________________________________|

% 7.- Evaluation of matrix  "f" for _____ 
% the new values of "x" and "y"         %|
                      % In this step     |
f1 =ff1(x,y);                  % f1 = 13 |
f2 =ff2(x,y);                  % f2 =  3 |
% Definition of matrix "f"               |
f = [f1;f2];                            %|
% _______________________________________|

% 8.- Jacobian ______________________________
% For the Jacobian,                          |
% J =[df1/dx, df1/dy; df2/dx, df2/dy]        |
% the derivatives are evaluated              |
%df1x =dff1x(x,y)    % in the first step     |                    %|
df1x = dff1x(x,y);   % df1x = 8|             |
df1y =dff1y(x,y);    % df1y = 2              |              %|
df2x =dff2x(x,y);    % df2x = 1              |              %|
df2y =dff2y(x,y);    % df2y = 4              |                       %|
% Jaccobian rows                             |
d1 = [df1x,df1y];   % [df1/dx, df1/dy]      %|
d2 = [df2x,df2y];   % [df2/dx, df2/dy]      %|
% Jacobian:                                  |
J =	[d1;d2];                                %|
% ___________________________________________|   

% 9.- Matrix backslash product J\f __________
Jinf = J\f;                                 %|
Jinf = Jinf';      % a column matrix        %|
% __________________________________________%|  

% 10.- New solution matrix___________________
% The new matrix is calculated              %|
New =	Old - Jinf;                         %|
% __________________________________________%|

% 11.- Evaluation of stop condition__________
% for a end or to continue iteration        %|
% First element of Jacobian                 %|
Jinf1= Jinf(1); abs1 = abs(Jinf1);          %|
% Second element of Jacobian                %|
Jinf2= Jinf(2); abs2 = abs(Jinf2);          %|
    if (abs1 < e) && (abs2 < e)             %|
      MaxDiff = max(abs1,abs2);             %|
      break                                 %|
    end                                     %|
% The "New" solution matrix is assigned     %|
% to the "Old" solution matrix              %|
  Old = New;                                %|
% __________________________________________%|

% 12.- End of for loop
end
% ==========================================

% 13.- Evaluacion of stop condition for the
% titel display 
x =New(1,1);y =New(1,2);
Jinf1= Jinf(1); abs1 = abs(Jinf1);
Jinf2= Jinf(2); abs2 = abs(Jinf2);
MaxDiff = max(abs1,abs2); 
if (i==n)
    disp('Maximum iterations number used :')
else
    disp('Solutions:')
end

% 14.- Display outputs
fprintf('x = %0.4f \n', x)
fprintf('y = %0.4f \n', y)
disp(' ')
if MaxDiff < e
    fprintf('Allow. iters: n = %0.0f \n', n)
    fprintf('Iterations: i = %0.0f \n', i) 
    fprintf('MaxDiff = %d \n', MaxDiff)
end

if MaxDiff > e 
    fprintf('Allowable iterations: i = %0.0f \n', n)
    fprintf('Iterations: i = %0.0f \n', i)
    disp('Maximum difference of both variables')
    disp('in the last two iterations: ')
    fprintf('dif = %d \n', MaxDiff)
    disp('prescribed maximum difference: ')
    fprintf('  e = %d \n', e)
    disp('Improve assumptions or increase the')
    disp('maximum allowable number of iterations')
    disp(' ')		 		
end

% Note 1 
% Problem statement
% Find the intersection points of the given
% circle (ff1) and hyperbola (ff2).
% This problem has two solutions, and the obtained
% solution depends on the initial assumptions.
% A graphication of both curves can be used to
% get an idea of the solution values. 
% In this example, the solutions are:
% (1.9319,0.5176)  and (-1.9319,-0.5176)
% If the assumptions are too bad, it will be 
% required a too large number of iterations  
% or the method will fail.

% Note 2
% Anonymous functions example: 
% The equation of a circle x^2 + y^2 - 4 = 0                                        %|
% is represented as anonymous function as 
% a function ff1, as follows:
% ff1 = @(x,y) x^2 + y^2 - 4
% The equation of the hyperbola  x*y=1
% is represented as anonymous function as 
% a function ff2, as follows:
% ff2 = @(x,y) x*y - 1      



