I am using DifferentialEquations.jl for solving an ODE, and it works great. However, I’d like to save intermediate outputs, and I can’t figure out how to do that.

Consider the example below:

using DifferentialEquations
α=1
u0=1/2
function f(t,not_used,u)
p = u^2 # how to save p??
return α*u
end
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
solve(prob)

How can I save p as a function of t?

This seems to be a pretty standard application of solving ODEs, consider the system below for instance

x_dot(t) = A * x(t) + B * u(t)
y(t) = C * x(t) + D * u(t)

Sure, I can do that. However, then I have to recompute all intermediate variables a second time. My application is fairly computationally intensive, so I’d very much prefer not to compute these intermediate variables twice.

It’s at most 1/8 of the cost, if y(t) >> x(t) in terms of size. That’s the reason why it’s not implemented in any ODE solver. I can show some ways to do it, but if you’re looking at this for efficiency you’re looking in the wrong place.

For the systems I’m interested in, typically size(y) is similar to size(x) – I’m not trying to save every single intermediate variable. Does your conclusion still apply then (or is the cost even more negligible)?

Out of curiosity, where does the 1/8 factor come from?

5th order methods like Tsit5() do 7 f evaluations per step, so if you do an entire extra one to compute the y part once more, then it’s at most 7+1 evaluations where you’re adding 1, so 1/8. If you’re only having to do half of that work, and if you’re outputting less points then the steps, and …, then that factor keeps dropping.