I have a a system of bivariate polynomials as follows:

$ E(u,v): e_2(u) v^2 + e_1(u) v + e_0(v) = 0 \\ F(u,v): f_2(u) v^2 + f_1(u) v + f_0(v) = 0$

where $e_n(u) = e_{n_2}u^2 + e_{n_1} u + e_{n_0}$ and similarly $f_n(u) = f_{n_2}u^2 + f_{n_1} u + f_{n_0}$

Thus far, I have solved via a Sylvester Matrix:

Let the Sylvester Matrix $\mathbf{G}$ be defined as follows: \begin{align} \mathbf{G} = \begin{bmatrix} e_{2}(u) & e_{1}(u) & e_{0}(u) & 0 \\ 0 & e_{2}(u) & e_{1}(u) & e_{0}(u) \\ f_{2}(u) & f_{1}(u) & f_{0}(u) & 0 \\ 0 & f_{2}(u) & f_{1}(u) & f_{0}(u) \end{bmatrix} \end{align}

Possible solutions for $u$ may be found by finding the determinant of $\mathbf{G}$, which will be a polynomial of degree 8, i.e.: \begin{align} G_u &= \det\left( \mathbf{G} \right) \\ &= \sum_{i=0}^{8} G_{u_i} u^i \end{align}

Real roots of $G_u$ may or may not be solutions of $E$ and $F$. To determine whether each real root, $u_i$ is a solution, substitute $u_i$ into $E$ and $F$ and solve for roots $v_{E_{i,j}}, j \in \{1, 2\}$ and $v_{F_{i,j}}, j \in \{1, 2\}$ respectively. Where there is a common real solution to $v$, i.e. \begin{equation} v_{E_{i,j}} = v_{F_{i,j}} \end{equation} then the pair $\left( u_i, v_{E_{i,j}} \right)$ is a solution.

The problem with this method is that is appears to be horribly numerically unstable. Ultimately, I needed to used MPFR and Eigen with a crazy number of bits (I think 500 or so) to get correct answers - double was nowhere near good enough.

In particular, fine variations of the coefficients of $G(u)$ caused by rounding or other numerical stability issues significantly moved the zeros of the polynomial, or in some instances, caused no real solutions at all.


In this instance, the real solutions are:

u = 283.118152022138386048081228332, v = 253.676766248666064015102940087
u = 283.531719376538811255222357595, v = 251.92153962772290946172986145
u = 305.922640195704152315990688064, v = 277.085127660836580117988944869
u = 307.305710607492739248674271296, v = 275.328781581725038016656158251

Where the original bivariate quadratic is:

e_2_2 = 2.96184510308460152308099163079e-05
e_2_1 = -0.0174923195188134962339249161373
e_2_0 = 2.58392910782947268103166838601

e_1_2 = -0.0156683206373072309167532302994
e_1_1 = 9.24288932470428612919593805961
e_1_0 = -1363.76987336669518990414286879

e_0_2 = 2.07344822818614405777141090841
e_0_1 = -1221.73998971149758113804476542
e_0_0 = 180057.979202842850602362010181

f_2_2 = 0.0057280951738335050897549908978
f_2_1 = -3.37584828933150293352939673891
f_2_0 = 497.643484709582863520217198114

f_1_2 = -3.02706498792933011045380441656
f_1_1 = 1781.95310671648910949906364811
f_1_0 = -262381.909585190823129099681101

f_0_2 = 400.16454659182969104246501941
f_0_1 = -235295.984529511565861184026453
f_0_0 = 34606200.151155688437914904126

And the $G(u)$ coefficients are (from $G_8 u^8 + ... + G_0$):


What techniques (or code ;) ) will help me out?

  • $\begingroup$ Aren't there methods to deal with these polynomials more directly? Algebraic geometry packages also implement methods to solve these sorts of systems. These problems are outside of my area of expertise, but one package I know of is Bertini; I'm sure there are others. $\endgroup$ Jun 25, 2012 at 23:44
  • $\begingroup$ @GeoffOxberry, the motivation for this problem is actually the intersection of three quadric surfaces. This particular problem is a sub-set of Xu's method. Others in the field are Lazard/Dupont, or methods derived from Levin. $\endgroup$
    – Damien
    Jun 26, 2012 at 3:02

1 Answer 1


To find all the solutions of a polynomial system of equations with the usual accuracy expectable from floating-point calculations, use a homotopy algorithm that exploits the polynomial structure. See, e.g.,

A Morgan and A. Sommese, Computing all solutions to polynomial systems using homotopy continuation, Applied Mathematics and Computation 24 (1987), 115-138.

and there is an associated software package called HOMPACK: (source: http://www.netlib.org/hompack/ ; paper: http://dl.acm.org/citation.cfm?id=279235 )

Alternatively, use a branch and bound package such as ICOS https://sites.google.com/site/ylebbah/icos . This gives guarantees even in the presence of rounding errors, and you can add inequality constraints to restrict to a desired domain.

  • $\begingroup$ Did you have a particular paper in mind? $\endgroup$
    – Damien
    Jun 26, 2012 at 22:27
  • $\begingroup$ I added more details to my answer. $\endgroup$ Jun 27, 2012 at 14:04

Your Answer

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

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