My code is running so long and never gives solution.

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]

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]

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]

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?

Answer 1

[Jan]

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.