# Distinguishing between the solutions of a system of equations

Let's define a sample list:

Clear["Global*"];
grd = {{0.2, 0.2}, {0.4, 0.8}, {0.8, 0.1}, {0.8, 0.46}, {1, 1}};


and then a couple of trivial equations:

r1 = Sqrt[(x + μ)^2 + y^2];
r2 = Sqrt[(x + μ - 1)^2 + y^2];
Ω = (q1*(1 - μ))/r1 + (q2*μ)/r2 + ω^2/2*(x^2 + y^2);
Ωx = D[Ω, x];
Ωy = D[Ω, y];

μ = 0.2;
ω = 0.25;


Then a module for finding roots of system of equations:

Options[FindRoots2D] = {WorkingPrecision -> 20, MaxRecursion -> 20};

FindRoots2D[funcs_, {x_, a_, b_}, {y_, c_, d_}, opts___] :=
Module[{fZero, seeds, signs, fy},
fy = Compile[{x, y}, Evaluate[funcs[]]];
fZero =
Cases[Normal[
ContourPlot[
funcs[] == 0, {x, a - (b - a)/97, b + (b - a)/103}, {y,
c - (d - c)/98, d + (d - c)/102},
Evaluate[FilterRules[{opts}, Options[ContourPlot]]]]],
Line[z_] :> z, Infinity];
seeds = Flatten[((signs = Sign[Apply[fy, #1, {1}]];
#1[[1 +
Flatten[
Position[Rest[signs*RotateRight[signs]], -1]]]]) &) /@
fZero, 1];
If[seeds == {}, {},
Select[Union[({x, y} /.
FindRoot[{funcs[],
funcs[]}, {x, #1[]}, {y, #1[]},
Evaluate[FilterRules[{opts}, Options[FindRoot]]]] &) /@
seeds, SameTest -> (Norm[#1 - #2] < 10^(-6) &)],
a <= #1[] <= b && c <= #1[] <= d &]]]


We use this module in order to obtain a new list data:

data = (
Q1 = #[];
Q2 = #[];
q1 = Q1^3;
q2 = Q2^3;
pts = FindRoots2D[{Ωx, Ωy}, {x, -3, 3}, {y, -3, 3}, PlotPoints -> 200];
Npts = Length[pts];
{Q1, Q2, Npts}) & /@ grd;


The system of equations has either 3 or 5 solutions. For example when $Q1 = Q2 = 0.2$ we have

{{-0.613042, 3.72089*10^-20}, {0.3, 0.0631203}, {0.3, -0.0631203},
{0.302369, 8.6956*10^-17}, {0.969193, 4.35961*10^-19}}


5 solutions, 3 on the $x$ axis and 2 on the $(x,y)$ plane.

I want the following:

Automatically distinguish between the roots. In particular:

• the root with the smallest $x$ value ---> let's call it $xL2$
• the root with highest $x$ value ---> let's call it $xL3$
• the root between the above ones ---> let's call it $xL1$
• the root with non zero $y$ value ---> this root is $(x,y) = (xL4,yL4)$

So the data list should contain the following information: {Q1, Q2, Npts, xL1, xL2, xL3, xL4, yL4}. In the case where there are only 3 roots the list should contain only {Q1, Q2, Npts, xL1, xL2, xL3}, without Null in the position of $xL4$ and $yL4$.

data = (Q1 = #[];
Q2 = #[];
q1 = Q1^3;
q2 = Q2^3;
pts =
Chop@FindRoots2D[{Ωx, Ωy}, {x, -3,
3}, {y, -3, 3}, PlotPoints -> 200];
Npts = Length[pts];

xL1 =
If[Npts == 5,
First@Select[pts,
Min@pts[[All, 1]] < #[] < Max@pts[[All, 1]] &&
Min@pts[[All, 2]] < #[] < Max@pts[[All, 2]] &],
Select[pts,
Min@pts[[All, 1]] < #[] < Max@pts[[All, 1]] &]];
xL2 = First@Select[pts, #[] == Min@pts[[All, 1]] &];
xL3 = First@Select[pts, #[] == Max@pts[[All, 1]] &];
xL4 =
If[Npts == 5, First@Select[pts, #[] == Max@pts[[All, 2]] &],
Nothing];
xL5 =
If[Npts == 5, First@Select[pts, #[] == Min@pts[[All, 2]] &],
Nothing];

{Q1, Q2, Npts, xL1, xL2, xL3, xL4, xL5}) & /@ grd;


Then

data A pre-Nothing version of Nothing may be Unevaluated@Sequence[].

You might do better with an approach that uses NSolve to find all solutions, given a specific set of values for the parameters {q1,q2}. I show one way to do this below.

The idea is to first make the system into explicit polynomials by clearing denominators and substituting new variables and polynomials to define the radicals. Also since they appear in denominators, we add new variables and defining relations that force them not to vanish (standard tactic: use var*varrecip-1). I remark that the use of relations defining the radicals can introduce "parasite" solutions, so we'll need to check for that at the end.

grd = Rationalize[{{0.2, 0.2}, {0.4, 0.8}, {0.8, 0.1}, {0.8,
0.46}, {1, 1}}];
r1 = Sqrt[(x + \[Mu])^2 + y^2];
r2 = Sqrt[(x + \[Mu] - 1)^2 + y^2];
\[Mu] = 1/5;
\[Omega] = 1/4;
\[CapitalOmega] = (q1*(1 - \[Mu]))/r1 + (q2*\[Mu])/
r2 + \[Omega]^2/2*(x^2 + y^2);
\[CapitalOmega]x = Numerator[Together[D[\[CapitalOmega], x]]];
\[CapitalOmega]y = Numerator[Together[D[\[CapitalOmega], y]]];

reprules = {Sqrt[16 - 40 x + 25 x^2 + 25 y^2] -> sqrt1,
Sqrt[1 + 10 x + 25 x^2 + 25 y^2] -> sqrt2};
polys = {\[CapitalOmega]x, \[CapitalOmega]y,
16 - 40 x + 25 x^2 + 25 y^2 - sqrt1^2,
1 + 10 x + 25 x^2 + 25 y^2 - sqrt2^2, sqrt1*sqrt1recip - 1,
sqrt2*sqrt2recip} /. reprules

(* Out= {-5120 q1 sqrt1 + 320 q2 sqrt2 - 12800 q1 sqrt1 x +
2800 q2 sqrt2 x + 16 sqrt1 sqrt2 x + 56000 q1 sqrt1 x^2 +
4000 q2 sqrt2 x^2 + 120 sqrt1 sqrt2 x^2 - 40000 q1 sqrt1 x^3 -
10000 q2 sqrt2 x^3 + 25 sqrt1 sqrt2 x^3 - 750 sqrt1 sqrt2 x^4 +
625 sqrt1 sqrt2 x^5 - 8000 q1 sqrt1 y^2 + 8000 q2 sqrt2 y^2 -
40000 q1 sqrt1 x y^2 - 10000 q2 sqrt2 x y^2 +
425 sqrt1 sqrt2 x y^2 - 750 sqrt1 sqrt2 x^2 y^2 +
1250 sqrt1 sqrt2 x^3 y^2 + 625 sqrt1 sqrt2 x y^4,
y (-25600 q1 sqrt1 - 400 q2 sqrt2 + 16 sqrt1 sqrt2 +
64000 q1 sqrt1 x - 4000 q2 sqrt2 x + 120 sqrt1 sqrt2 x -
40000 q1 sqrt1 x^2 - 10000 q2 sqrt2 x^2 + 25 sqrt1 sqrt2 x^2 -
750 sqrt1 sqrt2 x^3 + 625 sqrt1 sqrt2 x^4 - 40000 q1 sqrt1 y^2 -
10000 q2 sqrt2 y^2 + 425 sqrt1 sqrt2 y^2 -
750 sqrt1 sqrt2 x y^2 + 1250 sqrt1 sqrt2 x^2 y^2 +
625 sqrt1 sqrt2 y^4), 16 - sqrt1^2 - 40 x + 25 x^2 + 25 y^2,
1 - sqrt2^2 + 10 x + 25 x^2 + 25 y^2, -1 + sqrt1 sqrt1recip,
sqrt2 sqrt2recip} *)


For a given set of parameter values, numerically solve the system. I first eliminate the new radical and reciprocal variables, mostly to get rid of clutter.

AbsoluteTiming[
solns = NSolve[
GroebnerBasis[
polys /. Thread[{q1, q2} -> grd[]^3], {x, y}, {sqrt1, sqrt2,
sqrt1recip, sqrt2recip}, MonomialOrder -> EliminationOrder],
Reals]]

(* Out= {0.738787, {{x -> 0.969193, y -> 0.}, {x -> 0.957362,
y -> 0.}, {x -> 0.617768, y -> 0.}, {x -> -0.613042,
y -> 0.}, {x -> -0.606563, y -> 0.}, {x -> 0.514257,
y -> 0.}, {x -> 0.45252, y -> 0.}, {x -> 0.3,
y -> -0.0631203}, {x -> 0.302369, y -> 0.}, {x -> 0.3,
y -> 0.0631203}, {x -> -0.2, y -> 0.}}} *)


Now discard parasites. We accept as valid any solution that gives a residual of under 10^(-6). There are ore reliable ways to go about this but you get the idea.

validsolns =
Select[solns, (
Norm[{D[\[CapitalOmega], x], D[\[CapitalOmega], y]} /.
Thread[{q1, q2} -> grd[]^3] /. #] < 10^(-6)) &]

(* Out= {{x -> 0.969193, y -> 0.}, {x -> -0.613042, y -> 0.},
{x -> 0.3, y -> -0.0631203}, {x -> 0.302369, y -> 0.},
{x -> 0.3, y -> 0.0631203}} *)


This approach is likely to be more reliable for the stated purpose than using FindRoot repeatedly and hoping all valid solutions get found.

• I do not see how this can obtain the desired data` list with the information {Q1, Q2, Npts, xL1, xL2, xL3, xL4, yL4}. Oct 16 '16 at 17:38