I don't understand why this script is outputting imaginary values

The code is trying to solve 8 nonlinear functions for some fluids homework I am doing. fsolve runs but the output is wrong and I think it has something to do with the equations in the function. I can post the equations I want if that helps but something is wrong with one of the equations because I am getting complex values back.

format long g
xguess = [.0001,.0001,.0001,.0001,.0001,.0001,.0001,.0001];
f = @sub1;
xroot = fsolve(f,xguess);

Solver stopped prematurely.

fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 8.000000e+02.

xroot

xroot = 
  Columns 1 through 3

       6.97403850998521e-06 -  2.13433012297505e-06i       0.000100493295050067 -  1.50866907679873e-07i       1.13602962476239e-05 +  9.48591178840745e-07i

  Columns 4 through 6

       0.000101806929372439 +  1.50866907677379e-07i       -0.00709319523332905 +    0.0220307371181461i         0.0239772983255304 +    0.0107722101188694i

  Columns 7 through 8

           5.19445626759942 -      1.58970795366322i            12.690968180021 +      1.05980531707591i

froot = sub1(xroot); %check answer
froot

froot = 
  Columns 1 through 3

        0.00141471526413984 -  1.98018748310224e-14i       0.000628761226740271 -  4.28890114172002e-15i          -2.99979769977558 -  2.49371793627586e-18i

  Columns 4 through 6

       2.28476800103567e-10 +   7.1553479484986e-11i       1.90680093936635e-09 -  9.53681578153009e-13i        0.00122486904854746 -  2.62454502575338e-11i

  Columns 7 through 8

           4.59921007444615 -      6.43152450035871i           5.16900166768458 -       1.8039460299698i

function f = sub1(x)
f = x*0;
pipeDA = .3;
pipeDB = .45;
pipeLA = 500;
pipeLB = 800;
epsilon = .000045;
density = 720;
viscosity = .00029;
gravity = 9.8;
f(1) = (4*x(2)/(pi*pipeDA^2))-x(1);
f(2) = (4*x(4)/(pi*pipeDB^2))-x(3);
f(3) = x(2)+x(4)-3;
f(4) = (x(6)*(pipeLB/pipeDB)*(x(3)^2/(2*gravity)))-(x(5)*(pipeLA/pipeDA)*(x(1)^2/(2*gravity)));
f(5) = ((density*x(1)*pipeDA)/viscosity) - x(7);
f(6) = ((density*x(3)*pipeDB)/viscosity) - x(8);
f(7) = 1.737*log(.269*(epsilon/pipeDA)+1.257/(x(7)*sqrt(x(5))))+(1/sqrt(x(5)));
f(8) = 1.737*log(.269*(epsilon/pipeDB)+1.257/(x(8)*sqrt(x(6))))+(1/sqrt(x(6)));
end

Answer 1

[Walter Roberson]

There is one real solution, and three complex solutions.

format long g
syms x [1 8]
eqn = sub1(x)

eqn = 

