# Error control for NDSolve

I have a problem controlling the numerical error associated with the following non-linear ODE :

sol=NDSolve[{(R'[t])^2 + 2 R[t] R''[t] == -1, R[1]==1,R'[1]==2/3},R,{t,1,3}]


Once you solve this equation using the above code and then compute

R'[t]^2 + 2 R[t] R''[t] +1


you do not get exactly zero. If you plot this expression, you will see a highly oscillatory function, which I guess is an indication of poor error control.

I have tried increasing accuracy goal, precision goal, reducing the maximum step size, toying around with RungeKutta, but so far haven't had any luck.

• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory Tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. – bbgodfrey Jul 22 '15 at 1:33
• There are also the options StartingStepSize, WorkingPrecision to play with. The problem seems to be at the start. It usually does not help much to increase PrecisionGoal and AccuracyGoal unless you also increase WorkingPrecision. What error would be acceptable? Since NDSolve computes discrete steps, you should expect the error oscillate between the steps (but only a small amount). Having zero error throughout is not really an achievable goal. – Michael E2 Jul 22 '15 at 1:43
• Hi Mike. My goal is to keep error smaller than 10^-30 everywhere without making the computation painfully time consuming – Sina Jul 22 '15 at 18:46

The OP's -- oops, they're bbodfrey's -- pictures suggest the problem is with interpolation, as bbgodfrey also observed. Some of the problem can be ameliorated with the InterpolationOrder option.

In functions such as NDSolve, InterpolationOrder->All specifies that the interpolation order should be chosen to be the same as the order of the underlying solution method.

Another problem is that NDSolve seems to have trouble at the beginning of the integration. This sometimes can be dealt with the StartingStepSize option. It can also be dealt with by differentiating the differential equation.

To get a precision of 30 digits, WorkingPrecision usually needs to be at least twice as large or 60, and PrecisionGoal and AccuracyGoal need to be at least 30 -- a little higher is usually needed get strictly at least 30 digits. By default, these are set to half of WorkingPrecision. In some cases, the MaxStepSize option is needed to keep NDSolve from getting too ambitious and committing an "NDSolve::nderr: Error test failure".

The following three approaches give results with a residual error of less than 10^-30 throughout the interval of integration.

{rsol} = NDSolve[
{(R'[t])^2 + 2 R[t] R''[t] == -1, R[1] == 1, R'[1] == 2/3},
R, {t, 1, 3},
Method -> "ExplicitRungeKutta", WorkingPrecision -> 61,
InterpolationOrder -> All];

diffeq = (R'[t])^2 + 2 R[t] R''[t] == -1;
{rsol} = NDSolve[
{D[diffeq, t],
{R[1], R'[1], R''[1]} ==
({R[t], R'[t], R''[t]} /. First@Solve[diffeq, R''[t]] /. {R[t] -> 1, R'[t] -> 2/3})},
R, {t, 1, 3},
MaxStepSize -> 3*^-4, WorkingPrecision -> 62,
InterpolationOrder -> All];

{rsol} = NDSolve[
{(R'[t])^2 + 2 R[t] R''[t] == -1, R[1] == 1, R'[1] == 2/3},
R, {t, 1, 3},
StartingStepSize -> 1*^-8, MaxStepSize -> 1*^-4,
PrecisionGoal -> 33, AccuracyGoal -> 33, WorkingPrecision -> 70,
MaxSteps -> 2*^5, InterpolationOrder -> All];


The first one ("ExplicitRungeKutta") produces the smallest (memory use) solution (2.5s, 3.4MB). The second (differentiating the ODE) is fastest with excellent precision control (1.5s, 43MB). And the third, well, it works, and after it gets started is the most accurate (3.4s, 89MB).

Plot[
(R'[t])^2 + 2 R[t] R''[t] + 1 /. rsol // RealExponent // Evaluate,
{t, 1, 3}, PlotRange -> All, WorkingPrecision -> 70]


The Runge-Kutta solution.

Differentiating the differential equation.

Controlled step size.

You might notice that the largest error is near the beginning (in all cases). The interesting stuff around t == 2.5 is where the first derivative is zero.

The problem appears to be associated with interpolation by the InterpolatingFunction produced by NDSolve rather than by NDSolve itself. For instance, with

sol = NDSolve[{(R'[t])^2 + 2 R[t] R''[t] == -1, R[1] == 1,
R'[1] == 2/3}, {R}, {t, 1, 1.2}, WorkingPrecision -> 100];


a Plot of R and its first two derivatives yields near t = 0

Plot[{R[t], R'[t], R''[t]} /. sol, {t, 1, 1.00001}]


shows that R''[t] is ill behaved there. Magnifying the R-axis scale gives

Yet, plotting the data points in InterpolatingFunction indicates that they are well behaved there.

ListPlot[Transpose[{sol[[1, 1, 2]][[3, 1, 1 ;; 118]], sol[[1, 1, 2]][[4, 1 ;; 118, 3]]}]]


In a comment below, Guess who it is suggested specifying the Method. It turns out that Method -> "ExplicitRungeKutta"works well, although the other three principal methods offer little improvement.

sol = NDSolve[{(R'[t])^2 + 2 R[t] R''[t] == -1, R[1] == 1, R'[1] == 2/3},
{R}, {t, 1, 3}, WorkingPrecision -> 30, Method -> "ExplicitRungeKutta"];
Plot[Evaluate[{R[t], R'[t], R''[t]} /. sol], {t, 1, 3}]


The residual error is small.

Plot[(R'[t]^2 + 2 R[t] R''[t] + 1) /. sol, {t, 1, 3}, PlotRange -> All]


and can be made smaller yet by increasing WorkingPrecision.

• What happens if you use a different Method, like "ExplicitRungeKutta" or "Adams"? – J. M.'s technical difficulties Jul 22 '15 at 4:14
• @Guesswhoitis. Thanks for the suggestion. "ExplicitRungeKutta" works very well, although "Adams" and the other methods offer little improvement. – bbgodfrey Jul 22 '15 at 4:58
• Ah, thought so. NDSolve[] uses multistep methods by default (Method -> "LSODA"), so if "Adams" didn't work, then it confirms multistep methods cannot deal with this ODE. Since Runge-Kutta worked, I believe Bulirsch-Stoer (Method -> "Extrapolation") ought to as well. – J. M.'s technical difficulties Jul 22 '15 at 5:02
• @Guesswhoitis. The residual is an order magnitude larger with Method -> "Extrapolation" but certainly produces satisfactory results. Thanks, again. – bbgodfrey Jul 22 '15 at 5:06
• The error is still somewhat large. Its hard to interpret the result of what I was gonna do with it. – Sina Jul 22 '15 at 18:39