# How to increase NDSolve accuracy for 2nd order ODE?

My attempt to NDSolve a 2nd order nonlinear ODE

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

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

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[] /. 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 == -v};

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[] /. 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.

• It's always good to add a sentence on what outcome you expect. Aug 30, 2019 at 5:54
• You might consider changing the Method used for integrating your ODE; I've found that for really stringent precision requirements, Method -> 'StiffnessSwitching" (which uses Bulirsch-Stoer) does remarkably well. Nov 29, 2019 at 22:27

Use asymptotics at the initial value r == 0 to define u'' 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, u == u0},
u[r], {r, 0, 2}], r^2]
] ClearAll[r, u, u0];

(*
* SHOOTING METHOD - Find initial conditions for BVP
*)
psol = ParametricNDSolveValue[{ode2, u' == 0, u == u0},
u, {r, 0, 1}, {u0}, WorkingPrecision -> 32];
bvpic = FindRoot[psol[u0] + psol[u0]',
{u0, 0.00019111, 0.0002111}, (* prec = 111 so I can be lazy when changing WP *)
WorkingPrecision -> 32, Method -> "Brent"]
(*
* BVP SOLUTION
*)
sol = NDSolveValue[{ode2, u' == 0, u == u0} /. bvpic,
u, {r, 0, 1},
WorkingPrecision -> 32, InterpolationOrder -> All];
sol + sol'
(* 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.