Solve equation without symbolic toolbox

Hi. I am trying to solve "0.0178508*x^2 - 1/(0.694548*(1/x)^0.9 + 0.000136807)^5.44=0" this equation which is a simplified equation derived from Reynold's number equation for head loss.

syms x

solve(0.0178508*x^2 - 1/(0.694548*(1/x)^0.9 + 0.000136807)^5.44==0,x);

I tried this code. It's not giving me any solution for symbolically. Can anyone help me to solve this equation for x. Thanks.

FYI.

This is the actual equation.

hL=19.5, L=18, g=9.8, D=0.0889, Ɛ=0.000045, v=1.37, V̇=?

Answer 1

[Walter Roberson]

digits(50)
syms x
F = 0.0178508*x^2 - 1/(0.694548*(1/x)^0.9 + 0.000136807)^5.44;
sol = vpasolve(F, x, [1e10 1e15])

sol = 

242078482462.87555530675540801383196081859252378323

subs(F, x, sol)

ans = 

0.0000000000000000000000000000000000037615819226313200254999569191111861690197297816707

That is non-zero. Is there actually a root there?

subs(F, x, sol-1e-25)

ans = 

We see that there is a sign change in the function close to the location found, which lends weight to the possibility of a root. But is it an actual root or is there a singularity involved?

If you examine the 0.0178508*x^2 term, you can see that there is no discontinuity in that, at least over the reals

If you examine the (1/x) term, then you can see that would have a discontinuity at 0 but not at any other real value. So in itself 0.694548*(1/x)^0.9 would not have a discontinuity. But the overall term is 1/(0.694548*(1/x)^0.9 + 0.000136807)^5.44 : could it have a discontinuity? Well, if we think about positive real x, 0.694548*(1/x)^0.9 is positive real, and add a positive real constant to that gives you a positive real, so (0.694548*(1/x)^0.9 + 0.000136807) is a positive real. Could the ^5.44 part make it zero? It could if the base expression were less than eps(realmin)^(1/5.44) but with 0.694548*(1/x)^0.9 being non-negative for finite real x and with 0.000136807 being added, the 0.694548*(1/x)^0.9 + 0.000136807 expression cannot be as small as 1e-60 . So (0.694548*(1/x)^0.9 + 0.000136807)^5.44 cannot be 0 for positive real x, so that term cannot form a singularity either.

By continuity, we can therefore see that there must indeed be an actual zero for the function near the location indicated.

[Torsten]

And without symbolic toolbox:

format long
f = @(x)0.0178508*x.^2 - 1./(0.694548./x.^0.9 + 0.000136807).^5.44

f = function_handle with value:
    @(x)0.0178508*x.^2-1./(0.694548./x.^0.9+0.000136807).^5.44

sol = fzero(f,[1e10 1e15])

sol = 
     2.420784824628759e+11

[S M Inzamul Islam]

Hi. THanks for you answer. That is a root. I have added the actual equation.

[S M Inzamul Islam]

Also the possible solution I am expecting ranges from 1.97-2.4.

[Walter Roberson]

hL=19.5; L=18; g=9.8; D=0.0889; epsilon=0.000045; v=1.37;
hL_eqn = @(V) 1.07 * (V.^2.*L)./(g.*D.^5) .* (log(epsilon ./ (3.7*D) + 4.62.*(v.*D./V).^0.9).^(-2))

hL_eqn = function_handle with value:
    @(V)1.07*(V.^2.*L)./(g.*D.^5).*(log(epsilon./(3.7*D)+4.62.*(v.*D./V).^0.9).^(-2))

to_solve = @(V) hL - hL_eqn(V);

syms V
fplot(to_solve(V), [0 0.03] )

format long g
sol = fzero(to_solve, [0 0.03])

sol = 
         0.022609084060871

The equation has a discontinuity around 0.1 and after that does not rise as high as 0. There is no solution near 1.97 to 2.4 ... not even close to one.

[S M Inzamul Islam]

Hi Walter. Thank you so much. This answer is correct. 1.97-2.4 is the final answer. Sorry i forgot to mention that.

[John D'Errico]

A problem in all of this analysis is that the coefficients/exponents are given only to as few as 3 significant digits.

Answer 2

[John D'Errico]

Do you understand this is not a simple quadtatic polynomial? In fact, it is equivalent to a rather high order polynomial. And is is known (Abel-Ruffini, dating back to 1799) that a polynomial equation of degree higher than 4 has no algebraic solution (except for some trivial cases.) And your simply written problem is equivalent to a polynomial of degree over 100.

Anyway, I'm not sure what you were expecting to see as a result. You gave MATLAB no place to return anything into a variable, but then tyou eneded the line with a semi-colon, so nothing would be displayed on the command line.

syms x
F = 0.0178508*x^2 - 1/(0.694548*(1/x)^0.9 + 0.000136807)^5.44;
xsol = vpasolve(F==0,x)

xsol = 

Note that vpasolve will return ONE solution, but as I said, there may be literally hundreds of solutions. Why do I say that? You are raising things to the power 5.44. One can theoretically turn that into a polynomial, with only integer exponents. But now the exponents will be at least 100, and here, surely more. Most of those solutions will be complex valued of course.

Does a real solution exist at all? Perhaps not. It does not appear to cross zero, at least for positive x.

fplot(F,[0,1])

[Walter Roberson]

But now the exponents will be at least 100

5.44 = 136/25 so you could potentially get away with "only" degree 25. Well, at least for that part: the ^0.9 also needs to be taken into account. So degree 50 perhaps.

[S M Inzamul Islam]

I guess the question I am asking what code will give me a possible value!!!!

[John D'Errico]

I showed you code that returns a "value" It is a complex value. But I also point out it is entirely possible no real solution exists.

[S M Inzamul Islam]

Thanks @John D'Errico