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} /. 
(* 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]. *)
  • 2
    $\begingroup$ Can you provide a working code to illustrate your problem? $\endgroup$
    – codebpr
    Feb 25 at 15:01
  • $\begingroup$ NDSolve needs to process eqn symbolically before integrating it numerically, but ?NumericQ prevents doing it. So, replace eqn[t_?NumericQ] by eqn[t_] . $\endgroup$
    – bbgodfrey
    Feb 25 at 20:01

1 Answer 1


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

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.

  • $\begingroup$ 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 ... $\endgroup$ Feb 25 at 18:39
  • $\begingroup$ @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}] $\endgroup$
    – Goofy
    Feb 25 at 22:13
  • $\begingroup$ Related: (78641), (110950) $\endgroup$
    – Goofy
    Feb 25 at 22:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.