Coupled PDE with nonlinear coefficients: Changing step size to fix large errors?

Question

After help from user @xzczd I was successful in getting Mathematica to numerically solve my PDE (as you can see in this question).

When I give realistic parameters for my differential equation, I see that Mathematica is seeing some errors:

NDSolve::eerr: Warning: scaled local spatial error estimate of 123.96384315484316` at t = 0.000015` in the direction of independent variable z is much greater than the prescribed error tolerance. Grid spacing with 25 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options.

I am interested in solving this set of PDEs without having errors. It suggests using "MaxStepSize" and "MinPoints" but how do I know what to choose for these?

Here is my code:

s = NDSolve[{Derivative[1, 0][p11][t, z] == -3000000*I*Ep[t, z]*(p13[t, z] - p31[t, z]), Derivative[1, 0][p12][t, z] == 
     -100*p12[t, z] + (1/2)*I*(-60000000*p13[t, z] + 6000000*Ep[t, z]*p32[t, z]), Derivative[1, 0][p13][t, z] == 
     -3000000*p13[t, z] - (1/2)*I*(60000000*p12[t, z] + 6000000*Ep[t, z]*(p11[t, z] - p33[t, z])), Derivative[1, 0][p21][t, z] == 
     -3000000*p21[t, z] - (1/2)*I*(6000000*Ep[t, z]*p23[t, z] - 60000000*p31[t, z]), Derivative[1, 0][p22][t, z] == -30000000*I*(p23[t, z] - p32[t, z]), 
    Derivative[1, 0][p23][t, z] == (-(1/2))*I*(6000000*Ep[t, z]*p21[t, z] + 60000000*p22[t, z] - 6000000*I*p23[t, z] - 60000000*p33[t, z]), 
    Derivative[1, 0][p31][t, z] == (1/2)*I*(60000000*p21[t, z] + 6000000*I*p31[t, z] + 6000000*Ep[t, z]*(p11[t, z] - p33[t, z])), 
    Derivative[1, 0][p32][t, z] == (1/2)*I*(6000000*Ep[t, z]*p12[t, z] + 60000000*p22[t, z] + 6000000*I*p32[t, z] - 60000000*p33[t, z]), 
    Derivative[1, 0][p33][t, z] == (1/2)*I*(60000000*p23[t, z] + 6000000*Ep[t, z]*(p13[t, z] - p31[t, z]) - 60000000*p32[t, z] + 12000000*I*p33[t, z]), 
    Derivative[0, 1][Ep][t, z] + Derivative[1, 0][Ep][t, z]/300000000 == (0. + 150.*I)*p12[t, z], p11[0, z] == 1, p12[0, z] == 0, p13[0, z] == 0, p21[0, z] == 0, p22[0, z] == 0, p23[0, z] == 0, 
    p31[0, z] == 0, p32[0, z] == 0, p33[0, z] == 0, Ep[t, 0] == -3.3546262790251186*^-8 + 0.0001/E^(500000000000*(-(1/250000) + t)^2), Ep[0, z] == 0}, 
   {p11[t, z], p12[t, z], p13[t, z], p21[t, z], p22[t, z], p23[t, z], p31[t, z], p32[t, z], p33[t, z], Ep[t, z]}, {z, 0, 0.1}, {t, 0, 15/10^6}, StartingStepSize -> 0.1]

Answer 1

Here "MaxStepSize" and "MinPoints" are sub-options of "TensorProductGrid" method. (MaxStepSize is also an option of NDSolve, but MaxStepSize in the warning isn't refer to it. ) These sub-options are discussed in the tutorial The Numerical Method of Lines. (Just search in the document, or check the online version. ) To eliminate the warning, you just need to add an option like

Method -> {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> 50, "MinPoints" -> 50}}

to NDSolve.

BTW, since the adjustion for these options is frequently needed when dealing with PDE solving, I myself have put the following function in my SystemOpen@"init.m" file:

mol[n : _Integer | {_Integer ..}, o_: "Pseudospectral"] := {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, "MinPoints" -> n, 
    "DifferenceOrder" -> o}}

