clear
format compact
echo on
%This program provides the solution to the shower/W.C.
%problem using Newton's method.

global r26pass %We make the WC resistance a global variable.

%We have the temperatures:
th=125;
tc=60;

%x(1) is the W.C. flow, x(2) is the hot water shower
%flow, and x(3) is the cold water shower flow.
pause

%We must specify the resistance of the WC valve.
echo off
for flag = 1:2,
if flag == 1,
 disp('We start with the valve closed:')
 r26=10^10
else
 disp('Now for an open valve')
 r26=10
end
r26pass=r26;  %We send the resistance of the WC to the function.
echo on
%We have the initial guess for the flow rates:
x=zeros(3,1); %We use a column vector for the solution.
x(1)=1;
x(2)=1;
x(3)=1;
%We calculate the pressure drops to see if they are
%zero already:
f=pdrop(x);

%If the value of f is sufficiently different from
%zero we go through the Newton's method iterations:

while max(abs(f))>.0001,
  %Now we must calculate the Jacobian using a finite
  %difference scheme.  We shall use a center difference
  %formula for greater accuracy.  Remember that we
  %have set x as a column vector.
J=zeros(3,3); %We initialize the Jacobian
  for i=1:3,
    e=zeros(3,1);
    e(i)=.00001;
    J(:,i)=(pdrop(x+x.*e)-pdrop(x-x.*e))/(2*e(i)*x(i));
  end;
  dx=J\f;
  x=x-dx;
  f=pdrop(x);
echo off
end
echo on
disp('The flow rates are:')
x
disp('and the temperature is:')
if x(3)>0,
  Ts = (x(2)*th + x(3)*tc)/(x(2) + x(3))
else,
  Ts = th
end
if flag == 1,
 Tsclosed=Ts;
else
 Tsopen=Ts;
end
pause
end
%The increase in the temperature is thus:
deltaT=Tsopen-Tsclosed
%You can show that this difference decreases with
%decreasing r12 and r13 by changing the values in
%the attached function pdrop.m.  It turns out that
%you need to reduce r12 rather than r13 (e.g., fix
%the cold water line).  If you fix the hot water
%supply without changing the cold, you actually make
%things worse!
pause
%We can get the solution alot more easily if we
%use the matlab non-linear equation solver 'fsolve':
x=fsolve('pdrop',[1;1;1])
Ts = (x(2)*th + x(3)*tc)/(x(2) + x(3))
%which yields the same result.  This actually uses
%the same algorithm, but with extra stuff added to
%make the Newton method more robust.
echo off
end



__________________________



function y=pdrop(x)
%This function returns as an array the equations
%whose roots we are trying to determine.  First we
%set up the resistances:

global r26pass %We take in the WC resistance as a global variable.

r12=20;
r13=20;

r24=20;
r34=5;

r45=10;

%this is the resistance of the W.C.:
r26=r26pass;

dptot=50;

%We define the flow rates in terms of elements of
%the vector x(i):

qt=x(1);
qh=x(2);
qc=x(3);

%We define the pressure drops along each path:

dp12= sign(qt+qc) * r12 * (qt + qc)^2;
dp13= sign(qh) * r13 * (qh)^2;
dp24= sign(qc) * r24 * (qc)^2;
dp26= sign(qt) * r26 * (qt)^2;
dp34= sign(qh) * r34 * (qh)^2;
dp45= sign(qh+qc) * r45 * (qh + qc)^2;

%And, finally, we sum the pressure drops to determine
%whether they match the total pressure drop:

y=zeros(3,1); %We initialize the output array as a column vector.
y(1)= dptot- (dp12 + dp26);
y(2)= dptot- (dp13 + dp34 + dp45);
y(3)= dptot- (dp12 + dp24 + dp45);