# NDSolve with InterpolatingFunction object

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]. *)


• Can you provide a working code to illustrate your problem? Feb 25 at 15:01
• NDSolve needs to process eqn symbolically before integrating it numerically, but ?NumericQ prevents doing it. So, replace eqn[t_?NumericQ]  by eqn[t_] . Feb 25 at 20:01

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.

• Thank you! That is bringing me closer. Unfortunately, my rhs is vector valued, and the Diff Eq. is set up like {{x''[t],y''[t]}==rhs[x[t],y[t],t], + iniital cond.}. This yields to an error: "There are more dependent variables, ... than equations, so the system is underdetermined", I think becuase NDSolve does not see that rhs yields two values? Calling rhs twice would be inelegant, and inefficient again ... Feb 25 at 18:39
• @StefanGillessen I think you're right. You can solve it as a vector-valued function NDSolve[{xy''[t] == rhs[xy[t], t], xy[0] == {1, 0}, xy'[0] == {0, 1}}, {xy}, {t, 0, 3}] Feb 25 at 22:13
• Related: (78641), (110950) Feb 25 at 22:23