Save additional quantities using DifferentialEquations.jl pkg

Specific Domains Modelling & Simulations
1

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
2

What about using a DEDataArray style?

1 Like
3

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.

4

It may not update all internal vectors.

5

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

2 Likes
6

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.

7

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

8

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)
9

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

10

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.

11

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

12

Yes, that works.

13

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

14

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

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

15

Is there a plan to resolve this in the near future?

16

I’ve come across the same problem using DEDataVector, same error as @moesphere. Except I also see the error when using Tsit5().

Inside my ode function I calculate an intermediate vector quantity which is used in computing the RHS of the ODE. I want to save this vector at each output time. Seems simple but struggling to figure out how to do it!

17

We have a new interface for this coming out soon.

4 Likes
18

Any update on when the new interface will be released? This feature would be really useful!

19

It exists but needs better docs. This week and the next are ergonomics and docs.

2 Likes
20

@TimKnab using ModelingToolkit, you can add observational equations and use them from the solution and they will be factored out of the solver process for full performance. Example:

https://mtk.sciml.ai/dev/tutorials/acausal_components/