7
$\begingroup$

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.

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

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}]

plot

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

This update addresses an issue raised by the OP in comments to this answer.

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}]

plot

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you very much for the answer. I have been learning a lot on this community. Your answers in particular are outstanding. $\endgroup$ – Laura K Nov 27 '17 at 16:39
  • $\begingroup$ 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. $\endgroup$ – Laura K Nov 27 '17 at 16:42
  • $\begingroup$ 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}}} $\endgroup$ – Laura K Nov 27 '17 at 20:38
  • $\begingroup$ @LauraK. What version of Mathematica are you using? I solved the system using V11.1.1. $\endgroup$ – m_goldberg Nov 27 '17 at 20:58
  • 1
    $\begingroup$ @LauraK. I have updated my answer to address the issue of eliminating solution where (x - 0.5)/(y - 0.5) is not well defined. $\endgroup$ – m_goldberg Nov 27 '17 at 23:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.