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?

  • 2
    $\begingroup$ 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. $\endgroup$
    – nicoguaro
    Feb 22, 2018 at 13:21
  • $\begingroup$ Regarding the use of SciPy, I just found the roots using fsolve. You can check the gist here: gist.github.com/nicoguaro/3c964ae8d8f11822b6c22768dc9b3716 $\endgroup$
    – nicoguaro
    Feb 22, 2018 at 14:07
  • 2
    $\begingroup$ 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. $\endgroup$
    – JaneFlo
    Feb 22, 2018 at 16:11
  • $\begingroup$ 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). $\endgroup$
    – JaneFlo
    Mar 2, 2018 at 13:18

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

  • $\begingroup$ 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 ...) $\endgroup$
    – JaneFlo
    Mar 2, 2018 at 13:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.