# Good practice when working with multiple solutions of a system of equations

I have a system of nonlinear equations. I want to do some comparative statistics and I want to plot the results from the statistical computations. I have two questions about good practice on manipulating the solution lists.

1. What is the best practice for selecting the results from solving a system of equations without calling Solve multiple times?
2. When more than one solution for the same variable is available, how do I select values from the according to a given criterion; for example, the highest value?

The best example I could come up is:

EQ1[c_] := x - 0.5 - ((0.8*(z - 0.5))/(y - 0.5 - 0.8*(z - 0.5)))*(y - 0.5);
EQ2[c_] := 1 - x - y - z;
EQ3[c_] := 0.8*(z - c)^2 + (y - 0.5)^2 - 0.8*(-c)^2
SYSTEM = {EQ1[c], EQ2[c], EQ3[c]};
Table[
{c, Solve[SYSTEM == 0 && 0 <= x <= 1 && 0 <= y <= 1 && 0 <= z <= 1, {x, y, z}]},
{c, 0.24, 0.28, 0.01}]

{{0.24, {{x -> 0.0380158, y -> 0.48276, z -> 0.479225},
{x -> 0.291508, y -> 0.326764, z -> 0.381728}}},
{0.25, {{x -> 0.309524, y -> 0.309524, z -> 0.380952},
{x -> ComplexInfinity, y -> 0.5, z -> 0.5}}},
{0.26, {{x -> 0.324093, y -> 0.294404, z -> 0.381503}}},
{0.27, {{x -> 0.33648, y -> 0.280628, z -> 0.382892}}},
{0.28, {{x -> 0.347344, y -> 0.267788, z -> 0.384868}}}}


Now that we have the example, I will go back to my questions:

1. What is the best practice for building a table such that I can easily plot the results afterwards? The way I would do, which is really inefficient, is:

Tablex =
Table[
{c, Solve[SYSTEM == 0 && 0 <= x <= 1 && 0 <= y <= 1 && 0 <= z <= 1, {x, y, z}][[1,1,2]]},
{c, 0.24, 0.28, 0.01}]
Tabley =
Table[
{c, Solve[SYSTEM == 0 && 0 <= x <= 1 && 0 <= y <= 1 && 0 <= z <= 1, {x, y, z}][[1,2,2]]},
{c, 0.24, 0.28, 0.01}]
Tablez =
Table[
{c, Solve[SYSTEM == 0 && 0 <= x <= 1 && 0 <= y <= 1 && 0 <= z <= 1, {x, y, z}][[1,3,2]]},
{c, 0.24, 0.28, 0.01}]

2. We can see that for c <= 0.25 we have two possible solutions. How can we create a rule that, for example, always looks for the highest value of the first variable? I can see the results are in ascending order, therefore, I could get the highest value for all c such that c <= 0.25 with

Solve[SYSTEM == 0 && 0 <= x <= 1 && 0 <= y <= 1 && 0 <= z <= 1, {x, y, z}][[2,1,2]]


But if the parameters change, this threshold will change and so I would like my code to pick itself the highest x in the case of multiple solutions.

• How do you want the case x -> ComplexInfinity to be handled? Nov 27 '17 at 2:00
• First: thank you for the edits. We can ignore the ComplexInfinity case, m_goldberg. Nov 27 '17 at 16:32

Let's start with the results you show for your example and see if we can't exact the values you want from that list.

data =
{{0.24, {{x -> 0.0380158, y -> 0.48276, z -> 0.479225},
{x -> 0.291508, y -> 0.326764, z -> 0.381728}}},
{0.25, {{x -> 0.309524, y -> 0.309524, z -> 0.380952},
{x -> ComplexInfinity, y -> 0.5, z -> 0.5}}},
{0.26, {{x -> 0.324093, y -> 0.294404, z -> 0.381503}}},
{0.27, {{x -> 0.33648, y -> 0.280628, z -> 0.382892}}},
{0.28, {{x -> 0.347344, y -> 0.267788, z -> 0.384868}}}}

