How is the starting point chosen in GNU Scientific Library multiroot finder?

I'm using an implementation of the GNU Scientific Library's multiroot finder to solve for the unknowns (x and y) in the following system of non-linear equations:

I'm a bit confused, however, about the "starting point":

Solve(const double *x, int maxIter = 0, double absTol = 0, double relTol = 0) Find the root starting from the point X; Use the number of iteration and tolerance if given otherwise use default parameter values which can be defined by the static method SetDefault

How is the starting point chosen?

In one dimensional polynomial problem, this is a well-studied domain of Real-root_isolation. In two dimensions its less well know.

First note we can transform the equations into polynomials, take the first one

sqrt[(x-x1)^2 + (y-y1)^2] + s (t2 -t1) = sqrt[(x-x1)^2 + (y-y1)^2]

square both sides

[(x-x1)^2 + (y-y1)^2] + 2 sqrt[(x-x1)^2 + (y-y1)^2][s (t2 -t1)] + [s (t2 -t1)]^2 = [(x-x1)^2 + (y-y1)^2]

rearrange

2 sqrt[(x-x1)^2 + (y-y1)^2][s (t2 -t1)] = [(x-x1)^2 + (y-y1)^2] - [(x-x1)^2 + (y-y1)^2] - [s (t2 -t1)]^2

square again

4 [(x-x1)^2 + (y-y1)^2] [s (t2 -t1)]^2 = ([(x-x1)^2 + (y-y1)^2] - [(x-x1)^2 + (y-y1)^2] - [s (t2 -t1)]^2)^2

giving us a polynomial.(note)

Once we have a set of polynomial there are some algebraic techniques you can use.

A technique I use is the convert the polynomials into Bernstein polynomials. First, rescale your domain so both variables lie in [0,1]. If m, n are the degrees of the two polynomials then a polynomial can be expressed as the sum

sum i=0..m sum j=0..n b_ij mCi nCj x^i (1-x)^(n-i) y^j (1-y)^(n-j)

Where the mCi, nCj are binomial coefficients.

These polynomials have a nice convexity property. If the coefficients b_ij are all positive then the value of the polynomial will be positive for all x,y in [0,1]. Similar if the coefficients are all negative.

This allows you to say that "there is no solution in this region".

So a strategy to solve the problem might be:

1. pick the largest region the solutions might occur in.
2. Subdivide this into a set of smaller regions
3. For each region calculate the Bernstein polynomials for each of your equations
4. Examine the coefficients of the Bernstein polynomials, if they are all positive (or all negative) reject that region
5. You should now have rejected a large part of the domain. Start the multiroot finder using a point in each remaining region.

If you want details of how to construct Berstein polynomials you can read my paper at A new method for drawing Algebraic Surfaces

Note: by squaring we actually get more solutions that we want. In the initial problem we want the principle sqareroot, i.e. +sqrt(A), there are also solutions using the other root -sqrt(A). This makes the problem a bit easier rather than four solutions which we would get from the intersection of two hyperbolae, we just have the intersections of one branch so one or two solutions to the problem.

---

For your problem, there is quite a nice simple way to get one of the starting points.

Assume s=0. Then each equation will give the set of points that are equidistant from two points. This is a line, the perpendicular bisector of the line segment joining the points, then simply find the point of intersection of the three perpendicular bisectors. This will be the centre of the Circumscribed_circle. This might be enough for the algorithm. Even simpler is to simply take the average of the three points.