type SatelliteData <: AbstractVector
x # Vector of continuous values
discrete_field2 # Should type these
Now just make implement the array interface using
getindex(A::SatelliteData,i::Int...) = A.x[i....]
setindex!(A::SatelliteData,x,i::Int...) = (A.x[i...] = x)
length(A) = length(A.x)
Now you can write
# `s` is the input SatelliteData
# `ds` is the SatelliteData for the change in the continuous values
# Write the ODE here
# Discrete values are available via `s.discrete_field1`
Now you just have to choose how to handle the discrete changes. You can either do this via events are
tstop + dumb callback. Let’s do the
tstop way. Setup
const tstop = [ ...] # Values of t for which discrete changes occur
Now you make a dumb callback which checks if it’s a timepoint where discrete changes occur, and apply them:
sat_callback = @ode_callback begin
if t in tstop
# Apply discrete changes to `s` here
# s.discrete_field1 = 2
@ode_savevalues # make sure you save after!
and then just do the standard
prob = ODEProblem(f,s0,tspan)
solve(prob,callback = sat_callback, tstops = tstop)
and that should work (I might be missing a method in the array interface or something like that, but it should be clear from here how to do that).
This method will adaptively step as far as possible, unless there is a time in
tstop that is in the next interval, and then it will shorten to hit that point exactly, and after doing that continuous update it will apply the discrete update (and then save). This implementation using
t in tstop has a bad spot though since that check isn’t a trivial calculation, so you should replace it with something more related to the problem (i.e. if you
1:10000, instead just check if
t is the same as its rounded value. Of course, the optimization here depends on the problem). You may then want to play with
save_timeseries=false depending on your problem. Check the documentation for solver options there.
Another way, without using tstops, is to just allow the discrete change at the end of each timestep with no adaptivity. This means you just get rid of the conditional in the callback (i.e. have it use the update directly), and pass in
dt = ..... You can even do this without discrete changes every step, and just check
t in the callback as above. Adaptive timestepping should usually be better though.
These will work with any algorithm in OrdinaryDiffEq.jl (maybe not the Rosenbrock23? I’d have to see how it handles the Jacobian. It may “just work” there but there’s an optimization in the solver I can see that I’d want to do). After standardizing callbacks, that should extend to Sundials.jl and, if you interested in stochastic equations, StochsaticDiffEq.jl.
Again, I just typed this down and did test what I wrote, so excuse me if there’s a typo but it should be pretty simple from here.
So that’s one quick way to do it. The implementation time for this should be very simple: you just write the
f (which you probably already have written). There is one downside though in that the continuous variables are not named. If you only need the continuous variables to be named inside of
f, a quick way to handle this would be to use ParameterizedFunctions.jl to define the
For example, for the Lotka-Volterra equations:
f = @ode_def LotkaVolterra begin
dx = a*x - b*x*y
dy = -c*y + d*x*y
end a=>1.5 b=>1 c=3 d=1
It ends up defining the function on the indices of the input, so it would name
y and everything will work as long as the input type (as determined by the initial condition s0) is indexable (which it is given the interface we put on it).Also, this will calculate analytical functions for things like the Jacobian, so that could help speed things up (currently not used… but it will be!).
Now let’s say you want to actually want to make it all index as fields? I’ll put that in a different response because it’s a neat trick, but not truly necessary.