I have an ODE where n states evolve in time. I’d like to calculate this once, and then interpolate it at different time points to speed up other calculations. However, I can’t get DataInterpolations:LinearInterpolation to work on anything more than a single variable evolving in time.

I am posting a single-state version (which works), and a multiple-state version (which doesn’t) – is there an easy trick to make the second one work, that is more elegant than serially using n different interp functions and querying them in a for-loop?

Pkg.add("DataInterpolations") # for LinearInterpolation
using DataInterpolations
# Interpolation of a single variable (n=1 evolving state) -- this works
u = rand(5)
t = 0:3
interp = LinearInterpolation(u,t);
interp.u
interp.t # times corresponding to the u points
interp(2.5) # Gives the linear interpolation value at t=2.5
# Attempt to interpolate an ODE output with n=4 evolving states
# Setting up an example ODE output (i.e., the states at 4 timepoints)
ode_samples_t0 = [0.0 0.0 0.0 0.0]
ode_samples_t1 = [0.1 0.2 0.3 0.4]
ode_samples_t2 = [0.1 0.2 0.3 0.4]*2
ode_samples_t3 = [0.1 0.2 0.3 0.4]*3
ode_samples = [ode_samples_t0, ode_samples_t1, ode_samples_t2, ode_samples_t3]
u = ode_samples
# Setting up the interpolation seems to work
t = 0:3
interp = LinearInterpolation(u,t);
interp.u
interp.t # times corresponding to the u points
# But, running it does not work:
interp(2.5)
# ERROR: MethodError: objects of type
# LinearInterpolation{Array{Array{Float64,2},1},Array{Int64,1},true,Array{Float64,2}}
# are not callable
# Use square brackets [] for indexing an Array.

If you solve the ode using Ordinarydiffeq or Differentialequations.jl, the solution object has built in interpolation of the same order as the integrator. See the docs for how to use it, essentially call it as a function

However, do you know if it is really doing a simple & fast interpolation, or is it re-calculating the whole ODE? Re-calculating may get slow once I’m up at n=1000 or whatever, and I just want to be able to very quickly query this ODE at different timepoints, for input into another set of ODEs. (The ODEs I’m dealing with result in pretty smooth & nonstiff curves so a pretty simple interpolator should work fine.)

It doesn’t recalculate anything. You can build the interpolation using the internal steps of the Runge-Kutta method, or through a Hermite interpolation using the two values and the two end derivatives. There’s also a few interpolations in the Rosenbrock methods that use internal step behavior as well. But essentially, the k_i of the methods, which in some sense have “f-calls embedded in them”, are what are you used to give a high order continuous function. Using those values it’s just a search (for the values at a given time point), and a linear combination.