# Convergence Rate of NDSolve by Increasing the Spatial Grids

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.

• You're missing values for the parameters xMin, xMax, tMin, tMax, which prevents people from understanding your problem or giving concrete advice. – Michael E2 Oct 15 '20 at 15:05
• I've added in the question. – Nam Nguyen Oct 15 '20 at 15:59

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}  *)

• Wonderful idea of the Compiled function of cA. I really love it. I made a test between the cA, and the cACompile function. Thread[cA[Range[xMin, xMax, 0.00001], 0.5]]; // Timing and fhis, cACompile /@ Thread[List[Range[xMin, xMax, 0.00001], 0.5]]; // Timing. And I get {0.03125, Null} for the cA, and {0.703125, Null} for the Compiled version. So why the compiled version is slower? Thank you very much. – Nam Nguyen Oct 19 '20 at 9:04
• @NamNguyen The issue with the original problem is at t == 0. It prevents the computation for errorList to be done as a packed array. The compiled function gets around that with using If, and If is much faster inside Compile than in general. OTOH, if you avoid t == 0 as in your comment, Mma uses the vectorized Erf[]... – Michael E2 Oct 19 '20 at 13:01
• Vectorized functions in the MKL will generally be faster on arrays (on Intel machines, at least) than anything else. Note that the call to cA in Thread[cA[Range[xMin, xMax, 0.00001], 0.5]] results in a vectorized call to the arithmetic functions as well as Erf[]. Range produces a packed array. And the Thread does nothing in this case. -- In your second, compiled example in the comment, there are three things I think slow things down. First is Thread, which does a somewhat expensive memory operation. Second are the two If tests in caCompiled. Third is that it's not vectorized. – Michael E2 Oct 19 '20 at 13:10
• Thank you very much. I've another question on L2-norm. Please read the updated question. – Nam Nguyen Oct 19 '20 at 15:08

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}

• I'm increasing the number of the spatial grid, so the cell size h->0, right?. The temporal t grid is controlled by ODA (adaptive) based on the Accuracy Goal. So my question is, how to make the error go to zero. So as your opinion, I will try to decrease the temporal step size too, is it right? – Nam Nguyen Oct 15 '20 at 11:17
• So then, do you refine the time grid, too? After all, you use MaxSteps -> Automatic and and fixed values for AccuracyGoal  and PrecisionGoal . That looks like "no". – Henrik Schumacher Oct 15 '20 at 11:19
• I will update the question, by refining both space-grid, and time-grid. – Nam Nguyen Oct 15 '20 at 11:21
• Oh, and I see that you compute the relative error pointwise as $||(u - v)/v||_{L^1}$. So there is probaly also quite some catastrophic cancellation going on... Better try to use something like $||u - v||/||u||$ with some suitable norm (e.g., the $L^2$-norm). – Henrik Schumacher Oct 15 '20 at 11:41
• @NamNguyen I feel that the discontinuity in the boundary condition, u[x, tMin] == If[x == xMin, 0, 1], probably limits the convergence, but I'm not very familiar with PDEs. Related: mathematica.stackexchange.com/questions/129426/… – Michael E2 Oct 15 '20 at 15:05