You can add observed variables, which are just scalings of other state variables, and the system will automatically eliminate them from the solve. You then just do sol[r_scaled]
or whatnot and it’ll give you those values in the solution.
Hmm, I’m not sure how to implement that…Do you think you could illustrate this with a simple example? I’ve checked this page for more info but this particular topic doesn’t have any code…
I’m not sure how to add them back after the fact, but you can do this:
using ModelingToolkit, DifferentialEquations
using Plots
@parameters t β γ N
@variables S(t) I(t) R(t) s(t) i(t) r(t) x y
D = Differential(t)
SIR_eqs = [D(s) ~ -β*s*i,
D(i) ~ β*s*i - γ*i,
D(r) ~ γ*i,
0 ~ s*N - S,
0 ~ i*N - I,
0 ~ r*N - R,
]
SIR = ODESystem(SIR_eqs, name = :SIR)
# After this, the system has only 3 states [s,i,r]
SIR = structural_simplify(SIR)
u0 = [s => 0.99, i => 0.01, r => 0.00]
pars = [β => 0.20, γ => 0.15, N=>1e3]
prob = ODEProblem(SIR, u0, (0.0, 200.0), pars)
sol = solve(prob; saveat=1:1:200)
To plot scaled
plot(sol; vars=[s,i,r])
Or to plot unscaled
plot(sol; vars=[S,I,R])
@contradict Thanks, this helps a lot! Let me see if I understand correctly: by saying
SIR = structural_simplify(SIR)
the variables s,i,r
are the variables which are actually being integrated, whereas as the variables S,I,R are the “observables” which are calculated from the solution for s(t),i(t),r(t)?
Yes, in this case that is what happens. structural_simplify
is kind of a magic box that applies several tools to try to make a simpler system to solve. As such, it isn’t very controllable so you sometimes have to play with the input equations if you want a particular shape. The fun starts here if you want to see what it is doing.
Awesome, this looks to be exactly what I originally wanted!
Another question: How does structural_simplify
“know” that s,i,r are the observables? Is it the 0
on the LHS of the equations? Or the fact that the s,i,r equations don’t involve derivatives? Or did we just happen to get lucky in this case and structural_simplify
somehow inferred what the observables would be?
…I’ll check out the source code in more detail (though based on a quick glance, it seems a bit over my head, currently ).
It uses a graph algorithm to identify the minimal equation to solve and solves that, generating explicit representations of the other equations from the elimination equations.
Hmm, very interesting!
I have some follow-up questions on this. I think I’ll probably start a new thread though since it doesn’t really have to do with my original title.
Also, would it be ok if I made a PR for the documentation to include the change of variables example with the SIR equations?
That would be nice.