Testing the accuracy of numerically computed derivatives

Question

I am calculating approximate derivatives by using NDSolve`FiniteDifferenceDerivative, so this works:

 Subscript[Der, i_][yyy_] := 
   Module[{xx}, 
     xx = Length[yyy];
     NDSolve`FiniteDifferenceDerivative[
       Derivative[i], 
       N[yyy], 
       DifferenceOrder -> 2] @ "DifferentiationMatrix"
         // Normal // Developer`ToPackedArray // SparseArray];

xi = 1.;
xf = -1;
yy = 100; 
xgrid = Table[xi + i (xf - xi/yy), {i, 0, yy}];

(Der1 = Subscript[Der, 1][xgrid]) // MatrixForm;

numerical = Der1.Exp[-xgrid^2];
exact = -2*xgrid*Exp[-xgrid^2];

diff = numerical - exact; 
diffError = yy^2*diff

ListLinePlot[yy^2 Abs[diff]]

I want to show my solution is accurate by demonstrating that the difference between the numerical solution and the exact solution goes to zero as $\mathtt{yy}^{-2}$. For this I want to plot $\mathtt{yy}^2 |\mathrm{numerical} - \mathrm{exact}|$ for different values of $\mathtt{yy}$ but am not sure how to do this.

The code gives reasonable values for the differences, though I am not sure how to plot them for different $\mathtt{yy}$ values.

I obtained the follow plot from the code shown above.

Answer 1

xi = -1.;
xf = 1;
xgrid[yy_] := Range[xi, xf, (xf - xi)/yy]

Subscript[Der, i_][yyy_] := Module[{xx}, xx = Length[yyy];
 NDSolve`FiniteDifferenceDerivative[Derivative[i], N[yyy], 
    DifferenceOrder -> 2]@"DifferentiationMatrix" // Normal // 
 Developer`ToPackedArray // SparseArray];

Der1[yy_] := Subscript[Der, 1][xgrid[yy]]
numerical[yy_] := Der1[yy].Exp[-xgrid[yy]^2]
exact[yy_] := -2*xgrid[yy]*Exp[-xgrid[yy]^2]
diff[yy_] := numerical[yy] - exact[yy]

yyvals = {100, 300, 1000};

ListLinePlot[
 Table[Transpose[{xgrid[yy], yy^2 Abs[diff[yy]]}], {yy, yyvals}], 
 PlotRange -> All, PlotLegends -> StringTemplate["yy = ``"] /@ yyvals]

Max[diff[100]] / Max[diff[1000]] = 99.9756

This means the error ~ 1/yy^2. For better see this scaling low one can use logarithmic scale:

ListLinePlot[
 Table[Transpose[{xgrid[yy], Abs[diff[yy]]}], {yy, yyvals}], 
 PlotRange -> All, PlotLegends -> StringTemplate["yy = ``"] /@ yyvals,
 ScalingFunctions -> "Log", Frame -> True]

NonlinearModelFit[Table[{yy, Max[diff[yy]]}, {yy, 100, 10000, 500}], 
 a + b/x^2, {a, b}, x]

Answer 2

GeneralUtilities`BenchmarkPlot will fit data (by design, timings) to a number of models. Here's a hack to show that 1/error ~ yy^2 is the best model of the ones it tests.

der[i_][yyy_] := Module[{xx}, xx = Length[yyy];
   NDSolve`FiniteDifferenceDerivative[Derivative[i], N[yyy], 
        DifferenceOrder -> 2]@"DifferentiationMatrix" // Normal // 
     Developer`ToPackedArray // SparseArray];

ClearAll[err];
err[p_: Infinity][yy_] := Module[
   {xgrid, Der1, numerical, exact},
   xgrid = Subdivide[-1., 1., yy];
   (Der1 = der[1][xgrid]);
   numerical = Der1.Exp[-xgrid^2];
   exact = -2*xgrid*Exp[-xgrid^2];
   Norm[numerical - exact, p]/yy^(1/p)
   ];

(* error measured by the scaled infinity and two norms *)
Needs@"GeneralUtilities`";
Block[{RepeatedTiming = List}, (* hack timing to be the reciprocal error *)
 BenchmarkPlot[{1/err[Infinity][#] &, 1/err[2][#] &}, # &, 2^Range@12,
   "IncludeFits" -> True]
 ]