2

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:

enter image description here

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?

10GeV
  • 453
  • 2
  • 14
  • 1
    This is a problem that many numerical methods have. They work by improving the current estimate at each step. How does one get started? Usually that is outside the method itself. One plausible approach is to try random starting points, another is to solve a simpler version of the problem and then use that as your starting point for the full problem. Many other heuristics are possible. Good luck and have fun. – Robert Dodier Oct 03 '20 at 01:57
  • 1
    The third equation is redundant, since it follows from the first two. Each equation is the locus of points with constant difference of distances to two fixed points, known to be a [hyperbola](https://en.wikipedia.org/wiki/Hyperbola#Definition_of_a_hyperbola_as_locus_of_point). So you are looking for the intersection points of two hyperbolas, which can have up to 4 distinct solutions. You'll need some additional information about the equations to decide which point(s) you want to find, and that can hint to where to start. The problem can also be solved analytically, though it's not pretty. – dxiv Oct 03 '20 at 03:34
  • @dxiv A hyperbola gives the locus where the absolute difference of the distances to two points is a constant. i.e |d1 - d2| = k. Here we want the signed distance. This just gives one branch of the hyperbola. Rather that 4 solutions we expect a single solution to the problem. – Salix alba Oct 03 '20 at 09:11
  • I’m voting to close this question because choosing a starting point for an iterative solver isn't a programming problem. The answer already posted also strongly suggests that this is a mathematical problem. – High Performance Mark Oct 03 '20 at 11:26
  • 1
    @Salixalba You are right about the signedness of the difference here, but that still leaves up to two solutions in general. – dxiv Oct 03 '20 at 18:11
  • @KeithMadison Did you mean `s(t1-t3)` in the third equation? – dxiv Oct 10 '20 at 18:18

1 Answers1

2

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.

Salix alba
  • 7,536
  • 2
  • 32
  • 38