# NDSolve cannot consistently handle a simple 2D reaction-diffusion PDE with spiky (and random) initial conditions

I am still fighting with a modified version of the 2D reaction-diffusion PDE that I have already posted in a previous question ( Reaction-diffusion PDE with NDSolve: either very slow or very inaccurate ). Now I use a slightly different set of initial conditions, and I cannot convince NDSolve to do the task for me properly and within acceptable time. This is the system I use now:

    SetDirectory[NotebookDirectory[]];

Needs["DifferentialEquationsNDSolveProblems"]
Needs["DifferentialEquationsNDSolveUtilities"]
Needs["FunctionApproximations"]
a = AbsoluteTime[];

rs = 0.1;
rn = 0.1;
Diff = 1;
Difn = 0.001;

sigmas0 = 2;
sigmans0 = 1;

s0 = 1000.;           (* Initial biomass of Starters *)
cs =   10;            (* Number of starter colonies *)
ns0 = .2;             (* Initial biomass of Non-starters *)
cns = 10;             (* Number of non-starter colonies *)

time = 500;
size = 50;
rep = 1;

k = 2 cs + 1;
l = 2 cs + cns;

rnd = RandomReal[2 size, 2 cs + 2 cns];
rnd = rnd - size;

soln = NDSolve[{
D[S[t, x, y], t] == -rs S[t, x, y],
D[F[t, x, y], t] ==
Diff (D[F[t, x, y], x, x] + D[F[t, x, y], y, y]) +
rs S[t, x, y] - rn F[t, x, y] NS[t, x, y],
D[NS[t, x, y], t] ==
Difn (D[NS[t, x, y], x, x] + D[NS[t, x, y], y, y]) +
rn F[t, x, y] NS[t, x, y],

S[t, -size, y] == S[t, size, y], S[t, x, -size] == S[t, x, size],
F[t, -size, y] == F[t, size, y], F[t, x, -size] == F[t, x, size],
NS[t, -size, y] == NS[t, size, y],
NS[t, x, -size] == NS[t, x, size],

S[0, x, y] ==
Sum[s0/cs Exp[-((x - rnd[[n]])^2/(2 sigmas0^2) + (y -
rnd[[cs + n]])^2/(2 sigmas0^2))], {n, 1, cs}],
F[0, x, y] == 0.,
NS[0, x, y] ==
Sum[ns0/cns Exp[-((x - rnd[[nn]])^2/(2 sigmans0^2) + (y -
rnd[[cns + nn]])^2/(2 sigmans0^2))], {nn, k, l}]},

{S, F, NS}, {t, 0, time}, {x, -size, size}, {y, -size, size},
Method -> {"StiffnessSwitching", "NonstiffTest" -> False}];


Since the culprit is most probably the spiky initial condition both in S[0,x,y] and in NS[0,x,y], I have tried the "StiffnessSwitching" method in this attempt, but even then the kernel stops with the error message

NDSolve::ndsz: At t == 28.71100997366924, step size is effectively zero; singularity or stiff system suspected. >>


Singularity is out of the question (at least in principle it is), as the initial conditions are the sums of Gaussian functions (for S and NS) and the constant zero function (for F).

Does anyone have an idea how to make this work consistently (i.e., regardless of the actual parameters defining the initial conditions)? I have also tried setting the MaxStepSize option to 0.1 - then I got no error message, but no result either, for more than 10 hours...

• Hi Tamás! Could you give a quick background on the biology of your system and what the variables represent? Mar 31 '16 at 16:57
• Hi Chris, sure I can. This is a model of bacterial "cannibalism": One batch of bacterial strains (called "non-starters" feed on the corpses of another ("starter") strain that is distributed in patches within a solid medium (like agar). The lysed "starter" patches emit nutrients into the medium. The nutrients diffuse within the medium and are used by the (also randomly distributed) "non-starter" strain for its population growth. Apr 1 '16 at 7:25
• S=starters, F=food, NS=non-starters Apr 1 '16 at 18:51
• Following the documentation on the method of lines, I tried tweaking various options through Method -> {"PDEDiscretization" -> {"MethodOfLines", "TemporalVariable" -> t, "SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 100, "MaxPoints" -> 100, "DifferenceOrder" -> 4}}}. I couldn't find anything that helped with your default parameter set. But setting Difn=0 made it work very fast, even with default options. I bet this will be insignificantly different from Difn=0.001 in a 100x100 domain over 500 time steps. You can rely on Diff=1 to bring the food to the bacteria! Apr 3 '16 at 14:13
• Thank you, Chris! You are absolutely right: in fact the bacteria cannot really move within the semi-solid matrix. Yet, the NS colonies will definitely increase in size as their population increases - that's what I intended to model with their marginal diffusion rate. Since this aspect of the dynamics is not really essential, I may simply use somewhat larger initial patches (larger spatial variance for the initial NS patches) and zero diffusion rate for them. I will try that. Apr 3 '16 at 21:00

NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent.


So there is a problem from the very beginning. What can I recommend? We set

sigmas0 = 1/2;sigmans0 = 1/2;rnd = RandomReal[{-size/3, size/3}, 2 cs + 2 cns];


Accordingly, remove the following line

rnd = rnd - size;


After that, we have nice pictures not for every random distribution, but for many of them

    Table[Plot3D[F[t, x, y] /. soln, {x, -size, size}, {y, -size, size},
PlotRange -> All, ColorFunction -> Hue, Mesh -> None,
AxesLabel -> {"x", "y", ""}, PlotLabel -> t], {t, .2*time,
time, .2*time}]

Table[Plot3D[NS[t, x, y] /. soln, {x, -size, size}, {y, -size, size},
PlotRange -> All, ColorFunction -> Hue, Mesh -> None,
AxesLabel -> {"x", "y", ""}, PlotLabel -> t], {t, 0, time, .2*time}] 