Convergence Rate of NDSolve by Increasing the Spatial Grids

Question

I have a very simple PDE equation, with an analytical solution (exact solution). And I want to play with NDSolve and increasing the number of Spatial Grids.

Here is the exact solution:(Analytical Solution)

cA[x_, t_] := Erf[x/(2 Sqrt[t])]

Here is how I use the NDSolve:

solveUniformGrid[nPoint_, order_, xMin_, xMax_, pGoal_] := 
 Block[{xgrid, sol},
  xgrid = Join[Range[xMin, xMax, xMax/(nPoint - 1)], {xMax}];
  sol = NDSolve[
    {
     D[u[x, t], t] == D[u[x, t], x, x],
     u[xMax, t] == 1,
     u[x, tMin] == If[x == xMin, 0, 1],
     u[xMin, t] == 0
     }, u, {x, xMin, xMax}, {t, tMin, tMax},
    MaxSteps -> Infinity,
    InterpolationOrder -> Automatic,
    AccuracyGoal -> pGoal,
    PrecisionGoal -> pGoal,
    Method -> {
      "MethodOfLines", 
      "SpatialDiscretization" ->
       {
        "TensorProductGrid", 
         "DifferenceOrder" -> order,
        "Coordinates" -> {xgrid}
        },
      "DifferentiateBoundaryConditions" -> Automatic
      }];
  sol = First[u /. sol]
  ]

With xMin = tMin = 0, and xMax = 10; tMax = 1;

