# Testing the accuracy of numerically computed derivatives

I am calculating approximate derivatives by using NDSolveFiniteDifferenceDerivative, so this works:

 Subscript[Der, i_][yyy_] :=
Module[{xx},
xx = Length[yyy];
NDSolveFiniteDifferenceDerivative[
Derivative[i],
N[yyy],
DifferenceOrder -> 2] @ "DifferentiationMatrix"
// Normal // DeveloperToPackedArray // 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.

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

Subscript[Der, i_][yyy_] := Module[{xx}, xx = Length[yyy];
NDSolveFiniteDifferenceDerivative[Derivative[i], N[yyy],
DifferenceOrder -> 2]@"DifferentiationMatrix" // Normal //
DeveloperToPackedArray // 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]


• Thank for your answer, please note it should be plotting in the vertical axis yy^2 and note yy . Commented Nov 17, 2019 at 16:37
• Yes, I missed that in plot, please see my edited answer.
– Alx
Commented Nov 17, 2019 at 23:23

GeneralUtilitiesBenchmarkPlot 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];
NDSolveFiniteDifferenceDerivative[Derivative[i], N[yyy],
DifferenceOrder -> 2]@"DifferentiationMatrix" // Normal //
DeveloperToPackedArray // 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]
]
`