matlab code

function [Q] = RootFinder_P2(dens, visc, rough, eps) g = 9.81; Q = [10 1 1 1 1 2 2.2 6 6 2.3 15]/60; D = [405 200 200 200 355 355 250 305 305 305 305]*10^-3; L = [1250 500 400 500 400 600 1100 1100 1250 1000 1000]; Delta = 0.01; t = 0; while t == 0 F = flowRate(Q); t = termCrit(F); J = jacob(F,Q); deltaX = (-J\F)'; Q = Q + deltaX; end function [f] = flowRate(Q) f = [Q(3)+Q(7)+Q(8)-1000/3600 -Q(6)-Q(8)+700/3600 -Q(1)+Q(5)+Q(6)+800/3600 Q(1)+Q(2)-2000/3600 -Q(4)-Q(5)-Q(7)+1300/3600 -node(Q(8),D(8),L(8))+node(Q(6),D(6),L(6))-node(Q(5),D(5),L(5))+node(Q(7),D(7),L(7)) node(Q(5),D(5),L(5))+node(Q(1),D(1),L(1))-node(Q(2),D(2),L(2))-node(Q(4),D(4),L(4)) -node(Q(7),D(7),L(7))+node(Q(4),D(4),L(4))+node(Q(3),D(3),L(3)) ]; end %jacobian function [J] = jacob(f,Q) %find number for matrix size max = size(f); J = zeros(max); for i =1:max for j = 1:max Q(j) = Q(j) + Delta; f = flowRate(Q); J(i,j) = f(i); Q(j) = Q(j) - Delta; f = flowRate(Q); J(i,j) = (J(i,j) - f(i))/Delta; end end end %test termination criterion function [t] = termCrit(f) total = 0; for i = 1:size(f) total = total + f(i).^2; end Fu = sqrt(total); if Fu > eps t = 0; else t = 1; end end function [n] = node(Q,D,L) fric = frictFac(Q,D); n = fric * (8*L)/(pi^2*g*D^5)*Q*abs(Q); end function [f] = frictFac(Q,D) fe = -1.8*log10(((rough/D)/3.7)^1.11 +6.9/reynolds(Q,D)); f = (1/fe)^2; end function [Re] = reynolds(Q,D) Re = (4*abs(Q)*dens)/(pi*D*visc); end end