Save additional quantities using DifferentialEquations.jl pkg

Within my ODE function additional quantities are calculated which I would like to save to extract them later for e.g. inspection or plotting.

There is the approach with the SavingCallback; however, the save_func(u, t, integrator) only has access to the states.

Another option might be to use other callbacks like ContinousCallback.

My idea was to extend the differential equations with additional states and setting the corresponding differential equations to zero. It works I guess, but there might be better ways? One shortcoming is that I then have to supply initial values for the additional quantities. Any other ideas?

2 Likes

What about using a DEDataArray style?

1 Like

Thanks for the fast answer!

I am using that for input on the right hand side of the differential equations, e.g. control variables.

But why not using it on the LHS as well for saving values like u.f1 = 1.0.

It may not update all internal vectors.

You can also get access to your ODE function in the SavingCallback as integrator.f.f (if I remember correctly).

2 Likes

Iā€™m using a PeriodicCallback for control simulation and just implemented your du[4] = 0.0 idea. In my controller, on every tick, I update the state variable with integrator.u[4] = var_of_interest. Seems to work great.

Could you pose the model as a DAE and use a DAE solver?

Considering the DEDataArray approach:

The minimal example shown below throws an error as an object of type Array{T,2} for the field x of SimType needs to be assigned somewhere in the code.

It works if the inner constructor of SimType vectorizes the x argument, i.e. SimType(x,f1) = new{eltype(x)}(vec(x),f1).

using DifferentialEquations
mutable struct SimType{T} <: DEDataVector{T}
    x::Array{T,1}
    f1::T
end

function f(du,u,p,t)
    du[1] = -0.5*u[1]
    du[2] = -0.5*u[2]
    u.f1 = 1.
end

u0 = SimType([10.0;10.0], 0.0)
prob = ODEProblem(f,u0,(0.0,10.0))

sol = solve(prob)

If you want to use a matrix, you should type it as a matrix, i.e.

mutable struct SimType{T} <: DEDataMatrix{T}
    x::Array{T,2}
    f1::T
end

It has the same semantics as Array

It is not a matrix though, but a vector, see the example which is taken from the documentation. Somewhere in the DifferentialEquations.jl pkg it changes to a matrix while solving.

If I remove the u.f1 = 1. in the ode function, it works fine.

It makes a matrix for the stiff ODE solver. Does choosing Tsit5() work?

Yes, that works.

What does false .* u .* u' generate when u is a SimType?

ERROR: MethodError: no method matching SimType(::Array{Float64,2}, ::Float64)

Same error when you try to solve the ode with a stiff solver.