Solving ill-posed non-linear equations numerically in python/SymPy

Question

I'm trying to get a solution by running the code below.

Python just "hangs" and won't find a numeric solution. I can use an app on my phone (Desmos) to graph the functions and find a numeric solution easily, 0.024. Does python have limitations when solving for 2 decimal places?

import sympy

x = sympy.symbols('x')
e_1 = x**-0.5
e_2 = -2*sympy.log(0.0001*3.7**-1*0.05**-1+2.51*350000**-1*x**-0.5, 10)
sol = sympy.solve(e_2 - e_1, x, 0.024)
num = float(sol[0])
print(num)

nowhere in your code are the values of `f_x` and `g_x` defined. — President James K. Polk, Sep 22 '22 at 11:23
What is ```f_x``` and ```g_x```? If you mean ```e_1``` and ```e_2``` instead, the equation is simply not analytically solvable. — BernieD, Sep 22 '22 at 11:24
A numerical solution is retrieved with "nsolve" instead of "solve". — Michael Butscher, Sep 22 '22 at 11:26
`diameter` is not defined. Please actually test the code in a fresh python process before posting. — Oscar Benjamin, Sep 22 '22 at 11:57
Why do you use the ** operator? Awful to read, inefficient to execute. Just divide by the value. — duffymo, Sep 22 '22 at 19:48

score 4 · Answer 1 · answered Sep 22 '22 at 12:08

Usually, nsolve is the SymPy tool used to numerically solve an equation (or a system of equations). However, I wasn't able to use: it kept raising errors. The problem is that your function is defined on a very small region, and the zero is very close to the boundary:

So, in this case we can try numerical root finding techniques. Again, some tools might fail, but after a few tries I've found that bisect works fine:

from sympy import *
from scipy.optimize import root, brentq, bisect
x = symbols('x')

# you didn't provide the diameter, so I've computed it based on your expected solution
d = 1.56843878182221
e_1 = x**-0.5
e_2 = -2 * log(0.00013 * 7-1*d-1+2.51350000**-1*x**-0.5, 10)

# convert the symbolic expression to a numerical function
ff = lambdify(x, e_1 - e_2)
root, output = bisect(ff, 0.023, 0.025, full_output=True)
print(root)
# 0.024000000001862646
print(output)
#      converged: True
#           flag: 'converged'
# function_calls: 32
#     iterations: 30
#           root: 0.024000000001862646

smichr · Answer 2 · 2022-09-22T22:10:35.913

The fixed point method is a great one to use for situations like this. Or at least the principles of transforming the equation into a compatible form can benefit standard solvers by providing a less ill-behaved form of the function.

You have an ill-defined equation in the form y - g(y) where y = 1/sqrt(x). So let's get the inverse of g (call it G) so we can solve G(y) - G(g(y)) = G(y) - y instead.

>>> g = e_2.subs(1/x**.5, y)
>>> d = Dummy()
>>> G = solve(g - d, y)[0].subs(d, y)
>>> nsolve(G - y, 6)
6.45497224367903
>>> solve(1/x**.5 - _, dict=True)
{x: 0.024}

The process of rearranging an equation f(x) into form x - g(x) could probably use a built-in method in SymPy, but it's not too hard to do this by replacing each x with a dummy variable, solving for it, and then replacing the dummy symbols with x again. Different g will be more favorable for finding different roots as can be seen in the example below where the purple dashed line is good for finding the root near 1 while the solid blue is better near the smaller root.

Here is a possibility for a "fixed point form" function:

def fixedpoint_Eqs(eq, x=None):
    """rearrange to give eq in form x = g(x)"""
    f = eq.free_symbols
    fp = []
    if x is None:
        assert len(f) == 1, 'must specify x in this case'
        x = list(f)[0]
    Xeq = eq.replace(lambda _:_ == x, lambda _:Dummy())
    X = Xeq.free_symbols - f
    reps = {xi: x for xi in X}
    for xi in X:
        try:
            g = solve(Xeq, xi, dict=True)
            if len(g) != 1:
                raise NotImplementedError
            fp.append(Eq(x, g[0][xi].xreplace(reps)))
        except NotImplementedError:
            pass
    return fp

>>> fixedpoint_Eqs(x+exp(x)+1/x-5)
Eq(x, -1/(x + exp(x) - 5))
Eq(x, -exp(x) + 5 - 1/x)
Eq(x, log(-x + 5 - 1/x))

Solving ill-posed non-linear equations numerically in python/SymPy

2 Answers2