# Can Mathematica provide a reliable estimate of the numerical error from NDSolve?

In the Details section of the Mathematica documentation for PrecisionGoal, one is told that

Even though you may specify PrecisionGoal -> n, the results you get may sometimes have much less than n‐digit precision.

and that

With PrecisionGoal -> p and AccuracyGoal -> a, the Wolfram Language attempts to make the numerical error in a result of size $$x$$ be less than $$10^{-a}+10^{-p} \left| x\right|$$.

According to these, after setting, for example, WorkingPrecision -> 100, PrecisionGoal -> 23, AccuracyGoal -> 50, then NDSolve may give an output with a precision of, say, $$19$$ significant digits.

How can one fix this? If not, how can one be aware of the fact that Mathematica output is less precise than WorkingPrecision -> 100, PrecisionGoal -> 23, AccuracyGoal -> 50? Can Mathematica be induced to give some estimate of the resulting numerical error?

• You can calculate the error by substituting a numerical solution in the equation on the grid. – Alex Trounev Nov 29 '19 at 13:27
• @Alex, that only works if the resulting interpolating function has high enough order for the PDE at hand. A better way to do this is to use manufactured solutions. – user21 Nov 29 '19 at 15:43
• @user21 OK! Why not? Use what you need. In any case, to determine the error, you need to use the equation and solution. Or are there other options?:) – Alex Trounev Nov 29 '19 at 15:55
• @AlexTrounev, if you want I can send you a detailed explanation I have that will be published in V12.1. ruebenkoATwolfram.com. If not I can post a link here after 12.1 is released. – user21 Nov 29 '19 at 16:00
• @AlexTrounev There's some discussion of it here in an answer to another of the OP's questions. – Michael E2 Nov 29 '19 at 16:27

Options[MonitorMethod] = {Method -> Automatic,
"MonitorFunction" ->
Function[{h, state, meth},
Print[{"H" -> h, "SD" -> state@"SolutionData"}]]};
MonitorMethod /:
NDSolveInitializeMethod[MonitorMethod, stepmode_, sd_, rhs_,
state_, OptionsPattern[MonitorMethod]] := Module[{submethod, mf},
mf = OptionValue["MonitorFunction"];
submethod = OptionValue[Method];
submethod =
NDSolveInitializeSubmethod[MonitorMethod, submethod, stepmode,
sd, rhs, state];
MonitorMethod[submethod, mf]];
MonitorMethod[submethod_, mf_]["Step"[f_, h_, sd_, state_]] :=
Module[{res},
res = NDSolveInvokeMethod[submethod, f, h, sd, state];
Return[\$Failed]]; (* submethod not valid for monitoring *)
mf[h, state, submethod];
If[SameQ[res[[-1]], submethod], res[[-1]] = Null,
res[[-1]] = MonitorMethod[res[[-1]], mf]];
res];
MonitorMethod[___]["StepInput"] = {"Function"[All], "H",
"SolutionData", "StateData"};
MonitorMethod[___]["StepOutput"] = {"H", "SD", "MethodObject"};
MonitorMethod[submethod_, ___][prop_] := submethod[prop];


If the Method implements the "StepError" method, it will return the (scaled) step error estimate. (The only way to know the true error is to know the true solution and compare.) By "scaled," Mathematica means $$\text{scaled error} = {|\text{error}| \over 10^{-\text{ag}} + 10^{-\text{pg}} | x |} \,,$$ which will be between 0 and 1 when the $$\text{error}$$ satisfies the AccuracyGoal $$\text{ag}$$ and the PrecisionGoal $$\text{pg}$$.

The MonitorMethod takes a "MonitorFunction" option, which should be a function of the form

Function[{h, state, meth}, <...body...>]


where h is the step size, state is the NDSolveStateData object, and meth is the Method object of the submethod.

Example use:

