how to decrease time needed to find numerical solution of a single nonlinear equation [sympy]

Question

I have a function that draws a line based on input values (which are all floats) and finds a numerical value for B_slv using the sympy function 'nsolve'. The value it finds is correct, but the calculation time is too high as the number of function calls is also high. Is there a way to improve the speed of this calculation? Sacrificing precision for speed is not an issue in many cases.

def cord_line_function(x1,y1,x2,y2, H, w, accuracy, B_old):

    from sympy import symbols, nsolve, cosh, solve
    from numpy import linspace, sinh
    from numpy import cosh as npcosh
    # Solve for B . Start with the initial estimation B_old in the iterative process:
    if not(isnan(B_old)):
        B_estimate = B_old
    else:
        B_estimate = -0.5*(x2-x1)
        B     = symbols('B')
        B_slv = float(nsolve((H/w * (cosh(w*(x1+B)/H)-1) - y1) - (H/w * 
    (cosh(w*(x2+B)/H)-1) - y2), B, B_estimate, quick = True))
    # return a function based on the found values:
    x = linspace(float(x1), float(x2), int((x2-x1)*accuracy))
    y = H/w * (npcosh(w*(x+B_slv)/H)-1)

    # since the first point of the graph is set as 0, we know that y_real = y - y[0]
    y = y - y[0]

    # The cord length would be given by L:
    # L = Σ sqrt(dx^2 + dy^2) = ∫sqrt(1 + (dy/dx)^2)dx
    # Since y = H/w * (cosh(w*(x+B)/H)-1) - A
    # dy/dx = sinh(w*(x+B)/H)
    # gives: L =      ∫sqrt(1+(sinh(w*(x+B)/H))^2) dx
    #          =      ∫sqrt(1+(sinh(w*u    /H))^2) du      with u = B+x   ;     du = dx
    #          = H/w  ∫sqrt(1+(sinh(s)        )^2) ds      with s=w*u/H   ;     ds = w/H * du
    #          = H/w  ∫sqrt(   cosh(s)         ^2) ds
    #          = H/w  ∫        cosh(s)             ds
    #          = H/w  *        sinh(s) + Constant
    s1 = w*(x1+B_slv)/H
    s2 = w*(x2+B_slv)/H
    length = H/w * sinh(s2) - H/w * sinh(s1)

    return [x, y, length, B_slv]

Answer 1

This particular equation can be solved symbolically by SymPy:

In [56]: H, w, x1, x2, y1, y2, B = symbols('H, w, x1, x2, y1, y2, B')

In [57]: eq = (H/w * (cosh(w*(x1+B)/H)-1) - y1) - (H/w * (cosh(w*(x2+B)/H)-1) - y2)

In [58]: eq
Out[58]: 
  ⎛    ⎛w⋅(B + x₁)⎞    ⎞     ⎛    ⎛w⋅(B + x₂)⎞    ⎞          
H⋅⎜cosh⎜──────────⎟ - 1⎟   H⋅⎜cosh⎜──────────⎟ - 1⎟          
  ⎝    ⎝    H     ⎠    ⎠     ⎝    ⎝    H     ⎠    ⎠          
──────────────────────── - ──────────────────────── - y₁ + y₂
           w                          w                      

In [59]: eq.rewrite(exp).expand()
Out[59]: 
   B⋅w  w⋅x₁      B⋅w  w⋅x₂      -B⋅w   -w⋅x₂       -B⋅w   -w⋅x₁           
   ───  ────      ───  ────      ─────  ──────      ─────  ──────          
    H    H         H    H          H      H           H      H             
H⋅ℯ   ⋅ℯ       H⋅ℯ   ⋅ℯ       H⋅ℯ     ⋅ℯ         H⋅ℯ     ⋅ℯ                
──────────── - ──────────── - ──────────────── + ──────────────── - y₁ + y₂
    2⋅w            2⋅w              2⋅w                2⋅w                 

In [60]: s1, s2 = solve(eq.rewrite(exp).expand(), B)

In [61]: print(s1)
H*log(w*(y1 - y2)/(H*(exp(w*x1/H) - exp(w*x2/H))) - sqrt((H**2*exp(2*w*x1/H) - 2*H**2*exp(w*x1/H)*exp(w*x2/H) + H**2*exp(2*w*x2/H) + w**2*y1**2*exp(w*x1/H)*exp(w*x2/H) - 2*w**2*y1*y2*exp(w*x1/H)*exp(w*x2/H) + w**2*y2**2*exp(w*x1/H)*exp(w*x2/H))*exp(w*x1/H)*exp(w*x2/H))*exp(-w*x1/H)*exp(-w*x2/H)/(H*(exp(w*x1/H) - exp(w*x2/H))))/w

