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:


    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...

  • $\begingroup$ Hi Tamás! Could you give a quick background on the biology of your system and what the variables represent? $\endgroup$
    – Chris K
    Mar 31 '16 at 16:57
  • $\begingroup$ 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. $\endgroup$ Apr 1 '16 at 7:25
  • $\begingroup$ S=starters, F=food, NS=non-starters $\endgroup$ Apr 1 '16 at 18:51
  • $\begingroup$ 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! $\endgroup$
    – Chris K
    Apr 3 '16 at 14:13
  • $\begingroup$ 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. $\endgroup$ 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}]

fig 1


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.