Trying to create a new list by selecting values with a set condition

I am trying to plot 'data' in function of m. To achieve this I tried to get 501 datapoints using divide and then running the equations for each m. This however outputs two solutions, of which the ones with negative numbers are not possible. I have removed the negative numbers but would now like to select from each row the solution without a zero or previously a negative number.

ClearAll[m, L, w, alpha, sigma, H, F]
L = 2000; (* Queen daily laying rate*)
w = 27000;(* Rate at which eclosion approaches L as N gets large*)
alpha = 0.25;(*maximum rate that hive bees will become foragers when \
there are no foragers in the colony*)
sigma = 3/4;(*social inhibition factor*)
(*H = hive bees, F= foraging bees, m= death rate *)

mvector = Subdivide[0.01, 0.4, 500];

(* getting all the data points *)
data = Table[N[#] & /@ Solve[{
0 == L*(H + F)/(w + H + F) - H*(alpha - sigma*(F/(H + F))),
0 == H*(alpha - sigma*(F/(H + F))) - m*F}, {H, F}], {m, mvector}];
data = data /. {x_?Negative -> 0};


Extra: I don't know if the method is the most efficient, I have very little experience. The rest of the code where I plot both solutions can be found below:

(* Calculating N=H+F of the first row -- not sure how to approach the \
two rows *)
newData = Table[{data[[i, 1]][[1]][[2]] + data[[i, 1]][[2]][[2]]}, {i, 1,501}];
newData2 = Table[{data[[i, 2]][[1]][[2]] + data[[i, 2]][[2]][[2]]}, {i, 1,501}];
(* Flatten and removal of negative populations *)

Tble = Flatten[newData];
Tble2 = Flatten[newData2];
Tble = Tble /. {x_?Negative -> 0};
Tble2 = Tble2 /. {x_?Negative -> 0};

Result = Transpose[{mvector, Tble}];
Result2 = Transpose[{mvector, Tble}];

(* Plot the List *)
ListPlot[{Result}]

ListPlot[{Result2}]

• It would probably be more programmatically logical to find the lists/rows/columns/entries with the zeros or negatives or whatever you don’t want, then select what you do want from the remaining lists. Commented May 12, 2020 at 16:18

1 Answer

Your equations can be solved analytically.

ClearAll[m, L, w, alpha, sigma, H, F];
L = 2000;(*Queen daily laying rate*)
w = 27000;(*Rate at which eclosion approaches L as N gets large*)
alpha =
1/4;(*maximum rate that hive bees will become foragers when there are no \
foragers in the colony*)
sigma =
3/4;(*social inhibition factor*)(*H=hive bees,F=foraging bees,m=death rate*)


Include the constraints that H > 0, F > 0, m > 0

soln = Solve[{
0 == L*(H + F)/(w + H + F) - H*(alpha - sigma*(F/(H + F))),
0 == H*(alpha - sigma*(F/(H + F))) - m*F,
H > 0, F > 0, m > 0}, {H, F}] // FullSimplify;

H /. soln


The function is identical in both intervals

SameQ @@ (First /@ (H /. soln))

(* True *)


The function is defined at the intervals' shared boundary

First /@ (H /. soln) /. {m -> 2/9}

(* {8000, 8000} *)


The solution for H is then

solH[m_] =
ConditionalExpression[(H /. soln)[[1, 1]], 0 < m < (2/57)*(6 + Sqrt[17])]


Similarly for F

The function is identical in both intervals

SameQ @@ (First /@ (F /. soln))

(* True *)


The function is defined at the intervals' shared boundary

First /@ (F /. soln) /. {m -> 2/9} // Simplify

(* {1000 (-13 + Sqrt[241]), 1000 (-13 + Sqrt[241])} *)


The solution for F is then

solF[m_] =
ConditionalExpression[(F /. soln)[[1, 1]], 0 < m < (2/57)*(6 + Sqrt[17])]


Plot[{solH[m], solF[m]}, {m, 0.01, 0.36},
PlotLegends -> Placed[{H, F}, {0.5, 0.5}]]


LogPlot[{solH[m], solF[m]}, {m, 0.01, 0.36},
PlotLegends -> Placed[{H, F}, {0.5, 0.3}]]


• Thank you very much, this is much better than my approach.
– Luca
Commented May 12, 2020 at 22:42