# My code is running so long and never gives solution.

3 views (last 30 days)
Heather Holzer on 23 Aug 2019
Answered: Jan on 27 Aug 2019
clear;
close all;
g=32.2; nu=1.21e-5; %ft^2/s
QC=6;
QD=8;
QE=11;
QA=QC+QD+QE;
%Roughness for concrete e=0.001 ft
%rr=e/D
L1=600; D1=1.5; rr1=0.00067; Re1=1.7e6; Q1=11;
L2=600; D2=1.5; rr2=0.00067; Re2=1.7e6; Q2=14;
L3=800; D3=1.25; rr3=0.0008; Re3=1.7e6; Q3=7;
L4=800; D4=1.25; rr4=0.0008; Re4=1.7e6; Q4=7;
L5=400; D5=1; rr5=0.0001; Re5=1.7e6; Q5=4;
L6=400; D6=1; rr6=0.0001; Re6=1.7e6; Q6=1;
L7=800; D7=1.25; rr7=0.0008; Re7=1.7e6; Q7=5;
L8=400; D8=1; rr8=0.0001; Re8=1.7e6; Q8=1;
L9=400; D9=1; rr9=0.0001; Re9=1.7e6; Q9=4;
L10=400; D10=1; rr10=0.0001; Re10=1.7e6; Q10=4;
fF1=FrictionFactor(rr1, Re1);
fF2=FrictionFactor(rr2, Re2);
fF3=FrictionFactor(rr3, Re3);
fF4=FrictionFactor(rr4, Re4);
fF5=FrictionFactor(rr5, Re5);
fF6=FrictionFactor(rr6, Re6);
fF7=FrictionFactor(rr7, Re7);
fF8=FrictionFactor(rr8, Re8);
fF9=FrictionFactor(rr9, Re9);
fF10=FrictionFactor(rr10, Re10);
QV=[Q1; Q2; Q3; Q4; Q5; Q6; Q7; Q8; Q9; QA];
fFV=[fF1, fF2, fF3, fF4, fF5, fF6, fF7, fF8, fF9, fF10]';
eF=1;
cnf=0;
tolQ=QA*1e-6;
while eF>1e-6
cnf=cnf+1;
eQ=tolQ+1;
K1=fF1*L1*8/g/D1^5/pi^2;
K2=fF2*L2*8/g/D2^5/pi^2;
K3=fF3*L3*8/g/D3^5/pi^2;
K4=fF4*L4*8/g/D4^5/pi^2;
K5=fF5*L5*8/g/D5^5/pi^2;
K6=fF6*L6*8/g/D6^5/pi^2;
K7=fF7*L7*8/g/D7^5/pi^2;
K8=fF8*L8*8/g/D8^5/pi^2;
K9=fF9*L9*8/g/D9^5/pi^2;
K10=fF10*L10*8/g/D10^5/pi^2;
cnQ=0;
while eQ>tolQ
cnQ=cnQ+1;
F1=QV(1)+QV(2)-QA;
F2=-QV(1)+QV(3)+QV(4);
F3=QC+QV(6)-QV(2)-QV(5);
F4=QD+QV(8)-QV(4);
F5=QE-QV(6)-QV(10);
F6=QV(7)+QV(5)-QV(3);
F7=QV(9)+QV(10)-QV(8);
F8=K1*QV(1)*abs(QV(1))+K3*QV(3)*abs(QV(3))+K5*QV(5)*abs(QV(5))-K2*QV(2)*abs(QV(2));
F9=K4*QV(4)*abs(QV(4))+K8*QV(8)*abs(QV(8))+K9*QV(9)*abs(QV(9))-K7*QV(7)*abs(QV(7))-K3*QV(3)*abs(QV(3));
F10=K7*QV(7)*abs(QV(7))-K9*QV(9)*abs(QV(9))+K10*QV(10)*abs(QV(10))-K6*QV(6)*abs(QV(6))-K5*QV(5)*abs(QV(5));
Jac=[1, 1, 0, 0, 0, 0, 0, 0, 0, 0;
-1, 0, 1, 1, 0, 0, 0, 0, 0, 0;
0, -1, 0, 0, -1, 1, 0, 0, 0, 0;
0, 0, 0, -1, 0, 0, 0, 1, 0, 0;
0, 0, 0, 0, 0, -1, 0, 0, 0, -1;
0, 0, -1, 0, 1, 0, 1, 0, 0, 0;
0, 0, 0, 0, 0, 0, 0, -1, 1, 1;
2*K1*abs(QV(1)), -2*K2*abs(QV(2)), 2*K3*abs(QV(3)), 0, 2*K5*abs(QV(5)), 0, 0, 0, 0, 0;
0, 0, 0, 2*K4*abs(QV(4)), 0, 0, -2*K7*abs(QV(7)), 2*K8*abs(QV(8)), 2*K9*abs(QV(9)), 0;
0, 0, 0, 0, -2*K5*abs(QV(5)), -2*K6*abs(QV(6)), 2*K7*abs(QV(7)), 0, -2*K9*abs(QV(9)), 2*K10*abs(QV(10))];
b=[-F1; -F2; -F3; -F4; -F5; -F6; -F7; -F8; -F9; -F10];
dQV=Jac\b;
QVn=QV+dQV;
eQ=max(abs(QVn-QV));
tolQ=mean(abs(QV))/1000;
QV=QVn;
end
Re1=4*QV(1)/nu/D1/pi;
Re2=4*QV(2)/nu/D2/pi;
Re3=4*QV(3)/nu/D3/pi;
Re4=4*QV(4)/nu/D4/pi;
Re5=4*QV(5)/nu/D5/pi;
Re6=4*QV(6)/nu/D6/pi;
Re7=4*QV(7)/nu/D7/pi;
Re8=4*QV(8)/nu/D8/pi;
Re9=4*QV(9)/nu/D9/pi;
Re10=4*QV(10)/nu/D10/pi;
fF1n=FrictionFactor(rr1, Re1);
fF2n=FrictionFactor(rr2, Re2);
fF3n=FrictionFactor(rr3, Re3);
fF4n=FrictionFactor(rr4, Re4);
fF5n=FrictionFactor(rr5, Re5);
fF6n=FrictionFactor(rr6, Re6);
fF7n=FrictionFactor(rr7, Re7);
fF8n=FrictionFactor(rr8, Re8);
fF9n=FrictionFactor(rr9, Re9);
fF10n=FrictionFactor(rr10, Re10);
eF=max(abs([fF1, fF2, fF3, fF4, fF5, fF6, fF7, fF8, fF9, fF10]-[fF1n, fF2n, fF3n, fF4n, fF5n, fF6n, fF7n, fF8n, fF9n, fF10n]));
fF1=fF1n;
fF2=fF2n;
fF3=fF3n;
fF4=fF4n;
fF5=fF5n;
fF6=fF6n;
fF7=fF7n;
fF8=fF8n;
fF9=fF9n;
fF10=fF10n;
fFV=[fF1, fF2, fF3, fF4, fF5, fF6,fF7, fF8, fF9, fF10]';
end
function [f, e]=FrictionFactor(rr, Re)
e=1; f=0.01;
while e>1e-6
fn=(-0.5/log10(rr/3.7+2.51/Re/sqrt(f)))^2;
e=abs(fn-f);
f=fn;
end
end

