Solving a system of quadratic equations in Python

Question

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?

Answer 1

Why don't you use regular Newton? Your system is simple enough that you can find its closed-form Jacobian and write your own Newton solver. If you just need one solution which is close to a given starting point (like you wrote on MO), then it rates to work pretty well.

(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$.)