char(eqn.')

ans = '[(400*x2)/(9*pi) - x1; (1600*x4)/(81*pi) - x3; x2 + x4 - 3; (40000*x3^2*x6)/441 - (12500*x1^2*x5)/147; (21600000*x1)/29 - x7; (32400000*x3)/29 - x8; (1737*log(1257/(1000*x5^(1/2)*x7) + 1488652246748361/36893488147419103232))/1000 + 1/x5^(1/2); (1737*log(1257/(1000*x6^(1/2)*x8) + 496217415582787/18446744073709551616))/1000 + 1/x6^(1/2)]'

solv = vpasolve(eqn)

solv = struct with fields:
    x1: 12.995077674207773655531170462597
    x2: 0.91856791246769726078808622233061
    x3: 13.087217992724881133113111008922
    x4: 2.0814320875323027392119137776694
    x5: 0.0032728747813005144831856166371721
    x6: 0.0030252672961245925107603753927112
    x7: 9679092.3366513072744645959307618
    x8: 14621581.48152710857630568264445

partsol = solve(eqn(1:end-1), x(2:end))

Warning: Solutions are only valid under certain conditions. To include parameters and conditions in the solution, specify the 'ReturnConditions' value as 'true'.

partsol = struct with fields:
    x2: [4×1 sym]
    x3: [4×1 sym]
    x4: [4×1 sym]
    x5: [4×1 sym]
    x6: [4×1 sym]
    x7: [4×1 sym]
    x8: [4×1 sym]

lasteqn = subs(eqn(end), partsol);
residue = lasteqn;


%{
tiledlayout('flow')
X = linspace(0, 40);
for K = 1 : length(lasteqn)
    F = matlabFunction(residue(K));
    Y = F(X);
    nexttile; plot(X, real(Y), X, imag(Y)); title(string(K));
end
%}

solparts = subs(x(2:end), partsol);

for K = 1 : length(lasteqn)
    lastsol(K) = vpasolve(lasteqn(K), 10);
    var = symvar(lasteqn(K));
    if isempty(var)
        fprintf('no var for eqn %d\n', K);
        vpa(lasteqn(K), 10)
    elseif length(var) > 1
        fprintf('multiple var for eqn %d\n', K)
        vpa(lasteqn(K), 10)
    else
        fullsol(K,:) = [lastsol(K), subs(solparts(K,:), var, lastsol(K))];
    end
end

double(fullsol)

ans = 
  Columns 1 through 3

           12.9950776742078 +                     0i          0.918567912467697 +                     0i           13.0872179927249 +                     0i
           14.0282884446127 -      5.43288418574407i          0.991601278200763 -     0.384027953529099i           12.6280132058782 +      2.41461519366403i
           14.0282884446127 +      5.43288418574407i          0.991601278200763 +     0.384027953529099i           12.6280132058782 -      2.41461519366403i
            16.736275670398 -      9.83592544438072i           1.18301761562796 -     0.695260600139776i           11.4244633277514 +      4.37152241972476i

  Columns 4 through 6

            2.0814320875323 +                     0i        0.00327287478130051 +                     0i        0.00302526729612459 +                     0i
           2.00839872179924 +     0.384027953529099i        0.00103781043600269 +   0.00211622344437779i        0.00302528961896786 -  6.03670368035582e-06i
           2.00839872179924 -     0.384027953529099i        0.00103781043600269 -   0.00211622344437779i        0.00302528961896786 +  6.03670368035582e-06i
           1.81698238437204 +     0.695260600139776i      -0.000278096916218121 +   0.00125077411029544i        0.00302535657597038 -   1.2072598835405e-05i

  Columns 7 through 8

           9679092.33665131 +                     0i           14621581.4815271 +                     0i
            10448656.220815 -      4046562.01420937i           14108538.8920846 +      2697708.00947291i
            10448656.220815 +      4046562.01420937i           14108538.8920846 -      2697708.00947291i
           12465639.8096758 -      7326068.60684909i           12763883.1661775 +      4884045.73789939i

function f = sub1(x)
f = x*0;
pipeDA = .3;
pipeDB = .45;
pipeLA = 500;
pipeLB = 800;
epsilon = .000045;
density = 720;
viscosity = .00029;
gravity = 9.8;
f(1) = (4*x(2)/(pi*pipeDA^2))-x(1);
f(2) = (4*x(4)/(pi*pipeDB^2))-x(3);
f(3) = x(2)+x(4)-3;
f(4) = (x(6)*(pipeLB/pipeDB)*(x(3)^2/(2*gravity)))-(x(5)*(pipeLA/pipeDA)*(x(1)^2/(2*gravity)));
f(5) = ((density*x(1)*pipeDA)/viscosity) - x(7);
f(6) = ((density*x(3)*pipeDB)/viscosity) - x(8);
f(7) = 1.737*log(.269*(epsilon/pipeDA)+1.257/(x(7)*sqrt(x(5))))+(1/sqrt(x(5)));
f(8) = 1.737*log(.269*(epsilon/pipeDB)+1.257/(x(8)*sqrt(x(6))))+(1/sqrt(x(6)));
end

Answer 2

[John D'Errico]

Hint: the log of a negative number is imaginary.

Hint: The sqrt of a negative number is imaginary.

Do you take logs in there? (Yes.) Do you take square roots? (Yes.)

Does fsolve assure that none of the numbers returned will not cause a problem? (No.)