2

I'm having trouble using Sympy to solve an equation. When I run the code, for example:

print(correction(10))

I expect it to print a number which is f. Instead it gives me ERROR: execution aborted.

def correction(r):

    from sympy import cosh, log, exp, symbols, solve
    f = symbols('f')

    def equation():
        return cosh(((r - 1.0)/(r + 1.0))*(log(2.0)/f)) - 0.5*exp(log(2.0)/f)

    correction = solve(equation(),f)
    return correction

What is the problem?

user
  • 5,370
  • 8
  • 47
  • 75
Emma Smith
  • 381
  • 1
  • 4
  • 9
  • Sound like you are running into some time or memory limit. If I try this with `r != 0`, this just sits, consuming CPU with growing memory use for the backtracker. Or in conclusion: `sympy` searches itself to death. – dhke Dec 14 '16 at 17:21

1 Answers1

4

Your equation is highly non-linear and my guess is that a closed-form solution cannot be found. That's why sympy.solve fails. The only option you have is to solve the equation numerically. Sympy offers the nsolve function for this purpose, which, as typical in numerical solvers, requires an estimate of the solution.

import sympy as sp
r, f = sp.symbols('r, f')
expr = sp.cosh(((r - 1)/(r + 1))*(sp.log(2)/f)) - sp.Rational(1,2)*sp.exp(sp.log(2)/f)

sol = sp.nsolve(expr.subs({r:10}), f, 0.5)
print(sol)

0.699259455239414

Stelios
  • 5,271
  • 1
  • 18
  • 32
  • Clever, as always. And how did you choose 1/2? – Bill Bell Dec 14 '16 at 21:08
  • 1
    @BillBell Simply plotted the function and spotted a zero near that point :) – Stelios Dec 14 '16 at 21:21
  • 1
    Can you please explain what the `sp.Rational` and `expr.subs({r:10})` parts are for? I haven't come across these before. – Emma Smith Dec 15 '16 at 13:12
  • 1
    @EmmaSmith `sp.Rational(a,b)` defines a rational number `a/b` in Sympy. The `subs({r:10})` method substitutes every instance of symbol `r` in `expr` with `10`. Please see the Sympy documentation for further details. Note, however, that these elements of the code are not essential for the answer. – Stelios Dec 15 '16 at 13:29