# How can I speed up NDSolve

I have written a simulation that solve a pde, I used a method I found on this site to speed up the calculation but it is still slow. Is there any way to speed up the prosses?

L = 25;(*length of square*)
T = 70;(*Time integration*)
pts = 1000;
f = 0.2; (*The infiltration contrast between vegetated and bare-soil \
areas*)
A = 40.0; (*The maximal infiltration rate, obtained in densely \
vegetated areas*)
n = 30.0; (*the soil-water evaporation rate*)
R = 30.0; (*the shading effect*)
E1 = 2.0; (*root to shoot ratio of B1*)
E2 = 0.0;(*root to shoot ratio of B2*)
M1 = 0.75; (*the mortality rate of B1*)
M2 = 7.5; (*the mortality rate of B2*)
K1 = 5.0; (*the maximum standing biomass of B1*)
K2 = 0.5; (*the maximum standing biomass of B2*)
Λ1 = 0.05; (*growth rate of B1*)
Λ2 = 1.0; (*growth rate of B2*)
Γ1 = 10.0; (*the water uptake rate of B1*)
Γ2 = 15.0; (*the water uptake rate of B2*)
q = 5.0; (*reference biomass value beyond which he effect of \
vegetation on infiltration becomes significant*)
Y1 = 1.0; (*relative contribution of the woody specie (B1)*)
Y2 = 100; (*relative contribution of the herbaceous specie (B2)*)
DB1 = 0.1; (*the biomass expansion rate of B1*)
DB2 = 0.05; (*the biomass expansion rate of B2*)
DW = 3.0; (*the soil-water diffusion*)
DH = 1.0; (*the rate of surface-water spread*)
i = A (Y1 B1[t, x, y] + Y2 B2[t, x, y] + q f)/(
Y1 B1[t, x, y] + Y2 B2[t, x, y] + q );
P = 1.7;

(*system of nonlinear PDE*)
pde = {D[B1[t, x, y],
t] == Λ1 W[t, x, y] B1[t, x,
y] (1 - B1[t, x, y]/K1) (1 + E1 B1[t, x, y])^2 -
M1 B1[t, x, y] +
DB1 (D[B1[t, x, y], x, x] + D[B1[t, x, y], y, y]),
D[B2[t, x, y],
t] == Λ2 W[t, x, y] B2[t, x,
y] (1 - B2[t, x, y]/K2) (1 + E2 B2[t, x, y])^2 -
M2 B2[t, x, y] +
DB2 (D[B2[t, x, y], x, x] + D[B2[t, x, y], y, y]),
D[W[t, x, y], t] ==
i H[t, x, y] - (n W[t, x, y])/(
1 + R B1[t, x, y]/K1 ) - Γ1 W[t, x, y] B1[t, x,
y] (1 + E1 B1[t, x, y])^2 - Γ2 W[t, x, y] B2[t,
x, y] (1 + E2 B2[t, x, y])^2 +
DW (D[W[t, x, y], x, x] + D[W[t, x, y], y, y]),
D[H[t, x, y], t] ==
P - i H[t, x, y] +
DH (D[H[t, x, y]^2, x, x] + D[H[t, x, y]^2, y, y])};
(*Newman boundary condition*)
bc = {(D[B1[t, x, y], x] /. x -> -L) ==
0, (D[B1[t, x, y], x] /. x -> L) ==
0, (D[B1[t, x, y], y] /. y -> -L) ==
0, (D[B1[t, x, y], y] /. y -> L) ==
0, (D[B2[t, x, y], x] /. x -> -L) ==
0, (D[B2[t, x, y], x] /. x -> L) ==
0, (D[B2[t, x, y], y] /. y -> -L) ==
0, (D[B2[t, x, y], y] /. y -> L) ==
0, (D[W[t, x, y], x] /. x -> -L) ==
0, (D[W[t, x, y], x] /. x -> L) ==
0, (D[W[t, x, y], y] /. y -> -L) ==
0, (D[W[t, x, y], y] /. y -> L) == 0,
(D[H[t, x, y], x] /. x -> -L) == 0, (D[H[t, x, y], x] /. x -> L) ==
0, (D[H[t, x, y], y] /. y -> -L) ==
0, (D[H[t, x, y], y] /. y -> L) == 0};

ic = {B1[0, x, y] == If[y < 0, 1, 0],
B2[0, x, y] == If[Abs[y] < 1, 1, 0],
W[0, x, y] == If[y < 0, 9, 12],
H[0, x, y] == If[y < 0, 3, 5]};

eqns = Flatten@{pde, bc, ic};

sol = NDSolve[eqns, {B1, B2, W, H}, {t, 0, T}, {x, -L, L}, {y, -L, L},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> pts, "MaxPoints" -> pts},
Method -> {"IDA", "ImplicitSolver" -> {"GMRES"}}}];


maybe using multiple CPU?

