Skip to content

Instantly share code, notes, and snippets.

What would you like to do?
Solve a system of quadratic equations
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
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,
if not close_in(sols, sol, tol):
print(np.round(np.asarray(sols), 3))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.