How does SymPy's solve() find solutions to nonlinear systems?

user_14073182 0
Viewed 189

The SymPy symbolic mathematics library for Python (see) offers a solvers module for solving both singular equations and systems of equations. An example of usage follows:

import sympy as sym
x,y,z = sym.symbols('x,y,z')
c1 = sym.Symbol('c1')
eq1 = sym.Eq(2*x**2+y+z,1)
eq2 = sym.Eq(x+2*y+z,c1)
eq3 = sym.Eq(-2*x+y,-z)
result = sym.solve([eq1,eq2,eq3],(x,y,z))
print(result)

I'm using solve() after unsuccessfully trying a breadth of custom methods to solve the following nonlinear system. It works great, however I cannot find any information on how SymPy actually solves these. After trying for such time to implement my own solver (just as an interesting exercise), I'm very interested in learning how this is accomplished.

            (x - x_i) ** 2 + (y - y_i) ** 2 + (z - z_i) ** 2 - (299792458 * (t_i - t)) ** 2
            (x - x_j) ** 2 + (y - y_j) ** 2 + (z - z_j) ** 2 - (299792458 * (t_j - t)) ** 2
            (x - x_k) ** 2 + (y - y_k) ** 2 + (z - z_k) ** 2 - (299792458 * (t_k - t)) ** 2
            (x - x_m) ** 2 + (y - y_m) ** 2 + (z - z_m) ** 2 - (299792458 * (t_m - t)) ** 2

Where x, y, z and t are unknown, and I'm aiming to find x,y,z

solver sympy python
1 Answers

The example systems you show are polynomial in the unknowns. The general framework for solving a polynomial system is given by the Gröbner basis: https://en.wikipedia.org/wiki/Gr%C3%B6bner_basis

You can compute this in sympy:

In [14]: ps = [e.rewrite(Add) for e in (eq1, eq2, eq3)]                                                                           

In [15]: groebner(ps, [x, y, z])                                                                                                  
Out[15]: 
             ⎛⎡                                      2              2                     ⎤                                  ⎞
GroebnerBasis⎝⎣-c₁ + 5⋅x - z, -2⋅c₁ + 5⋅y + 3⋅z, 2⋅c₁  + 10⋅c₁ + 2⋅z  + z⋅(4⋅c₁ + 10) - 25⎦, x, y, z, domain=ℤ[c₁], order=lex⎠

The idea is that this presents a form in which it is hopefully easiest to work through solving the equations one by one, analogous to the back-substitution step of Gaussian elimination. Here in the basis we see that there are 3 polynomials and one of them only involves z. Since it is quadratic we can apply the formula to solve for two possible values of z. Substituting that into the other two polynomials gives separate linear equations for x and y.

Computing the Gröbner basis can be extremely slow for larger systems. As the degree of the polynomials or the number of unknowns increases it also gets very difficult to solve anything using any method but there are certainly faster methods for special cases.

Note that sympy's polys module is huge and has a lot of code that all comes together to make it possible to compute a Gröbner basis and also to find the roots of a polynomial.

That being said I think that there are improvements that could be made in solve for these kinds of systems. Ideally solve should first solve the linear equations and eliminate as many variables as possible. It's also potentially more efficient to use polynomial resultants: https://en.wikipedia.org/wiki/Resultant

Also worth noting is the Abel Ruffini theorem which gives a strong limit on what any analytic solver for polynomial equations can achieve: https://en.wikipedia.org/wiki/Abel%E2%80%93Ruffini_theorem

user_94509910
Related