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 If
condition 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.
xMin
,xMax
,tMin
,tMax
, which prevents people from understanding your problem or giving concrete advice. $\endgroup$ – Michael E2 Oct 15 '20 at 15:05