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!
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$