MATLAB Answers

0

My code is running so long and never gives solution.

Asked by Heather Holzer on 23 Aug 2019
Latest activity Answered by 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

  3 Comments

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.
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);
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?

Sign in to comment.

Products


Release

R2018b

1 Answer

Answer by 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.

  0 Comments

Sign in to comment.