用Python求解非线性方程组

Question

我可以在Python中用参数来解决一组非线性方程吗？有没有示例或者教程？我在Maple中可以很轻松地做到这一点，但我这个特定系统的表达式比较复杂，复制过来很麻烦。

举个例子：

sigma*(y-x) = 0
x*(rho-z)-y = 0
x*y-beta*z = 0

你应该能得到这些解：

[[x = 0, y = 0, z = 0], [x = sqrt(beta*rho-beta), y = sqrt(beta*rho-beta), z = rho-1],
[x = -sqrt(beta*rho-beta), y = -sqrt(beta*rho-beta), z = rho-1]]

我之所以问这个，是因为我有一个很大的非线性常微分方程组。我想找出固定点（这个是可以做到的，之前在Maple中做过，但结果很复杂且难看）。我想根据这些固定点创建更多的表达式，然后使用SciPy中的优化工具包。我更希望能全部在Python中完成，而不是来回翻译，因为这样效率低下，而且容易出错。

Answer 1

SymPy 可能会有帮助；我不太确定它在解决非线性方程方面的效果如何：http://scipy-lectures.github.io/advanced/sympy.html#id23

你应该能够运行类似下面的代码（这是从上面链接的例子中来的）：

from sympy import *
x = Symbol('x')
y = Symbol('y')
z = Symbol('z')
beta = Symbol('beta')
rho = Symbol('rho')
sigma = Symbol('sigma')

solve([sigma*(y-x), x*(rho-z)-y, x*y-beta*z], [x, y, z])

我还没测试过它是否有效（我在这台机器上没有安装它）。

Answer 2

警告 我是Sage的开发者，所以我可能会有点偏见。

我不知道怎么用纯Python来做到这一点，但我推荐你使用Sage系统，它的界面是用Python写的（实际上，它的命令行是特别配置过的IPython），可以做到这样的事情：

+--------------------------------------------------------------------+
| Sage Version 5.10, Release Date: 2013-06-17                        |
| Type "notebook()" for the browser-based notebook interface.        |
| Type "help()" for help.                                            |
+--------------------------------------------------------------------+
sage: var("sigma y x rho beta z")
(sigma, y, x, rho, beta, z)
sage: sys = [sigma*(y-x), x*(rho-z)-y, x*y-beta*z]
sage: solve(sys, x, y, z)
[[x == sqrt(beta*rho - beta), y == (beta*rho - beta)/(sqrt(rho - 1)*sqrt(beta)), z == rho - 1], [x == -sqrt(beta*rho - beta), y == -(beta*rho - beta)/(sqrt(rho - 1)*sqrt(beta)), z == rho - 1], [x == 0, y == 0, z == 0]]

通常这样使用会更简单：

sage: solve(sys, x, y, z, solution_dict=True)
[{z: rho - 1,
  x: sqrt(beta*rho - beta),
  y: (beta*rho - beta)/(sqrt(rho - 1)*sqrt(beta))},
 {z: rho - 1,
  x: -sqrt(beta*rho - beta),
  y: -(beta*rho - beta)/(sqrt(rho - 1)*sqrt(beta))},
 {z: 0, x: 0, y: 0}]

主要的缺点是，Sage是一个完整的数学软件包，它自带了自己的Python解释器（还有很多其他用不同语言写的东西，包括C/C++、Cython、lisp、fortran），如果你想用自己的解释器，安装起来会非常麻烦。

好消息是，Scipy已经和Sage一起打包好了。

Answer 3

重申一下@Russ的回答，这个可以很简单地用 sympy 来实现。比如说：

In [1]: import sympy as sp
In [2]: x, y, z = sp.symbols('x, y, z')
In [3]: rho, sigma, beta = sp.symbols('rho, sigma, beta')
In [4]: f1 = sigma * (y - x)
In [5]: f2 = x * (rho - z) - y
In [6]: f3 = x * y - beta * z
In [7]: sp.solvers.solve((f1, f2, f3), (x, y, z))
Out[7]: 
[(0, 0, 0),
 (-sqrt(beta*rho - beta), -sqrt(beta*(rho - 1)), rho - 1),
 (sqrt(beta*rho - beta), sqrt(beta*(rho - 1)), rho - 1)]

这里的输出格式是三个可能的元组，表示 (x, y, z) 的可能值。