Solving system of PDEs with NDSolve

Question

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}];

Thank you for your help!

Answer 1

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.