Numerically stable real solution(s) to a system of bivariate quadratics

Question

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.

Example:

In this instance, the real solutions are:

Solutions: 
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$):

3.96025792973163717278031784693e-13
-9.33800975227555950436364008271e-10
9.63557996744215866033970248179e-07
-0.000568300516844107153566292037895
0.209542890052501016806589147067
-49.4609340406527708387460620101
7298.7082158003621329697479114
-615609.87547890187406435557892
22722530.2050891103731763098963

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

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.