1
$\begingroup$

Consider a (non-linear) system of equations with two parameters x and y and multiple variables z1, z2, z3, z4

 eqs = {
   -0.001 z1 + (0.001 + 0.1 z1) (1 - z1 - z2 - z3 - z4) - 0.06 z1 (z3 + z4) == 0,
   0.001 z1 - 0.011 z2 + (0.05 + x (1 - y)) z3 - 0.6 z2 (z3 + z4) == 0,
   (1 - z1 - z2 - z3 - z4) (0.001 + 0.6 z3 + 0.6 z4) - (0.05 + x) z4 == 0,
   0.001 z2 - (0.06 + x (1 - y)) z3 + 0.06 z1 (z3 + z4) + 0.6 z2 (z3 + z4) == 0
   };

The system has as many equations as variables (4) and can be solved given values for the parameters

NSolve[Join[eqs /. {x -> .6, y -> .5}, {z1 >= 0, z2 >= 0, z3 >= 0, z4 >= 0}], {z1, z2, z3, z4}]

Out: {{z1 -> 0.897951, z2 -> 0.0915341, z3 -> 0.000372715, z4 -> 0.0000192353}}

Sometimes it has more than one solution

NSolve[Join[eqs /. {x -> .7, y -> .8}, {z1 >= 0, z2 >= 0, z3 >= 0, z4 >= 0}], {z1, z2, z3, z4}]

Out:
{
 {z1 -> 0.894114, z2 -> 0.0942418, z3 -> 0.00107801, z4 -> 0.0000233464},
 {z1 -> 0.274034, z2 -> 0.256676, z3 -> 0.249387, z4 -> 0.0417326},
 {z1 -> 0.636264, z2 -> 0.250085, z3 -> 0.0626664, z4 -> 0.00259097}
}

Now I want to show the solutions for z1 in the x,y plane. The following works ok

list = Flatten[Table[{x, y, z1 /. #} & /@ NSolve[Join[eqs, {z1 >= 0, z2 >= 0, z3 >= 0, z4 >= 0}], {z1, z2, z3, z4}], {x, 0, 1, .025}, {y, 0, 1, .025}], 2];
ListPointPlot3D[list, ColorFunction -> "Rainbow", AxesLabel -> {"x", "y", "z1"}]

enter image description here

But I find it a bit difficult to interpret. I was hoping to show the solution as continuous planes, as in the help file for ContourPlot3D

enter image description here

However I don't know how to use ContourPlot3D because of the other variables in the system z2, z3, z4. My first idea was to reduce the system to a single equation with just z1. That didn't work though because the system is non-linear and Solve didn't find a solution.

Any ideas?

Thanks!

$\endgroup$

1 Answer 1

3
$\begingroup$

Edit

The surface still have some holes.

eqs = {-0.001 z1 + (0.001 + 0.1 z1) (1 - z1 - z2 - z3 - z4) - 
     0.06 z1 (z3 + z4) == 0, 
   0.001 z1 - 0.011 z2 + (0.05 + x (1 - y)) z3 - 0.6 z2 (z3 + z4) == 
    0, (1 - z1 - z2 - z3 - z4) (0.001 + 0.6 z3 + 0.6 z4) - (0.05 + 
        x) z4 == 0, 
   0.001 z2 - (0.06 + x (1 - y)) z3 + 0.06 z1 (z3 + z4) + 
     0.6 z2 (z3 + z4) == 0};
eqs = Rationalize[eqs, 0];
data = Flatten[
   Table[{x, y, z1} /. 
     NSolve[Join[
       eqs, {z1 >= 0, z2 >= 0, z3 >= 0, z4 >= 0}], {z1}, {z2, z3, 
       z4}], {x, 0, 1, .01}, {y, 0, 1, .01}], 2];
surf = ListSurfacePlot3D[data, MaxPlotPoints -> 60, 
  ColorFunction -> "Rainbow", PerformanceGoal -> "Quality", 
  ViewPoint -> {-2.60, -1.74, 1.27}]

enter image description here

Original

The solution have several branch.

eqs = {-0.001 z1 + (0.001 + 0.1 z1) (1 - z1 - z2 - z3 - z4) - 
     0.06 z1 (z3 + z4) == 0, 
   0.001 z1 - 0.011 z2 + (0.05 + x (1 - y)) z3 - 0.6 z2 (z3 + z4) == 
    0, (1 - z1 - z2 - z3 - z4) (0.001 + 0.6 z3 + 0.6 z4) - (0.05 + 
        x) z4 == 0, 
   0.001 z2 - (0.06 + x (1 - y)) z3 + 0.06 z1 (z3 + z4) + 
     0.6 z2 (z3 + z4) == 0};
eqs = Rationalize[eqs, 0];
sol = Solve[eqs, {z1}, {z2, z3, z4}];
Plot3D[z1 /. sol// Evaluate, {x, 0, 1}, {y, 0, 1},MaxRecursion -> 4, ColorFunction -> "Rainbow", 
 AxesLabel -> {"x", "y", "z1"}, BoxRatios -> {1, 1, 1}]

enter image description here

$\endgroup$
2
  • $\begingroup$ Thank you!! I managed to get rid of the holes in the edited solution by increasing the number of datapoints, now it's perfect. $\endgroup$
    – tukan
    Jan 17, 2022 at 17:04
  • $\begingroup$ Still, I'm curious about the original solution. Is there a way to show with Plot3D only the relevant solutions (all variables z1, z2, z3 and z4 being real and non-negative)? $\endgroup$
    – tukan
    Jan 17, 2022 at 17:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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