# Mathematica 12: avoid crashing when using NDSolve on a large domain

I have a set of coupled PDEs, that I'd like to solve numerically using Mathematica (11.3 or 12.0). The problem is, that I need to solve in on a large domain to make sure that there are no boundary effects. Mathematica seems to have problems with such calculationg (probably just due to my RAM). So I decided to split the problem into small time intervals and just export parts of the solution and start with the last solution as new initial conditions. Now, Mathematica keeps crashin in this calculation, e.g., after 30 runs of the loop. Why so? Is there any way to avoid that?

Here is the code:

$HistoryLength = 0; (* save *) (*PDEs*) pde11 := D[pp[t, x], t] == 1.*Laplacian[pp[t, x], {x}] + pp[t, x]*(1 - c11*pp[t, x] - z[t, x]/(1 + pp[t, x]^2)); pde21 := D[z[t, x], t] == 1.*Laplacian[z[t, x], {x}] + z[t, x]*(eps*pp[t, x]/(1 + pp[t, x]^2) - m); (*Initial conditions*) lo = 7498; hi = 7502; domlen = 15000; ic11[x_] := Which[x > lo && x < hi, 6, True, 0]; ic21[x_] := Which[x < hi && x > lo, 0.5, True, 1/c11]; eps = 1.4434; m = 0.3; c11 = 0.1732; tfin = 30; For[i = 0, i <= IntegerPart[6000/30], i++, Print[i]; sol1d = NDSolve[{pde11, pde21, z[0, x] == ic11[x], pp[0, x] == ic21[x]}, {pp, z}, {t, 0, tfin}, {x, 0, domlen}, MaxStepSize -> 0.1]; resultsForExport = {}; For[j = 0, j < tfin, j = j + 0.1, resultsForExport = Append[resultsForExport, Evaluate[z[j, 7500]] /. sol1d]; ]; resultsForExport = Flatten[resultsForExport]; Export["largedomain" <> ToString[i] <> ".dat", resultsForExport]; ic11[x_] := sol1d[[1, 2, 2]][tfin, x]; ic21[x_] := sol1d[[1, 1, 2]][tfin, x]; ]  I guesss my exporing may be not too elegant - sorry for that. I've ran this several times now (takes quite long) and it keeps crashing between $$i=20$$ and $$i=30$$. Any thoughts/help/workarounds/comments are appreciated. ## 1 Answer This problem can be solved in 0.03 seconds on my laptop. It is only necessary to add boundary conditions and get rid of Which[] (replaced by Piecewise[]) (*PDEs*)pde11 = D[pp[t, x], t] == 1.*Laplacian[pp[t, x], {x}] + pp[t, x]*(1 - c11*pp[t, x] - z[t, x]/(1 + pp[t, x]^2)); pde21 = D[z[t, x], t] == 1.*Laplacian[z[t, x], {x}] + z[t, x]*(eps*pp[t, x]/(1 + pp[t, x]^2) - m); bc = {z[t, 0] == ic11[0], z[t, domlen] == ic11[domlen], pp[t, 0] == ic21[0], pp[t, domlen] == ic21[domlen]}; ic = {z[0, x] == ic11[x], pp[0, x] == ic21[x]};(*Initial conditions*) lo = 7498; hi = 7502; domlen = 15000; ic11[x_] := Piecewise[{{6, lo < x < hi}, {0, True}}]; ic21[x_] := Piecewise[{{1/2, lo < x < hi}, {1/c11, True}}]; eps = 1.4434; m = 0.3; c11 = 0.1732; tfin = 6000; sol1d = NDSolve[{pde11, pde21, bc, ic}, {pp, z}, {t, 0, tfin}, {x, 0, domlen}];  Figure 1 shows the functions on a large and small scale in time. It is seen that time is spent mostly on useless calculations. {Plot3D[pp[t, x] /. sol1d, {t, 0, tfin}, {x, 0, domlen}, Mesh -> None, PlotRange -> All], Plot3D[z[t, x] /. sol1d, {t, 0, tfin}, {x, 0, domlen}, Mesh -> None, PlotRange -> All]} {Plot3D[pp[t, x] /. sol1d, {t, 0, 10}, {x, 0, domlen}, Mesh -> None, PlotRange -> All], Plot3D[z[t, x] /. sol1d, {t, 0, 10}, {x, 0, domlen}, Mesh -> None, PlotRange -> All]}  • If you plot your solution until$t=50$, you will see that there are numerical uncertainties, i.e., negative values emerge which by looking at the equations cannot happen. I tried to include MaxStepSize in your method which makes calculations substantially longer. – gumpel Sep 15 at 5:01 • @gumpel You see that nothing happens at t> 100. Therefore, you can limit the time as tfin = 100. – Alex Trounev Sep 15 at 5:08 • If I reduce MaxStepSize, there is something happening at$t>100\$. – gumpel Sep 15 at 5:09
• @gumpel Then you need to smooth the initial data. Such initial data are not suitable for fine calculations. – Alex Trounev Sep 15 at 5:14
• However, even if I take smooth initial conditions, I need a smaller step size. – gumpel Sep 15 at 5:45