Search code examples
pythonsympysubstitutioninverse

Inverting a function using sympy and evaluating the inverted function gives me a wrong answer

I am trying to invert a function and evaluate the inverted function. Basically, I have the equation:

y = 1 / f * ln(1+c*x) / x

where c and f are numerical constants. In python, I have the following code to find x as a function of y and evaluate x for y = 1.5:

import sympy as sp
import math

x, y_prime = sp.symbols("x y_prime", positive=True)
c = 10
f = math.log(1+c) - c/(1+c)
y = 1 / f * sp.log(1+c*x) / x
eqn = sp.Eq(y_prime, y)
sol_x = sp.solve(eqn, x, rational=False)[0]
print(sol_x.subs(y_prime, 1.5))

The output is get is x=0, which I know is wrong. Also, I tried this in Mathematica, where I got x = 1.120173048953996, which I know is correct because sticking in this value of x in the expression for y gives y = 1.5. I am wondering how I can get the correct answer in python. Any help is appreciated.

Solution 
  • from sympy import *
var("f, c, x, y")
eq = Eq(y,  1 / f * log(1+c*x) / x)
sol = solve(eq, x)
print(sol[0])
# out: -LambertW(-f*y*exp(-f*y/c)/c)/(f*y) - 1/c

    As you can see, there is the Lambert W function, which has two branches:

    enter image description here

    The solution returned by SymPy is represented by the upper branch (the blue line on the plot). For this case, the lower (orange) branch is the correct one.

    So, let's modify the solution:

    mod_sol = sol[0].replace(lambda t: isinstance(t, LambertW), lambda t: LambertW(*t.args, -1))
print(mod_sol)
# out: -LambertW(-f*y*exp(-f*y/c)/c, -1)/(f*y) - 1/c

    Now we can evaluate it to find x:

    import math
c_val = 10
f_val = math.log(1+c_val) - c_val/(1+c_val)

print(mod_sol.subs({f: f_val, c: c_val, y: 1.5}).n())
# out: 1.120173048954

    EDIT to satisfy comment:

    I know what the end goal of the code sol[0].replace(lambda t: isinstance(t, LambertW), lambda t: LambertW(*t.args, -1)) is, but I'm struggling to understand how the code goes about to achieve that.

    replace is an "advanced" method that allows to perform pattern-matching operations. In this example, I passed to replace the following arguments:

    • a query in the form of a lambda function: lambda t: isinstance(t, LambertW). SymPy is going to apply this query to every object contained in the symbolic expression; the t argument represents the current object sympy is looking at. The query returns True if a match is found, otherwise False. For example, SymPy start looking at sol[0] (please, look above to see how this expression is composed):

      • At the beginning, t is going to be sol[0]. Is sol[0] a Lambert W function? No, because sol[0] is an addition of two terms, one of which contains the Lambert W function. So, let's move on.
      • Then, t is going to be -1/c: is it a Lambert W function? clearly not, let's move on.
      • t is going to be -LambertW(-f*y*exp(-f*y/c)/c)/(f*y). Is it a Lambert W function? No, it is a multiplication, in which one of the terms is a Lambert W function. So, let's move on.
      • Eventually, t is going to be LambertW(-f*y*exp(-f*y/c)/c). Is this a Lambert W function? Yes, so the query function returns True.

    • a function that returns a replacement object: lambda t: LambertW(*t.args, -1)). If the query function returns True for some object, SymPy will take that object and pass it to the replacement function. Specifically, here I'm creating a new Lambert W function with the same arguments of the original term t, but I'm asking it to consider the lower branch with the -1.

    In the future, if I come across the LambertW function, how do I know which branch gives the correct solution?

    This was the first time I came across the Lambert W function. I looked at Wikipedia, saw the two branches, and I understood that there must be a limiting value somewhere. So I plotted sol[0] and mod_sol (see the chart below). As you can see, up to y ~ 6.5 the lower branch provides the correct results, after that the upper branch must be considered.

    enter image description here