Since the example in the documentation is broken, how do I solve a non-linear system of equations numerically in SymPy?

Question

The documentation on nonlinsolve gives this example:

from sympy.core.symbol import symbols
from sympy.solvers.solveset import nonlinsolve
x, y, z = symbols('x, y, z', real=True)
nonlinsolve([x*y - 1, 4*x**2 + y**2 - 5], [x, y])
{(-1, -1), (-1/2, -2), (1/2, 2), (1, 1)}

but even in the live shell on their website, that throws an error:

>>> from sympy.solvers.solveset import nonlinsolve
Traceback (most recent call last):
  File "<string>", line 1, in <module>
ImportError: cannot import name nonlinsolve

How can I use nonlinsolve to solve a system of equations numerically? I know I can use ufuncify to convert the equations into a system that scipy.optimize.fsolve can solve, but I would rather avoid those couple of lines of boilerplate and just use SymPy directly.

According to the SymPy documentation on solve, using solve is not recommended. For nonlinear systems of equations, the documentation recommends sympy.solvers.solveset.nonlinsolve, which is what I'm trying to use here.

Answer 1

If you want to solve the systems numerically, use nsolve. It requires an initial guess for the solution (there are also many options you can pass to use different solvers, see http://docs.sympy.org/latest/modules/solvers/solvers.html#sympy.solvers.solvers.nsolve and http://mpmath.org/doc/current/calculus/optimization.html).

In [1]: nsolve([x*y - 1, 4*x**2 + y**2 - 5], [x, y], [1, 1])
Out[1]:
matrix(
[['1.0'],
 ['1.0']])

For symbolic solutions, I would recommend using the old solve, until nonlinsolve matures.