3

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.

ellipse314
  • 47
  • 3
  • sympy seems to have issues when being combined with other modules like `math`. When using `sp.log` instead of `math.log` your code gives another result than 0, but still not the one you are looking for. – MagnusO_O Sep 04 '22 at 10:07
  • I made a slight change to my code. Since `f` is a float, I used the `nsimplify` function when I write down the expression for `y` as follows: `y = sp.nsimplify(1 / f * sp.log(1+c*x) / x)`. `nsimplify` expresses all the floats in fractional form where the numerator and denominator are integers. After about 3 minutes of computation, the print statement outputted: `0.0666666666666667*(110*LambertW(0.109071590447994 - 0.0109071590447994*log(285311670611)) - 15.0 + 1.5*log(285311670611))/(10 - log(285311670611))`. How do I ask python to evaluate this expression and give me a single numerical value? – ellipse314 Sep 04 '22 at 10:59
  • I tried to manually evaluate the above expression using `scipy`'s `lambertw` function but still got an answer very close to `0` – ellipse314 Sep 04 '22 at 11:10

1 Answers1

1
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

Davide_sd
  • 10,578
  • 3
  • 18
  • 30
  • Thanks very much for your answer. I am a little new to Python and have a couple of questions. 1) 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. With regards to the `replace` function, I am more used to explicitly specifying a function like `sin` or `cos` as its arguments. Would you be able to explain that part of your code. 2) In the future, if I come across the `LambertW` function, how do I know which branch gives the correct solution? – ellipse314 Sep 05 '22 at 12:17