Does mathematica cache the interpolating function from ParametricNDSolve?

Question

I want to evaluate the solution of a system of non-linear ODEs using ParametricNDSolve. The output of ParametricNDSolve is a ParametricFunction object. Let's call it $f_\theta(t)$, where $\theta$ is a list of all parameters used in ParametricNDSolve.

For a given choice of the parameters $\hat\theta$, the ParametricFunction becomes a InterpolatingFunction object. I want to evaluate this function at multiple points $\{t_1,t_2,...,t_n\}$, which is, obtaining the value $f_{\hat\theta}(t_i)$.

Does Mathematica fully solve the system of equations each time I call $f_{\hat\theta}(t_i)$ for a different $i$? Or does it cache the InterpolatingFunction after the first call (let's say, $f_{\hat\theta}(t_1)$) and uses it to obtain the value of $f_{\hat\theta}$ at $t_2,...,t_n$?

Answer 1

To extend this Q&A beyond one that is "easily found in the documentation," where some basic facts about "ParametricCaching" are given, I'll add some undocumented suboptions to Method -> {"ParametricCaching" -> {Automatic, ...opts...}}:

sol = ParametricNDSolveValue[{x''[t] + a x[t] == 0, x[0] == b, 
    x'[0] == 0}, x, {t, 0, 1000}, {a, b}, MaxSteps -> ∞];
sol[[-1, -1, -3]]
(*
{"Cache" -> True, "CacheTableLength" -> 19, "CacheTableWidth" -> 7, 
 "CacheKeyMaxBytes" -> 1000000, "CacheResultMaxBytes" -> 1000000, 
 "KeyComparison" -> None, "ResultComparison" -> LessEqual}
*)

For instance, we can use the timings to verify when a cached solution is reused. This confirms that the "TableCache*" suboptions seem to set the dimensions of the cache table:

sol1 = ParametricNDSolveValue[{x''[t] + a x[t] == 0, x[0] == b, 
    x'[0] == 0}, x, {t, 0, 1000}, {a, b}, MaxSteps -> ∞, 
   Method -> {"ParametricCaching" -> {Automatic, 
       "CacheTableLength" -> 1, "CacheTableWidth" -> 1}}]; 
sol1[1, 1][1] // AbsoluteTiming
sol1[1, 1][2] // AbsoluteTiming
sol1[1, 2][2] // AbsoluteTiming
sol1[1, 1][2] // AbsoluteTiming
(*
  {0.017357, 0.540302}
  {0.000014, -0.416147}
  {0.020204, -0.832294}
  {0.01881, -0.416147}   <-- N.B. Recomputed
*)

sol2 = ParametricNDSolveValue[{x''[t] + a x[t] == 0, x[0] == b, 
   x'[0] == 0}, x, {t, 0, 1000}, {a, b}, MaxSteps -> ∞, 
  Method -> {"ParametricCaching" -> {Automatic, 
      "CacheTableLength" -> 1, "CacheTableWidth" -> 2}}]; 
sol2[1, 1][1] // AbsoluteTiming
sol2[1, 1][2] // AbsoluteTiming
sol2[1, 2][2] // AbsoluteTiming
sol2[1, 1][2] // AbsoluteTiming
(*
  {0.019945, 0.540302}
  {0.000017, -0.416147}
  {0.01882, -0.832294}
  {0.000022, -0.416147}  <-- N.B. Reused
*)

The memory-related options should be self-explanatory and useful. We can test "CacheResultMaxBytes" but I'm not sure how to test "CacheKeyMaxBytes".

(* bytes of a result *)
sol1[1, 1] // ByteCount
(*  305704  *)

sol3 = ParametricNDSolveValue[{x''[t] + a x[t] == 0, x[0] == b, 
   x'[0] == 0}, x, {t, 0, 1000}, {a, b}, MaxSteps -> ∞, 
  Method -> {"ParametricCaching" -> {Automatic, 
      "CacheResultMaxBytes" -> 400000}}]; 
sol3[1, 1][1] // AbsoluteTiming
sol3[1, 1][2] // AbsoluteTiming
(*
  {0.015709, 0.540302}
  {0.000013, -0.416147}  <-- N.B. Reused
*)

sol4 = ParametricNDSolveValue[{x''[t] + a x[t] == 0, x[0] == b, 
   x'[0] == 0}, x, {t, 0, 1000}, {a, b}, MaxSteps -> ∞, 
  Method -> {"ParametricCaching" -> {Automatic, 
      "CacheResultMaxBytes" -> 1000}}];
sol4[1, 1][1] // AbsoluteTiming
sol4[1, 1][2] // AbsoluteTiming
(*
  {0.018927, 0.540302}
  {0.022862, -0.416147}  <-- N.B. Recomputed
*)

I don't about keys and their function in ParametricNDSolve caching, nor how "ResultComparison" is used.

(Note it's important to rerun ParametricNDSolveValue each time you want to test to reset the cache table.)