Star Strider on 24 Aug 2019
Note that ‘eF’ never changes and always appears to be:
eF = 0.018448303231474
apparently because:
eF=max(abs([fF1, fF2, fF3, fF4, fF5, fF6, fF7, fF8, fF9, fF10]-[fF1n, fF2n, fF3n, fF4n, fF5n, fF6n, fF7n, fF8n, fF9n, fF10n]));
I have no idea what you are doing, so I will defer to you to sort that.
Heather Holzer on 24 Aug 2019
Thank you. Im trying to find new QV values. i think this makes more sense for that issue.
fF1=fF1n+eF;
However, this is where it stops after running for 10 min.
K8=fF8*L8*8/g/D8^5/pi^2;
K9=fF9*L9*8/g/D9^5/pi^2;
K10=fF10*L10*8/g/D10^5/pi^2;
cnQ=0;
while eQ>tolQ
cnQ=cnQ+1; %<---------
F1=QV(1)+QV(2)-QA;
F2=-QV(1)+QV(3)+QV(4);
F3=QC+QV(6)-QV(2)-QV(5);
Devineni Aslesha on 26 Aug 2019
It is mentioned in the comment that fF1=fF1n+eF; makes more sense for the issue while finding new QV values. However, this equation is not there in the provided code. Is the code supposed to be this way?

Jan on 27 Aug 2019
To accelerate your code, you can store the results of e.g. D1^5/pi^2 in a variable. This avoids 20 expensive power operations in the outer loop.
Matlab is optimized for vector and matrix operations. I assume the pile of different variables is less efficient.
Is there any hint, that the inner loop converges? In single typo can prevent eQ>tolQ from beeing reachable.