This is the official documentation of the solveset module in solvers. It contains the frequently asked questions about our new module to solve equations.


For a beginner-friendly guide focused on solving common types of equations, refer to Solve Equations.

What’s wrong with solve():

SymPy already has a pretty powerful solve function. But it has some deficiencies. For example:

  1. It doesn’t have a consistent output for various types of solutions It needs to return a lot of types of solutions consistently:

    • Single solution : \(x = 1\)

    • Multiple solutions: \(x^2 = 1\)

    • No Solution: \(x^2 + 1 = 0 ; x \in \mathbb{R}\)

    • Interval of solution: \(\lfloor x \rfloor = 0\)

    • Infinitely many solutions: \(\sin(x) = 0\)

    • Multivariate functions with point solutions: \(x^2 + y^2 = 0\)

    • Multivariate functions with non-point solution: \(x^2 + y^2 = 1\)

    • System of equations: \(x + y = 1\) and \(x - y = 0\)

    • Relational: \(x > 0\)

    • And the most important case: “We don’t Know”

  2. The input API has a lot of parameters and it can be difficult to use.

  3. There are cases like finding the maxima and minima of function using critical points where it is important to know if it has returned all the solutions. solve does not guarantee this.

Why Solveset?

  • solveset has an alternative consistent input and output interface: solveset returns a set object and a set object takes care of all types of output. For cases where it does not “know” all the solutions a ConditionSet with a partial solution is returned. For input it only takes the equation, the variables to solve for and the optional argument domain over which the equation is to be solved.

  • solveset can return infinitely many solutions. For example solving for \(\sin{(x)} = 0\) returns \(\{2 n \pi | n \in \mathbb{Z}\} \cup \{2 n \pi + \pi | n \in \mathbb{Z}\}\), whereas solve only returns \([0, \pi]\).

  • There is a clear code level and interface level separation between solvers for equations in the complex domain and the real domain. For example solving \(e^x = 1\) when \(x\) is to be solved in the complex domain, returns the set of all solutions, that is \(\{2 n i \pi | n \in \mathbb{Z}\}\), whereas if \(x\) is to be solved in the real domain then only \(\{0\}\) is returned.

Why do we use Sets as an output type?

SymPy has a well developed sets module, which can represent most of the set containers in mathematics such as:

  • FiniteSet

    Represents a finite set of discrete numbers.

  • Interval

    Represents a real interval as a set.

  • ProductSet

    Represents a Cartesian product of sets.

  • ImageSet

    Represents the image of a set under a mathematical function

    >>> from sympy import ImageSet, S, Lambda
    >>> from import x
    >>> squares = ImageSet(Lambda(x, x**2), S.Naturals)  # {x**2 for x in N}
    >>> 4 in squares
  • ComplexRegion

    Represents the set of all complex numbers in a region in the Argand plane.

  • ConditionSet

    Represents the set of elements, which satisfies a given condition.

Also, the predefined set classes such as:

  • Naturals, $mathbb{N}

    Represents the natural numbers (or counting numbers), which are all positive integers starting from 1.

  • Naturals0, \(\mathbb{N_0}\)

    Represents the whole numbers, which are all the non-negative integers, inclusive of 0.

  • Integers, \(\mathbb{Z}\)

    Represents all integers: positive, negative and zero.

  • Reals, \(\mathbb{R}\)

    Represents the set of all real numbers.

  • Complexes, \(\mathbb{C}\)

    Represents the set of all complex numbers.

  • EmptySet, \(\emptyset\)

    Represents the empty set.

The above six sets are available as Singletons, like S.Integers.

It is capable of most of the set operations in mathematics:

  • Union

  • Intersection

  • Complement

  • SymmetricDifference

The main reason for using sets as output to solvers is that it can consistently represent many types of solutions. For the single variable case it can represent:

  • No solution (by the empty set).

  • Finitely many solutions (by FiniteSet).

  • Infinitely many solutions, both countably and uncountably infinite solutions (using the ImageSet module).

  • Interval

  • There can also be bizarre solutions to equations like the set of rational numbers.

No other Python object (list, dictionary, generator, Python sets) provides the flexibility of mathematical sets which our sets module tries to emulate. The second reason to use sets is that they are close to the entities which mathematicians deal with and it makes it easier to reason about them. Set objects conform to Pythonic conventions when possible, i.e., x in A and for i in A both work when they can be computed. Another advantage of using objects closer to mathematical entities is that the user won’t have to “learn” our representation and she can have her expectations transferred from her mathematical experience.

For the multivariate case we represent solutions as a set of points in a n-dimensional space and a point is represented by a FiniteSet of ordered tuples, which is a point in \(\mathbb{R}^n\) or \(\mathbb{C}^n\).

Please note that, the general FiniteSet is unordered, but a FiniteSet with a tuple as its only argument becomes ordered, since a tuple is ordered. So the order in the tuple is mapped to a pre-defined order of variables while returning solutions.

For example:

>>> from sympy import FiniteSet
>>> FiniteSet(1, 2, 3)   # Unordered
{1, 2, 3}
>>> FiniteSet((1, 2, 3))  # Ordered
{(1, 2, 3)}

Why not use dicts as output?

Dictionary are easy to deal with programmatically but mathematically they are not very precise and use of them can quickly lead to inconsistency and a lot of confusion. For example:

  • There are a lot of cases where we don’t know the complete solution and we may like to output a partial solution, consider the equation \(fg = 0\). The solution of this equation is the union of the solution of the following two equations: \(f = 0\), \(g = 0\). Let’s say that we are able to solve \(f = 0\) but solving \(g = 0\) isn’t supported yet. In this case we cannot represent partial solution of the given equation \(fg = 0\) using dicts. This problem is solved with sets using a ConditionSet object:

    \(sol_f \cup \{x | x ∊ \mathbb{R} ∧ g = 0\}\), where \(sol_f\) is the solution of the equation \(f = 0\).

  • Using a dict may lead to surprising results like:

    • solve(Eq(x**2, 1), x) != solve(Eq(y**2, 1), y)

      Mathematically, this doesn’t make sense. Using FiniteSet here solves the problem.

  • It also cannot represent solutions for equations like \(|x| < 1\), which is a disk of radius 1 in the Argand Plane. This problem is solved using complex sets implemented as ComplexRegion.

