Constraining Function to be positive with NDSolve, Allee effect

Question

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!

Answer 1

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

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