Dixon Resultant Solution of Nonlinear Systems of Equations in Geodesy

Paláncz (B^1)., Awange (J^2)., Zaletnyik (P^3). and Grafarend (E^4).

Abstract
Dixon's resultant, a powerful mathematical tool is here proposed as an alternative to Groebner basis, Sylvester Resultant, Sturmfels Resultant or Multipolynomial Resultant for solving nonlinear system of equations in geodesy. We demonstrate its power in solving  intersection, GPS ranging and C7 Conformal transformation problems. Its overraiding advantage compared to the other closed form solutions listed above is its small size which enables faster solution in-case of relatively small systems with 3-4 variables. Geodetic users uncomfortable with lengthy expressions of Groebner basis or Multipolynomial Resultant may find it useful.

1. Introduction

Solutions of nonlinear system of equations in geodesy has received wide coverage as evidenced in the works of  Bancroft (1985) , Grafarend and Schaffrin (1993), Singer et al. (1993), Kleusberg (1994, 2003), Lichtenegger (1995), and Awange (2003). The development of Computer Algebraic System CAS to explain algorithms, methods or computation techniques is providing the possibility of exploiting methods which were not achievable in the past largely due the the poor state of computers. Thanks to the development in the computer industry, most of CAS related approaches are now being unlocked. Due to Computer Algebra Systems (CAS) like Mathematica or Maple  the usage of these methods could be very efficient. Although, most of the scientists believe, even today, that CAS language is a programming language, however this is not so. CAS can be used like "live'' mathematics for creating, proving as well evaluating formulas and expressions in numeric or symbolic form, see for example Freeman (1994), Helton and Merino (1998), Bellomo et al. (2000) and Romeny (2003).  

Applications of CAS in Geosciences, Geodesy and Geoinformatics have recently started appearing in scientific literature therefore warranting special consideration. For instance, Haneberg (2004) presents Mathematica notebooks providing instant computational facility for practitioners. Awange and Grafarend (2005) on the other hand have presented several closed form algorithms which are applicable to wide range of geodetic and geoinformatic problems. The advantages of these closed forms solution are well known in geodesy. As an extension to the previous works, this contribution, realising the difficulties of the previous closed form solutions of Groebner basis and Multipolynomial resultants, i.e., lengthy expressions and longer computational time, proposes an alternative approach based on Dixon's resultant. It is hoped that with the addition of computational tool for closed form solution, geodetic community will have a wider range of choice to meet their needs. In this paper we try to prove the utility of one of the above mentioned methods in geodetical problems , the so called Dixon resultant which is not very known in the geodetical society, but is a very powerful tool. First a general review is given about Dixon resultant and after the efficiency of the method is demonstrated in geodetical problems. In the first example, considering GPS ranging we demonstrate how to use Dixon resultant. In the second example, solving intersection problem Dixon resultant is compared with other resultant method, like Strumfels multipolynomial resultants techique (MR). In the third example, Dixon method is compared with the Groebner basis solution, in case of C7 conformal transformation problem. Examples were running on hp workstation xw 4100 with XP operation system, 3 MHz P4 Intel processor and 1 MB RAM.

2. Dixon's Resultant: Basic Concepts

Consider a system of polynomials
ℱ ={ (Subscript[f, 0]),(Subscript[f, 1]),...,(Subscript[f, d])} 
where
(Subscript[f, i])= (∑_(α ∈  (Subscript[𝒜, i])))(Subscript[c, i, α]) (x^α)
and
(x^α) = (Subscript[((Subscript[x, 1])^α), 1]) (Subscript[((Subscript[x, 2])^α), 2]) ...  
          (Subscript[((Subscript[x, 1])^α), d])
for each   i = 0,...,d,  and  ( (Subscript[c, i, α])) are parameters.
The goal is to derive a condition on parameters   ( (Subscript[c, i,  α])) such that the polynomial system ℱ = 0 has a solution. We can consider this problem as the elimination of variables from the polynomial system. Elimination theory tells us that such a condition exists for large family of polynomial systems, and is called the resultant of the polynomial system. Since the number of equations is more than the number of variables, in general, for arbitrary values of (Subscript[c, i,  α]) the polynomial system does not have any solution.
If there is only one such polynomial represent this family, then its defining equation is called the resultant.

One way to compute the resultant of a given polynomial system is to construct a matrix with property that whenever the polynomial system has a solution, such a matrix has a deficient rank, thereby implying that the determinant of any maximal minor is a multiple of the resultant.
A simple way to construct a resultant matrix is to use the dialytic method, i.e. multiply each polynomial with a finite set of monomials, and rewrite the resulting system in matrix form. We call such matrix the dialytic matrix. This alone, however, does not guarantee that a matrix so constructed is a resultant matrix.
Note that such matrices are usually quite sparse: matrix entries are either zero or coefficients of the polynomials in the original system. Good examples of resultant dialytic matrices are Sylvester (Sylvester, 1853) for univariate case, and Macaulay (Macaulay, 1916) as well as Newton sparse  matrices of Canny and Emiris (2000) for the multivariate case; they all differ only in the selection of multiplier monomial sets.

In contrast to dialytic matrices, the Dixon matrix is dense since its entries are determinants of the coefficients of the polynomials in the original system. It has the advantage of being an order of magnitude smaller in comparison to a dialytic matrix, which important as the computation of the determinant of a matrix symbolic entries is sensitive to its size. The Dixon matrix is constructed through the computation of the Dixon polynomial, which is expressed in matrix form.

