Solving a system of quadratic equations in Python

I'd like to solve numerically a system of quadratic equations:

$A_{11}x_1+A_{12}x_2+A_{13}x_3+B_{12}x_1x_2+B_{13}x_1x_3=C_1$ $A_{21}x_1+A_{22}x_2+A_{23}x_3+B_{21}x_2x_1+B_{23}x_2x_3=C_2$ $A_{31}x_1+A_{32}x_2+A_{33}x_3+B_{31}x_3x_1+B_{32}x_3x_2=C_3$

where $A_{ij}$, $B_{ij}$, $C_i$ are real numbers and $x_i$ the variables (here only three but I have in mind to solve systems up to 50 variables).

I've tried Newton-Krylov method (implemented in scipy) but it does not seem to be stable. I'm now trying to use routines implemented in sympy but I face a problem illustrated by this example:

from sympy import *
from sympy.abc import x,y,z
Equations = [x**2 + y + z - 1,
x + y**2 + z - 1,
x + y + z**2 - 1]


If we use solve_poly system we obtain

solve_poly_system(Equations, x, y, z)
[(0, 0, 1), (0, 1, 0), (1, 0, 0), (-1 + sqrt(2), -1 + sqrt(2), -1 + sqrt(2)), (-sqrt(2) - 1, -sqrt(2) - 1, -sqrt(2) - 1)]


and, if we use solve_triangulated we obtain

solve_triangulated(Equations, x, y, z)
[(0, 0, 1), (0, 1, 0), (1, 0, 0)]


Some solutions are not found by solve_triangulated. What is the reason? What kind of solutions can we find with solve_triangulated?

• If your are planning to solve for 50 variables, an analytical solution is probably not the way to go. Regarding the differences you see between solve_poly_system and solve_triangulated you should ask in a SymPy channel like the gitter chat room or the mailing list. – nicoguaro Feb 22 '18 at 13:21
• Regarding the use of SciPy, I just found the roots using fsolve. You can check the gist here: gist.github.com/nicoguaro/3c964ae8d8f11822b6c22768dc9b3716 – nicoguaro Feb 22 '18 at 14:07
• Thanks a lot! I've just made some tests with random systems and I'm surprised to see that fsolve is much more stable than newton_krylov (which often cannot find a solution and raise an error even for small systems). Could be due to the fact that I don't know how to parametrize newton_krylov. I will continue with fsolve, it can deal apparently with (quadratic) non linear systems with 100 variables in a few minutes. – JaneFlo Feb 22 '18 at 16:11
• After many tests, it seems that scipy.optimize.root with method=lm and explicit jacobian in input is the best solver for my specific problem (quadratic non linear systems with a few dozens of equations). – JaneFlo Mar 2 '18 at 13:18

(As I wrote on MO, I guess that there can be up to $2^{\text{number of variables}}$ real solutions, so finding all of them is going to be problematic for $n\approx 50$.)
• Thanks for the suggestion. I've implemented it but I face the problem of convergence/divergence of the Newton method. There are some solutions which are indeed found very quickly (the algorithm converges) and there are other solutions which cannot be found even if I give in input an initial starting point which is very close to the solution (because the algorithm diverges). Since scipy.optimize.root works well for my problem, I will continue with it (don't think I'll be able to do better with my limited computational skills ...) – JaneFlo Mar 2 '18 at 13:24