In [62]: print(s2)
H*log(w*(y1 - y2)/(H*(exp(w*x1/H) - exp(w*x2/H))) + sqrt((H**2*exp(2*w*x1/H) - 2*H**2*exp(w*x1/H)*exp(w*x2/H) + H**2*exp(2*w*x2/H) + w**2*y1**2*exp(w*x1/H)*exp(w*x2/H) - 2*w**2*y1*y2*exp(w*x1/H)*exp(w*x2/H) + w**2*y2**2*exp(w*x1/H)*exp(w*x2/H))*exp(w*x1/H)*exp(w*x2/H))*exp(-w*x1/H)*exp(-w*x2/H)/(H*(exp(w*x1/H) - exp(w*x2/H))))/w

You can just copy paste those solution expressions into your code or otherwise use lambdify to evaluate the expressions efficiently.

Here solve finds two solutions. I am not sure how the two different solutions relate to your original problem but presumably one of them is the one that you are interested in.

Answer 2

OP has stated that the input values are not sympy symbols, but rather standard python floats.

As such, and building on Oscar Benjamin's answer, we can state that in order to speed up the code, you should solve your equation once to get a closed-form solution, and then turn this closed-form solution into a function which can be numerically calculated (rather than symbolically).

In code:

from sympy import symbols, solve, cosh, exp, lambdify,
H, w, x1, x2, y1, y2, B = symbols('H, w, x1, x2, y1, y2, B')
eq = (H/w * (cosh(w*(x1+B)/H)-1) - y1) - (H/w * (cosh(w*(x2+B)/H)-1) - y2)
sol = solve(eq.rewrite(exp).expand(), B)
# this will produce a numerically evaluated function from the solution
numerical_function = lambdify((H, w, x1, x2, y1, y2), sol)

... 
# call the created function later as you would a normal function
B_slv_1, B_slv_2 = numerical_function(H, w, x1, x2, y1, y2)

In fact, if we're smart about it, we can make it even faster by specifying the correct module. By default, lambdaify is will use numpy, which is optimised for vector-like computations. However, here we're dealing with scalar inputs, so the native python math module would be faster! Thankfully, we can easily tell lambdaify to switch to the math module using the modules keyword.

To illustrate the speedup we get, we can use the timeit module :

from timeit import timeit
setup = '''
from sympy import symbols, solve, nsolve, cosh, exp, Eq, lambdify, nsolve

def calculation_function_factory(module):
    H, w, x1, x2, y1, y2, B = symbols('H, w, x1, x2, y1, y2, B')
    eq: Eq = (H/w * (cosh(w*(x1+B)/H)-1) - y1) - \
        (H/w * (cosh(w*(x2+B)/H)-1) - y2)

    s1, s2 = solve(eq.rewrite(exp).expand(), B)
    return lambdify((H, w, x1, x2, y1, y2), s1, modules=module)


calculation_mt = calculation_function_factory("math")
calculation_np = calculation_function_factory("numpy")


def symbolic_solve(H, w, x1, x2, y1, y2):
    B = symbols('B')
    B_estimate = -0.5*(x2-x1)
    return float(nsolve(
        (H/w * (cosh(w*(x1+B)/H)-1) - y1) - (H/w * (cosh(w*(x2+B)/H)-1) - y2),
        B,
        B_estimate,
        quick=True))

'''

calc_m = "calculation_mt(1,2,3,4,5,6)"
calc_n = "calculation_np(1,2,3,4,5,6)"
solv = "symbolic_solve(1,2,3,4,5,6)"

ITER = 1000

t_solve = timeit(solv, setup=setup, number=ITER)
print(t_solve) # >>> 2.0169183000107296 seconds for ITER iterations
t_calc_n = timeit(calc_n, setup=setup, number=ITER)
print(t_calc_n) # >>> 0.013456200002110563 seconds for ITER iterations
t_calc_m = timeit(calc_m, setup=setup, number=ITER)
print(t_calc_m) # >>> 0.00293279999459628 seconds for ITER iterations