3
$\begingroup$

I am trying to numerically solve the steady-state behaviour of the Allee effect (with one spatial dimension, diffusion):

$\frac{\partial n}{\partial t} = D \frac{\partial^2 n }{\partial x^2} + rn\left(1-\frac{n}{K}\right)\left(1-\frac{a}{n}\right)$, $0 < a < K$,

where $D$ is the diffusivity, $r$ is the growth rate, $K$ is the carrying capacity, $a$ is the threshold population, and with boundary conditions:

$n(0,t)=n(L,t)=0$,

using NDSolve. However, I run into some issues with the population, $n$, becoming negative, when I use the parameters $D = 0.1$, $L = 50.0$, $r \approx 0.008$, $K = 150$, and $a = 100$; with initial population distribution $n_0\left(x\right) = 15 \cos \left( \frac{2 \pi}{L} x \right) + 175$. I'm aware that $n=0$ gives division by zero, and that $n_0(x)$ does not satisfy the boundary conditions, but I wanted to try this anyways. Sometime around $t=1800$ there will be negative values for $n$ near the boundaries, $x = 0.0$ and $x = 50.0$.

A sample of my code is here:

Diffusivity = 0.1;
L = 50.0;
r = 0.008;

ic = {n[x,0] == 15*Cos[(2*Pi*x)/L]+175,n[0,t] == 0 , n[L,t] == 0};
eqn = D[n[x, t], t] - Diffusivity*D[n[x, t], x, x] - 
   r*n[x, t]*(1 - (1/150)*n[x, t])*(1 - 100/n[x, t]);
s = NDSolve[{eqn == 0}~Join~ic,n,{x,0,L},{t,0,100}];

and for plotting solutions at different times ($t=0$ here):

f0[x_] = Evaluate[n[x,t] /.s /.t −> 0];

I read a similar posting here, but that was for a slightly different situation.

Is there any way I can add the condition $n\left(x,t\right) \geq 0$ into my already existing code?

This is a problem I had for course which I completed recently. So, I'm really just curious how to go about this problem with Mathematica, and maybe the solution to this problem will help someone else in the future.

Thanks in advanced for any suggestions!

$\endgroup$
  • 1
    $\begingroup$ While there is no universal equation for an Allee effect, I think the one you're using is quite problematic because dn/dt needs to equal zero when there is no population. This is surely the cause of your negative population size, which a good population model will never produce. Compare with the formulation in the wikipedia article. $\endgroup$ – Chris K Apr 26 '16 at 21:57
  • $\begingroup$ Okay, thank you Chris! $\endgroup$ – Canada709 Apr 26 '16 at 22:28
4
$\begingroup$

As mentioned by Chiris K in the comment above, the unexpected behavior of the solution is just the nature of the not-that-proper model, so I think the best solution is to modify the model a bit i.e. adding a small number eps to denominator of $\frac{a}{n}$:

Diffusivity = 0.1; 
L = 50.; 
r = 0.008; 
ic = {n[x, 0] == 15 Cos[(2 Pi x)/L] + 175, n[0, t] == 0, n[L, t] == 0}; 
eqn = D[n[x, t], t] - Diffusivity D[n[x, t], x, x] - 
       r n[x, t] (1 - (1/150) n[x, t]) (1 - 100/(n[x, t] + eps)); 
s = ParametricNDSolveValue[Join[{eqn == 0}, ic], n, {x, 0, L}, 
       {t, 0, 1800}, eps]; 

smodified = s[1]    
soriginal = s[0]

Plot3D[smodified[x, t], {x, 0, L}, {t, 0, 1800}]

Mathematica graphics

As one can see, the modified model i.e. smodified is always nonnegative. Comparing it to the original one i.e. soriginal, we see their behaviors are almost the same before soriginal is out of whack:

Manipulate[Plot[{soriginal[x, t], smodified[x, t]}, {x, 0, L}, 
  PlotRange -> {-50, 190}], {t, 0, 766}]

enter image description here

$\endgroup$
  • $\begingroup$ Sorry that I disappeared, after I read and thought about Chris K's answer I accepted that I was doing something wrong and started focusing elsewhere. Your solution looks quite nice, thank you!! Really like the Manipulate demonstration. Hope it has helped others! $\endgroup$ – Canada709 Sep 20 '18 at 20: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.