How to solve a pair of nonlinear equations using Python?
Solution 1
for numerical solution, you can use fsolve:
http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html#scipy.optimize.fsolve
from scipy.optimize import fsolve
import math
def equations(p):
x, y = p
return (x+y**2-4, math.exp(x) + x*y - 3)
x, y = fsolve(equations, (1, 1))
print equations((x, y))
Solution 2
If you prefer sympy you can use nsolve.
>>> nsolve([x+y**2-4, exp(x)+x*y-3], [x, y], [1, 1])
[0.620344523485226]
[1.83838393066159]
The first argument is a list of equations, the second is list of variables and the third is an initial guess.
Solution 3
Short answer: use fsolve
As mentioned in other answers the simplest solution to the particular problem you have posed is to use something like fsolve
:
from scipy.optimize import fsolve
from math import exp
def equations(vars):
x, y = vars
eq1 = x+y**2-4
eq2 = exp(x) + x*y - 3
return [eq1, eq2]
x, y = fsolve(equations, (1, 1))
print(x, y)
Output:
0.6203445234801195 1.8383839306750887
Analytic solutions?
You say how to "solve" but there are different kinds of solution. Since you mention SymPy I should point out the biggest difference between what this could mean which is between analytic and numeric solutions. The particular example you have given is one that does not have an (easy) analytic solution but other systems of nonlinear equations do. When there are readily available analytic solutions SymPY can often find them for you:
from sympy import *
x, y = symbols('x, y')
eq1 = Eq(x+y**2, 4)
eq2 = Eq(x**2 + y, 4)
sol = solve([eq1, eq2], [x, y])
Output:
⎡⎛ ⎛ 5 √17⎞ ⎛3 √17⎞ √17 1⎞ ⎛ ⎛ 5 √17⎞ ⎛3 √17⎞ 1 √17⎞ ⎛ ⎛ 3 √13⎞ ⎛√13 5⎞ 1 √13⎞ ⎛ ⎛5 √13⎞ ⎛ √13 3⎞ 1 √13⎞⎤
⎢⎜-⎜- ─ - ───⎟⋅⎜─ - ───⎟, - ─── - ─⎟, ⎜-⎜- ─ + ───⎟⋅⎜─ + ───⎟, - ─ + ───⎟, ⎜-⎜- ─ + ───⎟⋅⎜─── + ─⎟, ─ + ───⎟, ⎜-⎜─ - ───⎟⋅⎜- ─── - ─⎟, ─ - ───⎟⎥
⎣⎝ ⎝ 2 2 ⎠ ⎝2 2 ⎠ 2 2⎠ ⎝ ⎝ 2 2 ⎠ ⎝2 2 ⎠ 2 2 ⎠ ⎝ ⎝ 2 2 ⎠ ⎝ 2 2⎠ 2 2 ⎠ ⎝ ⎝2 2 ⎠ ⎝ 2 2⎠ 2 2 ⎠⎦
Note that in this example SymPy finds all solutions and does not need to be given an initial estimate.
You can evaluate these solutions numerically with evalf
:
soln = [tuple(v.evalf() for v in s) for s in sol]
[(-2.56155281280883, -2.56155281280883), (1.56155281280883, 1.56155281280883), (-1.30277563773199, 2.30277563773199), (2.30277563773199, -1.30277563773199)]
Precision of numeric solutions
However most systems of nonlinear equations will not have a suitable analytic solution so using SymPy as above is great when it works but not generally applicable. That is why we end up looking for numeric solutions even though with numeric solutions: 1) We have no guarantee that we have found all solutions or the "right" solution when there are many. 2) We have to provide an initial guess which isn't always easy.
Having accepted that we want numeric solutions something like fsolve
will normally do all you need. For this kind of problem SymPy will probably be much slower but it can offer something else which is finding the (numeric) solutions more precisely:
from sympy import *
x, y = symbols('x, y')
nsolve([Eq(x+y**2, 4), Eq(exp(x)+x*y, 3)], [x, y], [1, 1])
⎡0.620344523485226⎤
⎢ ⎥
⎣1.83838393066159 ⎦
With greater precision:
nsolve([Eq(x+y**2, 4), Eq(exp(x)+x*y, 3)], [x, y], [1, 1], prec=50)
⎡0.62034452348522585617392716579154399314071550594401⎤
⎢ ⎥
⎣ 1.838383930661594459049793153371142549403114879699 ⎦
Solution 4
Try this one, I assure you that it will work perfectly.
import scipy.optimize as opt
from numpy import exp
import timeit
st1 = timeit.default_timer()
def f(variables) :
(x,y) = variables
first_eq = x + y**2 -4
second_eq = exp(x) + x*y - 3
return [first_eq, second_eq]
solution = opt.fsolve(f, (0.1,1) )
print(solution)
st2 = timeit.default_timer()
print("RUN TIME : {0}".format(st2-st1))
->
[ 0.62034452 1.83838393]
RUN TIME : 0.0009331008900937708
FYI. as mentioned above, you can also use 'Broyden's approximation' by replacing 'fsolve' with 'broyden1'. It works. I did it.
I don't know exactly how Broyden's approximation works, but it took 0.02 s.
And I recommend you do not use Sympy's functions <- convenient indeed, but in terms of speed, it's quite slow. You will see.
Solution 5
An alternative to fsolve
is root
:
import numpy as np
from scipy.optimize import root
def your_funcs(X):
x, y = X
# all RHS have to be 0
f = [x + y**2 - 4,
np.exp(x) + x * y - 3]
return f
sol = root(your_funcs, [1.0, 1.0])
print(sol.x)
This will print
[0.62034452 1.83838393]
If you then check
print(your_funcs(sol.x))
you obtain
[4.4508396968012676e-11, -1.0512035686360832e-11]
confirming that the solution is correct.
Related videos on Youtube
AIB
A programmer, researcher. Celebrating 10million users :) SOreadytohelp
Updated on July 05, 2022Comments
-
AIB almost 2 years
What's the (best) way to solve a pair of non linear equations using Python. (Numpy, Scipy or Sympy)
eg:
- x+y^2 = 4
- e^x+ xy = 3
A code snippet which solves the above pair will be great
-
AIB over 12 yearsyea I know that..I wish to do it in python, because I want to do it repetitively for different sets of equations
-
Blender over 12 yearsYou can
import sage
from any Python script. -
endolith almost 10 yearssage does it by being a wrapper for sympy and maxima, so you could just use those directly.
-
aquirdturtle about 7 yearsWhy is this being down-voted? Seems like a legitimate way to proceed...
-
Sander Heinsalu about 6 yearsI get the error 'name y is not defined' with the code in this answer.
-
Krastanov about 6 years@SanderHeinsalu, just follow what the error message is saying. If "name y is not defined", define it (python can not magically know what you want undefined variables to be). For instance here you want y to be a symbol object you can use to build bigger symbolic objects:
y = Symbol('symbol_name_string')
. Probably you want to keep the same symbol name, soy = Symbol('y')
. -
Edoardo Serra almost 4 yearsDon't know why this isn't the most voted answer, however, is there a way to convert the analytical solutions given by SymPy to a list of approximate numerical values? As of my understanding, the only way to find all of the solutions is via the analytical method, but having those solutions converted could be very useful.
-
Oscar Benjamin almost 4 yearsYou can numerically evaluate any sympy expression that does not have free symbols using
expr.evalf()
: docs.sympy.org/latest/modules/evalf.html -
Oscar Benjamin almost 4 yearsI've added an example with evalf
-
Edoardo Serra almost 4 yearsThank you very much!
-
Edoardo Serra almost 4 yearsYet another question I'm sorry. How can I separate the real from the complex solutions, for example by displaying only the real ones?
-
eric almost 4 years@Edoardo
np.real(x)
andnp.imag(x)
-
Andrew over 2 yearsI get (4.4508396968012676e-11, -1.0512035686360832e-11) as an answer, but this does not work: x+y^2 = 4 != 4.4508396968012676e-11+(-1.0512035686360832e-11)**2 = 4.4508396968123175e-11. equations() returns (0,0) according to what was entered and the original question, so apparently these two small #s are its attempt at that? Also what is the "1, 1" , and where does it come from. Just trying to understand... thanks.
-
creanion over 2 yearsFor sympy, you use the
sympy.real
andsympy.imag
functions, too. (In numpy, I'd just acess .real and .imag as attributes on the relevant values!) -
Oscar Benjamin over 2 yearsIf the symbols are declared as
symbols('x, y', real=True)
then solve will try to return only real solutions. Alternatively you can filter the solutions with something likeReals & set(solve(z**4 - 1, z))
or[s for s in sol if s.is_real]
. Note that sometimes it isn't easy to verify whether or not a solution expression is real. -
billjoie over 2 yearsAn advantage of
root
overfsolve
is the easy ability to specify the resolution method and tolerance. -
Jed over 2 years@Andrew, the output of
equations((x, y))
is the result ofx + y ** 2 - 4
andmath.exp(x) + x * y - 3
. This shows you that the 2 formulas that were set to 0 in the functionequations
are now 0 with the values found forx
andy
. If youprint((x, y))
you'll get the solutions you're looking for.