1
$\begingroup$

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.

accuracytest

$\endgroup$
2
$\begingroup$
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]

enter image description here

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]

enter image description here

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

Blockquote

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank for your answer, please note it should be plotting in the vertical axis yy^2 and note yy . $\endgroup$ – RiPHX Nov 17 '19 at 16:37
  • $\begingroup$ Yes, I missed that in plot, please see my edited answer. $\endgroup$ – Alx Nov 17 '19 at 23:23
2
$\begingroup$

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

enter image description here

| improve this answer | |
$\endgroup$

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.