And then I calculate the Relative Error between the Analytical, and the Numerical solution: (The average of relative error in the domain of u[x,t]

calulateError[analytical_, numerical_, xMin_, xMax_, tMin_, tMax_, 
  nx_, nt_] :=
 Block[{xgrid, tgrid, errorList, percenterror},
  xgrid = Join[Range[xMin, xMax, xMax/(nx - 1)], {xMax}];
  tgrid = Join[Range[tMin, tMax, tMax/(nt - 1)], {tMax}];
  errorList = 
   Quiet[Abs[analytical[xgrid, #] - numerical[xgrid, #]]/
       analytical[xgrid, #] & /@ tgrid];
  errorList = 
   errorList /. {ComplexInfinity -> 0., Indeterminate -> 0.};
  percenterror = Mean[Flatten@errorList]
  ]

Now, I will play with the Grid Refinement, I will increase the number of spatial grid.

nPointList = {10, 20, 50, 100, 200};
solList1 = solveUniformGrid[#, 4, xMin, xMax, 8] & /@ nPointList;

And I plot the spatial error as a function of the number of grid points.

errorConcList1 = 
  calulateError[cA, #, xMin, xMax, tMin, tMax, 100, 50] & /@ solList1;

ListLogLogPlot[Transpose[{nPointList, errorConcList1}], 
 Joined -> True, Mesh -> All, Frame -> True, PlotRange -> All]

However, when I check the Convergence Rate, there is something I didn't understand. when I increase the number of grids, the Relative Error is stuck at a level. (The relative error is only 0.01).

In general, when we increase the grid points, the relative error will decrease further. Can someone explain this? Thank you

Update 01

It is the way to calculate the Relative Error between the analytical, and the numerical function that I didn't do it well. With the help of Henrik, by using the L2-Norm, I've got the correct convergence rate.

Update 02 - Improvement and Question about L2-Norm As you can see, the cA (analytical solution) is which is undefined at t == 0. So @Michael E2 has a very nice solution to compile and add the Ifcondition here.

Here is the 3 analytical solutions:, cA is the original analytical solution, cACompile is the compiled version by Michael E2, and cAImprove is just the non-compiled solution with an If condition to avoid the underfined problem at t==0.

cA[x_, t_] := Erf[x/(2 Sqrt[t])]

cACompile = Compile[{{xt, _Real, 1}},(*call:cA[{x,t}]*)
   Module[{x = First[xt], t = Last[xt]},
    If[x == 0,
     0.,
     If[t == 0,
      1.,
      Erf[x/(2 Sqrt[t])]
      ]]],
   RuntimeAttributes -> {Listable}, Parallelization -> True];

cAImprove[x_, t_] := If[x == 0,
  0.,
  If[t == 0,
   1.,
   Erf[x/(2 Sqrt[t])]
   ]]

I made a performance test on 1000 000 grid points

Thread[cA[Range[xMin, xMax, 0.00001], 0.5]]; // Timing
Thread[cAImprove[Range[xMin, xMax, 0.00001], 0.5]]; // Timing
cACompile /@ Thread[List[Range[xMin, xMax, 0.00001], 0.5]]; // Timing

And I get:

{0.03125, Null}

{0., Null}

{0.6875, Null}

My 1st question is:

Why the cAImprove with an If injected is faster than the original cA? It should be slower, right?

My 2nd question is:

How to obtain an L2-Norm between the two functions (exact, and approximation) for all of the domaine {xgrid,tgrid}?

Here is the 3 definitions of L2-Norm so far:

(* Integrate of L2-Norm on domain by Henrik Schumacher - Rather SLOW *)
globalIntegrateL2[anal_, num_] := 
 Divide[Sqrt[
   NIntegrate[
    Abs[anal[x, t] - num[x, t]]^2, {x, xMin, xMax}, {t, tMin, tMax}]],
   Sqrt[NIntegrate[
    Abs[anal[x, t]]^2, {x, xMin, xMax}, {t, tMin, tMax}]]]

(* L2 Norm on Grid by Michael E2 *)
traprule[yy_, xx_] := 
  Fold[#2.MovingAverage[#, 2] &, yy, Differences /@ xx];
globalGridL2[anal_, num_] := With[
  {
   xt = num@"Coordinates",
   exact = Apply[anal, num@"Grid", {2}],(*exact values on grid*)
   approx = num@"ValuesOnGrid"
   },(*computed solution on grid*)
  Divide @@ {traprule[(approx - exact)^2, xt] // Sqrt, 
    traprule[exact^2, xt] // Sqrt}]

(* L2 Norm on Grid by myself *)
globalGridL2Own[anal_, num_] := With[
  {
   exact = Apply[anal, num@"Grid", {2}],(*exact values on grid*)
   approx = num@"ValuesOnGrid"
   },
  Divide[Norm[approx - exact, 2], Norm[exact, 2]]
  ]

Here is the convergence rate of the 3 error functions. I don't know which one is correct.

Answer 1

I wish to add two things, separating the error of the discrete solution computed by NDSolve from the interpolating error between the interpolation grid, and speeding up the computation of the norm. The interpolation error is noticeable, but not significant in the $L^2$ norm, although it is probably the explanation of the OP's original difficulty with the point-wise-relative 1-norm.

One of the problems with speed is cA, which is undefined at t == 0. A different definition, especially with Compile speeds up computation.

ClearAll[cA];
cA = Compile[{{xt, _Real, 1}}, (* call: cA[{x, t}] *)
   Module[{x = First[xt], t = Last[xt]},
    If[x == 0,
     0.,
     If[t == 0,
      1.,
      Erf[x/(2 Sqrt[t])]
      ]]],
   RuntimeAttributes -> {Listable}, Parallelization -> True];

Some parameters. I memoized the solutions so I could play with them without recomputing them. It's unneeded but some of the rest of the code assumes calling sol[nx] won't be slow.

pGoal = 8;
xMin = tMin = 0;
xMax = 10;
tMax = 1;
order = 4;

nxList = {25, 100, 400, 1600, 6400, 25600}; (* discretization sequence *)

ClearAll[sol];
mem : sol[nx_] := With[{xgrid = Subdivide[N@xMin, xMax, nx]},
   mem = NDSolveValue[
     {D[u[x, t], t] == D[u[x, t], x, x],
      u[xMax, t] == 1, u[x, tMin] == If[x == xMin, 0, 1], 
      u[xMin, t] == 0},
     u, {x, xMin, xMax}, {t, tMin, tMax},
     MaxSteps -> Automatic, InterpolationOrder -> Automatic, 
     PrecisionGoal -> pGoal,
     Method -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"TensorProductGrid", 
         "DifferenceOrder" -> order, "Coordinates" -> {xgrid}}, 
       "DifferentiateBoundaryConditions" -> Automatic}]
   ];

Do[sol[nx], {nx, nxList}] (* optional: pre-compute & memoize solutions *)

We compute the integral for the $L^2$ norm from the solution grid used by NDSolve. The values we need, except for "Grid", are stored in the solution and can be obtained from the InterpolatingFunction; the "Grid" is computed efficiently from the "Coordinates". The value of the "Coordinates" has the form xx = {{x0, x1,..., xj}, {t0, t1,..., tk}}, that is, a list of the x-grid and t-grid.

traprule[yy_, xx_] := 
  Fold[#2.MovingAverage[#, 2] &, yy, Differences /@ xx];

Table[With[{
   xt = sol[nx]@"Coordinates",
   exact = cA@ sol[nx]@"Grid",        (* exact values on grid *)
   approx = sol[nx]@"ValuesOnGrid"},  (* computed solution on grid *)
  Divide @@ {
    traprule[(approx - exact)^2, xt] // Sqrt,
    traprule[exact^2, xt] // Sqrt
    }
  ],
 {nx, {25, 100, 400, 1600, 6400, 25600}}]
ListLogPlot[%, Joined -> True]

(*  {0.00202437, 0.000244795, 0.0000493161, 0.0000394941, 0.000039159, 0.0000393847}  *)

Answer 2

For a parabolic PDE as yours, the a priori error estimates are typically of the form $\mathrm{error} \sim (h^k + \tau^{k/2})$ for some $k \geq 0$ that depends on both the method and the norm in which you measure the error. Here $\tau >0$ is the time step size (i.e. total time divided by number of spatial grids) and $h>0$ is the maximal cell size in the spatial grid. The error goes to $0$ only if both $h \to 0$ and $\tau \to 0$, but you let go only $h \to 0$ while you keep the choice of $\tau$ to NDSolve by using MaxSteps -> Automatic. It it might not converge to $0$.

Edit:

After a some redious reverse engineering, I do not understand what the problem is. I get a nicely decaying sequence of relative $L^2$-errors from the following:

cA[x_, t_] := Erf[x/(2 Sqrt[t])]
pGoal = 8;
xMin = tMin = 0;
xMax = 10;
tMax = 1;
order = 4;

Table[
 xgrid = Subdivide[N@xMin, xMax, nx];
 sol = NDSolveValue[{
    D[u[x, t], t] == D[u[x, t], x, x],
    u[xMax, t] == 1,
    u[x, tMin] == If[x == xMin, 0, 1],
    u[xMin, t] == 0
    },
   u,
   {x, xMin, xMax}, {t, tMin, tMax},
   MaxSteps -> Automatic,
   InterpolationOrder -> Automatic,
   PrecisionGoal -> pGoal,
   Method -> {
     "MethodOfLines",
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "DifferenceOrder" -> order, "Coordinates" -> {xgrid}},
     "DifferentiateBoundaryConditions" -> Automatic}
   ];
 Divide[
  Sqrt[NIntegrate[Abs[sol[x, t] - cA[x, t]]^2, {x, xMin, xMax}, {t, tMin, tMax}]],
  Sqrt[NIntegrate[Abs[cA[x, t]]^2, {x, xMin, xMax}, {t, tMin, tMax}]]
  ],
 {nx, {100, 200, 400, 800}}]

{0.000645067, 0.000229846, 0.0000890372, 0.0000484648}