I am trying to solve a system of coupled PDEs with zero-flux boundary conditions on a large domain. I have two problems:
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.
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!