With this function, we just need

Method -> mol[50, 4]

to achieve the same option adjustion. (The default setting of DifferenceOrder is 4. )

---

The setting above is tested in v9.0.1. At least in v12.0.1 and v12.1.1, one need an additional

MaxStepSize -> {5 10^-7, Automatic}

to eliminate the warning.

Here 5 10^-7 is the max step size setting in t direction, and Automatic in z direction. (It's set to Automatic because we've already controlled the grid size in z direction using sub-option of TensorProductGrid. ) Notice this MaxStepSize is the option for NDSolve.

This looks like a backslide in my view.

Answer 2

In addition to xzczd answer we can use next code been tested with versions 12, 12.1:

var = {p11, p12, p13, p21, p22, p23, p31, p32, p33, Ep};
sol = NDSolve[{Derivative[1, 0][p11][t, z] == -3000000*I*
      Ep[t, z]*(p13[t, z] - p31[t, z]), 
    Derivative[1, 0][p12][t, 
      z] == -100*p12[t, z] + (1/2)*
       I*(-60000000*p13[t, z] + 6000000*Ep[t, z]*p32[t, z]), 
    Derivative[1, 0][p13][t, 
      z] == -3000000*p13[t, z] - (1/2)*
       I*(60000000*p12[t, z] + 
         6000000*Ep[t, z]*(p11[t, z] - p33[t, z])), 
    Derivative[1, 0][p21][t, 
      z] == -3000000*p21[t, z] - (1/2)*
       I*(6000000*Ep[t, z]*p23[t, z] - 60000000*p31[t, z]), 
    Derivative[1, 0][p22][t, z] == -30000000*
      I*(p23[t, z] - p32[t, z]), 
    Derivative[1, 0][p23][t, z] == (-(1/2))*
      I*(6000000*Ep[t, z]*p21[t, z] + 60000000*p22[t, z] - 
        6000000*I*p23[t, z] - 60000000*p33[t, z]), 
    Derivative[1, 0][p31][t, z] == (1/2)*
      I*(60000000*p21[t, z] + 6000000*I*p31[t, z] + 
        6000000*Ep[t, z]*(p11[t, z] - p33[t, z])), 
    Derivative[1, 0][p32][t, z] == (1/2)*
      I*(6000000*Ep[t, z]*p12[t, z] + 60000000*p22[t, z] + 
        6000000*I*p32[t, z] - 60000000*p33[t, z]), 
    Derivative[1, 0][p33][t, z] == (1/2)*
      I*(60000000*p23[t, z] + 
        6000000*Ep[t, z]*(p13[t, z] - p31[t, z]) - 
        60000000*p32[t, z] + 12000000*I*p33[t, z]), 
    Derivative[0, 1][Ep][t, z] + 
      Derivative[1, 0][Ep][t, z]/300000000 == (0. + 150.*I)*p12[t, z],
     p11[0, z] == 1, p12[0, z] == 0, p13[0, z] == 0, p21[0, z] == 0, 
    p22[0, z] == 0, p23[0, z] == 0, p31[0, z] == 0, p32[0, z] == 0, 
    p33[0, z] == 0, 
    Ep[t, 0] == -3.3546262790251186*^-8 + 
      0.0001/E^(500000000000*(-(1/250000) + t)^2), Ep[0, z] == 0}, 
   var, {z, 0, 0.1}, {t, 0, 15/10^6}, 
   Method -> {"PDEDiscretization" -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"TensorProductGrid", 
         "MinPoints" -> 81, "MaxPoints" -> 81, 
         "DifferenceOrder" -> "Pseudospectral"}}}]; 

Visualization

Table[Plot3D[
  Evaluate[Abs[var[[i]][z, t] /. sol]], {z, 0, 0.1}, {t, 0, 15/10^6}, 
  ColorFunction -> Hue, Mesh -> None, PlotLabel -> var[[i]], 
  AxesLabel -> Automatic, PlotRange -> All, Boxed -> False, 
  PlotTheme -> "Scientific"], {i, Length[var]}]