How to increase NDSolve accuracy for 2nd order ODE?

Question

My attempt to NDSolve a 2nd order nonlinear ODE

    rmin = 10^(-40); (* as close to 0 as possible*)

    bc = {u'[rmin] == 0, u'[1] == -u[1]}; 

    ode = r*u''[r] + 2*u'[r] + r*(Pi/64)*Exp[-128*r^2]*(1 - u[r]^5) == 0;


     s = NDSolve[{ode, bc}, u, {r, rmin, 1}, WorkingPrecision -> 70, 
        AccuracyGoal -> 20];

resulted in an accuracy of about $10^{-12}$ as can be seen from the plot

Plot[ode[[1]] /. s, {r, rmin, 1}, PlotRange -> All]

An alternative approach ( from this post)

ClearAll[s, u, v, rmin]

rmin = 10^(-40);

defv = u'[r] == v[r];

odev = r*v'[r] + 2 v[r] + r*(Pi/64)*Exp[-128*r^2]*(1 - u[r]^5) == 0;

bcv = {v[rmin] == 0, u[1] == -v[1]};

AbsoluteTiming[
 s = NDSolve[{defv, odev, bcv}, {u, v}, {r, rmin, 1}, 
    StartingStepSize -> 1*^-8, MaxStepSize -> 1*^-4, 
    PrecisionGoal -> 33, AccuracyGoal -> 33, WorkingPrecision -> 70, 
    MaxSteps -> 2*^5, InterpolationOrder -> All];]

resulted in an increased accuracy $10^{-19}$

Plot[odev[[1]] /. s, {r, rmin, 1}, PlotRange -> All

However computation time also increased significantly (about 1 min).

Is this the best accuracy Mathematica can achieve?

Motivation: The above boundary value problem corresponds to the construction of initial data for General Relativity free evolution (see Okawa, Cardoso, Pani, Phys.Rev.D, 90, 104032 (2014), eq.23 ).

Here is an attempt to address the full problem.

This kind of evolution ivolves only PDE's with respect to time and can in princple be performed in Mathematica.

However Mathematica probably lacks some tips and tricks used by Numerical GR algorithms in which case it will be more prone to crash due to numerical error on the initial data surface.

Answer 1

Use asymptotics at the initial value r == 0 to define u''[0] for ode2:

ode = r*u''[r] + 2*u'[r] + r*(Pi/64)*Exp[-128*r^2]*(1 - u[r]^5) == 0;
ode2 = u''[r] == Piecewise[{
    {u''[r] /. First@Solve[ode, u''[r]], r != 0}
    },
   2 Coefficient[
     AsymptoticDSolveValue[{1/64 E^(-128 r^2) \[Pi] r (1 - u[r]^5) + 
         2 u'[r] + r u''[r] == 0, u'[0] == 0, u[0] == u0}, 
      u[r], {r, 0, 2}], r^2]
   ]

ClearAll[r, u, u0];

(*
 * SHOOTING METHOD - Find initial conditions for BVP
 *)
psol = ParametricNDSolveValue[{ode2, u'[0] == 0, u[0] == u0}, 
   u, {r, 0, 1}, {u0}, WorkingPrecision -> 32];
bvpic = FindRoot[psol[u0][1] + psol[u0]'[1],
 {u0, 0.00019`111, 0.0002`111}, (* prec = 111 so I can be lazy when changing WP *)
 WorkingPrecision -> 32, Method -> "Brent"]
(*
 * BVP SOLUTION
 *)
sol = NDSolveValue[{ode2, u'[0] == 0, u[0] == u0} /. bvpic,
   u, {r, 0, 1},
   WorkingPrecision -> 32, InterpolationOrder -> All];
sol[1] + sol'[1]
(* omitted some underflow/overflow errors that don't seem to matter
  {u0 -> 0.00019174759848606314091019577338064}
  0.*10^-37
*)

Solution:

ListLinePlot@sol

Interpolating error as measured by the differential equation (high points):

Plot[ode2 /. Equal -> Subtract /. u -> sol // RealExponent, {r, 0, 1},
  PlotRange -> Automatic, WorkingPrecision -> 32]

The step errors as measured by ode2 should be zero plus rounding error, since that is how step methods work. They can be seen in the low points of the above graph. A common way to estimate the accuracy of a method is to run it at two different working precisions and measure the difference. If we rerun the above at MachinePrecision to produce a solutions solMP and at WorkingPrecision -> 40 to produce sol40, we get the following evidence of convergence:

Plot[{solMP[r] - sol[r], sol[r] - sol40[r]} // RealExponent,
 {r, 0, 1}, PlotRange -> Automatic, WorkingPrecision -> 40]

So solMP has an accuracy better than 10^-8, and sol has an accuracy better than 10^-16.