# Problem with ParametricNDSolve

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?

• Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. – bbgodfrey Feb 18 '15 at 3:10

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


• I believe setting PrecisionGoal is enough for this problem – Dr. belisarius Feb 18 '15 at 3:31
• Thanks! i tried that before, but for some reason i thought it didn't work. – Ivactheseeker Feb 18 '15 at 12:33

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.

• Thanks for the answer! So the step size is a function of the domain of the independent variable? I couldn't figure out how changing tmin affects the whole integration bacause, as you said, the numerical integration starts at s=0 and goes backwards. – Ivactheseeker Feb 19 '15 at 16:21
• @Ivactheseeker Yes, it's a function of the domain among other things. NDSolve tests a few values at the start of the integration interval and adjusts the step size according to its error estimate. In this case, it jumps to MaxStepFraction. For tmin = -50, it tries the step -5; but that lands where a'[s] is changing, and the error estimate is too large. So it tries a step of -5/4 = -1.25, which it uses until it needs to adjust the step size again. With tmin = -100, the step of -10 jumps over the variation and the error is estimated to be small; so it keeps on going. – Michael E2 Feb 19 '15 at 18:29
• @Michael_E2 Ok, thanks again!. – Ivactheseeker Feb 20 '15 at 11:14