Search code examples
pythonpython-2.7mathematical-optimizationpolynomial-math

Solving system of nonlinear equations with python


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.


Solution

  • 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).