Input API of solveset

solveset has simpler input API, unlike solve. It takes a maximum of three arguments:

solveset(equation, variable=None, domain=S.Complexes)


The equation to solve.


The variable for which the equation is to be solved.


The domain in which the equation is to be solved.

solveset removes the flags argument of solve, which had made the input API more complicated and output API inconsistent.

What is this domain argument about?

Solveset is designed to be independent of the assumptions on the variable being solved for and instead, uses the domain argument to decide the solver to dispatch the equation to, namely solveset_real or solveset_complex. It’s unlike the old solve which considers the assumption on the variable.

>>> from sympy import solveset, S
>>> from import x
>>> solveset(x**2 + 1, x) # domain=S.Complexes is default
{-I, I}
>>> solveset(x**2 + 1, x, domain=S.Reals)

What are the general methods employed by solveset to solve an equation?

Solveset uses various methods to solve an equation, here is a brief overview of the methodology:

  • The domain argument is first considered to know the domain in which the user is interested to get the solution.

  • If the given function is a relational (>=, <=, >, <), and the domain is real, then solve_univariate_inequality and solutions are returned. Solving for complex solutions of inequalities, like \(x^2 < 0\) is not yet supported.

  • Based on the domain, the equation is dispatched to one of the two functions solveset_real or solveset_complex, which solves the given equation in the complex or real domain, respectively.

  • If the given expression is a product of two or more functions, like say \(gh = 0\), then the solution to the given equation is the Union of the solution of the equations \(g = 0\) and \(h = 0\), if and only if both \(g\) and \(h\) are finite for a finite input. So, the solution is built up recursively.

  • If the function is trigonometric or hyperbolic, the function _solve_real_trig is called, which solves it by converting it to complex exponential form.

  • The function is now checked if there is any instance of a Piecewise expression, if it is, then it’s converted to explicit expression and set pairs and then solved recursively.

  • The respective solver now tries to invert the equation using the routines invert_real and invert_complex. These routines are based on the concept of mathematical inverse (though not exactly). It reduces the real/complex valued equation \(f(x) = y\) to a set of equations: \(\{g(x) = h_1(y), g(x) = h_2(y), ..., g(x) = h_n(y) \}\) where \(g(x)\) is a simpler function than \(f(x)\). There is some work needed to be done in this to find invert of more complex expressions.

  • After the invert, the equations are checked for radical or Abs (Modulus), then the method _solve_radical tries to simplify the radical, by removing it using techniques like squaring, cubing etc, and _solve_abs solves nested Modulus by considering the positive and negative variants, iteratively.

  • If none of the above method is successful, then methods of polynomial is used as follows:

    • The method to solve the rational function, _solve_as_rational, is called. Based on the domain, the respective poly solver _solve_as_poly_real or _solve_as_poly_complex is called to solve f as a polynomial.

    • The underlying method _solve_as_poly solves the equation using polynomial techniques if it’s already a polynomial equation or, with a change of variables, can be made so.

  • The final solution set returned by solveset is the intersection of the set of solutions found above and the input domain.

How do we manipulate and return an infinite solution?

  • In the real domain, we use our ImageSet class in the sets module to return infinite solutions. ImageSet is an image of a set under a mathematical function. For example, to represent the solution of the equation \(\sin{(x)} = 0\), we can use the ImageSet as:

    >>> from sympy import ImageSet, Lambda, pi, S, Dummy, pprint
    >>> n = Dummy('n')
    >>> pprint(ImageSet(Lambda(n, 2*pi*n), S.Integers), use_unicode=True)
    {2⋅n⋅π │ n ∊ ℤ}

    Where n is a dummy variable. It is basically the image of the set of integers under the function \(2\pi n\).

  • In the complex domain, we use complex sets, which are implemented as the ComplexRegion class in the sets module, to represent infinite solution in the Argand plane. For example to represent the solution of the equation \(|z| = 1\), which is a unit circle, we can use the ComplexRegion as:

    >>> from sympy import ComplexRegion, FiniteSet, Interval, pi, pprint
    >>> pprint(ComplexRegion(FiniteSet(1)*Interval(0, 2*pi), polar=True), use_unicode=True)
    {r⋅(ⅈ⋅sin(θ) + cos(θ)) │ r, θ ∊ {1} × [0, 2⋅π)}

    Where the FiniteSet in the ProductSet is the range of the value of \(r\), which is the radius of the circle and the Interval is the range of \(\theta\), the angle from the \(x\) axis representing a unit circle in the Argand plane.

    Note: We also have non-polar form notation for representing solution in rectangular form. For example, to represent first two quadrants in the Argand plane, we can write the ComplexRegion as:

    >>> from sympy import ComplexRegion, Interval, pi, oo, pprint
    >>> pprint(ComplexRegion(Interval(-oo, oo)*Interval(0, oo)), use_unicode=True)
    {x + y⋅ⅈ │ x, y ∊ (-∞, ∞) × [0, ∞)}

    where the Intervals are the range of \(x\) and \(y\) for the set of complex numbers \(x + iy\).

How does solveset ensure that it is not returning any wrong solution?

Solvers in a Computer Algebra System are based on heuristic algorithms, so it’s usually very hard to ensure 100% percent correctness, in every possible case. However there are still a lot of cases where we can ensure correctness. Solveset tries to verify correctness wherever it can. For example:

Consider the equation \(|x| = n\). A naive method to solve this equation would return {-n, n} as its solution, which is not correct since {-n, n} can be its solution if and only if n is positive. Solveset returns this information as well to ensure correctness.

