Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
"""Solvers of systems of polynomial equations. """
from __future__ import print_function, division
from sympy.core import S from sympy.polys import Poly, groebner, roots from sympy.polys.polytools import parallel_poly_from_expr from sympy.polys.polyerrors import (ComputationFailed, PolificationFailed, CoercionFailed) from sympy.simplify import rcollect from sympy.utilities import default_sort_key, postfixes
class SolveFailed(Exception): """Raised when solver's conditions weren't met. """
def solve_poly_system(seq, *gens, **args): """ Solve a system of polynomial equations.
Examples ========
>>> from sympy import solve_poly_system >>> from sympy.abc 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))]
""" except PolificationFailed as exc: raise ComputationFailed('solve_poly_system', len(seq), exc)
except SolveFailed: pass
def solve_biquadratic(f, g, opt): """Solve a system of two bivariate quadratic polynomial equations.
Examples ========
>>> from sympy.polys import Options, Poly >>> from sympy.abc import x, y >>> from sympy.solvers.polysys import solve_biquadratic >>> NewOption = Options((x, y), {'domain': 'ZZ'})
>>> a = Poly(y**2 - 4 + x, y, x, domain='ZZ') >>> b = Poly(y*2 + 3*x - 7, y, x, domain='ZZ') >>> solve_biquadratic(a, b, NewOption) [(1/3, 3), (41/27, 11/9)]
>>> a = Poly(y + x**2 - 3, y, x, domain='ZZ') >>> b = Poly(-y + x - 4, y, x, domain='ZZ') >>> solve_biquadratic(a, b, NewOption) [(-sqrt(29)/2 + 7/2, -sqrt(29)/2 - 1/2), (sqrt(29)/2 + 7/2, -1/2 + \ sqrt(29)/2)] """
return None
raise SolveFailed
def solve_generic(polys, opt): """ Solve a generic system of polynomial equations.
Returns all possible solutions over C[x_1, x_2, ..., x_m] of a set F = { f_1, f_2, ..., f_n } of polynomial equations, using Groebner basis approach. For now only zero-dimensional systems are supported, which means F can have at most a finite number of solutions.
The algorithm works by the fact that, supposing G is the basis of F with respect to an elimination order (here lexicographic order is used), G and F generate the same ideal, they have the same set of solutions. By the elimination property, if G is a reduced, zero-dimensional Groebner basis, then there exists an univariate polynomial in G (in its last variable). This can be solved by computing its roots. Substituting all computed roots for the last (eliminated) variable in other elements of G, new polynomial system is generated. Applying the above procedure recursively, a finite number of solutions can be found.
The ability of finding all solutions by this procedure depends on the root finding algorithms. If no solutions were found, it means only that roots() failed, but the system is solvable. To overcome this difficulty use numerical algorithms instead.
References ==========
.. [Buchberger01] B. Buchberger, Groebner Bases: A Short Introduction for Systems Theorists, In: R. Moreno-Diaz, B. Buchberger, J.L. Freire, Proceedings of EUROCAST'01, February, 2001
.. [Cox97] D. Cox, J. Little, D. O'Shea, Ideals, Varieties and Algorithms, Springer, Second Edition, 1997, pp. 112
Examples ========
>>> from sympy.polys import Poly, Options >>> from sympy.solvers.polysys import solve_generic >>> from sympy.abc import x, y >>> NewOption = Options((x, y), {'domain': 'ZZ'})
>>> a = Poly(x - y + 5, x, y, domain='ZZ') >>> b = Poly(x + y - 3, x, y, domain='ZZ') >>> solve_generic([a, b], NewOption) [(-1, 4)]
>>> a = Poly(x - 2*y + 5, x, y, domain='ZZ') >>> b = Poly(2*x - y - 3, x, y, domain='ZZ') >>> solve_generic([a, b], NewOption) [(11/3, 13/3)]
>>> a = Poly(x**2 + y, x, y, domain='ZZ') >>> b = Poly(x + y*4, x, y, domain='ZZ') >>> solve_generic([a, b], NewOption) [(0, 0), (1/4, -1/16)] """ """Returns True if 'f' is univariate in its last variable. """
"""Replace generator with a root so that the result is nice. """
p = p.expand(deep=False)
"""Recursively solves reduced polynomial systems. """
if not entry: return [] else: return None
else: raise NotImplementedError("only zero-dimensional systems supported (finite number of solutions)")
return []
return [ (zero,) for zero in zeros ]
except CoercionFailed: raise NotImplementedError
else: return None
def solve_triangulated(polys, *gens, **args): """ 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.
Examples ========
>>> from sympy.solvers.polysys import solve_triangulated >>> from sympy.abc 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)]
References ==========
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
""" G = groebner(polys, gens, polys=True) G = list(reversed(G))
domain = args.get('domain')
if domain is not None: for i, g in enumerate(G): G[i] = g.set_domain(domain)
f, G = G[0].ltrim(-1), G[1:] dom = f.get_domain()
zeros = f.ground_roots() solutions = set([])
for zero in zeros: solutions.add(((zero,), dom))
var_seq = reversed(gens[:-1]) vars_seq = postfixes(gens[1:])
for var, vars in zip(var_seq, vars_seq): _solutions = set([])
for values, dom in solutions: H, mapping = [], list(zip(vars, values))
for g in G: _vars = (var,) + vars
if g.has_only_gens(*_vars) and g.degree(var) != 0: h = g.ltrim(var).eval(dict(mapping))
if g.degree(var) == h.degree(): H.append(h)
p = min(H, key=lambda h: h.degree()) zeros = p.ground_roots()
for zero in zeros: if not zero.is_Rational: dom_zero = dom.algebraic_field(zero) else: dom_zero = dom
_solutions.add(((zero,) + values, dom_zero))
solutions = _solutions
solutions = list(solutions)
for i, (solution, _) in enumerate(solutions): solutions[i] = solution
return sorted(solutions, key=default_sort_key) |