The solvers module in Diofant implements methods for solving equations.

Algebraic equations

Use solve() to solve algebraic equations. We suppose all equations are equaled to 0, so solving x**2 == 1 translates into the following code:

>>> from diofant.solvers import solve
>>> from diofant import Symbol
>>> x = Symbol('x')
>>> solve(x**2 - 1, x)
[-1, 1]

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

diofant.solvers.solvers.solve(f, *symbols, **flags)[source]

Algebraically solves equations and systems of equations.

Currently supported are:
  • polynomial,
  • transcendental
  • piecewise combinations of the above
  • systems of linear and polynomial equations
  • systems containing relational expressions.

Input is formed as:

  • f
    • a single Expr or Poly that must be zero,
    • an Equality
    • a Relational expression or boolean
    • iterable of one or more of the above
  • symbols (object(s) to solve for) specified as
    • none given (other non-numeric objects will be used)
    • single symbol
    • denested list of symbols e.g. solve(f, x, y)
    • ordered iterable of symbols e.g. solve(f, [x, y])
  • flags
    ‘dict’=True (default is False)

    return list (perhaps empty) of solution mappings

    ‘set’=True (default is False)

    return list of symbols and set of tuple(s) of solution(s)

    ‘exclude=[] (default)’

    don’t try to solve for any of the free symbols in exclude; if expressions are given, the free symbols in them will be extracted automatically.

    ‘check=True (default)’

    If False, don’t do any testing of solutions. This can be useful if one wants to include solutions that make any denominator zero.

    ‘numerical=True (default)’

    do a fast numerical check if f has only one symbol.

    ‘minimal=True (default is False)’

    a very fast, minimal testing.

    ‘warn=True (default is False)’

    show a warning if checksol() could not conclude.

    ‘simplify=True (default)’

    simplify all but polynomials of order 3 or greater before returning them and (if check is not False) use the general simplify function on the solutions and the expression obtained when they are substituted into the function which should be zero

    ‘force=True (default is False)’

    make positive all symbols without assumptions regarding sign.

    ‘rational=True (default)’

    recast Floats as Rational; if this option is not used, the system containing floats may fail to solve because of issues with polys. If rational=None, Floats will be recast as rationals but the answer will be recast as Floats. If the flag is False then nothing will be done to the Floats.

    ‘manual=True (default is False)’

    do not use the polys/matrix method to solve a system of equations, solve them one at a time as you might “manually”

    ‘implicit=True (default is False)’

    allows solve to return a solution for a pattern in terms of other functions that contain that pattern; this is only needed if the pattern is inside of some invertible function like cos, exp, ....

    ‘particular=True (default is False)’

    instructs solve to try to find a particular solution to a linear system with as many zeros as possible; this is very expensive

    ‘quick=True (default is False)’

    when using particular=True, use a fast heuristic instead to find a solution with many zeros (instead of using the very slow method guaranteed to find the largest number of zeros possible)

    ‘cubics=True (default)’

    return explicit solutions when cubic expressions are encountered

    ‘quartics=True (default)’

    return explicit solutions when quartic expressions are encountered

    ‘quintics=True (default)’

    return explicit solutions (if possible) when quintic expressions are encountered


solve() with check=True (default) will run through the symbol tags to elimate unwanted solutions. If no assumptions are included all possible solutions will be returned.

>>> from diofant import Symbol, solve
>>> x = Symbol("x")
>>> solve(x**2 - 1)
[-1, 1]

By using the positive tag only one solution will be returned:

>>> pos = Symbol("pos", positive=True)
>>> solve(pos**2 - 1)

Assumptions aren’t checked when \(solve()\) input involves relationals or bools.

When the solutions are checked, those that make any denominator zero are automatically excluded. If you do not want to exclude such solutions then use the check=False option:

>>> from diofant import sin, limit
>>> solve(sin(x)/x)  # 0 is excluded

If check=False then a solution to the numerator being zero is found: x = 0. In this case, this is a spurious solution since sin(x)/x has the well known limit (without dicontinuity) of 1 at x = 0:

>>> solve(sin(x)/x, check=False)
[0, pi]

In the following case, however, the limit exists and is equal to the the value of x = 0 that is excluded when check=True:

>>> eq = x**2*(1/x - z**2/x)
>>> solve(eq, x)
>>> solve(eq, x, check=False)
>>> limit(eq, x, 0, '-')
>>> limit(eq, x, 0, '+')

Disabling high-order, explicit solutions

