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

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?

• Welcome to Mathematica.SE! 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – user9660 Feb 10 '16 at 13:24
• A 2D PDE will probably be time-consuming. I notice that your two diffusion coefficients differ by five orders of magnitude, which probably contributes to the slowness. Could you make them more equal or get away with setting Difn=0? My preliminary running of your code shows that these speed up the NDSolve around 10X. – Chris K Feb 10 '16 at 13:47
• @Chris K Thanks a lot for your reply, Chris! Difn=0 is a possibility for now, but I need to test the results against varying Difn later, so I will have to keep it positive. Are there no options of NDSolve possibly resulting in smooth InterplolatingFunctions without dramatically increasing computing time? – Tamás Czárán Feb 10 '16 at 14:55
• @Chris K Spatial configuration is a key factor in my problem, and it needs to be as close to a real situation as possible, but I will probably have to downgrade to 1D if everything else fails :-) – Tamás Czárán Feb 10 '16 at 15:25
• Why do you have the restriction MaxStepSize -> .5. If you remove that it integrates much faster. Is the result not correct? – user21 Feb 10 '16 at 18:06

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.84. and 0.66064. 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!

• Thanks for the suggestion - I will try it now with the AdaptiveMonteCarlo algorithm. My worry in the first place (the reason I did not try it) was that the initial conditions for S and NS are very spiky Gaussians, which might not be hit by a random sampling algorithm, but that was just a gut feeling - you might well be right that I should use it instead. – Tamás Czárán Feb 15 '16 at 10:50
• Indeed, the single spike of the NS function - which remains spiky all along the simulation - is often missed by the AdaptiveMonteCarlo algorithm - roughly in a fifth of the cases the integral is of zero value which it definitely should not be. (I can't copy code here, unfortunately.) Probably there is a way to increase the MC sampling density, isn't there? – Tamás Czárán Feb 15 '16 at 14:23
• @Tamás Czárán Try to look into tutorial/NIntegrateIntegrationStrategies#65285686/Global Adaptive Monte Carlo and Quasi Monte Carlo Strategies. There are some examples including those with spikes. – Alexei Boulbitch Feb 15 '16 at 16:05
• Thank you, Alexei, the AdaptiveMonteCarlo method with the MinRecursion option set sufficiently high above zero seems to do the trick. The speed is not very impressive, but acceptable. I doubt I can expect much more than that with my 2D diffusion problem :-) Thanks again for your help! – Tamás Czárán Feb 16 '16 at 11:07

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["IntegrationNIntegrateUtilities"]

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

• Thanks a lot, Anton! I'm going to check the NIntegrateSamplingPoints option out, it looks like a promising approach indeed. I'm still baffled how your Mathematica could NDSolve the system without the MaxStepSize option set - mine (ver. 10.3) couldn't... But removing the t=0 solution is a good idea - there the problem is that F is indeed 0 all over the integration domain, so the complaint of NIntegrate is justified (even if I don't really see why this should be so - the integral of the constant zero function should be zero, without complaints). Thank you once again! – Tamás Czárán Feb 15 '16 at 11:06
• @TamásCzárán Integral estimation with sampling points and fitted polynomials has difficulties applying reliably stoping criteria based on precision and accuracy goals when the integral is close to zero. NIntegrate would not know that the integral is zero, may be the integral value is just small. Also, NIntegrate does try to evaluate several scenarios when the integral seems to be zero is it really zero. – Anton Antonov Feb 23 '16 at 14:00
• Yes, that may be the problem here: the integral is often very close to zero. Is there a way to tell NIntegrate to set it to zero after a certain number of evaluation attempts, and go on? Of course it should not miss the spikes in the meantime - maybe these two criteria cannot be met simultaneously... – Tamás Czárán Apr 21 '16 at 7:50