Changing ParametricNDSolveValue options for only part of the run

Question

Here is an example of the system that I am running,

eqns = {y'[t] == y[t] (x[t] - 1), x'[t] == x[t] (a - y[t]), x[0] == 1,y[0] == 2.7};
tmax = 100;
pfun = ParametricNDSolveValue[eqns, Integrate[y[t], {t, 0, tmax}] , {t, 0, tmax},{a}];
Plot[pfun[a], {a, 0, 1}, AxesLabel -> {"a", "y"}]

and for a particular range of the parameter, "a", I get numerical giberish as shown here:

To check what is going on in this region I looked at the timeseries of x and y for different values of "a". If you plot example timeseries for "a" values higher than the numerical noise you get nice oscillating curves

This was found using this code:

tmax = 100;
amin = 1;
amax = 3;
aint = 1;
eqns1[a_] := {y'[t]==y[t] (x[t] - 1), x'[t]==x[t] (a - y[t]), x[0]==1, y[0]==2.7};
pn1 = ParametricNDSolveValue[eqns1[a], x, {t, 0, tmax}, {a}]
pn2 = ParametricNDSolveValue[eqns1[a], y, {t, 0, tmax}, {a}]

plt1 =Plot[Evaluate[Table[pn1[a][t],{a,amin,amax,aint}]], {t,0,tmax},PlotRange-> All];
plt2 =Plot[Evaluate[Table[pn2[a][t],{a,amin,amax,aint}]], {t,0,tmax},PlotLegends->  
Range[amin, amax, aint], PlotRange -> All];
GraphicsRow[{plt1, plt2}]

However, to get the timeseries to plot for the a values around the noise, i.e. between 0 and 0.1, I had to add Workingprecision-> 20 and Method-> "LSODA" into pn1 and pn2 (into the ParametricNDSolveValue options). Which nicely renders the curves

So the problem is that these options are not necessary for all the "a" values when I plot pfun above. In fact, if I include these options into my real system (which is larger) the whole plot takes exponentially longer to finish. So I would like to turn on these options for only some "a" values and not others, when I make a plot of pfun. Or use another trick you might know.

I would greatly appreciate any suggestions on how to do this, or other comments. Thank you.

Answer 1

You can switch between the options by defining a function that calls one form of ParametricNDSolve if a is in a certain range, and otherwise calls the other form.

eqns = {y'[t] == y[t] (x[t] - 1), x'[t] == x[t] (a - y[t]), x[0] == 1,
    y[0] == 2.7`20};
tmax = 100;

pfun = ParametricNDSolveValue[eqns, Integrate[y[t], {t, 0, tmax}], {t, 0, tmax}, {a}];
pfun2 = ParametricNDSolveValue[eqns, Integrate[y[t], {t, 0, tmax}], {t, 0, tmax}, {a}, 
   WorkingPrecision -> 20, Method -> "LSODA"];

myp[a_] /; 0.04 < a < 0.1 := pfun2[a];
myp[a_] := pfun[a];

You still may have some work to do on your DE system, however:

Plot[myp[a], {a, 0, 1}, AxesLabel -> {"a", "y"}]

For example:

pfun2[0.062]

ParametricNDSolveValue::ndsz: At t$159261 == 83.68333550162991842139594420675849673339`20., step size is effectively zero; singularity or stiff system suspected. >>

InterpolatingFunction::dmvali: The integration endpoint 100 in dimension 1 lies outside the range of data in the interpolating function. Extrapolation will be used. >>

(* Out= 0.*10^63 *)