Error control for NDSolve

Question

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.

Answer 1

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.

From InterpolationOrder:

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.

Answer 2

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

Addendum

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.