Possible to adaptively shrink number of state variables/ODEs?

I’m solving an ODE system adaptively in time (eg RK23). At minimum, I’m solving for the evolution of 10 coupled state variables so I have 10 (coupled non-linear) ODEs that are solved from some initial time to final time. This is easy to do on its own. But along the way, I have other objects that “live and die” and that have their own 5 state variables each and contribute coupled terms to the RHS of the 10 main ODEs. They all start at the same initial time as my main object, but many may not survive until the overall final time – after some earlier time, all their state variables and ODEs will be zero because they’ve “died”. There can be tens or hundreds of these objects so 5*N ODEs. If most of these are “dead” and my state variable array has lots of unchanging zeros, I think it’s memory-inefficient to keep tracking those state variables / computing their ODEs.

Is there a technical name for this kind of ODE system problem and can Julia solvers deal with it? Essentially I think it implies that I need an ODE solver that is not only adaptive in time but also can shrink the size of the state variable array (after an object “dies”, I can delete its 5 zeroed-out state variables and stop tracking its ODE terms for memory/performance).

I guess using callbacks is an option? Or are Julia solvers capable of handling ODE systems with 100-1000 state variables, some of which will be zeroed-out / unevolving for most of the overall integration time?

My advice would be to just try it with the full state of the 10 + 5 * N variables, even though most will be zero most of the time. I’m not an expert on the different solvers and their individual problems when the state becomes really large and mostly zero (perhaps that could cause some issues, but I can’t imagine why right now), but from my own experience, ODEs of size 100-1000 are not a problem per se – it mostly comes down to optimizing the ODE function and choosing a proper solver.

For instance, you can easily get ODEs with 100s or 1000s of states when solving PDEs with the method of lines (discretizing the space dimension). There are comprehensive docs of the (de facto) standard library on solving ODEs

This example features an ODE (from discretizing a PDE) with 100x100x2 = 20_000 state variables.

Of course, the structure of your equations will play a big role in what methods are best for your problem, but from what you describe (the “living-and-dying” species only coupling to the 10 “core” variables) it sounds like there could be a lot of sparsity in your problem. This page of the docs might be helpful in that case. (Note that it also recommends to just try solve(prob) first and continue from there if the performance or accuracy are not satisfactory.


When I really don’t find a solver that can handle the problem in an acceptable amount of time I would start looking for more non-standard solutions. I have also thought about the kind of approach you describe (dynamically changing the state size), but never looked for nor heard of an approach like that.

Depending on what you are actually trying to simulate, perhaps there is also a different kind of representation which you can choose to make the problem more tractable – e.g. rewriting/approximating/coarse-graining the equations or going to agent-based simulations of some sort (but ODEs are usually the most straightforward in my experience).

It’s definitely possible, the documentation has examples of this. See for example the splitting cells model:

This uses resize(integrator, n) inside a callback to change the number of ODEs. But there’s other tools as well. There is indeed a whole section on resizing tools:

These days I’d say just use a vector and resize it around as needed. However, I previously used to work on a lot more complex structures for handling these kinds of hierarchical size-changing dynamical systems, so there’s a tool MultiScaleArrays.jl to be aware of.

Its tests showcase all sorts of cases where you’re adding cells and deleting cells as the differential equation is being solved:

You can then couple this with JumpProcesses.jl in order to bring in discrete-time Gillespie type processes for the birth-death process

https://docs.sciml.ai/JumpProcesses/stable/

See for example the hybrid model tutorial

https://docs.sciml.ai/JumpProcesses/stable/tutorials/jump_diffusion/

A lot of the ODE and SDE solvers in Julia were first built for these whole-cell dynamically sized problems.

1 Like