Solve a system of quadratic equations

quad_system_scipy.py

from __future__ import division, print_function
import numpy as np
from scipy.optimize import fsolve


def fun(X):
    x, y, z = X
    return [x**2 + y + z - 1,
               x + y**2 + z - 1,
               x + y + z**2 - 1]


def jac(X):
    x, y, z = X
    return np.array([
            [2*x, 1, 1],
            [1, 2*y, 1],
            [1, 1, 2*z]])


def close_in(points, point, tol=1e-6):
    nvals = len(points)
    dists = [np.allclose(point, points[cont], atol=tol)
             for cont in range(nvals)]
    return True in dists


npts = 20
tol = 1e-4
np.random.seed(seed=3)
initial_guesses = np.random.uniform(-10, 10, (npts, 3))
sols = []
for x0 in initial_guesses:
    sol, info, ler, msg = fsolve(fun, x0, fprime=jac, full_output=True,
                                 xtol=tol)
    if not close_in(sols, sol, tol):
        sols.append(sol)

print(np.round(np.asarray(sols), 3))

Greetings Nico, can you please direct me to a method for transforming a system of quadratic equations into this matrix form?

Thank you!

Greetings Nico, can you please direct me to a method for transforming a system of quadratic equations into this matrix form?

Since it is a system of equations you can write it in vector form

$$\mathbf{F}(\mathbf{x}) = (f_0(\mathbf{x}), f_1(\mathbf{x}), \cdots, f_n(\mathbf{x})), ,$$

and that's the required form for Newton's method.

The matrix used in the code is the Jacobian Matrix. You compute it by arranging the partial derivatives of the components of your vector function into a matrix. I'd suggest that you use SymPy for this purpose to avoid mistakes.