clear
echo on
% In this example, we will analyze a piping network and demonstrate
% how to solve for the flow rates in each of the pipes. During this
% process, we will make use of the Newton's method for finding roots
% of a system of nonlinear equations.
%
% OK. Let's get started. The piping network is given as follows:
%
%
%                   Q1       A  Q3, L3     C     
%                ----->----------->--------|------>------)      
%                            |             |             |
%                           \|/ Q2        \|/  Q4        | Q5
%                            |  L2         |   L4        | L5
%                            |             |             |
%                   Q1       |    Q3,L3    |             |
%                -----<-----------<--------|-------------)      
%                            B             D
%
%   The pressure drop across each of the pipes is given by kLQ^2 where
%   k is a constant and the same for all the pipes, L is the length of
%   the pipe. Given the flow rate Q1, we are to solve for the flow rates
%   Q2, Q3, Q4 and Q5. 
%
%   Since we have four unknowns we need four equations.

pause

%
%   Let's first write mass balance equations:
%  
%        Q1 = Q2 + Q3           (1)
%        Q3 = Q4 + Q5           (2)
%  
%   Next, let's make use of the fact that the pressure drop across two
%   points must be the same irrespective of the path we traverse between
%   the two points. This is analogous to the electrical circuits. 
%
%   Pressure drop across AB:
%
%     k* L2 * Q2^2 = k * L3 * Q3^2 + k * L4 * Q4^2 + k * L3 * Q3^2
%
%   Simplifying, we get
%
%   L2 * Q2^2 = 2* L3 * Q3^2 + L4 * Q4^2  (3)
% 
%   Similarly we can write for the points C and D.
%
%   L5 * Q5^2 = L4 * Q4^2                 (4)
%
%
%   We have our four equations (1)-(4). Note that the first two equations
%   are linear in flow rates while eqns. (3) and (4) and nonlinear, due to
%   the presence of quadratic terms
%
%   Having set up the problem, let's give some values and solves
%
%   Q1 = 1, L3 = L2 = L4 = 1 and L5 = 3
%
%
pause
%   We want to use Newton's method to solve this, so we need a script for
%   function evaluation and one for the Jacobian evaluation
%
%    J(x(k)) * \delta x(k) = -f(x(k))
%
%
%    f = [ Q1-Q2-Q3 ; Q3-Q4-Q5 ; L2 * Q2^2 - 2 * L3 * Q3^2 - L4 * Q4^2;
%          L5 * Q5^2 - L4 * Q4^2]
%
%
%   J = [ -1 -1 0 0 ; 0 -1 -1 -1; 2*L2*Q2, -2*L3*Q3, -2*L4*Q4 0 ;
%       0 0 2*L5*Q5 2*L4*Q4]
%
%   so we will write scripts to do this as well.
pause
%
% 

Q1 = 1
lengths = [ 1 1 1 3] ;

% We can use the anonymous function command to get the expressions of both
% the function we are solving and the Jacobian matrix.
fevaluate=@(x) [Q1-x(1)-x(2);x(2)-x(3)-x(4);...
     lengths(1)*x(1)^2-2*lengths(2)*x(2)^2-lengths(3)*x(3)^2;...
     lengths(4) * x(4)^2 - lengths(3)*x(3)^2];
 
jacobian=@(x) [ -1  -1   0    0;...
            0    1   -1   -1 ;...
            2*lengths(1)*x(1)  -2*lengths(2)*x(2)  -2*lengths(3)*x(3) 0;...
            0     0   -2*lengths(3)*x(3)  2*lengths(4)*x(4)];

x = [ 0.5 0.5 0.25 0.25]' ;
x0 = x; 
fout = fevaluate(x)  ;

n_iter = 0;
while sum(abs(fout))> 0.00001
    % calcualte the Jacobian
    jacobian_x = jacobian(x) ;
    % calculate the correction ;
    dx = jacobian_x\fout ;
    x = x - dx ;
    %recompute fout
    fout = fevaluate(x) ;
    n_iter = n_iter + 1;
    echo off
end
echo on

% the flow rates are
x

pause
% number of iterations is
n_iter

% use the inbuilt fsolve command
[x_fsolve, fval]=fsolve(fevaluate, x0)

% calculate the difference between the two solutions
norm(x-x_fsolve)