Return to Answer

deleted 101 characters in body

Source Link

edit approved Feb 21, 2020 at 4:33

5.8k
2
14
43

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*)\[CapitalLambda]1Λ1 = 0.05;(*growth rate of B1*)\[CapitalLambda]2Λ2 = \
1.0;(*growth rate of B2*)\[CapitalGamma]1Γ1 = 10.0;(*the water uptake \
rate of B1*)\[CapitalGamma]2Γ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] == \[CapitalLambda]1Λ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] == \[CapitalLambda]2Λ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) - \[CapitalGamma]1Γ1 W[t, x, y] B1[t, x, 
       y] (1 + E1 B1[t, x, y])^2 - \[CapitalGamma]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]}

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*)\[CapitalLambda]1 = 0.05;(*growth rate of B1*)\[CapitalLambda]2 = \
1.0;(*growth rate of B2*)\[CapitalGamma]1 = 10.0;(*the water uptake \
rate of B1*)\[CapitalGamma]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] == \[CapitalLambda]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] == \[CapitalLambda]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) - \[CapitalGamma]1 W[t, x, y] B1[t, x, 
       y] (1 + E1 B1[t, x, y])^2 - \[CapitalGamma]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]}

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

Source Link

created Feb 21, 2020 at 1:09

Alex Trounev

48.8k
3
51
115

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*)\[CapitalLambda]1 = 0.05;(*growth rate of B1*)\[CapitalLambda]2 = \
1.0;(*growth rate of B2*)\[CapitalGamma]1 = 10.0;(*the water uptake \
rate of B1*)\[CapitalGamma]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] == \[CapitalLambda]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] == \[CapitalLambda]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) - \[CapitalGamma]1 W[t, x, y] B1[t, x, 
       y] (1 + E1 B1[t, x, y])^2 - \[CapitalGamma]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]}