(the ecology tag related because this simulation is for a model that calculate the time and space growth of plants in desert like area.

• If no option is specified the integration is faster. In other words why do you specify pts=1000? Also, this does run on multiple CPUs. Feb 20, 2020 at 12:40
• no particular reason, just wanted high definition. and it does not run on multiple CPU, how can i make it so it does? Feb 20, 2020 at 12:46
• Works in parallel on my machine. Have you changed the setup on yours? Feb 20, 2020 at 15:17
• As @user21 suggests, requiring 1000 points is a big problem. Since you have a 2D system, I think that means that NDSolve discretizes the PDE into 1000*1000 ODEs. Not surprising that's going to be slow. Feb 20, 2020 at 16:26

There are no gradients in x and y at L> 10. Nothing happens when T> 10. The number pts can be taken 81. Therefore, we can take L=10, T=10, pts=81. The initial data can be expressed using UnitStep[]. Then we get

L = 10;(*length of square*)T = 10;(*Time integration*)pts = 81;
f = 0.2;(*The infiltration contrast between vegetated and bare-soil \
areas*)A = 40.0;(*The maximal infiltration rate,obtained in densely \
vegetated areas*)n = 30.0;(*the soil-water evaporation rate*)R = \
30.0;(*the shading effect*)E1 = 2.0;(*root to shoot ratio of B1*)E2 = \
0.0;(*root to shoot ratio of B2*)M1 = 0.75;(*the mortality rate of \
B1*)M2 = 7.5;(*the mortality rate of B2*)K1 = 5.0;(*the maximum \
standing biomass of B1*)K2 = 0.5;(*the maximum standing biomass of \
B2*)Λ1 = 0.05;(*growth rate of B1*)Λ2 = \
1.0;(*growth rate of B2*)Γ1 = 10.0;(*the water uptake \
rate of B1*)Γ2 = 15.0;(*the water uptake rate of B2*)q \
= 5.0;(*reference biomass value beyond which he effect of vegetation \
on infiltration becomes significant*)Y1 = 1.0;(*relative contribution \
of the woody specie (B1)*)Y2 = 100;(*relative contribution of the \
herbaceous specie (B2)*)DB1 = 0.1;(*the biomass expansion rate of \
B1*)DB2 = 0.05;(*the biomass expansion rate of B2*)DW = 3.0;(*the \
soil-water diffusion*)DH = 1.0;(*the rate of surface-water spread*)i =
A (Y1 B1[t, x, y] + Y2 B2[t, x, y] + q f)/(Y1 B1[t, x, y] +
Y2 B2[t, x, y] + q);
P = 1.7;
us[y_] := 2 (UnitStep[y] - 1/2)
(*system of nonlinear PDE*)
pde = {D[B1[t, x, y],
t] == Λ1 W[t, x, y] B1[t, x,
y] (1 - B1[t, x, y]/K1) (1 + E1 B1[t, x, y])^2 -
M1 B1[t, x, y] +
DB1 (D[B1[t, x, y], x, x] + D[B1[t, x, y], y, y]),
D[B2[t, x, y],
t] == Λ2 W[t, x, y] B2[t, x,
y] (1 - B2[t, x, y]/K2) (1 + E2 B2[t, x, y])^2 -
M2 B2[t, x, y] +
DB2 (D[B2[t, x, y], x, x] + D[B2[t, x, y], y, y]),
D[W[t, x, y], t] ==
i H[t, x,
y] - (n W[t, x, y])/(1 +
R B1[t, x, y]/K1) - Γ1 W[t, x, y] B1[t, x,
y] (1 + E1 B1[t, x, y])^2 - Γ2 W[t, x, y] B2[t,
x, y] (1 + E2 B2[t, x, y])^2 +
DW (D[W[t, x, y], x, x] + D[W[t, x, y], y, y]),
D[H[t, x, y], t] ==
P - i H[t, x, y] +
DH (D[H[t, x, y]^2, x, x] + D[H[t, x, y]^2, y, y])};
(*Newman boundary condition*)
bc = {(D[B1[t, x, y], x] /. x -> -L) ==
0, (D[B1[t, x, y], x] /. x -> L) ==
0, (D[B1[t, x, y], y] /. y -> -L) ==
0, (D[B1[t, x, y], y] /. y -> L) ==
0, (D[B2[t, x, y], x] /. x -> -L) ==
0, (D[B2[t, x, y], x] /. x -> L) ==
0, (D[B2[t, x, y], y] /. y -> -L) ==
0, (D[B2[t, x, y], y] /. y -> L) ==
0, (D[W[t, x, y], x] /. x -> -L) ==
0, (D[W[t, x, y], x] /. x -> L) ==
0, (D[W[t, x, y], y] /. y -> -L) ==
0, (D[W[t, x, y], y] /. y -> L) ==
0, (D[H[t, x, y], x] /. x -> -L) ==
0, (D[H[t, x, y], x] /. x -> L) ==
0, (D[H[t, x, y], y] /. y -> -L) ==
0, (D[H[t, x, y], y] /. y -> L) == 0};

ic = {B1[0, x, y] == -(us[y] - 1)/2,
B2[0, x, y] == (-us[y - 1] + 1) (us[y + 1] + 1)/4,
W[0, x, y] == 9 + (12 - 9) (us[y] + 1)/2,
H[0, x, y] == 3 + (5 - 3) (us[y] + 1)/2};

eqns = Flatten@{pde, bc, ic};
sol = NDSolve[eqns, {B1, B2, W, H}, {t, 0, T}, {x, -L, L}, {y, -L, L},
Method -> {"IndexReduction" -> Automatic,
"EquationSimplification" -> "Residual",
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> pts, "MaxPoints" -> pts,
"DifferenceOrder" -> 2}}}] // Quiet;

{DensityPlot[B1[T, x, y] /. sol, {x, -L, L}, {y, -L, L},
ColorFunction -> "Rainbow", PlotLegends -> Automatic],
DensityPlot[B2[T, x, y] /. sol, {x, -L, L}, {y, -L, L},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotRange -> All],
DensityPlot[W[T, x, y] /. sol, {x, -L, L}, {y, -L, L},
ColorFunction -> "Rainbow", PlotLegends -> Automatic],
DensityPlot[H[T, x, y] /. sol, {x, -L, L}, {y, -L, L},
ColorFunction -> "Rainbow", PlotLegends -> Automatic]}