Let (Subscript[π, i])((x^α)) = \(\(\(σ\_1\^α\_1 ... \) σ\_i\^α\_i\ x\_\(i + 1\)\^α\_\(i + 1\) ... \) 
            x\_d\^α\_d\ where\ i ∈ {0, 1,  ... , 
              d}\ and\ the\ σ\_i' s\ are\ new\ variables\); (Subscript[π, 0])((x^α))= (x^α) (Subscript[π, i])  is extended to polynomials in a natural way as

(Subscript[π, i]) (f ((Subscript[x, 1]),  ... , (Subscript[x, d])))  =  
      f ((Subscript[σ, 1]), ... (Subscript[σ, i]),  (Subscript[x, i + 1]),  ... , 
          (Subscript[x, d]))

Given a polynomial system ℱ = {(Subscript[f, 0]), (Subscript[f, 1]),...(Subscript[f, d])} define its Dixon polynomial as

θ((Subscript[f, 0]),(Subscript[f, 1]),...(Subscript[f, d]))=δ ((Subscript[x, 1]),  ... , 
            (Subscript[x, d]), (Subscript[σ, 1]), ... (Subscript[σ, d]))=(∏_(i = 1)^d)(1/((Subscript[σ, i]) - (Subscript[x, i])))Det(({{(Subscript[π, 0]) ((Subscript[f, 0])), (Subscript[π, 0]) ((Subscript[f, 1])), ., ., ., (Subscript[π, 0]) ((Subscript[f, d]))}, {(Subscript[π, 1]) ((Subscript[f, 0])), (Subscript[π, 1]) ((Subscript[f, 1])), ., ., ., (Subscript[π, 1]) ((Subscript[f, d]))}, {., ., ., ., ., .}, {., ., ., ., ., .}, {., ., ., ., ., .}, {(Subscript[π, d]) ((Subscript[f, 0])), (Subscript[π, d]) ((Subscript[f, 1])), ., ., ., (Subscript[π, d]) ((Subscript[f, d]))}}))

The order in which original variables in x are replaced by new variables in σ is significant in the sense that Dixon polynomial computed using two different variable orderings may be different.

A Dixon polynomial δ(x, σ) can be witten in bilinear form  as

\(δ \((x\_1,  ... , 
            x\_d, σ\_1, \(\(...\) \(σ\_d\)\))\) = Ξ\ Θ\ X\^T\)

where Ξ = ((σ^((Subscript[β, 1]))),...,(σ^((Subscript[β, k])))) and X = ((x^((Subscript[α, 1]))),...,(x^((Subscript[α, s])))) are row vectors. The k × s matrix Θ is called the Dixon matrix.

3. Construction of  Dixon Resultant

3.1 Cayley's Formulation of Bézout's Method

Let us recall here the Cayley's formulation Cayle (1865) of Bezout's method for solving two polynomial equations in univariate case. However that method is actually due to Euler!

Let us consider two univariate polynomial f(x) and g(x), and let

deg = max (degree(f),degree(g))

and let σ be an auxiliary variable. The quantity

δ (x, σ) =(1/(x - σ))det(({{f (x), g (x)}, {f (σ), g (σ)}}))=((f (x) g (σ) - 
              f (σ) g (x))/(x - σ))

is a symmetric polynomial in x and σ of deg-1 which is called Dixon polynomial of f and g. Every common zero of f and g is a zero of δ(x,σ) for all values of σ.

Example

Let us consider two univariate polynomials with parameter π,

f (x) = (x - 1) (x + 3) (x - 4) = 
        12 - 11 x - 2 (x^2) + (x^3)

g (x) = (x - π) (x + 4) = -4 π + 
          4 x - π x + (x^2)

Then the Dixon polynomial is

δ=(1/(x - σ))det((({{f (x), g (x)}, {f (σ), g (σ)}})))=-48 + 56 π - 12 x + 8 π x - 4 π (x^2) - 
          12 σ + 8 π σ + 3 x σ - 
          2 π x σ + 4 (x^2) σ - π (x^2) σ - 
          4 π (σ^2) + 4 x (σ^2) - π x (σ^2) + 
          (x^2) (σ^2)

Hence at a common zero, each coefficient of (σ^i) in δ is vanishing

(σ^2): -4 π + 4 x - π x + 
          (x^2)=0

(σ^1): -12 + 8 π + 3 x - 
          2 π x + 4 (x^2) - π (x^2)=0

(σ^0):-48 + 56 π - 12 x + 
          8 π x - 4 π (x^2)=0

It is a homogeneous system in variables  (x^0), (x^1) and (x^2).

(({{-4 π, 4 - π, 1}, {-12 + 8  π, 3 - 2  π, 4 - π}, {-48 + 56 π, -12 + 
                    8  π, -4 π}}))(({{(x^0)}, {(x^1)}, {(x^2)}}))=(({{0}, {0}, {0}}))

This system has non-trivial solutions if and only if its determinant, D is zero. D is called the Dixon resultant of f and g. The matrix of the system M is the Dixon matrix,

M =(({{-4 π, 4 - π, 1}, {-12 + 8 π, 3 - 2  π, 4 - π}, {-48 + 56 π, -12 + 
                8  π, -4 π}}))

and its determinant is

D = -480 + 440 π + 80 (π^2) - 40 (π^3)

The polynomials have common zeros if  D is vanishing. Really, solving the equation D = 0 we get

(Subscript[π, 1]) = -3, (Subscript[π, 2]) = 1 and  (Subscript[π, 3]) = 4

3.2 Dixon's Generalization of the Cayley- Bézout Method

Dixon generalized Cayley's approach to Bezout's method to systems of three polynomials equations in three unknowns, Dixon (1908). Let

f(x,y,z) = 0

g(x,y,z) = 0

h(x,y,z) = 0

Now the Dixon polynomial is defined by

δ (x, y, z, σ, ξ)=(1/((x - σ) (y - ξ)))det(({{f (x, y, z), g (x, y, z), h (x, y, 
                    z)}, {f (σ, y, z), g (σ, y, 
                    z), h (σ, y, z)}, {f (σ, ξ, z), g (σ, ξ, 
                    z), h (σ, ξ, z)}}))

Example

Let us consider three threevariate polynomials

f=(x^2) + (y^2) - 1

g=(x^2) + (z^2) - 1

h=(y^2) + (z^2) - 1

The Dixon polynomial is

δ(x,y,z,σ,ξ ) =(1/((x - σ) (y - ξ))) det(({{(x^2) + (y^2) - 1, (x^2) + (z^2) - 1, (y^2) + (z^2) - 1}, {(σ^2) + (y^2) - 1, (σ^2) + (z^2) - 1, (y^2) + 
                (z^2) - 1}, {(σ^2) + (ξ^2) - 1, (σ^2) + (z^2) - 
                1, (ξ^2) + (z^2) - 1}}))

=y z - 2 (x^2) y z + y ξ - 2 (x^2) y ξ + 
          z σ - 2 (x^2) z σ + ξ σ - 
          2 (x^2) ξ σ

Now we  eliminate the variables  y and z, so considering that

δ (x, y, 
          z, σ, ξ) = σ ξ (1 - 
              2 (x^2)) + ξ y (1 - 2 (x^2)) + σ z (1 - 
              2 (x^2)) + y z (1 - 2 (x^2))

and  system of equations is

({{., (y^0) (z^0), (y^1) (z^0), (y^0) 
            (z^1), (y^1) (z^1)}, {(σ^0) (ξ^0):, 0, 0, 0, 1 - 2 (x^2)}, {(σ^1) (ξ^0):, 0, 0, 1 - 2 (x^2), 0}, {(σ^0) (ξ^1):, 0, 1 - 2 (x^2), 0, 0}, {(σ^1) (ξ^1):, 1 - 2 (x^2), 0, 0, 0}})

Consequently the Dixon matrix

M=(({{0, 0, 0, 1 - 2 (x^2)}, {0, 0, 1 - 2 (x^2), 0}, {0, 1 - 2 (x^2), 0, 0}, {1 - 2 (x^2), 0, 0, 0}}))

and its determinant, the Dixon resultant

D=1 - 8 (x^2) + 24 (x^4) - 32 (x^6) + 16 (x^8)

Dixon proved that for three 2 degree polynomials, the vanishing of D is a necessary condition for the existence of a common zero. Furthermore, D is not identically zero. Dixon's method and proofs easily generalize to a system of n+1 generic ndegree polynomials in n unknowns. Recall, that a polynomial is generic if all its coefficients are independent parameters, unrelated to each other. A polynomial in n variables is ndegree if all powers to the maximum of each variable appear in it, see Kapur, Saxena and Yang (1994) for more details.

4. Application of Dixon Resultant

Dixon's method applies only to generic ndegree polynomials. If this condition fails then one can face the following problems:

1. The Dixon matrix may be singular

2. After we remove the rows  and columns containig zeros, the vanishing of the determinant of the Dixon matrix may not give a necessary condition for the existence of a common zero.

3. After we remove the zero rows and columns containig zeros, the Dixon matrix may not even be square. Hence, its determinant cannot be defined.

Kapur, Sexana and Yang addressed all three problems successfully, provided a certain precondition holds.
 Assuming that the column that corresponds to the monomial 1 = ((Subscript[x, 1])^0)((Subscript[x, 2])^0)... of the Dixon matrix is not linear combination of the remaining ones. If this precondition is true, then D = 0 is a necessary condition of common zeros. This theorem yields a simple algorithm for obtaining the necessary condition D = 0, which is called Kapur -Sexana -Yang Dixon resultant. The algorithm is the following:
 
 1. Compute the Dixon matrix M. If the precondition holds continue
 2. Row reduce M without scaling to row echelon form, M'.
 3. Compute the product D of the pivots of M'.

This algorithm was implemented into Mathematica by Nakos and Williams (1997), and  MathSource (2002).
Let us employ it for the last example,

<< Resultant`Dixon`

\(\(\(\(p = x\^2 + y\^2 - 1\) \)\(;\)\(\(q = 
        x\^2 + z\^2 - 1\) \)\(;\)\(\(r = 
        y\^2 + z\^2 - 1\) \)\(;\)\(
\)\)\)

DixonPolynomial[{p, q, r}, {y, z}, {σ, ξ}]

M = DixonMatrix[{p, q, r}, {y, z}, {σ, ξ}]; 
    MatrixForm[M]

DixonResultant[{p, q, r}, {y, z}, {σ, ξ}]

After seeing the general functioning scheme of Dixon resultant let us solve some tipical geodetic problems with it.

Clear["Global`*"]

Example 1. Ranging by GPS

Throughout history, position determination has been one of the most important tasks of mountaineers, pilots, sailor, civil engineers etc. In modern times, Global Positioning System (GPS) employing Global Navigation Satellite Systems (GNSS) provide an ultimate method to accomplish this task. If one has a hand held GPS receiver, the receiver measures the travel time of the signal transmitted from the satellites. Then this distance can be computed by multiplying the measured time by the speed of light in vacuum. The distance of the receiver from the i-th satellite, (Subscript[d, i]) is related to the unknown position of the receiver, {X,Y, Z}

\(d\_i = \@\(\((a\_i - X)\)\^2 + \((b\_i - Y)\)\^2 + \((c\_i - Z)\)\^2\) + \ \ ξ\)

where {(Subscript[a, i]),(Subscript[b, i]),(Subscript[c, i])},i =0.  .. 3 are the coordinates of the i-th satellite.

The distance is influenced also by the satellite and receiver' clock biases. The satellite clock biases can be modelled while the receiver' clock biases have to be considered as an unknown variable, ξ. This means, we have four unknowns, consequently we need four satellite signals as minimum observation. Let us employ (Subscript[x, 1]), (Subscript[x, 2]), (Subscript[x, 3]) and (Subscript[x, 4]) variables for the four unknowns, {X,Y, Z, ξ,} then our equations are, (Subscript[e, i]) = 0, i = 1...4, where

e1 = (((x1 - (Subscript[a, 0])))^2) + (((x2 - (Subscript[b, 0])))^2) + (((x3 - (Subscript[c, 0])))^2) - (((x4 - (Subscript[d, 0])))^2);

e2 = (((x1 - (Subscript[a, 1])))^2) + (((x2 - (Subscript[b, 1])))^2) + (((x3 - (Subscript[c, 1])))^2) - (((x4 - (Subscript[d, 1])))^2);

e3 = (((x1 - (Subscript[a, 2])))^2) + (((x2 - (Subscript[b, 2])))^2) + (((x3 - (Subscript[c, 2])))^2) - (((x4 - (Subscript[d, 2])))^2);

e4 = (((x1 - (Subscript[a, 3])))^2) + (((x2 - (Subscript[b, 3])))^2) + (((x3 - (Subscript[c, 3])))^2) - (((x4 - (Subscript[d, 3])))^2);

First, this equation  system of polynomials, will be transformed into  a system of linear equations and a quadratic equation. Let us expand and multiply by minus one, and arranged the original equations,

eqsL = {e1L, e2L, e3L, e4L} = 
        Map[-Sort[Expand[#]] &, {e1, e2, e3, e4}]

Subtruct the fourth equation from the other three ones,

q = 
      Table[eqsL[[i]] - eqsL[[4]], {i, 1, 3}] // Simplify

This is a system of three linear equations, and they can be written as,

(Subscript[g, 1]) = (Subscript[a, 0, 3]) x1 + (Subscript[b, 0, 3]) x2 + (Subscript[c, 0, 3]) 
            x3 + (Subscript[d, 3, 0]) x4 + (Subscript[e, 0, 3]);

(Subscript[g, 2]) = (Subscript[g, 1]) /. {1 -> 2, 0 -> 1}

(Subscript[g, 3]) = (Subscript[g, 1]) /. {1 -> 3, 0 -> 2}

The coefficients {(Subscript[a, i, 3]), (Subscript[b, i, 3]), 
          (Subscript[c, i, 3]), (Subscript[d, 3, i]), (Subscript[e, i, 3])}, i = 0.  ..  2, can be determined  as,

coeffs0 = 
        Table[{Coefficient[q[[i + 1]], x1], 
              Coefficient[q[[i + 1]], x2], 
              Coefficient[q[[i + 1]], x3], 
              Coefficient[q[[i + 1]], x4]} // Factor, {i, 0, 
            2}];

which are the coefficients of the variables {(Subscript[x, 1]),(Subscript[x, 2]), (Subscript[x, 3]), (Subscript[x, 4])}. The constant part is,

coeffs1 = 
        Table[{q[[i]] - coeffs0[[i]] . {x1, x2, x3, x4} // 
              Simplify}, {i, 1, 3}];

Therefore, all of the coefficients are,

coeffs = 
      Table[Union[coeffs0[[i]], coeffs1[[i]]], {i, 1, 3}]

Let us assign these coefficients to the linear system,

coeffsn = 
      Flatten[Table[
          Inner[#1 -> #2 &, {(Subscript[a, i, 3]), (Subscript[b, i, 3]), (Subscript[c, i, 3]), 
              (Subscript[d, 3, i]), (Subscript[e, i, 3])}, coeffs[[i + 1]], List], {i, 0, 
            2}]]

In addition, we take one of the nonlinear equations, let say the fourth one,

(Subscript[e, 4]) = (((x1 - (Subscript[a, 3])))^2) + (((x2 - (Subscript[b, 3])))^2) + (((x3 - (Subscript[c, 3])))^2) - (((x4 - (Subscript[d, 3])))^2);

Now, we shall solve the linear system for the variables {(Subscript[x, 1]),(Subscript[x, 2]), (Subscript[x, 3])}, with (Subscript[x, 4]) as parameter. It means, the relations (Subscript[x, 1])= g((Subscript[x, 4])), (Subscript[x, 2])= g((Subscript[x, 4])) and (Subscript[x, 3])= g((Subscript[x, 4])) will be computed. To do that, different methods can be employed. Applying Kapur -Sexana -Yang Dixon resultant   algorithm

drx1 = 
      DixonResultant[{(Subscript[g, 1]), (Subscript[g, 2]), (Subscript[g, 3])}, {x2, x3}, {u2, u3}]

This is a linear expression containing only  (Subscript[x, 1]) and (Subscript[x, 4]),

Exponent[drx1, {x1, x2, x3, x4}]

Then the solution  for (Subscript[x, 1])= g((Subscript[x, 4]))

solx1 = Solve[drx1 == 0, x1]

Similarly, for the two additional variables, (Subscript[x, 2])= g((Subscript[x, 4])) and (Subscript[x, 3])= g((Subscript[x, 4])),

drx2 = 
      DixonResultant[{(Subscript[g, 1]), (Subscript[g, 2]), (Subscript[g, 3])}, {x1, x3}, {u1, u3}]

Exponent[drx2, {x1, x2, x3, x4}]

solx2 = Solve[drx2 == 0, x2]

and

drx3 = 
      DixonResultant[{(Subscript[g, 1]), (Subscript[g, 2]), (Subscript[g, 3])}, {x1, x2}, {u1, u2}]

Exponent[drx3, {x1, x2, x3, x4}]

solx3 = Solve[drx3 == 0, x3]

After substitution, we have again a quadratic equation for (Subscript[x, 4]),

G = 
        (Subscript[e, 4]) /. {solx1[[1, 1]], solx2[[1, 1]], 
            solx3[[1, 1]]};

Exponent[G, x4, List]

The coefficents of the quadratic equation are,

h2d = Coefficient[G, (x4^2)];

h1d = Coefficient[G, x4];

h0d = Simplify[G - (h2d (x4^2) + h1d x4)];

Let us suppose that the observation data are, see Grafarend and Shan (2005),

data = {(Subscript[a, 0]) -> 1.483230866 (10^7), 
          (Subscript[a, 1]) -> -1.579985405 (10^7), 
          (Subscript[a, 2]) -> 1.98481891 (10^6), 
          (Subscript[a, 3]) -> -1.248027319 (10^7), 
(Subscript[b, 0]) -> -2.046671589 (10^7), (Subscript[b, 1]) -> -1.330112917 (10^7), 
          (Subscript[b, 2]) -> -1.186767296 (10^7), 
(Subscript[b, 3]) -> -2.338256053 (10^7), 
(Subscript[c, 0]) -> -7.42863475 (10^6), (Subscript[c, 1]) -> 1.713383824 (10^7), 
          (Subscript[c, 2]) -> 2.371692013 (10^7), 
(Subscript[c, 3]) -> 
            3.27847268 (10^6),  
(Subscript[d, 0]) -> 
            2.4310764064 (10^7), (Subscript[d, 1]) -> 2.2914600784 (10^7), 
          (Subscript[d, 2]) -> 2.0628809405 (10^7), 
(Subscript[d, 3]) -> 
            2.3422377972 (10^7)};

The numerical evaluations are following. The variable x4 is,

solx4 = Solve[h2d (x4^2) + h1d x4 + h0d == 0, x4] /. 
          coeffsn /. data

The other variables can be computed using the second value  of x4. Then x1 is

x1 /. solx1 /. coeffsn /. data /. solx4[[2]]

and x2

x2 /. solx2 /. coeffsn /. data /. solx4[[2]]

in addition x3

x3 /. solx3 /. coeffsn /. data /. solx4[[2]]

Example 2. Intersection

In intersection method, angular (direction) observations are taken like in resection method, but in this case measurement are taken from the known points ((Subscript[P, 1]),(Subscript[P, 2]), (Subscript[P, 3])) to the unknown point ((Subscript[P, 0])) while in resection the situation is just the opposite, the unknown point is occupied and directions are measured to the known points.

The observation equations are

((Subscript[x, 2])^2) = 
      ((Subscript[x, 1])^2) + ((Subscript[S, 12])^2) - 
        2 (Subscript[x, 1]) (Subscript[S, 12]) cos ((Subscript[ψ, 12]))
((Subscript[x, 3])^2) = 
      ((Subscript[x, 2])^2) + ((Subscript[S, 23])^2) - 
        2 (Subscript[x, 2]) (Subscript[S, 23]) cos ((Subscript[ψ, 23]))
((Subscript[x, 1])^2) = 
      ((Subscript[x, 3])^2) + ((Subscript[S, 31])^2) - 
        2 (Subscript[x, 3]) (Subscript[S, 31]) cos ((Subscript[ψ, 31]))

Employing the following substitution rule,

r1 = {a0 -> (S12^2), b0 -> (S23^2), c0 -> (S31^2), 
          b1 -> -2 S12 cosφ12, 
          b2 -> -2 S23 cosφ23, 
          b3 -> -2 S31 cosφ31};

the observation equations are

f1 = (x1^2) + b1 x1 - (x2^2) + a0 == 0;

f2 = (x2^2) + b2 x2 - (x3^2) + b0 == 0;

f3 = (x3^2) + b3 x3 - (x1^2) + c0 == 0;

The Dixon polynomial is

DixonPolynomial[{f1[[1]], f2[[1]], f3[[1]]}, {x2, 
            x3}, {X2, X3}] // Simplify // Timing

Let us compute the Dixon matrix

Timing[𝒟 = 
          DixonMatrix[{f1[[1]], f2[[1]], f3[[1]]}, {x2, x3}, {X2, 
              X3}];]

which is quite dense,

MatrixForm[𝒟]

In order to compare the Dixon matrix with that of  an other MR method, we solve the problem with Sturmfels method, Sturmfels (1994). First we should homogenized the equations (f1, f2, f3), namely hiding x1, we get

f11 = -(x2^2) + (a0 + b1 x1 + (x1^2)) (x4^2)

f21 = (x2^2) + b2 x2 x4 - (x3^2) + b0 (x4^2)

f31 = (x3^2) + b3 x3 x4 + (-(x1^2) + c0) (x4^2)

The corresponding coefficients are

{a11, a12, a13, a14, a15, a16} = 
      Coefficient[
        f11, {(x2^2), (x3^2), (x4^2),  x2 x3,  x2 x4, x3 x4}]

{a21, a22, a23, a24, a25, a26} = 
      Coefficient[
        f21, {(x2^2), (x3^2), (x4^2),  x2 x3,  x2 x4, x3 x4}]

{a31, a32, a33, a34, a35, a36} = 
      Coefficient[
        f31, {(x2^2), (x3^2), (x4^2),  x2 x3,  x2 x4, x3 x4}]

Now we compute the Jacobian determinant

J=Det[(({{(∂_x2) f11, (∂_x3) 
                    f11, (∂_x4) f11}, {(∂_x2) f21, (∂_x3) 
                    f21, (∂_x4) f21}, {(∂_x2) f31, (∂_x3) 
                    f31, (∂_x4) f31}}))]//Simplify

and its partial derivatives are

Jx2 = D[J, x2]

Jx3 = D[J, x3]

Jx4 = D[J, x4]

The corresponding coeffients

{b11, b12, b13, b14, b15, b16} = 
      Coefficient[
        Jx2, {(x2^2), (x3^2), (x4^2),  x2 x3,  x2 x4, x3 x4}]

Let us just check this result

(Jx2 - {b11, b12, b13, b14, b15, b16} . {(x2^2), (x3^2), (x4^2),  
              x2 x3,  x2 x4, x3 x4}) // Simplify

and

{b21, b22, b23, b24, b25, b26} = 
      Coefficient[
        Jx3, {(x2^2), (x3^2), (x4^2),  x2 x3,  x2 x4, x3 x4}]

{b31, b32, b33, b34, b35, b36} = 
      Coefficient[
        Jx4, {(x2^2), (x3^2), (x4^2),  x2 x3,  x2 x4, x3 x4}]

The final step in computing the resultant involves the computation of the determinant of a 6 × 6 matrix given by

𝒮=(({{a11, a12, a13, a14, a15, a16}, {a21, a22, a23, a24, a25, a26}, {a31, a32, a33, a34, a35, a36}, {b11, b12, b13, b14, b15, b16}, {b21, b22, b23, b24, b25, b26}, {b31, b32, b33, b34, b35, b36}}));

MatrixForm[𝒮]

If you compare the two matrices, Dixon has a size 4 × 4, and its density ( the number of nonzero elements/ the number of all elements) is 0.87, while the  Sturmfels matrix has a size 6 × 6 and its density is 0.639.

The Dixon resultant is

Timing[ℛ = 
          Det[𝒟] // Simplify;]

ℛ

while the Sturmfels resultant is,

Timing[ℛℯ𝓈\[ScriptTwo]\[ScriptTwo]\[ScriptTwo] = Det[𝒮] // Simplify;]

ℛℯ𝓈\[ScriptTwo]\[ScriptTwo]\[ScriptTwo]

However, using directly Kapur -Sexana -Yang Dixon resultant to eliminate the variables x2 and x3

Timing[drx1 = 
            DixonResultant[{f1[[1]], f2[[1]], f3[[1]]}, {x2, 
                  x3}, {X2, X3}] // Simplify;]

we get

drx1

which is the same as the classical Dixon resultant

Simplify[drx1 - ℛ]

But the computation time is considerably longer!

It is a monomial of variable (Subscript[x, 1])

Exponent[drx1, x1]

for the Sturmfels resultant

Exponent[ℛℯ𝓈\[ScriptTwo]\[ScriptTwo]\[ScriptTwo], x1]

The numerical values for the example, see Awange and Grafarend (2005),

NumericalValues = {S12 -> 785.619, S23 -> 755.868, S31 -> 1082.873, 
          cosφ12 -> 0.8456221`, 
          cosφ23 -> 0.7108359, 
          cosφ31 -> 0.9752445};

The solutions are

solx1 = 
      NSolve[(drx1 /. r1 /. NumericalValues) == 0, 
        x1]

Let us select the positive ones

s1 = Select[solx1, (First[#][[2]] > 0) &]

for the Sturmfels resultant

solx1 = 
      NSolve[(ℛℯ𝓈\[ScriptTwo]\[ScriptTwo]\[ScriptTwo] /. r1 /. NumericalValues) == 0, 
        x1]

Let us select the positive ones

s1 = Select[solx1, (First[#][[2]] > 0) &]

The result is the same!

We can compute the other two variables similarly.

drx2 = 
      DixonResultant[{f1[[1]], f2[[1]], f3[[1]]}, {x1, x3}, {X1, 
            X3}] // Simplify

Exponent[drx2, x2]

solx2 = 
      NSolve[(drx2 /. r1 /. NumericalValues) == 0, 
        x2]

s2 = Select[solx2, (First[#][[2]] > 0) &]

drx3 = 
      DixonResultant[{f1[[1]], f2[[1]], f3[[1]]}, {x1, x2}, {X1, 
            X2}] // Simplify

Exponent[drx3, x3]

solx3 = 
      NSolve[(drx3 /. r1 /. NumericalValues) == 0, 
        x3]

s3 = Select[solx3, (First[#][[2]] > 0) &]

Now, in order to find the unique solution, every possible combination should be considered. However the lists of the different variables can be different, therefore

ni = Length[s1];
nj = Length[s2];
nk = Length[s3];

Then the combinations

hu = 
      Partition[
        Table[{s1[[i]], s2[[j]], s3[[k]]}, {i, 1, ni}, {j, 1, 
              nj}, {k, 1, nk}] // Flatten, 3]

All of them should be substituted into the observation equations

tr = 
      Map[Abs[{f1[[1]], f2[[1]], f3[[1]]} /. r1 /. 
                NumericalValues /. #] &, hu]

The proper solution is the very triplet, which gives practically zero for all equations. Therefore the sum of the remainders of  every equations should be also small,

yu = Map[Apply[Plus, #] &, tr]

We take solution belonging to the triplet, which gives the smallest sum of the remainders,

hu[[First[Position[yu, Min[yu]]]]]

Example 3. 7-parameter Datum Transformation

The (Subscript[C, 7])(3) Helmert transformation uses three translation parameters ((Subscript[X, 0]),  (Subscript[Y, 0]), (Subscript[Z, 0])), three rotation parameters (α, β, γ) in R rotation matrix) and a scale parameter  (s) to transform one set of coordinates in a given system into another coordinate system.

(({{(Subscript[x, i])}, {(Subscript[y, i])}, {(Subscript[z, i])}}))=s R (({{(Subscript[X, i])}, {(Subscript[Y, i])}, {(Subscript[Z, i])}}))+(({{(Subscript[X, 0])}, {(Subscript[Y, 0])}, {(Subscript[Z, 0])}}))

The rotation matrix can be exprsessed with the skew-symmetric matrix,

S=(({{0, -c, b}, {c, 0, -a}, {-b, a, 0}}))

Then the rotation matrix is

R = ((((Subscript[I, 3]) - S))^(-1)) ((Subscript[I, 3]) + S)

where (Subscript[I, 3]) is a 3×3  identity matrix.

S=(({{0, -c, b}, {c, 0, -a}, {-b, a, 0}}));(Subscript[I, 3]) = IdentityMatrix[3];

R = Inverse[((Subscript[I, 3]) - S)]  . ((Subscript[I, 3]) + S) // Simplify;  
    MatrixForm[R]

for which the next restriction is true:

R (R^★)  =  (Subscript[I, 3])

R . Transpose[R] // Simplify // MatrixForm

For the determination of the 7 parameters of the transformation (a, b, c, (Subscript[X, 0]), (Subscript[Y, 0]), (Subscript[Z, 0]), s) we need minimally (at least) 3 points with known coordinates in both coordinate systems. 
Expressing the rotation matrix with the skew-symmetric matrix, the nonlinear  system to be solved for ( (Subscript[f, i]) =0 ) determining the 9 parameters is:

For i = 1 point

(({{(Subscript[f, 1])}, {(Subscript[f, 2])}, {(Subscript[f, 3])}}))=(((Subscript[I, 3]) - S).(({{(Subscript[x, i])}, {(Subscript[y, i])}, {(Subscript[z, i])}}))-s ((Subscript[I, 3]) + S).(({{(Subscript[X, i])}, {(Subscript[Y, i])}, {(Subscript[Z, i])}}))-((Subscript[I, 3]) - S).(({{(Subscript[X, 0])}, {(Subscript[Y, 0])}, {(Subscript[Z, 0])}})))/.i -> 1//Expand;MatrixForm[(({{(Subscript[f, 1])}, {(Subscript[f, 2])}, {(Subscript[f, 3])}}))]

For  i = 2 point

(({{(Subscript[f, 4])}, {(Subscript[f, 5])}, {(Subscript[f, 6])}}))=(((Subscript[I, 3]) - S).(({{(Subscript[x, i])}, {(Subscript[y, i])}, {(Subscript[z, i])}}))-s ((Subscript[I, 3]) + S).(({{(Subscript[X, i])}, {(Subscript[Y, i])}, {(Subscript[Z, i])}}))-((Subscript[I, 3]) - S).(({{(Subscript[X, 0])}, {(Subscript[Y, 0])}, {(Subscript[Z, 0])}})))/.i -> 2//Expand;MatrixForm[(({{(Subscript[f, 4])}, {(Subscript[f, 5])}, {(Subscript[f, 6])}}))]

For  i = 3 point

(({{(Subscript[f, 7])}, {(Subscript[f, 8])}, {(Subscript[f, 9])}}))=(((Subscript[I, 3]) - S).(({{(Subscript[x, i])}, {(Subscript[y, i])}, {(Subscript[z, i])}}))-s ((Subscript[I, 3]) + S).(({{(Subscript[X, i])}, {(Subscript[Y, i])}, {(Subscript[Z, i])}}))-((Subscript[I, 3]) - S).(({{(Subscript[X, 0])}, {(Subscript[Y, 0])}, {(Subscript[Z, 0])}})))/.i -> 3//Expand;MatrixForm[(({{(Subscript[f, 7])}, {(Subscript[f, 8])}, {(Subscript[f, 9])}}))]

or

sys = Table[(Subscript[f, i]), {i, 1, 9}]

In equation form

eqs9 = Map[# == 0 &, sys]

Let us consider the numerical values of  3 Hungarian points in the system of ETRS89 (\(TraditionalForm\`\(\(x\_1\)\(,\)\(\ \)\(y\_1\)\(,\)\(\ \)\(z\_1\)\(,\)\(\ \)\(...\)\(\ \)\(,\)\(\ \)\(z\_3\)\(\ \)\)\)) and in the local Hungarian system HD72 (Hungarian Datum 1972) (\(TraditionalForm\`\(X\_! \), \ Y\_1, \ Z\_1, \  ... , \ Z\_3\)).

numericalValues = {(Subscript[x, 1]) -> 4171409.677, 
          (Subscript[x, 2]) -> 4146957.889, (Subscript[x, 3]) -> 3955632.880, 
          (Subscript[y, 1]) -> 1470823.777, (Subscript[y, 2]) -> 1277033.850, 
          (Subscript[y, 3]) -> 1611863.197, 
(Subscript[z, 1]) -> 
            4580140.907, (Subscript[z, 2]) -> 4659439.264, (Subscript[z, 3]) -> 4720991.316, 
          (Subscript[X, 1]) -> 4171352.311, (Subscript[X, 2]) -> 4146901.301, 
          (Subscript[X, 3]) -> 3955575.649, 
(Subscript[Y, 1]) -> 
            1470893.887, (Subscript[Y, 2]) -> 1277104.509, (Subscript[Y, 3]) -> 1611933.124, 
          (Subscript[Z, 1]) -> 4580150.178, (Subscript[Z, 2]) -> 4659448.287, 
          (Subscript[Z, 3]) ->  4721000.952};

eq = Table[{(Subscript["f", i]), sys[[i]]}, {i, 1, 9}]; 
    TableForm[eq]

We can reduce the equation system by subtracting the equations. In this way the translation parameters can be eliminated.
The elimination can be done according to this Table below, see the equation system.

Table. 1

({{(Subscript[X, 0]), (Subscript[Y, 0]), (Subscript[Z, 0])}, {(Subscript[f, 1]) - (Subscript[f, 4]), (Subscript[f, 2]) - (Subscript[f, 5]), (Subscript[f, 3]) - (Subscript[f, 6])}, {(Subscript[f, 1]) - (Subscript[f, 7]), (Subscript[f, 2]) - (Subscript[f, 8]), (Subscript[f, 3]) - (Subscript[f, 9])}, {(Subscript[f, 4]) - (Subscript[f, 7]), (Subscript[f, 5]) - (Subscript[f, 8]), (Subscript[f, 6]) - (Subscript[f, 9])}})

The number of the possible sequences is

Binomial[9, 7]

The possible choices are the following sequences:

seq = 
      Partition[Map[# &, Flatten[Subsets[Range[9], {7}]]], 7]

For example, let us consider the 31st sequence

neq = seq[[31]]

One can select  the following eliminations, see the Table 1. below,

Table 1. Possible equation-pairs

({{(Subscript[X, 0]), (Subscript[Y, 0]), (Subscript[Z, 0])}, {(Subscript[f, 1]) - (Subscript[f, 4]), (Subscript[f, 2]) - (Subscript[f, 5]), (Subscript[f, 3]) - (Subscript[f, 6])}, {(Subscript[f, 1]) - (Subscript[f, 7]), (Subscript[f, 2]) - (Subscript[f, 8]), (Subscript[f, 3]) - (Subscript[f, 9])}, {(Subscript[f, 4]) - (Subscript[f, 7]), (Subscript[f, 5]) - (Subscript[f, 8]), (Subscript[f, 6]) - (Subscript[f, 9])}})

After selecting a sequence, see above, we should take those eliminations from the table, which indices cover the sequence, namely ( 2, 3, 4, 5, 6, 8, 9), see eliminations with green backround in the table

r1 = (Subscript[f, 2]) - (Subscript[f, 8]);

r2 = (Subscript[f, 2]) - (Subscript[f, 5]);

r3 = (Subscript[f, 3]) - (Subscript[f, 9]);

r4 = (Subscript[f, 6]) - (Subscript[f, 9]);

In the remaining 4 equations will remain only 4 unknown parameters (a, b, c,  s). The second order nonlinarity is caused by the parameter s. Therefore if the value of s is known then the system will be linear for a, b, c.

sysR = {r1, r2, r3, r4}

Now we can solve the original system, eliminating the variables  a,  b and c to get a polynomial for s

Timing[drs = 
          DixonResultant[{r1, r2, r3, r4}, {a, b, 
                c}, {(Subscript[σ, 1]), (Subscript[σ, 2]), (Subscript[σ, 3])}] // 
            Simplify;]

drs

This is a quadratic equation for s

Exponent[drs, {a, b, c, s}]

The coefficients of these equations are

h1s = Coefficient[drs, s];

h2s = Coefficient[drs, (s^2)];

h3s = Coefficient[drs, (s^3)];

h4s = Coefficient[drs, (s^4)];

and

h0s = 
      Simplify[drs - (h1s s + h2s (s^2) + h3s (s^3) + 
              h4s (s^4))]

Let us compare the result of this solution with that of the direct numerical one. The quadratic polynomial for s,

qps = (h0s + s h1s + (s^2) h2s + (s^3) h3s + (s^4) h4s) /. 
        numericalValues

sols = NSolve[qps == 0, s]

The solution must be positive real

s0 = 
      Select[s /. sols, Im[#] == 0 ∧ # > 0 &] // First

In case of the value of s is known the solution of the system sysR is easy because its linearity. We have more equations than uknown, therefore pseudo inverse solution is proper.

<< LinearAlgebra`MatrixManipulation`

Ab = 
        LinearEquationsToMatrices[{r1, r2, r3, r4} /. 
            Join[{s -> s0}, numericalValues], {a, b, c}];

AA = Ab[[1]]; bb = Ab[[2]];

{a0, b0, c0} = PseudoInverse[AA] . bb // Simplify

In order to get the result in proper form, it is necessary transform these parameters into Cardan angles

Parametrization of the rotation matrix by Cardan angles

R (α, β, γ)=(({{cos β cos γ, \(cos\ β\ sin\ γ\), -sinβ}, {sin α sin β cos γ - 
                cos α sin γ, \(sin\ α\ sin\ β\ sin\ γ + 
                cos\ α\ cos\ γ\), \(sin\ αcos\ β\)}, {cos α sin β cos γ + 
                sin α sin γ

              , cos α sin β sin γ - 
                sin α cos γ, \(cos\ α\ cos\ β\)}}))

The cardan angles can be obtained from rotation matrix through:

\(tan\ α\  = \(r\_23\/r\_33\  => \ Λ\_Γ = arctan\ r\_23\/r\_33\)\)

\(tan\ β\  = \(\(-r\_31\)\/\@\(r\_11\^2 + r\_12\^2\) => ϕ\_Γ = arctan\ \(-r\_31\)\/\@\(r\_11\^2 + r\_12\^2\)\), \ 
    or\ β = \(-arcsin\)\ r\_13\)

\(tan\ γ = \(r\_12\/r\_11 => Ξ\_Γ = arctan\ r\_12\/r\_11\)\)

Consequently

R0 = R /. {a -> a0, b -> b0, c -> c0};

R0 // MatrixForm

α0 = 
        ArcTan[R0[[3, 3]], R0[[2, 3]]]*(180/π)*3600;

β0 = -ArcSin[R0[[1, 3]]]*(180/π)*3600;

γ0 = 
        ArcTan[R0[[1, 1]], R0[[1, 2]]]*(180/π)*3600;

{α0,  β0,  γ0}

To determine the translation parameters we can use the least squares method, because we have 7 linear equations with 3 variables. The equivalent 7 equations are:

neq

Coosing these:

sys7 = 
      Map[(sys[[#]] == 0) &, neq] /. 
        Join[numericalValues, {a -> a0, b -> b0, c -> c0, 
            s -> s0}]

Ab = 
        LinearEquationsToMatrices[sys7, {(Subscript[X, 0]), (Subscript[Y, 0]), (Subscript[Z, 0])}];

AA = Ab[[1]]; bb = Ab[[2]];

{X00, Y00, Z00} = PseudoInverse[AA] . bb

Now let us try to solve the problem with Groebner Basis. The built-   in functionGroebnerBasis can be also used to eliminate directly prespecified variables form the basis, Trott (2005). In this case, one may employ this function in similar way than the Dixon method does. It goes without saying,that now the function should be called three times with different arguments.

\(\(\(\(sysR // FullSimplify\) // TableForm\)\(
\)
    \)\)

For solving this system of equation with Groebner, it may helpful to introduce new variables. It is clear that the introduction of the following variables can simplify the description of the system, and consequently may reduces running time, see later.

newVarsB = {(Subscript[x, 2]) -> (Subscript[x, 1]) - x12, (Subscript[x, 3]) -> (Subscript[x, 1]) - x13, 
          (Subscript[x, 3]) -> (Subscript[x, 2]) - x23, (Subscript[y, 2]) -> (Subscript[y, 1]) - y12, 
          (Subscript[y, 3]) -> (Subscript[y, 1]) - y13, 
          (Subscript[y, 3]) -> (Subscript[y, 2]) - y23, 
(Subscript[z, 2]) -> 
            (Subscript[z, 1]) - z12, (Subscript[z, 3]) -> (Subscript[z, 1]) - z13, 
          (Subscript[z, 3]) -> (Subscript[z, 2]) - z23, 
(Subscript[X, 2]) -> 
            (Subscript[X, 1]) - X12, (Subscript[X, 3]) -> (Subscript[X, 1]) - X13, 
          (Subscript[X, 3]) -> (Subscript[X, 2]) - X23, 
(Subscript[Y, 2]) -> 
            (Subscript[Y, 1]) - Y12, (Subscript[Y, 3]) -> (Subscript[Y, 1]) - Y13, 
          (Subscript[Y, 3]) -> (Subscript[Y, 2]) - Y23, 
(Subscript[Z, 2]) -> 
            (Subscript[Z, 1]) - Z12, (Subscript[Z, 3]) -> (Subscript[Z, 1]) - Z13, 
          (Subscript[Z, 3]) -> (Subscript[Z, 2]) - Z23};

then our equation system is

sysRG = sysR /. newVarsB // FullSimplify // Expand

Employing reduced Groebner basis for eliminating a, b and c

Timing[{gbs} = 
          GroebnerBasis[sysRG, {s, a, b, c}, {a, b, c}];]

It means,  in that case the Reduced Groebner algorthm  is very slow even with these new variables. However, using different variable order we get the result quickly

Timing[{gbs} = 
          GroebnerBasis[sysRG, {s, b, c, a}, {b, c, a}];]

gbs

On the other hand employing the original system without introducing new variable, the running time is considerably higher even employing this new variable order

Timing[{gbs} = 
          GroebnerBasis[{r1, r2, r3, r4}, {s, b, c, a}, {b, c, 
              a}];]

We can solve with Dixon resultant also this simplified system of equation.

Timing[drs=DixonResultant[sysRG,{a, b, c},{(Subscript[σ, 1]), (Subscript[σ, 2]), (Subscript[σ, 3])}];]

drs

gbs + drs

Using Dixon or Grobner for solving the reduced system the solution will be the same.

The running time for reduced Gröbner basis algorithm is definitely changes with the order of variables to be eliminated, see Table I for a few chosen sequences. In case of six sequences from the 36 possible choices, these changes are enourmous! One of them was choosen as illustration above, see also the last column in Table I.

Table I.

Running times (sec) in case of different sequences and order of the variables to be eliminated
using reduced Groebner basis

({{., 1, 2, 3, 4, 6, 7, 8, 1, 2, 3, 4, 5, 8, 9, 1, 2, 3, 4, 
          5, 7, 9, 2, 3, 4, 5, 6, 8, 9}, {a, b, c, 0.21899999999999986`, 0.6879999999999882`, 0.9849999999999568`, 3017.69}, {a, c, b, 0.3600000000000003`, 35.92099999999999`, 0.6720000000000255`, 2601.75` }, {b, a, c, 0.25, 0.9220000000000255` , 0.7649999999999864`, 3172.4849999999997`}, {b, c, a, 163.76500000000001` , 0.5470000000000255`, 33.985000000000014`, 0.484}, {c, a, b, 0.7189999999999941`, 52.327999999999975`, 0.561999999999955` , 1831.42}, {c, b, a, 174.547`, 0.75` , 47.79699999999997`,  0.532}})

Groebner basis works very slowly or not at all with absolute coordinates (after eliminating the translation parameters) and hence need introducing new, relative coordinates. In Table II. we ilustrated the computation times of the above descibed example with and without using the new, relative variables with two different order of variables.

Table II.

Running times (sec) in case of one of the six critical sequences from the 36, using absolute coordinates and relative ones

({{Sequence, 2, 3, 4, 5, 6, 8, 9, 2, 3, 4, 5, 6, 8, 9}, {Order of variables, Original System, System with new  variables, 
 
            with relative coordinates}, {a, b, c, -, 3017.69}, {b, c, a, 1320.09, 0.484}})

The Dixon Resultant is faster than grobner basis, and very robust concerning the order of variables. This feature can be very important from practical point of view, because in case of Groebner basis the user should find the proper order via trial-error method, if exists any! For example in the above described case (2,3,4,5,6,8,9) only two orders can provide solution in a rational time from the six possible orders of the variables. An other fact that shows the robustness of Dixon that it works also without introducing new variables with  relative coordinates.

Conclusions

Dixon resultant seems very attractive tool for solving system of multivariate polynomial equations in geodesy. 
Comparing it to other Multipolynomial Resultant like Strumfels method, its advantage, the small size and density of the Dixon matrix and consequently the fast evaluation in case of relatively small systems (system with 3-4 variables). Comparing with Groebner basis the Dixon resultant is more robust  and gives a faster solution also. It can solve problems easier when Grobner basis is sensitive to the order of variables, and requires more work to find out the good order if exists any.  However we should mention also that there are problems also when Dixon resultant does not give a solution. For example in the 3D resection problem the determinant of the Dixon matrix is identically zero, therefore we get no information on how to compute the common roots. 
From  practical point of view, Dixon resultant can be successfully applied in many areas of geodesy, where symbolic, or symbolic - numeric calculation are necessary, i.e. intersection, ranging, datum transformation and so on. It can be a good alternative of the Groebner basis, or Strumfels method which happens to be slower than its counterpart, the Dixon resultant in some cases.

Acknowledgements

The third author wish to thank to the Hungarian Eotvos Fellowship for  the support of her visit at the Department of Geodesy and Geoinformatics of  Stuttgart University, in Germany,  during this time-period this work has been accomplished!

References

[1]Haneberg, W.C.:Computational Geosciences with Mathematics, Springer , Berlin, (2004)

[2] Awange J.L. and Grafarend E.W.: Solving Algebraic Computational Problems in Geodesy and Geoinformatics, Springer, Berlin, (2005)

[3] Freeman, J.A.: Simulating neural networks with Mathematica, Addison-Wesley, (1994)

[4] Bart M. ter Haar Romeny: Front-End Vision and Multi-Scale Image Analysis, Kluwer, Dordrecht,  (2003)

[5] Helton, J.W and Merino O.: Classical Control using a Methods, SIAM,  (1998)

[6] Bellomo, N., Preziosi, L. and Romano, A.: Mechanics and Dynamical Systems with Mathematica, Birkhäuser, (2000)

[7] Bancroft, S.: An algebraic solution of the GPS equations, IEEE Transaction on Aerospace and Electronic Systems AES-21, (1985), pp.
56-59.

[8] Grafarend EW. and Schaffrin, B.: Ausglechungsrechnung in Linearen Modellen, B.I. Wissenschaftsverlag, Mannheim,  (1993)

[9] Singer P., Ströbel D., Hördt R., Bahndorf J., Linkwitz K.: Direkte Lösung des räumlichen Bogenschnitts. Zeitschrift für Vermessungswesen, 124, (1993), S. 295 - 297.

[10] Kleusberg A.: Die direkte Lösung des räumlichen Hyperbelschnitts. Zeitschrift für Vermessungswesen, 119,  (1994), S. 188 - 192.

[11] Kleusberg  A.: Analytical GPS navigation solution. In Grafarend EW., Krumm FW. and Schwarze VS. (eds) Geodesy -the Challenge of the
3rd Millennium, Springer, Heidelberg, (2003), pp. 93 - 96.

[12] Lichtenegger H.: Eine direkte Lösung des räumlichen Bogenschnitts., Österreichische Zeitschrift für Vermessung und Geoinformation 83,  (1995), S. 224 - 226.

[13] Awange J.L.: Gröbner bases, multipolynomial resultants and Gauss-Jacobi combinatorial algorithms - adjustment of nonlinear GPS/LPS observations, Dissertation (D93), Geodätisches Institute der Universität Stuttgart, (2002).

[14] Grafarend EW. and Shan J.: Closed - form solution of the nonlinear  pseudoranging equations (GPS), Artificial Satellites, Planetary Geodesy 31,  (1996), pp.133 - 147.

[15] Sturmfels B.: Multigraded resultant of  Sylvester type. Journal of Algebra 163, (1994), pp.115-127.

[16] Dixon A.L.: The eliminant of three quantics in two independent variables. Proc. London Math.Soc. 6,  (1908),  pp. 468-478.

[17] Kapur D, Saxena T. and Yang L. (1994): Algebraic and geometric reasoning using Dixon resultants. In: ACM ISSAC 94. Oxford, England, July 1994, pp. 99-107.

[18] Bernstein, D.S.: Matrix Mathematics, Theory, Facts, and Formulas with Application to Linear Systems Theory, Princeton University Press, Princeton and Oxford, (2005), p.142.

[19] Nakos G. and Williams R.: Elimination with The Dixon resultant. Mathematica in Education and Research, Vol.6. No.2,  (1997), pp. 11- 21.

[20] Nakos G. and Williams R.: A fast algorithm implemented in Mathematica provides one-step elimination of a block of unknowns from a system of polynomial equations, (2002), http://wolfram.com/infocenter/Articles/2597

[21] Buchberger B. Gröbner bases: An Algorithmic method in Polynomial Ideal theory, Multidimensional Systems Theory, N.K. Bose, ed., D.Reidel Publ.Co. (1985)

[22] Trott M.: Mathematica Guide Book for Symbolic Computation, only electronic version avaiable yet, in Mathematica Guide Book for Programming or Graphics, Springer, New York,  (2005)

[23] Jacobi CGI.: Deformatione et proprietatibus determinantum, Crelle's Journal für die reine und angewandte Mathematik, Bd. 22, (1841)

[24]  Wellisch S.: Theorie und Praxis der Ausgleichungsrechnung. Bd.2, (1910)

[25] Hornoch AT: Über die Zurückführung der Methode der kleinsten Quadrate  auf das Prinzipe des arithmetischen Mittels. Zeitschrift für Vermessungswesen 38, .(1950), S. 13 - 18.

[26] Zavoti J, Jancso T.: The solution of the 7-parameter datum transformation problem with- and without the Gröbner basis, Acta Geod. Geoph. Hung., Vol. 41(1), pp. 87-100 (2006)

[27] Hofman-Wellenhof B, Lichtenegger H, Collins J.: Global Positioning System: Theory and practice, 5th Ed. Springer, Wien