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.

Lepidus

use sympy.solve([x*y - 1, 4*x**2 + y**2 - 5], [x, y]) – Sirreverence

@Sirreverence According to the Sympy documentation, using solve is not recommended. I'll add a statement about that to my question. – Lepidus

The docs are perhaps a bit aggressive about that. solveset is indeed going to replace solve, but there are still instances where it is necessary, such as nonlinear systems. – Rafat

Answer 1

It seems that this module is not included into 1.0.0 release but is available only in current head. Anyway, it may be good idea to check out current state of sympy repo as their releases are rather rare.

According to the docs, this is as simple as

git clone git://github.com/sympy/sympy.git
cd sympy
sudo python setupegg.py develop

(I had to use python3 instead of python as I have both Python 2 and Python 3 system-wide installations. You also possibly can use virtualenv to avoid sudo.)

After that, latest github version will be available with ordinary import.

import sympy
print(sympy.__version__)
# 1.0.1.dev

Now the example works:

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)}

P.S. If you are really interested in numeric solution of your equation, it seems that ufuncify + fsolve thing will suit you better. SymPy is a computer algebra system, so nonlinsolve will do its best to solve nonlinear system analytically, not numerically. And of course it can fail (as with any analytic solvers):

nonlinsolve([x**5 + x**2 + 1], [x])
# {(CRootOf(x**5 + x**2 + 1, 0),), (CRootOf(x**5 + x**2 + 1, 1),), (CRootOf(x**5 + x**2 + 1, 2),), (CRootOf(x**5 + x**2 + 1, 3),), (CRootOf(x**5 + x**2 + 1, 4),)}

Enact

For Windows users, go to github.com/sympy/sympy, download the zip, expand it in a convenient new folder, CD inside the folder that contains setup.py and then execute python setup.py install. Prepare for a lengthy wait whilst this process bumps and grinds its way to a conclusion. Then check that you can import sympy and that sympy.__version__ is what you expected. – Polychasium

Answer 2

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.

Rafat