>>> from sympy import symbols, S, pprint, solveset
>>> x, n = symbols('x, n')
>>> pprint(solveset(abs(x) - n, x, domain=S.Reals), use_unicode=True)
{x │ x ∊ {-n, n} ∧ (n ∈ [0, ∞))}

Though, there still a lot of work needs to be done in this regard.

Search based solver and step-by-step solution

Note: This is under Development.

After the introduction of ConditionSet, the solving of equations can be seen as set transformations. Here is an abstract view of the things we can do to solve equations.

  • Apply various set transformations on the given set.

  • Define a metric of the usability of solutions, or a notion of some solutions being better than others.

  • Different transformations would be the nodes of a tree.

  • Suitable searching techniques could be applied to get the best solution.

ConditionSet gives us the ability to represent unevaluated equations and inequalities in forms like \(\{x|f(x)=0; x \in S\}\) and \(\{x|f(x)>0; x \in S\}\) but a more powerful thing about ConditionSet is that it allows us to write the intermediate steps as set to set transformation. Some of the transformations are:

  • Composition: \(\{x|f(g(x))=0;x \in S\} \Rightarrow \{x|g(x)=y; x \in S, y \in \{z|f(z)=0; z \in S\}\}\)

  • Polynomial Solver: \(\{x | P(x) = 0;x \in S\} \Rightarrow \{x_1,x_2, ... ,x_n\} \cap S\),

    where \(x_i\) are roots of \(P(x)\).

  • Invert solver: \(\{x|f(x)=0;x \in S\} \Rightarrow \{g(0)| \text{ all g such that } f(g(x)) = x\}\)

  • logcombine: \(\{x| \log(f(x)) + \log(g(x));x \in S\}\)

    \(\Rightarrow \{x| \log(f(x).g(x)); x \in S\} \text{ if } f(x) > 0 \text{ and } g(x) > 0\) \(\Rightarrow \{x| \log(f(x)) + \log(g(x));x \in S\} \text{ otherwise}\)

  • product solve: \(\{x|f(x)g(x)=0; x \in S\}\)

    \(\Rightarrow \{x|f(x)=0; x \in S\} U \{x|g(x)=0; x \in S\}\) \(\text{ given } f(x) \text{ and } g(x) \text{ are bounded.}\) \(\Rightarrow \{x|f(x)g(x)=0; x \in S\}, \text{ otherwise}\)

Since the output type is same as the input type any composition of these transformations is also a valid transformation. And our aim is to find the right sequence of compositions (given the atoms) which transforms the given condition set to a set which is not a condition set i.e., FiniteSet, Interval, Set of Integers and their Union, Intersection, Complement or ImageSet. We can assign a cost function to each set, such that, the more desirable that form of set is to us, the less the value of the cost function. This way our problem is now reduced to finding the path from the initial ConditionSet to the lowest valued set on a graph where the atomic transformations forms the edges.

How do we deal with cases where only some of the solutions are known?

Creating a universal equation solver, which can solve each and every equation we encounter in mathematics is an ideal case for solvers in a Computer Algebra System. When cases which are not solved or can only be solved incompletely, a ConditionSet is used and acts as an unevaluated solveset object.

Note that, mathematically, finding a complete set of solutions for an equation is undecidable. See Richardson’s theorem.

ConditionSet is basically a Set of elements which satisfy a given condition. For example, to represent the solutions of the equation in the real domain:

\[(x^2 - 4)(\sin(x) + x)\]

We can represent it as:

\(\{-2, 2\} ∪ \{x | x \in \mathbb{R} ∧ x + \sin(x) = 0\}\)

What is the plan for solve and solveset?

There are still a few things solveset can’t do, which solve can, such as solving nonlinear multivariate & LambertW type equations. Hence, it’s not yet a perfect replacement for solve. As the algorithms in solveset mature, solveset may be able to be used within solve to replace some of its algorithms.

How are symbolic parameters handled in solveset?

Solveset is in its initial phase of development, so the symbolic parameters aren’t handled well for all the cases, but some work has been done in this regard to depict our ideology towards symbolic parameters. As an example, consider the solving of \(|x| = n\) for real \(x\), where \(n\) is a symbolic parameter. Solveset returns the value of \(x\) considering the domain of the symbolic parameter \(n\) as well:

\[([0, \infty) \cap \{n\}) \cup ((-\infty, 0] \cap \{-n\}).\]

This simply means \(n\) is the solution only when it belongs to the Interval \([0, \infty)\) and \(-n\) is the solution only when \(-n\) belongs to the Interval \((- \infty, 0]\).

There are other cases to address too, like solving \(2^x + (a - 2)\) for \(x\) where \(a\) is a symbolic parameter. As of now, It returns the solution as an intersection with \(\mathbb{R}\), which is trivial, as it doesn’t reveal the domain of \(a\) in the solution.

Recently, we have also implemented a function to find the domain of the expression in a FiniteSet (Intersection with the interval) in which it is not-empty. It is a useful addition for dealing with symbolic parameters. For example:

>>> from sympy import Symbol, FiniteSet, Interval, not_empty_in, sqrt, oo
>>> from import x
>>> not_empty_in(FiniteSet(x/2).intersect(Interval(0, 1)), x)
Interval(0, 2)
>>> not_empty_in(FiniteSet(x, x**2).intersect(Interval(1, 2)), x)
Union(Interval(1, 2), Interval(-sqrt(2), -1))


Solveset Module Reference

Use solveset() to solve equations or expressions (assumed to be equal to 0) for a single variable. Solving an equation like \(x^2 == 1\) can be done as follows:

>>> from sympy import solveset
>>> from sympy import Symbol, Eq
>>> x = Symbol('x')
>>> solveset(Eq(x**2, 1), x)
{-1, 1}

Or one may manually rewrite the equation as an expression equal to 0:

>>> solveset(x**2 - 1, x)
{-1, 1}

The first argument for solveset() is an expression (equal to zero) or an equation and the second argument is the symbol that we want to solve the equation for.

sympy.solvers.solveset.solveset(f, symbol=None, domain=Complexes)[source]

