Problem with ParametricNDSolve

Question

i have a strange problem with ParametricNDSolve. Given the sistem of differential euqations

fun = 
  ParametricNDSolveValue[
    {x'[s] == v[s], v'[s] == a[s], a'[s] == a[s] - 2 x[s] E^(-x[s]^2),
     x[0] == x0, v[0] == v0, a[0] == 0},
    {x, v, a}, {s, tmin, 0}, {x0, v0, tmin}];

ParametricNDSolve returns the correct solution only for certain values of tmin. For example, with

{x0, v0} = {10, 2}

it gives the correct solution only for tmin < -74, while for tmin > -74 it doesn't.

It's possible to visualize this with

Table[Plot[fun[10, 2, -50][[i]][t], {t, -10, 0}, PlotRange -> All], {i, 1, 3}]
Table[Plot[fun[10, 2, -100][[i]][t], {t, -10, 0}, PlotRange -> All], {i, 1, 3}]

and

Plot[fun[10, 2, t][[1]][-5], {t, -100, 0}]

where is possible to see the dependence of the solution calculated in t = -5 on the value of tmin. I tried to increase the WorkingPrecision but to no avail. Any ideas?

Answer 1

Use AccuracyGoaland PrecisionGoal

fun = ParametricNDSolveValue[{x'[s] == v[s], v'[s] == a[s], a'[s] == a[s] - 2 x[s] E^(-x[s]^2), 
                   x[0] == x0, v[0] == v0, a[0] == 0}, {x, v, a}, {s, tmin, 0}, {x0, v0, tmin}, 
                                                      AccuracyGoal -> 10, PrecisionGoal -> 10];

Table[Plot[fun[10, 2, -50][[i]][t], {t, -10, 0},  PlotRange -> All], {i, 1, 3}]
Table[Plot[fun[10, 2, -80][[i]][t], {t, -10, 0},  PlotRange -> All], {i, 1, 3}]

Answer 2

The underlying problem is with the step size, which is controlled by various options, including PrecisionGoal (as shown in belisarius's answer. It is also controlled by AccuracyGoal, MaxStepSize, MaxStepFraction and some others.

The default setting for MaxStepFraction is 1/10, which is not fine enough to stumble upon a little blip near the initial condition (at the beginning of the integration, which NDSolve calculates "backwards" from right to left). The following shows the step size to be 10.

Clear[fun]; 
fun = ParametricNDSolveValue[{
   x'[s] == v[s], v'[s] == a[s], a'[s] == a[s] - 2 x[s] E^(-x[s]^2),
   x[0] == x0, v[0] == v0, a[0] == 0},
  {x, v, a}, {s, tmin, 0}, {x0, v0, tmin}];

With[{f0 = fun[10, 2, -100]},
 GraphicsRow@
  Table[Plot[f0[[i]][t], {t, -20, 0}, PlotRange -> All, 
    Mesh -> f0[[i]]["Coordinates"], 
    MeshStyle -> {PointSize[Medium], Red}], {i, 1, 3}]]

Adding any one of the following lines results in an accurate solution for tmin = -100.

MaxStepSize -> 2         (* should work for all tmin *)
MaxStepFraction -> 1/20  (* relative to the size of tmin *)
AccuracyGoal -> 17, PrecisionGoal -> 0   (* both together *)
PrecisionGoal -> 8.3     (* or higher; > 8.5 is fairly robust *)

The problem is undoubtedly that the error estimate near the initial condition is so small that NDSolve feels it can take the maximum step. This passes over the interval of s where x is small, which is where the acceleration a[s] changes significantly. Another approach is to slow down the integration when x[s] gets small. If x[s] == 5, then the term 2 x[s] E^(-x[s]^2) from a[s] will be about 10^-10, which is perhaps a little small (the default AccuracyGoal and PrecisionGoal is about 8). But the starting step size is small enough that the solution is found.

Clear[fun];
fun = ParametricNDSolveValue[{
   x'[s] == v[s], v'[s] == a[s], a'[s] == a[s] - 2 x[s] E^(-x[s]^2),
   x[0] == x0, v[0] == v0, a[0] == 0,
   WhenEvent[x[s] == 5, "RestartIntegration"]},
  {x, v, a}, {s, tmin, 0}, {x0, v0, tmin}];

With[{f0 = fun[10, 2, -100]},
 GraphicsRow@
  Table[Plot[f0[[i]][t], {t, -10, 0}, PlotRange -> All, 
    Mesh -> f0[[i]]["Coordinates"], 
    MeshStyle -> {PointSize[Small], Red}], {i, 1, 3}]]

This answer basically explains what (I think) is going on with NDSolve. The solution proposed by belisarius seems like a good first stab at solution. My own knee-jerk reaction in such a case is to try WorkingPrecision -> 20, which increases AccuracyGoal and PrecisionGoal to 10, as in belisarius's solution, and also uses arbitrary precision reals. The setting slightly above MachinePrecision with the precision tracking gives me some feedback about the numerical stability of the computation of the solution at MachinePrecision.