由于文档中的示例已损坏,如何在 SymPy 中以数值方式求解非线性方程组?

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

documentation on nonlinsolve给出了这个例子:

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

但即使在他们网站上的直播 shell 中,也会引发错误:

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

如何使用 nonlinsolve 对方程组进行数值求解?我知道我可以使用 ufuncify 将方程式转换为 scipy.optimize.fsolve 可以求解的系统,但我宁愿避免使用那几行样板文件,而是直接使用 SymPy。

根据 SymPy documentation on solve, using solve is not recommended. For nonlinear systems of equations, the documentation recommends sympy.solvers.solveset.nonlinsolve,这就是我在这里尝试使用的内容。

似乎此模块未包含在 1.0.0 版本中,但仅在当前 head 中可用。无论如何,检查 sympy 回购的当前状态可能是个好主意,因为它们的版本相当罕见。

根据docs,这很简单

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

(我不得不使用 python3 而不是 python 因为我同时安装了 Python 2 和 Python 3 系统范围的安装。您也可以使用 virtualenv 以避免 sudo.)

之后,最新的github版本将与普通的import

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

现在示例有效:

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。如果您真的对方程的 numeric 解感兴趣,那么 ufuncify + fsolve 似乎更适合您。 SymPy 是一个计算机代数系统,因此 nonlinsolve 将尽最大努力 分析 而非数值求解非线性系统。当然它可能会失败(与任何解析求解器一样):

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

如果您想对系统进行数值求解,请使用 nsolve。它需要对解决方案进行初步猜测(您还可以传递许多选项来使用不同的求解器,请参阅 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']])

对于符号解决方案,我建议使用旧的 solve,直到 nonlinsolve 成熟。