Solves a given inequality or equation with set as output


f : Expr or a relational.

The target equation or inequality

symbol : Symbol

The variable for which the equation is solved

domain : Set

The domain over which the equation is solved



A set of values for \(symbol\) for which \(f\) is True or is equal to zero. An EmptySet is returned if \(f\) is False or nonzero. A ConditionSet is returned as unsolved object if algorithms to evaluate complete solution are not yet implemented.

solveset claims to be complete in the solution set that it returns.



The algorithms to solve inequalities in complex domain are not yet implemented.


The input is not valid.


It is a bug, please report to the github issue tracker.


Python interprets 0 and 1 as False and True, respectively, but in this function they refer to solutions of an expression. So 0 and 1 return the domain and EmptySet, respectively, while True and False return the opposite (as they are assumed to be solutions of relational expressions).


>>> from sympy import exp, sin, Symbol, pprint, S, Eq
>>> from sympy.solvers.solveset import solveset, solveset_real
  • The default domain is complex. Not specifying a domain will lead to the solving of the equation in the complex domain (and this is not affected by the assumptions on the symbol):

>>> x = Symbol('x')
>>> pprint(solveset(exp(x) - 1, x), use_unicode=False)
{2*n*I*pi | n in Integers}
>>> x = Symbol('x', real=True)
>>> pprint(solveset(exp(x) - 1, x), use_unicode=False)
{2*n*I*pi | n in Integers}
  • If you want to use solveset to solve the equation in the real domain, provide a real domain. (Using solveset_real does this automatically.)

>>> R = S.Reals
>>> x = Symbol('x')
>>> solveset(exp(x) - 1, x, R)
>>> solveset_real(exp(x) - 1, x)

The solution is unaffected by assumptions on the symbol:

>>> p = Symbol('p', positive=True)
>>> pprint(solveset(p**2 - 4))
{-2, 2}

When a ConditionSet is returned, symbols with assumptions that would alter the set are replaced with more generic symbols:

>>> i = Symbol('i', imaginary=True)
>>> solveset(Eq(i**2 + i*sin(i), 1), i, domain=S.Reals)
ConditionSet(_R, Eq(_R**2 + _R*sin(_R) - 1, 0), Reals)
  • Inequalities can be solved over the real domain only. Use of a complex domain leads to a NotImplementedError.

>>> solveset(exp(x) > 1, x, R), oo)

See also


solver for real domain


solver for complex domain

sympy.solvers.solveset.solveset_real(f, symbol)[source]
sympy.solvers.solveset.solveset_complex(f, symbol)[source]
sympy.solvers.solveset.invert_real(f_x, y, x)[source]

Inverts a real-valued function. Same as invert_complex(), but sets the domain to S.Reals before inverting.

sympy.solvers.solveset.invert_complex(f_x, y, x, domain=Complexes)[source]

Reduce the complex valued equation \(f(x) = y\) to a set of equations

\[\left\{g(x) = h_1(y),\ g(x) = h_2(y),\ \dots,\ g(x) = h_n(y) \right\}\]

where \(g(x)\) is a simpler function than \(f(x)\). The return value is a tuple \((g(x), \mathrm{set}_h)\), where \(g(x)\) is a function of \(x\) and \(\mathrm{set}_h\) is the set of function \(\left\{h_1(y), h_2(y), \dots, h_n(y)\right\}\). Here, \(y\) is not necessarily a symbol.

\(\mathrm{set}_h\) contains the functions, along with the information about the domain in which they are valid, through set operations. For instance, if \(y = |x| - n\) is inverted in the real domain, then \(\mathrm{set}_h\) is not simply \(\{-n, n\}\) as the nature of \(n\) is unknown; rather, it is:

$$ left(left[0, inftyright) cap left{nright}right) cup

left(left(-infty, 0right] cap left{- nright}right)$$

By default, the complex domain is used which means that inverting even seemingly simple functions like \(\exp(x)\) will give very different results from those obtained in the real domain. (In the case of \(\exp(x)\), the inversion via \(\log\) is multi-valued in the complex domain, having infinitely many branches.)

If you are working with real values only (or you are not sure which function to use) you should probably set the domain to S.Reals (or use invert_real which does that automatically).


>>> from sympy.solvers.solveset import invert_complex, invert_real
>>> from import x, y
>>> from sympy import exp

When does exp(x) == y?

>>> invert_complex(exp(x), y, x)
(x, ImageSet(Lambda(_n, I*(2*_n*pi + arg(y)) + log(Abs(y))), Integers))
>>> invert_real(exp(x), y, x)
(x, Intersection({log(y)}, Reals))

When does exp(x) == 1?

>>> invert_complex(exp(x), 1, x)
(x, ImageSet(Lambda(_n, 2*_n*I*pi), Integers))
>>> invert_real(exp(x), 1, x)
(x, {0})
sympy.solvers.solveset.domain_check(f, symbol, p)[source]

Returns False if point p is infinite or any subexpression of f is infinite or becomes so after replacing symbol with p. If none of these conditions is met then True will be returned.


>>> from sympy import Mul, oo
>>> from import x
>>> from sympy.solvers.solveset import domain_check
>>> g = 1/(1 + (1/(x + 1))**2)
>>> domain_check(g, x, -1)
>>> domain_check(x**2, x, 0)
>>> domain_check(1/x, x, oo)
  • The function relies on the assumption that the original form of the equation has not been changed by automatic simplification.

>>> domain_check(x/x, x, 0) # x/x is automatically simplified to 1
  • To deal with automatic evaluations use evaluate=False:

>>> domain_check(Mul(x, 1/x, evaluate=False), x, 0)
sympy.solvers.solveset.solvify(f, symbol, domain)[source]

Solves an equation using solveset and returns the solution in accordance with the \(solve\) output API.


We classify the output based on the type of solution returned by \(solveset\).



A ConditionSet is the input.

Solution | Output

FiniteSet | list

