How to solve (x+1)e^x=c(constant) in Python?

Question

I would like to solve the (x+1)e^x=c equation in Python.

The equation has been successfully solved by hand using lambert w functions as depicted in the figure below:

Using same steps, I would like to solve (x+1)e^x programmatically. I have coded it using the module SymPy as per the step shown in the figure above , but without success.

Is there any to solve these kinds of equations in Python?

import numpy as np
from sympy import *
n = symbols('n')
sigmao=0.06866
sigmas=0.142038295
theta=38.9
rad=(np.pi/180)*38.9076
cos=np.cos(rad)
sec=1/np.cos(rad)
out = (0.06*0.7781598455*n*(1-exp(-2*0.42*sec*n))+exp(-2*0.42*n*sec)*sigmas)/sigmao
#Apply diff for the above expression. 
fin=diff(out, n)
print(solve(fin,n))

Answer 1

Your expression is very numeric. As sympy's solve tries to find a perfect symbolic solution, sympy gets into troubles.

To find numeric solutions, sympy has nsolve (which allows sympy's expressions but behind the scenes calls mpmath's numeric solver). Unlike solve, here an initial guess is needed:

from sympy import symbols, exp, diff, nsolve, pi, cos

n = symbols('n')
sigmao = 0.06866
sigmas = 0.142038295
theta = 38.9076
rad = (pi / 180) * theta
sec = 1 / cos(rad)
out = (0.06 * 0.7781598455 * n * (1 - exp(-2 * 0.42 * sec * n)) + exp(-2 * 0.42 * n * sec) * sigmas) / sigmao
# Apply diff for the above expression.
fin = diff(out, n)

result = nsolve(fin, n, 1)
print(result, fin.subs(n, result).evalf())

Result: 1.05992379637846 -7.28565300819065e-17

Note that when working with numeric values, you should be very careful to use as many digits as possible to avoid accumulation of errors. Whenever you have an exact expression, it is recommended to leave that expression into the code, instead of replacing it with digits. (Usually, 64 bits or about 16 digits are used in calculations, but for intermediate calculations 80 bits can be taken into account).

To solve the original question with sympy:

from sympy import symbols, Eq, exp, solve

x = symbols('x')
solutions = solve(Eq((x + 1) * exp(x), 20))
for s in solutions:
    print(s.evalf())

Result: 1.92309074332181