Search code examples
pythonsympysolver

How to use SymPy to solve a complex equation with an approximate solution


For example,

>>> integrate(sqrt(sin(u)*sin(u)+1), (u, 0, b)).subs(b, 0.22).evalf()
0.221745186045595

But I want to know reversedly which b can get 0.221745186045595. So I write

>>> solve(integrate(sqrt(sin(u)*sin(u)+1), (u, 0, b)) - 0.221745186045595, b)
[]

I know we can not get a very precise solution, so my question is: How can we set SymPy's solve to do that with some tolerated precision?

The function sqrt(sin(u)*sin(u)+1) is just an example. If possible, it should be an unpredictable, user-input function.


Solution

  • This is not what SymPy is for. "Sym" in SymPy means Symbolic, as opposed to Numeric. You want numeric computations. Use SciPy quad and some root-finding routine like root or fsolve. For example:

    import numpy as np
    from scipy import integrate, optimize
    target = 0.221745186045595
    f = lambda u: np.sqrt(np.sin(u)**2 + 1)
    x = optimize.root(lambda b: integrate.quad(f, 0, b)[0] - target, 0).x
    

    returns x as array([0.22]).

    Parsing user input

    For turning user input into a callable function like f above, SymPy's lambdify can be used. Example:

    from sympy import sympify, lambdify     
    f_string = "sqrt(sin(u)**2+1)"      # user input 
    f_expr = sympify(f_string)
    sym = next(iter(f_expr.free_symbols))
    f = lambdify(sym, f_expr, "numpy")
    

    Here f_expr is a SymPy expression parsed from a string, sym is the SymPy symbol (the argument of the function), and f is a Python function created by lambdify. This f is then used as above.