When solving polynomial expressions, one might not want explicit solutions (which can be quite long). If the expression is univariate, RootOf instances will be returned instead:

>>> solve(x**3 - x + 1)
[-1/((-1/2 - sqrt(3)*I/2)*(3*sqrt(69)/2 + 27/2)**(1/3)) - (-1/2 -
sqrt(3)*I/2)*(3*sqrt(69)/2 + 27/2)**(1/3)/3, -(-1/2 +
sqrt(3)*I/2)*(3*sqrt(69)/2 + 27/2)**(1/3)/3 - 1/((-1/2 +
sqrt(3)*I/2)*(3*sqrt(69)/2 + 27/2)**(1/3)), -(3*sqrt(69)/2 +
27/2)**(1/3)/3 - 1/(3*sqrt(69)/2 + 27/2)**(1/3)]
>>> solve(x**3 - x + 1, cubics=False)
[RootOf(x**3 - x + 1, 0), RootOf(x**3 - x + 1, 1), RootOf(x**3 - x + 1, 2)]

If the expression is multivariate, no solution might be returned:

>>> solve(x**3 - x + a, x, cubics=False)

Sometimes solutions will be obtained even when a flag is False because the expression could be factored. In the following example, the equation can be factored as the product of a linear and a quadratic factor so explicit solutions (which did not require solving a cubic expression) are obtained:

>>> eq = x**3 + 3*x**2 + x - 1
>>> solve(eq, cubics=False)
[-1, -1 + sqrt(2), -sqrt(2) - 1]

Solving equations involving radicals