{cVals, solns} = Transpose[data];
{xVals, yVals, zVals} =
Transpose[({x, y, z} /. MaximalBy[#, #[[1, 2]] &] & /@ solns)[[All, 1]]];

xVals


{0.291508, ComplexInfinity, 0.324093, 0.33648, 0.347344}

 yVals


{0.326764, 0.5, 0.294404, 0.280628, 0.267788}

 zVals


{0.381728, 0.5, 0.381503, 0.382892, 0.384868}

The complex infinity value at c = .25 will cause trouble when plotting x against c. You need to come up with a policy to handle infinite results. But, ignoring that for the moment, here is a plot.

ListLinePlot[{xVals, yVals, zVals},
DataRange -> {.24, .28},
PlotLabels -> {x, y, z}] The main point of this answer is not in the details of the expression manipulations I used; rather, it is that is good practice to solve the system of equations once and manipulate the expression returned.

### Update

BTW, when I set up your equations in the way I think they should be set up, I get very different results from those you present.

EQ1 = x - 0.5 - ((0.8*(z - 0.5))/(y - 0.5 - 0.8*(z - 0.5)))*(y - 0.5) == 0;
EQ2 = 1 - x - y - z == 0;
EQ3[c_] := 0.8*(z - c)^2 + (y - 0.5)^2 - 0.8*c^2 == 0
SYSTEM = EQ1 && EQ2 && EQ3[c];
solns =
Quiet @ Table[
Solve[SYSTEM && 0 <= x <= 1 && 0 <= y <= 1 && 0 <= z <= 1, {x, y, z}] // Chop,
{c, 0.24, 0.28, 0.01}]

{{{x -> 0.0380158, y -> 0.48276, z -> 0.479225},
{x -> 0.291508, y -> 0.326764, z -> 0.381728},
{x -> 0.5, y -> 0.5, z -> 0},
{x -> 0.5, y -> 0.5, z -> 0}},
{{x -> 0.309524, y -> 0.309524, z -> 0.380952},
{x -> 0.5, y -> 0.5, z -> 0.5},
{x -> 0.5, y -> 0.5, z -> 0},
{x -> 0.5, y -> 0.5, z -> 0}},
{{x -> 0.324093, y -> 0.294404, z -> 0.381503},
{x -> 0.5, y -> 0.5, z -> 0},
{x -> 0.5, y -> 0.5, z -> 0}},
{{x -> 0.33648, y -> 0.280628, z -> 0.382892},
{x -> 0.5, y -> 0.5, z -> 0},
{x -> 0.5, y -> 0.5, z -> 0}},
{{x -> 0.347344,y -> 0.267788, z -> 0.384868},
{x -> 0.5, y -> 0.5, z -> 0},
{x -> 0.5, y -> 0.5, z -> 0}}}


### Update 2

To eliminate solutions where (x - 0.5)/(y - 0.5) is not well-defined, I will delete solutions where y -> 0.5 appears before selecting the solutions with maximal x value.

{xVals, yVals, zVals} =
Transpose[
({x, y, z} /. MaximalBy[#, #[[1, 2]] &] & /@
DeleteCases[solns, {_, y -> _?(# == .5 &), _}, {2}])[[All, 1]]];


Then

ListLinePlot[{xVals, yVals, zVals},
DataRange -> {.24, .28},
PlotLabels -> {x, y, z}] • Thank you very much for the answer. I have been learning a lot on this community. Your answers in particular are outstanding. Nov 27 '17 at 16:39
• Regarding your edit, what I am finding odd is that if I copy and paste your code, I still get only the first roots you showed there. I am not sure why that happens. Nov 27 '17 at 16:42
• I am impressed that I indeed cannot see the other roots like you showed in your Edit. Any idea of why? {{{x -> 0.0380158, y -> 0.48276, z -> 0.479225}, {x -> 0.291508, y -> 0.326764, z -> 0.381728}}, {{x -> 0.309524, y -> 0.309524, z -> 0.380952}, {x -> ComplexInfinity, y -> 0.5, z -> 0.5}}, {{x -> 0.324093, y -> 0.294404, z -> 0.381503}}, {{x -> 0.33648, y -> 0.280628, z -> 0.382892}}, {{x -> 0.347344, y -> 0.267788, z -> 0.384868}}} Nov 27 '17 at 20:38
• @LauraK. What version of Mathematica are you using? I solved the system using V11.1.1. Nov 27 '17 at 20:58
• @LauraK. I have updated my answer to address the issue of eliminating solution where (x - 0.5)/(y - 0.5) is not well defined. Nov 27 '17 at 23:29