clear
format compact
format short e
echo on
%This program demonstrates the use of Newton's
%method for solving systems of non-linear equations.
%We use as our example a simple piping network.  A
%piping network is very similar to a network of
%electrical resistors, with the pressure drop
%between junctions of pipes being equivalent to a
%voltage drop, and the resistance to the flow of
%water equivalent to the resistance of the connecting
%wires.  The principal difference is that the
%electrical network is linear in that the voltage
%drop is proportional to the current, while the
%piping network is non-linear.  At high Reynolds
%number, the pressure drop in a pipe is proportional
%to the -square- of the flow rate.
pause
%In this problem we will examine what happened to
%the flow rate in my shower when my daughter turned
%on the water to her bath downstairs.  The principal
%difficulty is that we had a couple of filters on
%our -very- iron laden water that provided a tremendous
%front end resistance to the water flow throughout
%our house.  I finally got around to reworking the
%water system a few years ago, but for a while it made
%a good example of how -not- to design a residential
%piping network.
pause
%OK, what does the network look like?  We basically
%have the following:
%
%                       ***(tub valve)**(3)
%                      *
%(1)***(filters)***(2)*
%                      *
%                       ***(shower valve)**(4)
%
%The nodes in this piping network are 1) the source
%(in this case our surge tank), 2) the tee after
%the softener, 3) Eleanor's tub, and 4) my shower.
%The resistances are R12: the filters, R23: the tub
%valve, and R24: the shower valve.
pause
%The total head loss between any paths connecting
%the same two nodes in a piping network must be the
%same.  In this case, this means that the loss from
%(1) to (3) and (1) to (4) must be equal to the
%total head provided by the system.  We shall call
%this ptot.  Recognizing that the resistance provided
%by each of the parts of the system is proportional
%to the square of the flow rate, we obtain the pair
%of equations:
%
%ptot-r12*(x(1)+x(2))^2-r23*x(1)^2=0
%and
%ptot-r12*(x(1)+x(2))^2-r24*x(2)^2-dh=0
%
%where x(1) is the flow rate into the tub and x(2)
%is the flow rate into the shower.  dh is the change
%in water pressure due to my shower head being about
%15 feet above Eleanor's tub.
pause
%
%Now we are ready to solve the problem.  First, let's
%put in a few numbers:
ptot=50; %psig
r12=20;  %psig/(gal/min)^2
r24=20;  %psig/(gal/min)^2
dh=8;    %psig
%
%We also need the resistance r23 of Eleanor's tub
%valve.  If the valve is open we take the resistance
%to be r23=5 psig/(gal/min)^2 and if it is closed, we
%let r23=1e10 (actually, this is not unrealistic since
%her tub valve leaks a little).  We first look at
%the closed valve case.  Thus:
r23=1e10;
pause
%We start with an initial guess for the flow rates:
x(1)=1;
x(2)=1;
v=[-.5,2,-.5,2];
figure(1)
plot(x(1),x(2),'o')
grid
xlabel('tub flow rate (gal/min)')
ylabel('shower flow rate (gal/min)')
axis(v)
axis(axis)
hold on
f(1)=ptot-r12*(x(1)+x(2))^2-r23*x(1)^2;
f(2)=ptot-r12*(x(1)+x(2))^2-r24*x(2)^2-dh;
while sum(abs(f))>.0001,
  %Now we define the jacobian matrix:
  j(1,1)=-2*r12*(x(1)+x(2))-2*r23*x(1);
  j(1,2)=-2*r12*(x(1)+x(2));
  j(2,1)=-2*r12*(x(1)+x(2));
  j(2,2)=-2*r12*(x(1)+x(2))-2*r24*x(2);
  %and we compute the change in our guess x:
  dx=j\f';
  %we update x:
  x=x-dx';
  %and we recompute f:
  f(1)=ptot-r12*(x(1)+x(2))^2-r23*x(1)^2;
  f(2)=ptot-r12*(x(1)+x(2))^2-r24*x(2)^2-dh;
  figure(1)
  plot(x(1),x(2),'o')
  pause(1)
  echo off
end
echo on
hold off
pause
%The flow rates were:
x
%The flow rate in the shower is thus a nice
%one gallon/minute, about right for a low flow
%rate shower head.  The leak in the tub is about
%a cup a day (actually, we put a cup under it
%to give fresh water to the cat).
pause
%You should also note that it took quite a while
%to iterate into the correct solution.  This was
%because the dependence on x(1) was much stronger
%than the dependence on x(2) due to the large
%resistance of the tub valve.  The Jacobian was 
%approaching singularity:
cond(j)
pause
%Now let's see what happens when the tub valve is
%open.  The resistance is now:
r23=5;  %psig/(gal/min)^2
%We again need the initial guess for the flow rates.
%We will just use the values that we had for a closed
%tub valve:
figure(1)
hold on
plot(x(1),x(2),'or')
pause
%
%We start the iterative procedure:
f(1)=ptot-r12*(x(1)+x(2))^2-r23*x(1)^2;
f(2)=ptot-r12*(x(1)+x(2))^2-r24*x(2)^2-dh;
while sum(abs(f))>.0001,
  %Now we define the jacobian matrix:
  j(1,1)=-2*r12*(x(1)+x(2))-2*r23*x(1);
  j(1,2)=-2*r12*(x(1)+x(2));
  j(2,1)=-2*r12*(x(1)+x(2));
  j(2,2)=-2*r12*(x(1)+x(2))-2*r24*x(2);
  %and we compute the change in our guess x:
  dx=j\f';
  %we update x:
  x=x-dx';
  %and we recompute f:
  f(1)=ptot-r12*(x(1)+x(2))^2-r23*x(1)^2;
  f(2)=ptot-r12*(x(1)+x(2))^2-r24*x(2)^2-dh;
  figure(1)
  plot(x(1),x(2),'or')
  pause(1)
  echo off
end
echo on
hold off
pause
%Note that the iteration converged much faster
%this time.  This is reflected in the small
%condition number of the Jacobian:
cond(j)
%which is three orders of magnitude smaller
%than the earlier value.
pause
format short
%The new flow rates are:
x
%As you can see, Eleanor's tub really affected my 
%shower's flow rate!  This was very inconvenient
%when you had shampoo in your hair!
echo off