Solving system of nonlinear equations with python

Question

Can I solve a system of nonlinear equations in terms of parameters in python? Is there a example or tutorial? I can do this easily in maple, but the expressions for my particular system are pretty big and copying them over is quite hard.

Example:

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

You should get the solutions:

[[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]]

The reason I ask: I have a large system of nonlinear ODEs. I want to solve for the fixed points (this is doable, it's been done in maple, but they are large and ugly). I want to create further expressions from the fixed points and then use optimisation package in scipy. I'd rather do it all in python than translate things back and forth since it is very inefficient and mistakes can be made.

Answer 1

Reiterating @Russ's answer, this can be easily accomplished in sympy. For example:

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

where the output format is 3 possible tuples of the possible values for (x, y, z).

Answer 2

Warning I'm a Sage developper, so I might not be neutral.

I don't know how to do that in pure Python, but I would recommend the Sage system whose interface is in Python (actually the command line is a specifically configured IPython) and which allows to do such thing:

+--------------------------------------------------------------------+
| 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]]

It is usually easier to use like this:

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

The main drawback is that Sage is a full distribution of math Software which ships its own Python interpreter (together with a huge bunch of other things written in many language including C/C++, Cython, lisp, fortran) and is notoriously hard to install if you want to use your own interpreter.

A good news for your problem is that Scipy is already shipped with sage.

Answer 3

SymPy might help; I don't know how good it is at solving non-linear equations: http://scipy-lectures.github.io/advanced/sympy.html#id23

You should be able to execute code like (following from theexamples linked above):

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

I haven't tested if it works (I don't have it handy to this machine).