7
$\begingroup$

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]

enter image description here

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

enter image description here

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.

$\endgroup$
  • 1
    $\begingroup$ It's always good to add a sentence on what outcome you expect. $\endgroup$ – user21 Aug 30 at 5:54
  • $\begingroup$ 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. $\endgroup$ – J. M. will be back soon Nov 29 at 22:27
5
$\begingroup$

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

Mathematica graphics

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

enter image description here

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

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

enter image description here

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]

enter image description here

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.