ODE solver save intermediate outputs (y = C*x + D*u)

General Usage
1

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)

How does one save y(t)?

2

Just calculate y(t) from the output. solve outputs a continuous function.

3

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.

4

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.

5

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?

6

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.