I have a problem with sympy log function

from sympy import symbols, solve, log, sqrt
# import sympy
from math import pi

f_T = 0.027

T = 60   # [F]
g = 32.2   # [ft/s^2]
H = 44  # [ft]
H_2 = 0     # [ft]
Q = 26     # [ft^3/s]
L = 1700   # [ft]
K_i = 0.5
K_l = 0.2
K_e = 1
D_i = 0.5     # [ft]
f_i = 0.02
print('D = ' + str(D_i) + ' [ft]')

epsilon = 0.0005
nu = 1.21*10**(-5)    # [ft^2/s]

for i in range(0, 50):

    V = symbols('V')
    D = symbols('D')
    f = symbols('f')
    A = (pi / 4) * D_i ** 2  # [ft^2] value of Area is obtained
    V = solve(A * V - Q, V)  # [ft/s] with Area velocity can be obtained
    Re = V[0] * D_i / nu    # with that velocity Reynolds is calculated
    print('V = ' + str(V[0]) + ' [ft/s]')
    print('Re = ' + str(Re))
    # y = epsilon / D_i
    # print('epsilon/D = ' + str(y))
    # f = float(input('Enter f: '))   # f is obtained from Moody Diagram
    # print('f = ' + str(f))
    D = solve(H_2 + ((V[0] ** 2) / (2 * g)) * ((1700 / D) * f_i + (4 * K_l + K_i + K_e)) - H, D)
    D_i = D[0]  # [ft]
    print('D = ' + str(D_i) + ' [ft]')

    f = solve(-2 * log((epsilon / D_i) / 3.7 + 2.51 / (Re * sqrt(f)), 10) - 1 / sqrt(f), f)
    f_i = f[0]
    print('f = ' + str(f_i))

    cond = input('Continue (y/n): ')

    if cond == 'n':
        break
    else:
        continue

When I run the code, it stops right where f is set to solve the equation:

f = solve(-2 * log((epsilon / D_i) / 3.7 + 2.51 / (Re * sqrt(f)), 10) - 1 / sqrt(f), f)

I know it has to do with the log() function of sympy, because if remove it, it works, but of course the result it is an incorrect answer.

What do you mean with “it stops”? Does it throw an exception? If yes, please show the traceback

I meant, it keeps running and never stops unless I stop it manually.

Hi,

I performed some testing on your code. It appears that the solve function is having an issue with the sqrt function having as its argument the f symbol variable. I modified your code slightly and brought it to its minimum requirements for testing purposes. Entering a value that is independent of the f variable inside the sqrt function (inside the log function) appears to alleviate the previously observed issues. The temporary variable for now is temp.

I would advise you to recheck the formula to determine if it is set up correctly or if the f symbol variable is needed in that particular location of the formula. If it is needed there, then you may have to re-write the formula into smaller expressions, for example, to get the same expected results.

from sympy import symbols, solve, log, sqrt
# import sympy
from math import pi

f_T = 0.027

T = 60     # [F]
g = 32.2   # [ft/s^2]
H = 44     # [ft]
H_2 = 0    # [ft]
Q = 26     # [ft^3/s]
L = 1700   # [ft]
K_i = 0.5
K_l = 0.2
K_e = 1
D_i = 0.5     # [ft]
f_i = 0.02
print('D = ' + str(D_i) + ' [ft]')

epsilon = 0.0005
nu = 1.21*10**(-5)    # [ft^2/s]


y = epsilon / D_i

for i in range(0, 5):

    V = symbols('V')
    D = symbols('D')
    f = symbols('f')
    A = (pi / 4) * D_i ** 2  # [ft^2] value of Area is obtained
    V = solve(A * V - Q, V)  # [ft/s] with Area velocity can be obtained
    Re = V[0] * D_i / nu    # with that velocity Reynolds is calculated
    temp = float(input('Enter f: '))   # f is obtained from Moody Diagram

    D = solve(H_2 + ((V[0] ** 2) / (2 * g)) * ((1700 / D) * f_i + (4 * K_l + K_i + K_e)) - H, D)
    D_i = D[0]  # [ft]
    f = solve((-2 * log((epsilon / D_i) / 3.7 + 2.51 / (Re * sqrt(temp)), 10)) - 1 / sqrt(f), f)
    print(f)



The equation has the form of those that can be solved using Lambert’s function. So, solve is computing a change of variable that turns the equation into the one defining Lambert’s function.

"""
Given an expression assumed to be in the form
    ``F(X, a..f) = a*log(b*X + c) + d*X + f = 0``
where X = g(x) and x = g^-1(X), return the Lambert solution,
    ``x = g^-1(-c/b + (a/d)*W(d/(a*b)*exp(c*d/a/b)*exp(-f/a)))``.
"""

In your case (a,b,c,d,f)=(-2, 458717838856613, -8499189896411140, log(10), 42*log(10))

It is computing the part exp(c*d/a/b)*exp(-f/a) as exp(c*d-b*f)/a/b)

In the process, it falls into evaluating

13882669564194443*log(10)

and for this it evaluates

10**13882669564194443

The “never stops” is the time that it is evaluating this.


Try replacing some/all of the variables that are parameters of the equation, with symbols. That should make it more manageable. Unfortunately, even then there are going to be 9287 solutions.
I stopped it when it was computing all 9287 roots of a polynomial of the form x^{9287}-C=0, for certain constant C. It is taking a lot of time trying some heuristics. Maybe it will finish.