Because of Diofant’s use of the principle root (issue sympy/sympy#8789), some solutions to radical equations will be missed unless check=False:

>>> from diofant import root
>>> eq = root(x**3 - 3*x**2, 3) + 1 - x
>>> solve(eq)
>>> solve(eq, check=False)

In the above example there is only a single solution to the equation. Other expressions will yield spurious roots which must be checked manually; roots which give a negative argument to odd-powered radicals will also need special checking:

>>> from diofant import real_root, Rational
>>> eq = root(x, 3) - root(x, 5) + Rational(1, 7)
>>> solve(eq)  # this gives 2 solutions but misses a 3rd
[RootOf(7*_p**5 - 7*_p**3 + 1, 1)**15,
RootOf(7*_p**5 - 7*_p**3 + 1, 2)**15]
>>> sol = solve(eq, check=False)
>>> [abs(eq.subs(x,i).n(2)) for i in sol]
[0.48, 0.e-110, 0.e-110, 0.052, 0.052]

The first solution is negative so real_root must be used to see that it satisfies the expression:

>>> abs(real_root(eq.subs(x, sol[0])).n(2))

If the roots of the equation are not real then more care will be necessary to find the roots, especially for higher order equations. Consider the following expression:

>>> expr = root(x, 3) - root(x, 5)

We will construct a known value for this expression at x = 3 by selecting the 1-th root for each radical:

>>> expr1 = root(x, 3, 1) - root(x, 5, 1)
>>> v = expr1.subs(x, -3)

The solve function is unable to find any exact roots to this equation:

>>> eq = Eq(expr, v); eq1 = Eq(expr1, v)
>>> solve(eq, check=False), solve(eq1, check=False)
([], [])

The function unrad, however, can be used to get a form of the equation for which numerical roots can be found:

>>> from diofant.solvers.solvers import unrad
>>> from diofant import nroots, pprint
>>> e, (p, cov) = unrad(eq)
>>> pvals = nroots(e)
>>> inversion = solve(cov, x)[0]
>>> xvals = [inversion.subs(p, i) for i in pvals]

Although eq or eq1 could have been used to find xvals, the solution can only be verified with expr1:

>>> z = expr - v
>>> [xi.n(chop=1e-9) for xi in xvals if abs(z.subs(x, xi).n()) < 1e-9]
>>> z1 = expr1 - v
>>> pprint([xi.n(chop=1e-9) for xi in xvals if abs(z1.subs(x, xi).n()) < 1e-9])


The output varies according to the input and can be seen by example:

>>> from diofant import solve, Poly, Eq, Function, exp, I, sqrt
>>> from import x, y, z, a, b
>>> f = Function('f')
  • boolean or univariate Relational

    >>> solve(x < 3)
    And(-oo < x, x < 3)
  • to always get a list of solution mappings, use flag dict=True

    >>> solve(x - 3, dict=True)
    [{x: 3}]
    >>> solve([x - 3, y - 1], dict=True) == [{x: 3, y: 1}]
  • to get a list of symbols and set of solution(s) use flag set=True

    >>> solve([x**2 - 3, y - 1], set=True) == ([x, y], {(-sqrt(3), 1), (sqrt(3), 1)})
  • single expression and single symbol that is in the expression

    >>> solve(x - y, x)
    >>> solve(x - 3, x)
    >>> solve(Eq(x, 3), x)
    >>> solve(Poly(x - 3), x)
    >>> solve(x**2 - y**2, x, set=True) == ([x], {(-y,), (y,)})
    >>> solve(x**4 - 1, x, set=True) == ([x], {(-1,), (1,), (-I,), (I,)})
  • single expression with no symbol that is in the expression

    >>> solve(3, x)
    >>> solve(x - 3, y)
  • single expression with no symbol given

    In this case, all free symbols will be selected as potential symbols to solve for. If the equation is univariate then a list of solutions is returned; otherwise – as is the case when symbols are given as an iterable of length > 1 – a list of mappings will be returned.

    >>> solve(x - 3)
    >>> solve(x**2 - y**2) == [{x: -y}, {x: y}]
    >>> solve(z**2*x**2 - z**2*y**2) == [{x: -y}, {x: y}, {z: 0}]
    >>> solve(z**2*x - z**2*y**2) == [{x: y**2}, {z: 0}]
  • when an object other than a Symbol is given as a symbol, it is isolated algebraically and an implicit solution may be obtained. This is mostly provided as a convenience to save one from replacing the object with a Symbol and solving for that Symbol. It will only work if the specified object can be replaced with a Symbol using the subs method.

    >>> solve(f(x) - x, f(x))
    >>> solve(f(x).diff(x) - f(x) - x, f(x).diff(x))
    [x + f(x)]
    >>> solve(f(x).diff(x) - f(x) - x, f(x))
    [-x + Derivative(f(x), x)]
    >>> solve(x + exp(x)**2, exp(x), set=True) == ([exp(x)], {(-sqrt(-x),), (sqrt(-x),)})
    >>> from diofant import Indexed, IndexedBase, Tuple, sqrt
    >>> A = IndexedBase('A')
    >>> eqs = Tuple(A[1] + A[2] - 3, A[1] - A[2] + 1)
    >>> solve(eqs, eqs.atoms(Indexed)) == {A[1]: 1, A[2]: 2}
    • To solve for a symbol implicitly, use ‘implicit=True’:

      >>> solve(x + exp(x), x)
      >>> solve(x + exp(x), x, implicit=True)
    • It is possible to solve for anything that can be targeted with subs:

      >>> solve(x + 2 + sqrt(3), x + 2)
      >>> solve((x + 2 + sqrt(3), x + 4 + y), y, x + 2) == {y: -2 + sqrt(3), x + 2: -sqrt(3)}
    • Nothing heroic is done in this implicit solving so you may end up with a symbol still in the solution:

      >>> eqs = (x*y + 3*y + sqrt(3), x + 4 + y)
      >>> solve(eqs, y, x + 2) == {y: -sqrt(3)/(x + 3), x + 2: (-2*x - 6 + sqrt(3))/(x + 3)}
      >>> solve(eqs, y*x, x) == {x: -y - 4, x*y: -3*y - sqrt(3)}
    • if you attempt to solve for a number remember that the number you have obtained does not necessarily mean that the value is equivalent to the expression obtained:

      >>> solve(sqrt(2) - 1, 1)
      >>> solve(x - y + 1, 1)  # /!\ -1 is targeted, too
      [x/(y - 1)]
      >>> [_.subs(z, -1) for _ in solve((x - y + 1).subs(-1, z), 1)]
      [-x + y]
    • To solve for a function within a derivative, use dsolve.

  • single expression and more than 1 symbol

    • when there is a linear solution

      >>> solve(x - y**2, x, y)
      [{x: y**2}]
      >>> solve(x**2 - y, x, y)
      [{y: x**2}]
    • when undetermined coefficients are identified

      • that are linear

        >>> solve((a + b)*x - b + 2, a, b) == {a: -2, b: 2}
      • that are nonlinear

        >>> solve((a + b)*x - b**2 + 2, a, b, set=True) == ([a, b], {(-sqrt(2), sqrt(2)), (sqrt(2), -sqrt(2))})
    • if there is no linear solution then the first successful attempt for a nonlinear solution will be returned

      >>> solve(x**2 - y**2, x, y) == [{x: -y}, {x: y}]
      >>> solve(x**2 - y**2/exp(x), x, y)
      [{x: 2*LambertW(y/2)}]
      >>> solve(x**2 - y**2/exp(x), y, x) == [{y: -x*sqrt(exp(x))}, {y: x*sqrt(exp(x))}]
  • iterable of one or more of the above

    • involving relationals or bools

      >>> solve([x < 3, x - 2])
      Eq(x, 2)
      >>> solve([x > 3, x - 2])
    • when the system is linear

      • with a solution

        >>> solve([x - 3], x)
        {x: 3}
        >>> solve((x + 5*y - 2, -3*x + 6*y - 15), x, y) == {x: -3, y: 1}
        >>> solve((x + 5*y - 2, -3*x + 6*y - 15), x, y, z) == {x: -3, y: 1}
        >>> solve((x + 5*y - 2, -3*x + 6*y - z), z, x, y) == {x: -5*y + 2, z: 21*y - 6}
      • without a solution

        >>> solve([x + 3, x - 3])
    • when the system is not linear

      >>> solve([x**2 + y -2, y**2 - 4], x, y, set=True) == ([x, y], {(-2, -2), (0, 2), (2, -2)})
    • if no symbols are given, all free symbols will be selected and a list of mappings returned

      >>> solve([x - 2, x**2 + y]) == [{x: 2, y: -4}]
      >>> solve([x - 2, x**2 + f(x)], {f(x), x}) == [{x: 2, f(x): -4}]
    • if any equation doesn’t depend on the symbol(s) given it will be eliminated from the equation set and an answer may be given implicitly in terms of variables that were not of interest

      >>> solve([x - y, y - 3], x)
      {x: y}
diofant.solvers.solvers.solve_linear(lhs, rhs=0, symbols=[], exclude=[])[source]

Return a tuple derived from f = lhs - rhs that is either:

(numerator, denominator) of f

If this comes back as (0, 1) it means that f is independent of the symbols in symbols, e.g:

y*cos(x)**2 + y*sin(x)**2 - y = y*(0) = 0
cos(x)**2 + sin(x)**2 = 1

If it comes back as (0, 0) there is no solution to the equation amongst the symbols given.

If the numerator is not zero then the function is guaranteed to be dependent on a symbol in symbols.


(symbol, solution) where symbol appears linearly in the numerator of f, is in symbols (if given) and is not in exclude (if given).

No simplification is done to f other than and mul=True expansion, so the solution will correspond strictly to a unique solution.


>>> from diofant.solvers.solvers import solve_linear
>>> from import x, y, z

These are linear in x and 1/x:

>>> solve_linear(x + y**2)
(x, -y**2)
>>> solve_linear(1/x - y**2)
(x, y**(-2))

When not linear in x or y then the numerator and denominator are returned.

>>> solve_linear(x**2/y**2 - 3)
(x**2 - 3*y**2, y**2)

If the numerator is a symbol then (0, 0) is returned if the solution for that symbol would have set any denominator to 0:

>>> solve_linear(1/(1/x - 2))
(0, 0)
>>> 1/(1/x)  # to Diofant, this looks like x ...
>>> solve_linear(1/(1/x))  # so a solution is given
(x, 0)

If x is allowed to cancel, then this appears linear, but this sort of cancellation is not done so the solution will always satisfy the original expression without causing a division by zero error.

>>> solve_linear(x**2*(1/x - z**2/x))
(x**2*(-z**2 + 1), x)

You can give a list of what you prefer for x candidates:

>>> solve_linear(x + y + z, symbols=[y])
(y, -x - z)

You can also indicate what variables you don’t want to consider:

>>> solve_linear(x + y + z, exclude=[x, z])
(y, -x - z)

If only x was excluded then a solution for y or z might be obtained.

diofant.solvers.solvers.solve_linear_system(system, *symbols, **flags)[source]

Solve system of linear equations.

Both under- and overdetermined systems are supported. The possible number of solutions is zero, one or infinite.


system : Matrix

Nx(M+1) matrix, which means it has to be in augmented form. This matrix will not be modified.

*symbols : list

List of M Symbol’s


solution: dict or None

Respectively, this procedure will return None or a dictionary with solutions. In the case of underdetermined systems, all arbitrary parameters are skipped. This may cause a situation in which an empty dictionary is returned. In that case, all symbols can be assigned arbitrary values.


>>> from diofant import Matrix, solve_linear_system
>>> from import x, y

Solve the following system:

   x + 4 y ==  2
-2 x +   y == 14
>>> system = Matrix(( (1, 4, 2), (-2, 1, 14)))
>>> solve_linear_system(system, x, y) == {x: -6, y: 2}

A degenerate system returns an empty dictionary.

>>> system = Matrix(( (0,0,0), (0,0,0) ))
>>> solve_linear_system(system, x, y)
diofant.solvers.solvers.solve_undetermined_coeffs(equ, coeffs, sym, **flags)[source]

Solve equation of a type p(x; a_1, ..., a_k) == q(x) where both p, q are univariate polynomials and f depends on k parameters. The result of this functions is a dictionary with symbolic values of those parameters with respect to coefficients in q.

This functions accepts both Equations class instances and ordinary Diofant expressions. Specification of parameters and variable is obligatory for efficiency and simplicity reason.

>>> from diofant import Eq, Rational
>>> from import a, b, c, x
>>> from diofant.solvers import solve_undetermined_coeffs
>>> solve_undetermined_coeffs(Eq(2*a*x + a+b, x), [a, b], x) == {a: Rational(1, 2), b: -Rational(1, 2)}
>>> solve_undetermined_coeffs(Eq(a*c*x + a+b, x), [a, b], x) == {a: 1/c, b: -1/c}
diofant.solvers.solvers.nsolve(*args, **kwargs)[source]

Solve a nonlinear equation system numerically:

nsolve(f, [args,] x0, modules=['mpmath'], **kwargs)

f is a vector function of symbolic expressions representing the system. args are the variables. If there is only one variable, this argument can be omitted. x0 is a starting vector close to a solution.

Use the modules keyword to specify which modules should be used to evaluate the function and the Jacobian matrix. Make sure to use a module that supports matrices. For more information on the syntax, please see the docstring of lambdify.

Overdetermined systems are supported.

>>> from diofant import Symbol, nsolve
>>> import diofant
>>> import mpmath
>>> = 15
>>> x1 = Symbol('x1')
>>> x2 = Symbol('x2')
>>> f1 = 3 * x1**2 - 2 * x2**2 - 1
>>> f2 = x1**2 - 2 * x1 + x2**2 + 2 * x2 - 8
>>> print(nsolve((f1, f2), (x1, x2), (-1, 1)))
[ 1.27844411169911]

For one-dimensional functions the syntax is simplified:

>>> from diofant import sin, nsolve
>>> from import x
>>> print(nsolve(sin(x), x, 2))
>>> print(nsolve(sin(x), 2))

mpmath.findroot is used, you can find there more extensive documentation, especially concerning keyword parameters and available solvers. Note, however, that this routine works only with the numerator of the function in the one-dimensional case, and for very steep functions near the root this may lead to a failure in the verification of the root. In this case you should use the flag \(verify=False\) and independently verify the solution.

>>> from diofant import cos, cosh
>>> from import i
>>> f = cos(x)*cosh(x) - 1
>>> nsolve(f, 3.14*100)
Traceback (most recent call last):
ValueError: Could not find root within given tolerance. (1.39267e+230 > 2.1684e-19)
>>> ans = nsolve(f, 3.14*100, verify=False); print(ans)
>>> f.subs(x, ans).n(2)
>>> (f/f.diff(x)).subs(x, ans).n(2)

One might safely skip the verification if bounds of the root are known and a bisection method is used:

>>> bounds = lambda i: (3.14*i, 3.14*(i + 1))
>>> print(nsolve(f, bounds(100), solver='bisect', verify=False))
diofant.solvers.solvers.check_assumptions(expr, **assumptions)[source]

Checks whether expression \(expr\) satisfies all assumptions.

\(assumptions\) is a dict of assumptions: {‘assumption’: True|False, ...}.


>>> from diofant import Symbol, pi, I, exp
>>> from diofant.solvers.solvers import check_assumptions
>>> check_assumptions(-5, integer=True)
>>> check_assumptions(pi, extended_real=True, integer=False)
>>> check_assumptions(pi, extended_real=True, negative=True)
>>> check_assumptions(exp(I*pi/7), extended_real=False)
>>> x = Symbol('x', extended_real=True, positive=True)
>>> check_assumptions(2*x + 1, extended_real=True, positive=True)
>>> check_assumptions(-2*x - 5, extended_real=True, positive=True)

\(None\) is returned if check_assumptions() could not conclude.

>>> check_assumptions(2*x - 1, extended_real=True, positive=True)
>>> z = Symbol('z')
>>> check_assumptions(z, extended_real=True)
diofant.solvers.solvers.checksol(f, symbol, sol=None, **flags)[source]

Checks whether sol is a solution of equation f == 0.

Input can be either a single symbol and corresponding value or a dictionary of symbols and values. When given as a dictionary and flag simplify=True, the values in the dictionary will be simplified. f can be a single equation or an iterable of equations. A solution must satisfy all equations in f to be considered valid; if a solution does not satisfy any equation, False is returned; if one or more checks are inconclusive (and none are False) then None is returned.


>>> from diofant import symbols
>>> from diofant.solvers import checksol
>>> x, y = symbols('x,y')
>>> checksol(x**4 - 1, x, 1)
>>> checksol(x**4 - 1, x, 0)
>>> checksol(x**2 + y**2 - 5**2, {x: 3, y: 4})

To check if an expression is zero using checksol, pass it as f and send an empty dictionary for symbol:

>>> checksol(x**2 + x - x*(x + 1), {})

None is returned if checksol() could not conclude.

‘numerical=True (default)’
do a fast numerical check if f has only one symbol.
‘minimal=True (default is False)’
a very fast, minimal testing.
‘warn=True (default is False)’
show a warning if checksol() could not conclude.
‘simplify=True (default)’
simplify solution before substituting into function and simplify the function before trying specific simplifications
‘force=True (default is False)’
make positive all symbols without assumptions regarding sign.
diofant.solvers.solvers.unrad(eq, *syms, **flags)[source]

Remove radicals with symbolic arguments and return (eq, cov), None or raise an error:

None is returned if there are no radicals to remove.

NotImplementedError is raised if there are radicals and they cannot be removed or if the relationship between the original symbols and the change of variable needed to rewrite the system as a polynomial cannot be solved.

Otherwise the tuple, (eq, cov), is returned where:

``eq``, ``cov``
    ``eq`` is an equation without radicals (in the symbol(s) of
    interest) whose solutions are a superset of the solutions to the
    original expression. ``eq`` might be re-written in terms of a new
    variable; the relationship to the original variables is given by
    ``cov`` which is a list containing ``v`` and ``v**p - b`` where
    ``p`` is the power needed to clear the radical and ``b`` is the
    radical now expressed as a polynomial in the symbols of interest.
    For example, for sqrt(2 - x) the tuple would be
    ``(c, c**2 - 2 + x)``. The solutions of ``eq`` will contain
    solutions to the original equation (if there are any).
an iterable of symbols which, if provided, will limit the focus of radical removal: only radicals with one or more of the symbols of interest will be cleared. All free symbols are used if syms is not set.

flags are used internally for communication during recursive calls. Two options are also recognized:

``take``, when defined, is interpreted as a single-argument function
that returns True if a given Pow should be handled.

Radicals can be removed from an expression if:

*   all bases of the radicals are the same; a change of variables is
    done in this case.
*   if all radicals appear in one term of the expression
*   there are only 4 terms with sqrt() factors or there are less than
    four terms having sqrt() factors
*   there are only two terms with radicals


>>> from diofant.solvers.solvers import unrad
>>> from import x
>>> from diofant import sqrt, Rational, root, real_roots, solve
>>> unrad(sqrt(x)*x**Rational(1, 3) + 2)
(x**5 - 64, [])
>>> unrad(sqrt(x) + root(x + 1, 3))
(x**3 - x**2 - 2*x - 1, [])
>>> eq = sqrt(x) + root(x, 3) - 2
>>> unrad(eq)
(_p**3 + _p**2 - 2, [_p, -x + _p**6])

Ordinary Differential equations (ODEs)

See ODE.

Partial Differential Equations (PDEs)

See PDE.

Deutils (Utilities for solving ODE’s and PDE’s)

diofant.solvers.deutils.ode_order(expr, func)[source]

Returns the order of a given differential equation with respect to func.

This function is implemented recursively.


>>> from diofant import Function
>>> from diofant.solvers.deutils import ode_order
>>> from import x
>>> f, g = map(Function, ['f', 'g'])
>>> ode_order(f(x).diff(x, 2) + f(x).diff(x)**2 +
... f(x).diff(x), f(x))
>>> ode_order(f(x).diff(x, 2) + g(x).diff(x, 3), f(x))
>>> ode_order(f(x).diff(x, 2) + g(x).diff(x, 3), g(x))

Recurrence Equations

diofant.solvers.recurr.rsolve(f, y, init=None)[source]

Solve univariate recurrence with rational coefficients.

Given \(k\)-th order linear recurrence \(\operatorname{L} y = f\), or equivalently:

\[a_{k}(n) y(n+k) + a_{k-1}(n) y(n+k-1) + \ldots + a_{0}(n) y(n) = f(n)\]

where \(a_{i}(n)\), for \(i=0, \ldots, k\), are polynomials or rational functions in \(n\), and \(f\) is a hypergeometric function or a sum of a fixed number of pairwise dissimilar hypergeometric terms in \(n\), finds all solutions or returns None, if none were found.

Initial conditions can be given as a dictionary in two forms:

  1. {   n_0  : v_0,   n_1  : v_1, ...,   n_m  : v_m }
  2. { y(n_0) : v_0, y(n_1) : v_1, ..., y(n_m) : v_m }

or as a list L of values:

L = [ v_0, v_1, ..., v_m ]

where L[i] = v_i, for \(i=0, \ldots, m\), maps to \(y(n_i)\).


Lets consider the following recurrence:

\[(n - 1) y(n + 2) - (n^2 + 3 n - 2) y(n + 1) + 2 n (n + 1) y(n) = 0\]
>>> from diofant import Function, rsolve
>>> from import n
>>> y = Function('y')
>>> f = (n - 1)*y(n + 2) - (n**2 + 3*n - 2)*y(n + 1) + 2*n*(n + 1)*y(n)
>>> rsolve(f, y(n))
2**n*C0 + C1*factorial(n)
>>> rsolve(f, y(n), { y(0):0, y(1):3 })
3*2**n - 3*factorial(n)
diofant.solvers.recurr.rsolve_poly(coeffs, f, n, **hints)[source]

Given linear recurrence operator \(\operatorname{L}\) of order \(k\) with polynomial coefficients and inhomogeneous equation \(\operatorname{L} y = f\), where \(f\) is a polynomial, we seek for all polynomial solutions over field \(K\) of characteristic zero.

The algorithm performs two basic steps:

  1. Compute degree \(N\) of the general polynomial solution.
  2. Find all polynomials of degree \(N\) or less of \(\operatorname{L} y = f\).

There are two methods for computing the polynomial solutions. If the degree bound is relatively small, i.e. it’s smaller than or equal to the order of the recurrence, then naive method of undetermined coefficients is being used. This gives system of algebraic equations with \(N+1\) unknowns.

In the other case, the algorithm performs transformation of the initial equation to an equivalent one, for which the system of algebraic equations has only \(r\) indeterminates. This method is quite sophisticated (in comparison with the naive one) and was invented together by Abramov, Bronstein and Petkovšek.

It is possible to generalize the algorithm implemented here to the case of linear q-difference and differential equations.

Lets say that we would like to compute \(m\)-th Bernoulli polynomial up to a constant. For this we can use \(b(n+1) - b(n) = m n^{m-1}\) recurrence, which has solution \(b(n) = B_m + C\). For example:

>>> from diofant import Symbol, rsolve_poly
>>> n = Symbol('n', integer=True)
>>> rsolve_poly([-1, 1], 4*n**3, n)
C0 + n**4 - 2*n**3 + n**2


[R518]S. A. Abramov, M. Bronstein and M. Petkovšek, On polynomial solutions of linear operator equations, in: T. Levelt, ed., Proc. ISSAC ‘95, ACM Press, New York, 1995, 290-296.
[R519]M. Petkovšek, Hypergeometric solutions of linear recurrences with polynomial coefficients, J. Symbolic Computation, 14 (1992), 243-264.
  1. Petkovšek, H. S. Wilf, D. Zeilberger, A = B, 1996.
diofant.solvers.recurr.rsolve_ratio(coeffs, f, n, **hints)[source]

Given linear recurrence operator \(\operatorname{L}\) of order \(k\) with polynomial coefficients and inhomogeneous equation \(\operatorname{L} y = f\), where \(f\) is a polynomial, we seek for all rational solutions over field \(K\) of characteristic zero.

This procedure accepts only polynomials, however if you are interested in solving recurrence with rational coefficients then use rsolve which will pre-process the given equation and run this procedure with polynomial arguments.

The algorithm performs two basic steps:

  1. Compute polynomial \(v(n)\) which can be used as universal denominator of any rational solution of equation \(\operatorname{L} y = f\).
  2. Construct new linear difference equation by substitution \(y(n) = u(n)/v(n)\) and solve it for \(u(n)\) finding all its polynomial solutions. Return None if none were found.

Algorithm implemented here is a revised version of the original Abramov’s algorithm, developed in 1989. The new approach is much simpler to implement and has better overall efficiency. This method can be easily adapted to q-difference equations case.

Besides finding rational solutions alone, this functions is an important part of Hyper algorithm were it is used to find particular solution of inhomogeneous part of a recurrence.

See also



[R521]S. A. Abramov, Rational solutions of linear difference and q-difference equations with polynomial coefficients, in: T. Levelt, ed., Proc. ISSAC ‘95, ACM Press, New York, 1995, 285-289


>>> from import x
>>> from diofant.solvers.recurr import rsolve_ratio
>>> rsolve_ratio([-2*x**3 + x**2 + 2*x - 1, 2*x**3 + x**2 - 6*x,
... - 2*x**3 - 11*x**2 - 18*x - 9, 2*x**3 + 13*x**2 + 22*x + 8], 0, x)
C2*(2*x - 3)/(2*(x**2 - 1))
diofant.solvers.recurr.rsolve_hyper(coeffs, f, n, **hints)[source]

Given linear recurrence operator \(\operatorname{L}\) of order \(k\) with polynomial coefficients and inhomogeneous equation \(\operatorname{L} y = f\) we seek for all hypergeometric solutions over field \(K\) of characteristic zero.

The inhomogeneous part can be either hypergeometric or a sum of a fixed number of pairwise dissimilar hypergeometric terms.

The algorithm performs three basic steps:

  1. Group together similar hypergeometric terms in the inhomogeneous part of \(\operatorname{L} y = f\), and find particular solution using Abramov’s algorithm.
  2. Compute generating set of \(\operatorname{L}\) and find basis in it, so that all solutions are linearly independent.
  3. Form final solution with the number of arbitrary constants equal to dimension of basis of \(\operatorname{L}\).

Term \(a(n)\) is hypergeometric if it is annihilated by first order linear difference equations with polynomial coefficients or, in simpler words, if consecutive term ratio is a rational function.

The output of this procedure is a linear combination of fixed number of hypergeometric terms. However the underlying method can generate larger class of solutions - D’Alembertian terms.

Note also that this method not only computes the kernel of the inhomogeneous equation, but also reduces in to a basis so that solutions generated by this procedure are linearly independent


[R522]M. Petkovšek, Hypergeometric solutions of linear recurrences with polynomial coefficients, J. Symbolic Computation, 14 (1992), 243-264.
  1. Petkovšek, H. S. Wilf, D. Zeilberger, A = B, 1996.


>>> from diofant.solvers import rsolve_hyper
>>> from import x
>>> rsolve_hyper([-1, -1, 1], 0, x)
C0*(1/2 + sqrt(5)/2)**x + C1*(-sqrt(5)/2 + 1/2)**x
>>> rsolve_hyper([-1, 1], 1 + x, x)
C0 + x*(x + 1)/2

Systems of Polynomial Equations

diofant.solvers.polysys.solve_poly_system(seq, *gens, **args)[source]

Solve a system of polynomial equations.


>>> from diofant import solve_poly_system
>>> from import x, y
>>> solve_poly_system([x*y - 2*y, 2*y**2 - x**2], x, y)
[(0, 0), (2, -sqrt(2)), (2, sqrt(2))]
diofant.solvers.polysys.solve_triangulated(polys, *gens, **args)[source]

Solve a polynomial system using Gianni-Kalkbrenner algorithm.

The algorithm proceeds by computing one Groebner basis in the ground domain and then by iteratively computing polynomial factorizations in appropriately constructed algebraic extensions of the ground domain.


1. Patrizia Gianni, Teo Mora, Algebraic Solution of System of Polynomial Equations using Groebner Bases, AAECC-5 on Applied Algebra, Algebraic Algorithms and Error-Correcting Codes, LNCS 356 247–257, 1989


>>> from diofant.solvers.polysys import solve_triangulated
>>> from import x, y, z
>>> F = [x**2 + y + z - 1, x + y**2 + z - 1, x + y + z**2 - 1]
>>> solve_triangulated(F, x, y, z)
[(0, 0, 1), (0, 1, 0), (1, 0, 0)]

Diophantine Equations (DEs)

See Diophantine


See Inequality Solvers