# Solving system of PDEs with NDSolve

I am trying to solve a system of coupled PDEs with zero-flux boundary conditions on a large domain. I have two problems: 1) Is there a possibility to use results of NDSolve as inititial conditions? As the results are given by an interpolating function, I guess this could be difficult. However, my RAM is not sufficient for long evaluations and evaluation using {t,0,10} and then {t,10,20} (for instance) does not seem to match.

2) If I use my initial conditions (see following code), I think the evaluation is wrong. A stationary bump in z emerges. I think the reason for this is the error NDSolve::mxsst: Using maximum number of grid points 100 allowed by the MaxPoints or MinStepSize options for independent variable x. How can I solve the system anyway?

The code:

(* Parameters *)
eps = 1.4434; m = 0.3; c11 = 0.1732;
(* PDEs *)
pde1 := D[pp[t, x, y], t] == 0.05*Laplacian[pp[t, x, y], {x, y}] + pp[t, x, y]*(1 - c11*pp[t, x, y] - z[t, x, y]/(1 + pp[t, x, y]^2));
pde2 := D[z[t, x, y], t] == 0.05*Laplacian[z[t, x, y], {x, y}] + z[t, x, y]*(eps*pp[t, x, y]/(1 + pp[t, x, y]^2) - m);
(* Initial conditions *)
ic1[x_, y_] := Which[Sqrt[(x - 50)^2 + (y - 50)^2] < 1, 6, True, 0];
ic2[x_, y_] :=
Which[Sqrt[(x - 50)^2 + (y - 50)^2] < 1, 0.5, True, 1/c11];
(* Numerical approximation using NDSolve with zero-flux boundary conditions*)
soln2d = NDSolve[{pde1, pde2,
(D[pp[t, x, y], x] /. x -> 0) == 0,
(D[pp[t, x, y], y] /. y -> 0) == 0,
(D[z[t, x, y], x] /. x -> 0) == 0,
(D[z[t, x, y], y] /. y -> 0) == 0,
(D[pp[t, x, y], x] /. x -> 100) == 0,
(D[pp[t, x, y], y] /. y -> 100) == 0,
(D[z[t, x, y], x] /. x -> 100) == 0,
(D[z[t, x, y], y] /. y -> 100) == 0,
z[0, x, y] == ic1[x, y],
pp[0, x, y] == ic2[x, y]},
{pp, z}, {t, 0, 500}, {x, 0, 100}, {y, 0, 100}];


• 1. " As the results are given by an interpolating function, I guess this could be difficult." No, there'll be no difficulty. The usage of InterpolatingFunction is just the same as analytic functions like Sin, Exp, etc. in this case. If you fails to, please show a specific example. 2. For unsmooth i.c. NDSolve will automatically choose too dense a grid, you need to manually control the spatial grid in this case using e.g. Method-> {"MethodOfLines", "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 100, "MinPoints" -> 100, "DifferenceOrder" -> 4}}. – xzczd Sep 3 '19 at 5:27

Your description on how to check the solution was a bit vague, so I did not do it. Try this, though:

(*Parameters*)
eps = 1.4434; m = 0.3; c11 = 0.1732;
(*PDEs*)
pde1 :=
D[pp[t, x, y], t] ==
0.05*Laplacian[pp[t, x, y], {x, y}] +
pp[t, x, y]*(1 - c11*pp[t, x, y] - z[t, x, y]/(1 + pp[t, x, y]^2));
pde2 := D[z[t, x, y], t] ==
0.05*Laplacian[z[t, x, y], {x, y}] +
z[t, x, y]*(eps*pp[t, x, y]/(1 + pp[t, x, y]^2) - m);
(*Initial conditions*)

ic1[x_, y_] := Which[Sqrt[(x - 50)^2 + (y - 50)^2] < 1, 6, True, 0];
ic2[x_, y_] :=
Which[Sqrt[(x - 50)^2 + (y - 50)^2] < 1, 0.5, True, 1/c11];
(*Numerical approximation using NDSolve with zero-flux boundary \
conditions*)
soln2d =
Monitor[NDSolveValue[{pde1, pde2, z[0, x, y] == ic1[x, y],
pp[0, x, y] == ic2[x, y]}, {pp, z}, {t, 0,
500}, {x, y} \[Element] Rectangle[{0, 0}, {100, 100}],
EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])]
, monitor]


This uses the finite element method to solve the PDEs.

To your firs question: Yes you can use an InterpolatingFunction as an initial condition, just like any other function. Being able to do so is a strong point of Mathematica.

• 1: OK, then I'll try this. Is there also the possibility to store this in a file and load it later as initial condition? I ask as Mathematica sometimes crashes. 2: Does the "[Element] Rectangle[..]" contain zero flux bc? I copied this code and I get the following error messages: CoefficientArrays::poly: pp$2407-0.05 (pp$2408+pp$2409)-pp (1-0.1732 pp-z/(1+pp^2)) is not a polynomial. NDSolveValue::femper: PDE parsing error of {pp$3017-0.05 (pp$3018+pp$3019)-pp (1-0.1732 pp-z/(1+pp^2)),-(-0.3+(1.4434 pp)/(1+pp^2)) z+z$3017-0.05 (z$3020+z\$3021)}. Inconsistent equation dimensions. – gumpel Sep 3 '19 at 7:55
• @gumpel, you can store an interpolating function just like any other expression. – user21 Sep 3 '19 at 7:59
• OK, thanks. Any idea regarding the second problem? – gumpel Sep 3 '19 at 8:01
• Or does the code you have given work for you? – gumpel Sep 3 '19 at 8:42
• @gumpel, the code runs without errors. If it does not for you: what version do you have? what issues does it show? – user21 Sep 3 '19 at 9:10