Why does the domain of an interpolating function not match the corresponding ParametricNDSolve solution range?

Question

Question

I had believed - and mostly observed - that when I use ParametricNDSolve the obtained ParametricFunctions given as the solution to some ODE or DAE will be translated into InterpolatingFunctions that have a Domain that corresponds to the solution range of the differential (algebraic) equations.

But that does not seem to be the case. So the following questions arise:

1. Why does this happen (e.g. in which cases)?
2. What will this mean for the reliability of results?
3. In case there are negative implications, what can be done about it?

Minimal Working Example

In the following example I have a very simple growth model for a stock, where the growth rate is a discrete variable that will be changed by WhenEvents at the times 2018, 2021, 2026. The simulation is to run from 2016 to 2027. In the example, I am using the options for a more elaborate simulation of a DAE.

parfuncs = Association @ ParametricNDSolve[
    {
        stock'[t] == stock[t] × growthRate[t],
        stock[2016] == 1,
        growthRate[2016] == gr0,
        WhenEvent[ t > 2018, growthRate[t] -> gr1],
        WhenEvent[ t > 2021, growthRate[t] -> gr2], 
        WhenEvent[ t > 2026, growthRate[t] -> gr3]  
    },
    { stock, growthRate },
    { t, 2016, 2027 },
    { gr0, gr1, gr2, gr3 },
    DiscreteVariables -> { growthRate },
    InterpolationOrder -> Automatic, (* alternatively: 1, 2, or All *)
    Method -> {
        "EquationSimplification" -> { Automatic, "SimplifySystem" -> False },
        "ParametricSensitivity" -> None,
        "IndexReduction" -> {
            "Pantelides",
            "IndexGoal" -> 1,
            "ConstraintMethod" -> Automatic
        }
    } 
];

intFunc = ReplaceAll[
    parfuncs[stock][gr0,gr1,gr2,gr3],
    {
      gr0 -> 0.01,
      gr1 -> 0.03,
      gr2 -> 0.05,
      gr3 -> 0.01
    }
]

As we can see, the Domain is {2020., 2030.}. I also note, that the option InterpolationOrder seems to be rather irrelevant (e.g. setting it to 1 or 2 does not change anything at all, the Order remains to be 3).

Edit:

I note, that when I shift the actual times to 0 (e.g. 2016 $\rightarrow$ 0, 2018 $\rightarrow$ 2, 2021 $\rightarrow$ 5, 2026 $\rightarrow$ 10, 2027 $\rightarrow$ 11), the Domain is given correctly as {0., 11.}, so there may be rounding issues involved.

Checking this hunch for my MWE reveals:

intFunc["Domain"]

{{2016., 2027.}}

So it is rounding and it certainly is quite confusing, as for certain parameters, the function will give a domain of {2020.,2020.}.

Answer 1

The issue is that the formatting routine BoxForm`ArrangeSummaryBox for InterpolatingFunction sets PrintPrecision -> 3, so you get the domain rounded to three digits.

foo = parfuncs[stock][0.01, 0.03, 0.05, 0.01]

MakeBoxes[#, StandardForm] &@foo /. 
  HoldPattern[PrintPrecision -> 3] -> PrintPrecision -> 4 //
 RawBoxes

You can also somewhat guess this is the problem if you copy and paste the domain from the summary form.