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

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. $$v_{E_{i,j}} = v_{F_{i,j}}$$ 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?

• 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. – Geoff Oxberry Jun 25 '12 at 23:44
• @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. – Damien Jun 26 '12 at 3:02