ImageSet, | list (if \(f\) is periodic) Union |

Union | list (with FiniteSet)

EmptySet | empty list

Others | None


>>> from sympy.solvers.solveset import solvify
>>> from import x
>>> from sympy import S, tan, sin, exp
>>> solvify(x**2 - 9, x, S.Reals)
[-3, 3]
>>> solvify(sin(x) - 1, x, S.Reals)
>>> solvify(tan(x), x, S.Reals)
>>> solvify(exp(x) - 1, x, S.Complexes)
>>> solvify(exp(x) - 1, x, S.Reals)
sympy.solvers.solveset.linear_eq_to_matrix(equations, *symbols)[source]

Converts a given System of Equations into Matrix form. Here \(equations\) must be a linear system of equations in \(symbols\). Element M[i, j] corresponds to the coefficient of the jth symbol in the ith equation.

The Matrix form corresponds to the augmented matrix form. For example:

\[4x + 2y + 3z = 1\]
\[3x + y + z = -6\]
\[2x + 4y + 9z = 2\]

This system will return \(A\) and \(b\) as:

$$ A = left[begin{array}{ccc}

4 & 2 & 3 \ 3 & 1 & 1 \ 2 & 4 & 9 end{array}right] b = left[begin{array}{c} 1 \ -6 \ 2 end{array}right] $$

The only simplification performed is to convert Eq(a, b) \(\Rightarrow a - b\).



The equations contain a nonlinear term.


The symbols are not given or are not unique.


>>> from sympy import linear_eq_to_matrix, symbols
>>> c, x, y, z = symbols('c, x, y, z')

The coefficients (numerical or symbolic) of the symbols will be returned as matrices:

>>> eqns = [c*x + z - 1 - c, y + z, x - y]
>>> A, b = linear_eq_to_matrix(eqns, [x, y, z])
>>> A
[c,  0, 1],
[0,  1, 1],
[1, -1, 0]])
>>> b
[c + 1],
[    0],
[    0]])

This routine does not simplify expressions and will raise an error if nonlinearity is encountered:

>>> eqns = [
...     (x**2 - 3*x)/(x - 3) - 3,
...     y**2 - 3*y - y*(y - 4) + x - 4]
>>> linear_eq_to_matrix(eqns, [x, y])
Traceback (most recent call last):
symbol-dependent term can be ignored using `strict=False`

Simplifying these equations will discard the removable singularity in the first and reveal the linear structure of the second:

>>> [e.simplify() for e in eqns]
[x - 3, x + y - 4]

Any such simplification needed to eliminate nonlinear terms must be done before calling this routine.

sympy.solvers.solveset.linsolve(system, *symbols)[source]

Solve system of \(N\) linear equations with \(M\) variables; both underdetermined and overdetermined systems are supported. The possible number of solutions is zero, one or infinite. Zero solutions throws a ValueError, whereas infinite solutions are represented parametrically in terms of the given symbols. For unique solution a FiniteSet of ordered tuples is returned.

All standard input formats are supported: For the given set of equations, the respective input types are given below:

\[3x + 2y - z = 1\]
\[2x - 2y + 4z = -2\]
\[2x - y + 2z = 0\]
  • Augmented matrix form, system given below:

$$ text{system} = left[{array}{cccc}

3 & 2 & -1 & 1\ 2 & -2 & 4 & -2\ 2 & -1 & 2 & 0 end{array}right] $$

system = Matrix([[3, 2, -1, 1], [2, -2, 4, -2], [2, -1, 2, 0]])
  • List of equations form

system  =  [3x + 2y - z - 1, 2x - 2y + 4z + 2, 2x - y + 2z]
  • Input \(A\) and \(b\) in matrix form (from \(Ax = b\)) are given as:

$$ A = left[begin{array}{ccc}

3 & 2 & -1 \ 2 & -2 & 4 \ 2 & -1 & 2 end{array}right] b = left[begin{array}{c} 1 \ -2 \ 0 end{array}right] $$

A = Matrix([[3, 2, -1], [2, -2, 4], [2, -1, 2]])
b = Matrix([[1], [-2], [0]])
system = (A, b)

Symbols can always be passed but are actually only needed when 1) a system of equations is being passed and 2) the system is passed as an underdetermined matrix and one wants to control the name of the free variables in the result. An error is raised if no symbols are used for case 1, but if no symbols are provided for case 2, internally generated symbols will be provided. When providing symbols for case 2, there should be at least as many symbols are there are columns in matrix A.

The algorithm used here is Gauss-Jordan elimination, which results, after elimination, in a row echelon form matrix.


A FiniteSet containing an ordered tuple of values for the

unknowns for which the \(system\) has a solution. (Wrapping

the tuple in FiniteSet is used to maintain a consistent

output format throughout solveset.)

Returns EmptySet, if the linear system is inconsistent.



The input is not valid. The symbols are not given.


>>> from sympy import Matrix, linsolve, symbols
>>> x, y, z = symbols("x, y, z")
>>> A = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 10]])
>>> b = Matrix([3, 6, 9])
>>> A
[1, 2,  3],
[4, 5,  6],
[7, 8, 10]])
>>> b
>>> linsolve((A, b), [x, y, z])
{(-1, 2, 0)}
  • Parametric Solution: In case the system is underdetermined, the function will return a parametric solution in terms of the given symbols. Those that are free will be returned unchanged. e.g. in the system below, \(z\) is returned as the solution for variable z; it can take on any value.

>>> A = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
>>> b = Matrix([3, 6, 9])
>>> linsolve((A, b), x, y, z)
{(z - 1, 2 - 2*z, z)}

If no symbols are given, internally generated symbols will be used. The tau0 in the third position indicates (as before) that the third variable – whatever it is named – can take on any value:

>>> linsolve((A, b))
{(tau0 - 1, 2 - 2*tau0, tau0)}
  • List of equations as input

