Testing the accuracy of numerically computed derivatives - alternative method

So I have asked this question earlier today: Testing the accuracy of numerically computed derivatives,

This method works well for DifferenceOrder methods but upon reading https://reference.wolfram.com/language/tutorial/NDSolveMethodOfLines.html, I believe it would be more accurate to use DifferenceOrder -> "Pseudospectral" with an appropriate grid.

Implementing the code as posted in my question but with the following changes:

Subscript[Der, i_][yyy_] :=
Module[{xx},
xx = Length[yyy];
NDSolveFiniteDifferenceDerivative[
Derivative[i],
N[yyy],
DifferenceOrder -> "Pseudospectral"] @ "DifferentiationMatrix"
// Normal // DeveloperToPackedArray // SparseArray];

xi = 1.;
xf = -1;
yy = 100;
xgrid = Table[(a+b)/2 + (a-b)/2)*Cos[i \pi/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 = Exp[yy]*diff


I get a very accurate result for the difference though,when I try to plot it, by applying the method prompted in the answer to my previous question, I get an error message:

NDSolveFiniteDifferenceDerivative::spc: The coordinates {-0.498973, -0.495895,-0.49078, -0.483647, -0.474528, -0.463458, -0.450484, -0.435659, -0.419044, -0.400707, <<29>>, 0.419044, 0.435659, 0.450484, 0.463458, 0.474528, 0.483647, 0.49078, 0.495895, 0.498973, 0.5} for dimension 1 do not have spacing proportional to the Chebyshev points -Cos[Range[0,48] Pi/48]. Nonperiodic pseudospectral derivatives are only computed for this spacing.

It seems to not like my xgrid choice though for when I implement my code the results are really great for the derivative, so it would be great to show this for different yy values. Thanks in advance.

• Did you try the grid it suggested in the error message? Nov 18 '19 at 1:10

As mentioned in referenced help page you need to construct a special grid:

Subscript[Der, i_][yyy_] := Module[{xx}, xx = Length[yyy];
NDSolveFiniteDifferenceDerivative[Derivative[i], N[yyy],
DifferenceOrder -> "Pseudospectral"]@"DifferentiationMatrix" //
Normal // DeveloperToPackedArray // SparseArray];

xi = -1.;
xf = 1;
xgrid[yy_] := xi + (xf - xi)/2 (1 - Cos[π Range[0, yy - 1]/(yy - 1)]);

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];


And now plotting:

yyvals = {100, 300, 1000, 2000};

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


• For stability at large values of yy, I would instead do xgrid[yy_] := Rescale[Haversine[π Range[0, yy - 1]/(yy - 1)], {0, 1}, {xi, xf}]`. Nov 18 '19 at 10:25