# Making NDSolve choosing by itself a spatial meshing leading to correct integrated results

In a previous post, I was asking how to force a fine spatial meshing in NDSolve using the FiniteElement method. The solution was

NDSolveValue[{op == 0, ics}, u, {t, 0, 2}, {x} \[Element] reg,
Method -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> 0.005}}]


giving the correct result below

instead of the incorrect one obtained (without any error or warning message) when using simply

NDSolveValue[{op == 0, ics}, u, {t, 0, 2}, {x} \[Element] reg]


But this solution is not satisfactory because in the general case, (1) it is not so obvious to detect that a solution is False and (2) it is long and boring to run two subsequent NDSolve integrations with two different meshings in order to check whether it gives twice the same result, and (3) the user can be tempted to use a too fine meshing and consequently too much time in order to be sure to get a correct result. One could expect NDSolve to have an optional algorithm that automatically refine the spatial meshing if needed.

1) Does such an option exist?

2) Is it compatible with the FiniteElement method?

Denis

• “...it is long and boring to run...” is one of the reasons society has invested in the development of computers and programming. What I’m getting at is that one could write a program to carry out adaptive refinement. I think ToElementMesh[] probably has the tools/features needed. It would probably not be slow on 1D problems (ODEs) but higher dimensional PDEs might be more challenging. I’ve not done it nor do I know if there are error estimation schemes other than what you describe above. Feb 14 '19 at 13:28
• It is necessary to publish all the code, not two pictures. Feb 14 '19 at 15:12
• Here's an example for the method of lines, which is easier than what I had in mind for FEM since a message is emitted when further reduction of step size is needed. It uses recursion instead of a loop, but I assumed it would converge after just a few iterations. mathematica.stackexchange.com/a/73848/4999 Feb 15 '19 at 3:58
• Sorry Alex, the code is in the previous post linked at the top. I did not want to repeat it here because the question is different. Feb 15 '19 at 13:08
• Thank you Michael. I'll see what I can do along this line, but note that in my example, no error or warning is generated by NDSolve when the integration happens to be completely wrong Feb 15 '19 at 13:24

No, adaptive mesh refinement is not implemented in version 11.3 for PDE solving with the FEM. NIntegrate with the FEM uses an adaptive mesh refinement.

Here is a brute-force solution consisting in integrating the problem with stepsizes divided by 2, 4, 8, 16... until the absolute difference between two consecutive solutions goes below a tolerance threshold:

reg = Line[{{0}, {1}}];
shape = Cos[16 π (x - 0.5)] D[0.125 Erf[(x - 0.5)/0.125], x];
op = D[u[t, x], {x, 2}] - D[u[t, x], {t, 2}];
ics = {u[0, x] == shape, Derivative[1, 0][u][0, x] == 0};
sol[maxCellMeasure_] := NDSolveValue[{op == 0, ics}, u, {t, 0, 2}, {x} \[Element] reg,
Method -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> maxCellMeasure_}}];

decreaseStep[{maxCellMeasure_, so_}, print_: False] :=
Block[{},
If[print, Print["maxCellMeasure=" <> ToString[maxCellMeasure/2]]];
{maxCellMeasure/2, sol[maxCellMeasure/2]}]

compare[sol1_, sol2_] :=
Block[{points = sol1[Grid][[-1]]},
( (sol2[Sequence @@ #] & /@ points)
- (sol1[Sequence @@ #] & /@ points)) // Abs // Max]

solu[initialMaxCellMeasure_, tol_, maxSteps_: 10, print_: False] :=
NestWhileList[decreaseStep[#, print] &,
{initialMaxCellMeasure, sol[initialMaxCellMeasure]},
(compare[#1[[2]], #2[[2]]] > tol)&,
2, maxSteps][[-1]]
solu0 = solu[0.05, 0.001, 10, True]


Any more elegant solution would be appreciated.

Denis