Reaction-diffusion PDE with NDSolve: either very slow or very inaccurate

Question

I am struggling to have Mathematica 10.3 solve a system of PDE's (with periodic boundary conditions and random initial conditions), but either it produces a set of very noisy InterpolatingFunction objects that are then nearly impossible to NIntegrate, or it fails to solve it altogether. My best try so far was this:

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

sigmas0 = 1;
sigmans0 = 1;

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

time = 500;
size = 50;

rnd = RandomReal[2 size, 2 cs];
rnd = rnd - size;
top = 0;

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] == ns0 Exp[-(x^2/(2 sigmans0^2) + y^2/(2 sigmans0^2))]},

  {S, F, NS}, {t, 0, time}, {x, -size, size}, {y, -size, size}, 
  MaxStepSize -> .5]

This seems to work (except for Mathematica warning me that the initial and the boundary conditions are inconsistent, but that cannot be helped with random initial and periodic boundary conditions, I guess):

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

However, trying to NIntegrate the resulting InterpolatingFunction objects:

FoodCurve = 
 Table[NIntegrate[
   F[t, x, y] /. soln, {x, -size, size}, {y, -size, size}, 
   Method -> "InterpolationPointsSubdivision", 
   WorkingPrecision -> 100, PrecisionGoal -> 10, 
   AccuracyGoal -> 6], {t, 0, 10}]

I get error messages related to ...

NIntegrate::slwcon: Numerical integration converging too slowly; 
suspect one of the following: singularity, value of the integration is 
0, highly oscillatory integrand, or WorkingPrecision too small. >>

The problem might be either that the result of the integration approaches 0 towards the end of the simulation, or that the integrand is singular at places, which is the result of insufficiently dense sampling, I guess. I have figured out from the manual that decreasing MaxStepSize in NDSolve could smooth the resultant InterpolatingFunction, but if I set it to 0.1 the kernel crashes after half a day of working on the PDF.

Is there another way to speed up NDSolve for this problem and/or to smooth the resulting InterpolatingFunction so that NIntegrate does not complain and does the job fast on it?

Thanks for your attention!

Answer 1

I only address your question about the integration. I tried your code for the integral for t=1. Indeed there was such a message:

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small

I think a good idea would be here to try other methods. For example the "AdaptiveMonteCarlo" yields the following:

    NIntegrate[
  F[t, x, y] /. soln /. t -> 1, {x, -size, size}, {y, -size, size}, 
  Method -> "AdaptiveMonteCarlo", WorkingPrecision -> 3, 
  PrecisionGoal -> 3, AccuracyGoal -> 3] // AbsoluteTiming

(*  {14.8197, {598.}}   *)

and no warnings. So you do not need to have this enormous WorkingPrecision. The "AdaptiveQuasiMonteCarlo"gives the following:

    NIntegrate[
  F[t, x, y] /. soln /. t -> 1, {x, -size, size}, {y, -size, size}, 
  Method -> "AdaptiveQuasiMonteCarlo", WorkingPrecision -> 4, 
  PrecisionGoal -> 4, AccuracyGoal -> 4] // AbsoluteTiming

NIntegrate::maxp: The integral failed to converge after 1000100 integrand evaluations. NIntegrate obtained 597.8`4. and 0.6606`4. for the integral and error estimates. >>

(* {87.3805, {597.8}}  *)

Here as you see I increased the Precision and Accuracy and there is a warning. So generally, I would try different methods from this list:

Have fun!

Answer 2

I computed the PDE system solutions without the MaxStep specification. (As suggested by user21.)

The numerical integration seems to be fast enough with the default NIntegrate options. If I use AccuracyGoal -> 6 the integration becomes 3 times faster. If I remove the integration at $t=0$ I do not get messages.

As Alexei Boulbitch suggested we can try different NIntegrate methods and different precision and accuracy goals.

For this we can use the functions of the package "NIntegrateUtilities" to compare the estimates and integration patterns. Here is an example:

Needs["Integration`NIntegrateUtilities`"]

Table[Labeled[
  NIntegrateSamplingPoints@
   NIntegrate[F[t, x, y] /. soln, {x, -size, size}, {y, -size, size}],
   Row[{"t=", t}]], {t, 0, 10}]

The following image demonstrates that the symbolic pre-processing -- used in the commands above -- improves (and speeds-up) the integration. (It is prevented in the command shown).