{sol, {errdata}} = Reap[
NDSolveValue[{x''[t] + x[t] == 0, x[0] == 1, x'[0] == 1},
x, {t, 0, 2},
Method -> {MonitorMethod,
"MonitorFunction" ->
Function[{h, state, meth},
Sow[meth@"StepError", "ScaledStepError"]]},
MaxStepFraction -> 1, WorkingPrecision -> 100,
PrecisionGoal -> 23, AccuracyGoal -> 50],
"ScaledStepError"];

GraphicsRow[{
ListLinePlot[Transpose@{Flatten@Rest@sol@"Grid", errdata},
Mesh -> All, PlotRange -> {0, 1}, PlotRangePadding -> Scaled[.05],
PlotLabel -> "Scaled error estimate"],
Show[
Plot[Sin[t] + Cos[t], {t, 0, 2}, PlotStyle -> Red],
ListLinePlot[sol, Mesh -> All],
PlotLabel -> "Steps on top of exact solution"]
}]


In our example, we know the exact solution, so we can check the actual error:

Block[{t = Flatten@sol@"Grid", data},
data = Transpose@{t, (Sin[t] + Cos[t] - sol[t])/(
10^-50 + 10^-23 Abs[Sin[t] + Cos[t]])};
ListLinePlot[data,
Epilog -> {PointSize@Medium, Tooltip[Point[#], N@Last@#] & /@ data},
PlotRange -> All, PlotLabel -> "Actual scaled error"]
]


Of course, when I'm this interested in the error, it's usually because I have reasons for wondering whether the error estimates, which are based on discrete approximations to a function assumed to have a certain smoothness, are unreliable.

• (+1), funny just a few days we internally talked about this and that it's not well known and now it pops up. Great answer. – user21 Nov 29 '19 at 15:46
• @user21 Thanks. Of course, I've rooted around NDSolve more than the average user. I noticed that MonitorMethod does not work with the submethods LSODA and Adams (maybe others?). Aside from wondering why, I thought it would be nice if the restrictions were documented. – Michael E2 Nov 29 '19 at 15:49
• It absolutely should be documented or the restrictions removed, unfortunately sw spreads us out too thinly over the bread to keep on top of all these... I'll file it as a bug/suggestion but don't hold your breath! Sorry about that. – user21 Nov 29 '19 at 15:53
• I filed it. We will see.... – user21 Nov 29 '19 at 15:56
• @user21 Thank you very much. Kind of too bad about the LSODA restriction, since, being the usual Automatic method, it is the method one would most often want to monitor. – Michael E2 Nov 29 '19 at 16:21

I'll add more details when I have my copy of Wagon's "Mathematica in Action" again, but as I mentioned in a comment, one possible way to check your solution would be to "integrate backwards", with initial conditions taken from your first call to NDSolve[]. Using Michael's example:

sol1 = NDSolveValue[{x''[t] + x[t] == 0, x[0] == 1, x'[0] == 1}, x, {t, 0, 2},
AccuracyGoal -> 50, MaxStepFraction -> 1,
PrecisionGoal -> 23, WorkingPrecision -> 100];


Here's the "backwards" version:

sol2 = NDSolveValue[{x''[t] + x[t] == 0, x[2] == sol1[2], x'[2] == sol1'[2]}, x, {t, 2, 0},
AccuracyGoal -> 50, MaxStepFraction -> 1, PrecisionGoal -> 23,
WorkingPrecision -> 95];


(note that I had to shrink the WorkingPrecision slightly)

Plot them together:

Plot[sol1[t] - sol2[t], {t, 0, 2}, PlotRange -> All, WorkingPrecision -> 60]


and we see that there is a discrepancy of at most $$2.5\times 10^{-15}$$.

Another possibility is to leverage ParametricNDSolve[] to see what happens if the initial conditions are perturbed. Adapting one of the examples from the docs to Michael's example:

pf = ParametricNDSolveValue[{x''[t] + x[t] == 0, x[0] == 1 + h, x'[0] == 1 + h}, x,
{t, 0, 2}, {h}, AccuracyGoal -> 50, MaxStepFraction -> 1,
PrecisionGoal -> 23, WorkingPrecision -> 100];

With[{h = 1*^-6},
Plot[Evaluate[pf[0][t] + {-h, 0, h} pf'[0][t]], {t, 0, 5},
PlotStyle -> {Gray, {ColorData[97, 1]}, Gray},
Filling -> {1 -> {3}}, WorkingPrecision -> 60]]
`

A very sensitive ODE would display a much wider band.

• (I'll prolly edit this answer later to put in a more interesting example, like Duffing's equation.) – J. M.'s technical difficulties Nov 29 '19 at 21:56