>>> Eqns = [3*x + 2*y - z - 1, 2*x - 2*y + 4*z + 2, - x + y/2 - z]
>>> linsolve(Eqns, x, y, z)
{(1, -2, -2)}
  • Augmented matrix as input

>>> aug = Matrix([[2, 1, 3, 1], [2, 6, 8, 3], [6, 8, 18, 5]])
>>> aug
[2, 1,  3, 1],
[2, 6,  8, 3],
[6, 8, 18, 5]])
>>> linsolve(aug, x, y, z)
{(3/10, 2/5, 0)}
  • Solve for symbolic coefficients

>>> a, b, c, d, e, f = symbols('a, b, c, d, e, f')
>>> eqns = [a*x + b*y - c, d*x + e*y - f]
>>> linsolve(eqns, x, y)
{((-b*f + c*e)/(a*e - b*d), (a*f - c*d)/(a*e - b*d))}
  • A degenerate system returns solution as set of given symbols.

>>> system = Matrix(([0, 0, 0], [0, 0, 0], [0, 0, 0]))
>>> linsolve(system, x, y)
{(x, y)}
  • For an empty system linsolve returns empty set

>>> linsolve([], x)
  • An error is raised if any nonlinearity is detected, even if it could be removed with expansion

>>> linsolve([x*(1/x - 1)], x)
Traceback (most recent call last):
NonlinearError: nonlinear term: 1/x
>>> linsolve([x*(y + 1)], x, y)
Traceback (most recent call last):
NonlinearError: nonlinear cross-term: x*(y + 1)
>>> linsolve([x**2 - 1], x)
Traceback (most recent call last):
NonlinearError: nonlinear term: x**2
sympy.solvers.solveset.nonlinsolve(system, *symbols)[source]

Solve system of \(N\) nonlinear equations with \(M\) variables, which means both under and overdetermined systems are supported. Positive dimensional system is also supported (A system with infinitely many solutions is said to be positive-dimensional). In a positive dimensional system the solution will be dependent on at least one symbol. Returns both real solution and complex solution (if they exist).


system : list of equations

The target system of equations

symbols : list of Symbols

symbols should be given as a sequence eg. list


A FiniteSet of ordered tuple of values of \(symbols\) for which the \(system\)

has solution. Order of values in the tuple is same as symbols present in

the parameter \(symbols\).

Please note that general FiniteSet is unordered, the solution

returned here is not simply a FiniteSet of solutions, rather it

is a FiniteSet of ordered tuple, i.e. the first and only

argument to FiniteSet is a tuple of solutions, which is

ordered, and, hence ,the returned solution is ordered.

Also note that solution could also have been returned as an ordered tuple,

FiniteSet is just a wrapper {} around the tuple. It has no other

significance except for the fact it is just used to maintain a consistent

output format throughout the solveset.

For the given set of equations, the respective input types

are given below:

\[xy - 1 = 0\]
\[4x^2 + y^2 - 5 = 0\]
system  = [x*y - 1, 4*x**2 + y**2 - 5]
symbols = [x, y]


The input is not valid. The symbols are not given.


The input symbols are not \(Symbol\) type.


>>> from sympy import symbols, 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)}
  1. Positive dimensional system and complements:

>>> from sympy import pprint
>>> from sympy.polys.polytools import is_zero_dimensional
>>> a, b, c, d = symbols('a, b, c, d', extended_real=True)
>>> eq1 =  a + b + c + d
>>> eq2 = a*b + b*c + c*d + d*a
>>> eq3 = a*b*c + b*c*d + c*d*a + d*a*b
>>> eq4 = a*b*c*d - 1
>>> system = [eq1, eq2, eq3, eq4]
>>> is_zero_dimensional(system)
>>> pprint(nonlinsolve(system, [a, b, c, d]), use_unicode=False)
  -1       1               1      -1
{(---, -d, -, {d} \ {0}), (-, -d, ---, {d} \ {0})}
   d       d               d       d
>>> nonlinsolve([(x+y)**2 - 4, x + y - 2], [x, y])
{(2 - y, y)}

2. If some of the equations are non-polynomial then \(nonlinsolve\) will call the substitution function and return real and complex solutions, if present.

>>> from sympy import exp, sin
>>> nonlinsolve([exp(x) - sin(y), y**2 - 4], [x, y])
{(ImageSet(Lambda(_n, I*(2*_n*pi + pi) + log(sin(2))), Integers), -2),
 (ImageSet(Lambda(_n, 2*_n*I*pi + log(sin(2))), Integers), 2)}

3. If system is non-linear polynomial and zero-dimensional then it returns both solution (real and complex solutions, if present) using solve_poly_system():

>>> from sympy import sqrt
>>> nonlinsolve([x**2 - 2*y**2 -2, x*y - 2], [x, y])
{(-2, -1), (2, 1), (-sqrt(2)*I, sqrt(2)*I), (sqrt(2)*I, -sqrt(2)*I)}

4. nonlinsolve can solve some linear (zero or positive dimensional) system (because it uses the sympy.polys.polytools.groebner() function to get the groebner basis and then uses the substitution function basis as the new \(system\)). But it is not recommended to solve linear system using nonlinsolve, because linsolve() is better for general linear systems.

>>> nonlinsolve([x + 2*y -z - 3, x - y - 4*z + 9, y + z - 4], [x, y, z])
{(3*z - 5, 4 - z, z)}

5. System having polynomial equations and only real solution is solved using solve_poly_system():

>>> e1 = sqrt(x**2 + y**2) - 10
>>> e2 = sqrt(y**2 + (-x + 10)**2) - 3
>>> nonlinsolve((e1, e2), (x, y))
{(191/20, -3*sqrt(391)/20), (191/20, 3*sqrt(391)/20)}
>>> nonlinsolve([x**2 + 2/y - 2, x + y - 3], [x, y])
{(1, 2), (1 - sqrt(5), 2 + sqrt(5)), (1 + sqrt(5), 2 - sqrt(5))}
>>> nonlinsolve([x**2 + 2/y - 2, x + y - 3], [y, x])
{(2, 1), (2 - sqrt(5), 1 + sqrt(5)), (2 + sqrt(5), 1 - sqrt(5))}

