Strange chattering in derivative of InterpolatingFunction returned by NDSolve using FEM

When using the Finite Element Method of NDSolve to solve a set of ODE, the resulting InterpolatingFunction showcases a strange chattering in its derivatives. This behavior does not occur when using the default method or Runge-Kutta options.

Let me illustrate the issue with a very simple problem.

DEFAULT CASE

A reference solution can be obtained with the following code

w2 = 6;
m = 2;
T = 2.0;
sol = NDSolveValue[{q'[t] == ζ[t], ζ'[t] + w2*Sin[q[t]] == 0,
q[0] == Pi/3., ζ[0] == 0}, {q, ζ}, {t, 0, T}]


Solutions are smooth. One can easily plot the time derivative of the first solution with

Plot[sol[[1]]'[t], {t, 0, T}]


FEM CASE

The same solutions can be obtained by using the Finite Element Method option of NDSolve:

solFEM = NDSolveValue[{q'[t] == ζ[t], ζ'[t] + w2*Sin[q[t]] == 0,
DirichletCondition[{q[t] == Pi/3., ζ[t] == 0}, t == 0]},
{q, ζ}, t ∈ Line[{{0}, {T}}], Method -> {"FiniteElement"}]


Of course, the solutions are also smooth in this case and the time derivative of the first solution can be plotted with

Plot[solFEM[[1]]'[t], {t, 0, T}]


We can compare both plots and see that the one coming from the Finite Element Method has some chattering to it:

This is of course just one basic example and the issue shown here tends to be more troublesome for more complicated situations.

So, why does this happen? And since I really need the FEM option for a specific scenario, how can this be avoided or treated?

EDIT 1

On appearance, an easy fix would have been to reduce the cell size of the FEM method. This can be achieved by setting the "MaxCellMeasure" size to something smaller (like 0.01 or 0.001 for example):

 solFEM =
NDSolveValue[{q'[t] == ζ[t], ζ'[t] + w2*Sin[q[t]] == 0,
DirichletCondition[{q[t] == Pi/3., ζ[t] == 0}, t == 0]},
{q, ζ}, t ∈ Line[{{0}, {T}}],
Method -> {"FiniteElement", "MeshOptions" -> {MaxCellMeasure -> 0.01}}]


However, this only "hides" the issue to the bare eye. One can still see the chattering effect when zooming around the minimum for example:

Why would this matter?

My problem with this issue has nothing to do with the curve appearance. On a practical level, if I want to determine, say, the local minimum (or maximum) near some region, then FindMinimum (or FindMaximum) will inevitably fail. I can easily find the local minimum near t->0.5 for the "default" solution:

FindMinimum[sol[[1]]'[t], {t, 0.5}]
(* {-2.44949, {t -> 0.688223}} *)


However, this is not possible with the FEM solution because chattering obviously affects FindMinimum:

FindMinimum[solFEM[[1]]'[t], {t, 0.5}]


••• FindMinimum: The line search decreased the step size to within the tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the function. You may need more than MachinePrecision digits of working precision to meet these tolerances.

(* {-2.1981, {t -> 0.5}} *)


Because of chattering, this simple operation becomes impossible even when directly feeding a very close initial guess:

FindMinimum[solFEM[[1]]'[t], {t, 0.688223}]


••• FindMinimum: Line search unable to find a sufficient decrease in the function value with MachinePrecision digit precision.

(* {-2.45211, {t -> 0.69}} *)


This is only one example of the consequences of that chattering effect.

I'm not sure about how to explain this behavior, but it can be circumvented by extracting the original data from the InterpolatingFunction[…], calculating the derivative based on finite difference method (FDM) and recreate an InterpolatingFunction[…]:

createD = With[{values = #["ValuesOnGrid"], grid = #["Grid"] // Flatten},
With[{order = Ordering@grid},
ListInterpolation[
NDSolveFiniteDifferenceDerivative[1, grid[[order]],
values[[order]]], {grid[[order]]}]]] &;
dq = createD@solFEM[[1]];

Plot[{sol[[1]]'[t], dq[t]} // Evaluate, {t, 0, T}, PlotStyle -> {Automatic, Dashed}]


FindMinimum[dq[t], {t, 0.5}]
(* {-2.44878, {t -> 0.682447}} *)

• This would have been too easy and did try that solution before. However, when you zoom-in around the minimum with Plot[solFEM[[1]]'[t], {t, 0.6, 0.8}], you can easily see that the chattering is still there, only less perceptible with the bare eye... Nov 18, 2022 at 15:57
• @Meclassic I rewrite my answer, have a look :) . Nov 18, 2022 at 17:05
• Thanks, this seems to work, however I am starting to believe that the chattering does come from the interpolation performed by the FEM itself. It would be interesting to know if that behavior can be modified Nov 18, 2022 at 19:45
• @Meclassic Not that simple. Just remove the NDSolveFiniteDifferenceDerivative in my code to rebuild a InterpolatingFunction for q using the traditional "Hermit" method and plot the derivative, you'll see the chattering is still there. Setting Method->"Spline" will remove the discontinuity, but this time the derivative is waving. Nov 19, 2022 at 2:03