NDSolve with InterpolatingFunction object

Question

Imagine some functions like x1[t], y1[t] given as InterpolatingFunction objects, for example as a result of solving an ODE with NDSolve.

Now I have an ODE for x[t], y[t], where (repeated) combinations like x1[t]/Sqrt[(x[t] - x1[t])^2+(y[t] - y1[t])^2] or y1[t]/Sqrt[(x[t] - x1[t])^2+(y[t] - y1[t])^2] appear.

This works fine, however it is inefficient, as for each call of x1[t] or y1[t] the InterpolatingFunction object is evaluated during the NDSolve. What would be much more efficient, is if one could for each time step evaluate x1 and y1 (i.e. each once) and plug in the result. I have not managed so far to make that efficient, as NDSolve does not accept a function of type f[t_?NumericQ, ...] when t is the independent time variable.

Any hints on how to speed up?

A minimum example how to show what is not working:

x1Rule = 
 NDSolve[{x1''[t] == -x1[t], x1[0] == 0, x1'[0] == 1}, 
   x1, {t, 0, 5}][[1]];
(* works fine *)

eqn[t_?NumericQ] := 
 Block[{x1func = 
    x1[t] /. x1Rule}, {x''[t] == -x[t] + 0.01*x1func + 0.02*x1func^2, 
   x[0] == 0.2, x'[0] == 0.5}
(* this function tries to be efficient in caching the value of x1[t] and re-use it twice *)

eqn2 = {x''[t] == -x[t] + 0.01*x1[t], x[0] == 0.2, x'[0] == 0.5} /. 
   x1Rule;
(* this is the same, but using the solution twice *)

NDSolve[eqn2, x, {t, 0, 3}]
(* works fine *)

NDSolve[eqn[t], x, {t, 0, 3}]
(* yields: NDSolve::deqn: Equation or list of equations expected instead of eqn[t] in the first argument eqn[t]. *)
 

Answer 1

Make the right-hand side of the ODE a function of the variables:

ClearAll[rhs];
rhs[x_?NumericQ, t_?NumericQ] := 
  With[{x1func = x1[t] /. First@x1Rule},
   -x + 0.01*x1func + 0.02*x1func^2];

NDSolve[{x''[t] == rhs[x[t], t], x[0] == 0.2, x'[0] == 0.5}
 , x, {t, 0, 3}]

You could apply First to x1Rule = First@NDSolve[...] instead of inside rhs if you want to save a few microseconds per evaluation.

Side tip: You might use InterpolationOrder -> All in the x1Rule call. It ameliorates the weak singularities at each step of the x1 solution, but it makes the interpolating function slightly longer to evaluate. It might speed things up overall by leading to fewer steps in the second ODE. Or it might not.