6. It is better to use symbols instead of trigonometric functions or Function. For example, replace \(\sin(x)\) with a symbol, replace \(f(x)\) with a symbol and so on. Get a solution from nonlinsolve and then use solveset() to get the value of \(x\).

How Nonlinsolve Is Better Than Old Solver _solve_system :

1. A positive dimensional system solver: nonlinsolve can return solution for positive dimensional system. It finds the Groebner Basis of the positive dimensional system(calling it as basis) then we can start solving equation(having least number of variable first in the basis) using solveset and substituting that solved solutions into other equation(of basis) to get solution in terms of minimum variables. Here the important thing is how we are substituting the known values and in which equations.

2. Real and complex solutions: nonlinsolve returns both real and complex solution. If all the equations in the system are polynomial then using solve_poly_system() both real and complex solution is returned. If all the equations in the system are not polynomial equation then goes to substitution method with this polynomial and non polynomial equation(s), to solve for unsolved variables. Here to solve for particular variable solveset_real and solveset_complex is used. For both real and complex solution _solve_using_known_values is used inside substitution (substitution will be called when any non-polynomial equation is present). If a solution is valid its general solution is added to the final result.

3. Complement and Intersection will be added: nonlinsolve maintains dict for complements and intersections. If solveset find complements or/and intersections with any interval or set during the execution of substitution function, then complement or/and intersection for that variable is added before returning final solution.


sympy.solvers.solveset._transolve(f, symbol, domain)[source]

Function to solve transcendental equations. It is a helper to solveset and should be used internally. _transolve currently supports the following class of equations:

  • Exponential equations

  • Logarithmic equations


f : Any transcendental equation that needs to be solved.

This needs to be an expression, which is assumed to be equal to 0.

symbol : The variable for which the equation is solved.

This needs to be of class Symbol.

domain : A set over which the equation is solved.

This needs to be of class Set.



A set of values for symbol for which f is equal to zero. An EmptySet is returned if f does not have solutions in respective domain. A ConditionSet is returned as unsolved object if algorithms to evaluate complete solution are not yet implemented.

How To Use _transolve

_transolve should not be used as an independent function, because it assumes that the equation (f) and the symbol comes from solveset and might have undergone a few modification(s). To use _transolve as an independent function the equation (f) and the symbol should be passed as they would have been by solveset.


>>> from sympy.solvers.solveset import _transolve as transolve
>>> from sympy.solvers.solvers import _tsolve as tsolve
>>> from sympy import symbols, S, pprint
>>> x = symbols('x', real=True) # assumption added
>>> transolve(5**(x - 3) - 3**(2*x + 1), x, S.Reals)
{-(log(3) + 3*log(5))/(-log(5) + 2*log(3))}

How _transolve Works

_transolve uses two types of helper functions to solve equations of a particular class:

Identifying helpers: To determine whether a given equation belongs to a certain class of equation or not. Returns either True or False.

Solving helpers: Once an equation is identified, a corresponding helper either solves the equation or returns a form of the equation that solveset might better be able to handle.

  • Philosophy behind the module

The purpose of _transolve is to take equations which are not already polynomial in their generator(s) and to either recast them as such through a valid transformation or to solve them outright. A pair of helper functions for each class of supported transcendental functions are employed for this purpose. One identifies the transcendental form of an equation and the other either solves it or recasts it into a tractable form that can be solved by solveset. For example, an equation in the form \(ab^{f(x)} - cd^{g(x)} = 0\) can be transformed to \(\log(a) + f(x)\log(b) - \log(c) - g(x)\log(d) = 0\) (under certain assumptions) and this can be solved with solveset if \(f(x)\) and \(g(x)\) are in polynomial form.

How _transolve Is Better Than _tsolve

  1. Better output

_transolve provides expressions in a more simplified form.

Consider a simple exponential equation

>>> f = 3**(2*x) - 2**(x + 3)
>>> pprint(transolve(f, x, S.Reals), use_unicode=False)
 -2*log(3) + log(2)
>>> pprint(tsolve(f, x), use_unicode=False)
     /   3     \
     | --------|
     | log(2/9)|
[-log\2         /]
  1. Extensible

The API of _transolve is designed such that it is easily extensible, i.e. the code that solves a given class of equations is encapsulated in a helper and not mixed in with the code of _transolve itself.

  1. Modular

_transolve is designed to be modular i.e, for every class of equation a separate helper for identification and solving is implemented. This makes it easy to change or modify any of the method implemented directly in the helpers without interfering with the actual structure of the API.

  1. Faster Computation

Solving equation via _transolve is much faster as compared to _tsolve. In solve, attempts are made computing every possibility to get the solutions. This series of attempts makes solving a bit slow. In _transolve, computation begins only after a particular type of equation is identified.

How To Add New Class Of Equations

Adding a new class of equation solver is a three-step procedure:

  • Identify the type of the equations

    Determine the type of the class of equations to which they belong: it could be of Add, Pow, etc. types. Separate internal functions are used for each type. Write identification and solving helpers and use them from within the routine for the given type of equation (after adding it, if necessary). Something like:

    def add_type(lhs, rhs, x):
        if _is_exponential(lhs, x):
            new_eq = _solve_exponential(lhs, rhs, x)
    rhs, lhs = eq.as_independent(x)
    if lhs.is_Add:
        result = add_type(lhs, rhs, x)
  • Define the identification helper.

  • Define the solving helper.

Apart from this, a few other things needs to be taken care while adding an equation solver:

  • Naming conventions: Name of the identification helper should be as _is_class where class will be the name or abbreviation of the class of equation. The solving helper will be named as _solve_class. For example: for exponential equations it becomes _is_exponential and _solve_expo.

  • The identifying helpers should take two input parameters, the equation to be checked and the variable for which a solution is being sought, while solving helpers would require an additional domain parameter.

  • Be sure to consider corner cases.

  • Add tests for each helper.

  • Add a docstring to your helper that describes the method implemented. The documentation of the helpers should identify:

    • the purpose of the helper,

    • the method used to identify and solve the equation,

    • a proof of correctness

    • the return values of the helpers

