2
$\begingroup$

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];
    NDSolve`FiniteDifferenceDerivative[
      Derivative[i], 
      N[yyy], 
      DifferenceOrder -> "Pseudospectral"] @ "DifferentiationMatrix"
         // Normal // Developer`ToPackedArray // 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:

NDSolve`FiniteDifferenceDerivative::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.

$\endgroup$
  • 1
    $\begingroup$ Did you try the grid it suggested in the error message? $\endgroup$ – Michael E2 Nov 18 at 1:10
6
$\begingroup$

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

Subscript[Der, i_][yyy_] := Module[{xx}, xx = Length[yyy];
 NDSolve`FiniteDifferenceDerivative[Derivative[i], N[yyy], 
    DifferenceOrder -> "Pseudospectral"]@"DifferentiationMatrix" //
   Normal // Developer`ToPackedArray // 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]    

enter image description here

$\endgroup$
  • 1
    $\begingroup$ 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}]. $\endgroup$ – J. M. will be back soon Nov 18 at 10:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.