0
$\begingroup$

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}]
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$ May 12, 2020 at 16:18

1 Answer 1

2
$\begingroup$

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

enter image description here

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

enter image description here

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

enter image description here

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

enter image description here

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

enter image description here

$\endgroup$
1
  • $\begingroup$ Thank you very much, this is much better than my approach. $\endgroup$
    – Luca
    May 12, 2020 at 22:42

Your Answer

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

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