sympy.solvers.solveset._is_exponential(f, symbol)[source]

Return True if one or more terms contain symbol only in exponents, else False.


f : Expr

The equation to be checked

symbol : Symbol

The variable in which the equation is checked


>>> from sympy import symbols, cos, exp
>>> from sympy.solvers.solveset import _is_exponential as check
>>> x, y = symbols('x y')
>>> check(y, y)
>>> check(x**y - 1, y)
>>> check(x**y*2**y - 1, y)
>>> check(exp(x + 3) + 3**x, x)
>>> check(cos(2**x), x)
  • Philosophy behind the helper

The function extracts each term of the equation and checks if it is of exponential form w.r.t symbol.

sympy.solvers.solveset._solve_exponential(lhs, rhs, symbol, domain)[source]

Helper function for solving (supported) exponential equations.

Exponential equations are the sum of (currently) at most two terms with one or both of them having a power with a symbol-dependent exponent.

For example

\[5^{2x + 3} - 5^{3x - 1}\]
\[4^{5 - 9x} - e^{2 - x}\]

lhs, rhs : Expr

The exponential equation to be solved, \(lhs = rhs\)

symbol : Symbol

The variable in which the equation is solved

domain : Set

A set over which the equation is solved.


A set of solutions satisfying the given equation.

A ConditionSet if the equation is unsolvable or

if the assumptions are not properly defined, in that case

a different style of ConditionSet is returned having the

solution(s) of the equation with the desired assumptions.


>>> from sympy.solvers.solveset import _solve_exponential as solve_expo
>>> from sympy import symbols, S
>>> x = symbols('x', real=True)
>>> a, b = symbols('a b')
>>> solve_expo(2**x + 3**x - 5**x, 0, x, S.Reals)  # not solvable
ConditionSet(x, Eq(2**x + 3**x - 5**x, 0), Reals)
>>> solve_expo(a**x - b**x, 0, x, S.Reals)  # solvable but incorrect assumptions
ConditionSet(x, (a > 0) & (b > 0), {0})
>>> solve_expo(3**(2*x) - 2**(x + 3), 0, x, S.Reals)
{-3*log(2)/(-2*log(3) + log(2))}
>>> solve_expo(2**x - 4**x, 0, x, S.Reals)
  • Proof of correctness of the method

The logarithm function is the inverse of the exponential function. The defining relation between exponentiation and logarithm is:

\[{\log_b x} = y \enspace if \enspace b^y = x\]

Therefore if we are given an equation with exponent terms, we can convert every term to its corresponding logarithmic form. This is achieved by taking logarithms and expanding the equation using logarithmic identities so that it can easily be handled by solveset.

For example:

\[3^{2x} = 2^{x + 3}\]

Taking log both sides will reduce the equation to

\[(2x)\log(3) = (x + 3)\log(2)\]

This form can be easily handed by solveset.

sympy.solvers.solveset._solve_logarithm(lhs, rhs, symbol, domain)[source]

Helper to solve logarithmic equations which are reducible to a single instance of \(\log\).

Logarithmic equations are (currently) the equations that contains \(\log\) terms which can be reduced to a single \(\log\) term or a constant using various logarithmic identities.

For example:

\[\log(x) + \log(x - 4)\]

can be reduced to:

\[\log(x(x - 4))\]

lhs, rhs : Expr

The logarithmic equation to be solved, \(lhs = rhs\)

symbol : Symbol

The variable in which the equation is solved

domain : Set

A set over which the equation is solved.


A set of solutions satisfying the given equation.

A ConditionSet if the equation is unsolvable.


>>> from sympy import symbols, log, S
>>> from sympy.solvers.solveset import _solve_logarithm as solve_log
>>> x = symbols('x')
>>> f = log(x - 3) + log(x + 3)
>>> solve_log(f, 0, x, S.Reals)
{-sqrt(10), sqrt(10)}
  • Proof of correctness

A logarithm is another way to write exponent and is defined by

\[{\log_b x} = y \enspace if \enspace b^y = x\]

When one side of the equation contains a single logarithm, the equation can be solved by rewriting the equation as an equivalent exponential equation as defined above. But if one side contains more than one logarithm, we need to use the properties of logarithm to condense it into a single logarithm.

Take for example

\[\log(2x) - 15 = 0\]

contains single logarithm, therefore we can directly rewrite it to exponential form as

\[x = \frac{e^{15}}{2}\]

But if the equation has more than one logarithm as

\[\log(x - 3) + \log(x + 3) = 0\]

we use logarithmic identities to convert it into a reduced form


\[\log(a) + \log(b) = \log(ab)\]

the equation becomes,

\[\log((x - 3)(x + 3))\]

This equation contains one logarithm and can be solved by rewriting to exponents.

sympy.solvers.solveset._is_logarithmic(f, symbol)[source]

Return True if the equation is in the form \(a\log(f(x)) + b\log(g(x)) + ... + c\) else False.


f : Expr

The equation to be checked

symbol : Symbol

The variable in which the equation is checked


True if the equation is logarithmic otherwise False.


>>> from sympy import symbols, tan, log
>>> from sympy.solvers.solveset import _is_logarithmic as check
>>> x, y = symbols('x y')
>>> check(log(x + 2) - log(x + 3), x)
>>> check(tan(log(2*x)), x)
>>> check(x*log(x), x)
>>> check(x + log(x), x)
>>> check(y + log(x), x)
  • Philosophy behind the helper

The function extracts each term and checks whether it is logarithmic w.r.t symbol.

Diophantine Equations (DEs)

See Diophantine


See Inequality Solvers

Ordinary Differential equations (ODEs)

See ODE.

Partial Differential Equations (PDEs)

See PDE.