Systems of differential equations

General Usage
1

Consider the set of equations and its solution. There are two variables: c and ce.

using DifferentialEquations

function diffusion!(du, u, p, t)
    α, β = p
    c, ce = u
    dc  = -α*c
    dce = -β*ce
    du .= [dc, dce]  # note the broadcast operator!
end

u0 = [1., 2.]
tspan = (0., 3.)
p = [1.0, 1.]
p1 = [1., 2., 8., 5.]  # α, β, D1, D2

prob = ODEProblem(coupled!, u0, tspan, p1)
sol = solve(prob, Tsit5())

The code runs as expected. Consider these two equations a model of a biological neuron (this is a minimal example). I would like to create a collection of N neurons that will interact with each other.
So I need an array of variables: c[1], …, c[N] and ce[1] to ce[N]. These equations will interact with each other. One way to handle this is with a “u” vector that has length 2*N, but then it becomes hard to track the variables. I would prefer to create a vector “c_vec” of size N, and “ce_vec” of size N, create a function that computes the system in terms of these higher-level constructs, and then somehow interfaces with ODEProblem and solve, from DifferentialEquation.jl .

Alternatively, there may be tools to handle networks of interacting nodes where each node is described by a set of ODEs (perhaps a reduced version of Brian2, found at Brian 2 documentation — Brian 2 2.5.1 documentation.

I could create a small example to demonstrate what I am trying to do. Thanks.

1 Like
2

You may not have to track the variables if you do, for example, cs = @view(u[1:2:end]) and ces = @view(u[2:2:end]). That’s one way to go about it - I was also interested in this problem, and had tried a graph-based description, but that proved to be too fiddly.

3

Interesting approach, thanks. I assume there is almost no performance penalty? I will try this for now. Perhaps I will get other replies. I was wondering if the ModelingToolkit might be of use. There might be other packages as well to provide a high-level interface that could feed into DifferentialEquations.jl.

4

@wsphillips and I were just talking about this today, building a high level tool for neural modeling which generates ModelingToolkit IR to then auto-GPU/distributed parallelize. I think this is a direction we will be going in the near future.

(In fact, if you know of the right grant to shoot for in this area, please share :slight_smile:)

2 Likes
5

I have no idea, Chris. But I am open to collaboration, at least as a way to learn. That is what I am interested in. Do you know about Brian2? Ultra-cool system. But Julia could improve on it.

6

You should have a look at Struct of Arrays for this. I am interested in Brian2 and a Julia version of it too!