Efficient way to simulate large system of ODEs

General Usage Performance
, ,
1

Edit:

So the way I proposed doing it here does work, I just made a mistake that @Oscar_Smith pointed out, so using the functions this way is fine.

Original Q:

Hi there!

I’m still a bit new to Julia and trying to figure out the best way of using DifferentialEquations.jl to set up a large system of ODEs for simulation. I have a large system of equations (over 50) that represent biomass flux between containers and I’m looking for a way to define that system via a set of functions, because each container is technically defined by the same equation, just with different varying parameters.

As an example with only two containers, imagine a different system:

\frac{dX_1}{dt} = r_1X_1(1-\frac{X_1}{K_1}) + X_1(a_{12}X_2)(1-\frac{X_1}{K_1}) - X_1(a_{21}X_2) - m_1X_1
\frac{dX_2}{dt} = r_2X_2(1-\frac{X_2}{K_2}) + X_2(a_{21}X_1)(1-\frac{X_2}{K_2}) - X_2(a_{12}X_1) - m_2X_2

I know that I can set that up in Julia with the DifferentialEquations.jl package like this:

r1 = 1.5
r2 = 1.9
a12 = 0.010975
a21 = 0.01
m1 = 1.4
m2 = 1.3
k1 = 5000
k2 = 2000

function eqns!(dx, x, p, t)
    dx[1] = r1*x[1]*(1-(x[1]/k1)) + x[1]*(a12*x[2])*(1-(x[1]/k1)) - 
            x[1]*(a21*x[2]) - m1*x[1]
    dx[2] = r2*x[2]*(1-(x[2]/k2)) + x[2]*(a21*x[2])*(1-(x[2]/k2)) - 
            x[2]*a12*x[1] - m2*x[2]
end

x0 = [100.0; 1000.0]
tspan = (0.0, 10.0)
prob = ODEProblem(eqns!, x0, tspan)
sol = solve(prob)
plot(sol)

Which gives the resulting plot:
image

But in this example I have to define every parameter manually and then write out each equation. I want a way to set each parameter programatically so that in a 50-container system, rather than manually setting a_{ij} manually for each combo of i and j where i \text{ and } j = \{x \space \epsilon \space \mathbb{Z} \space | \space 0 < x \leq 50 \} I would simply call some function get_a(i, j) and get that value for any i and j combo.

I’m thinking of this in the context of turning each term into a function so that for any X_i, the eq’n would be:
\frac{dX_i}{dt} =
indep_gain(i) + dep_gain(i) - dep_removal(i) - indep_removal(i)

For my example, I thought of doing that this way where I define a series of functions (this is just a small play example that I came up with) and then the functions can be called for any X_i:

 function indep_growth(i) 
    #define r
    if i == 1
        r = 1.5
        k = 5000
    else 
        r = 1.8
        k = 2000
    end 
    # return product of r and abundance of focal x 
    return (r*x[i](1-(x[i]/k)))
end 

function dep_growth(i)
    # determine if a12 or a21 should be used 
    if i == 1 # use a12 
        a = 010975
        k = 5000
        return x[i]*(a*x[2])*(1-(x[i]/k))
    else 
        a = 0.01
        k = 2000
        return x[i]*(a*x[1])*(1-(x[i]/k))
    end 
end 

function dep_death(i)
    if i == 1
        death = dep_growth(2)
    else 
        death = dep_growth(1)
    end 
    return death 
end 

function indep_death(i)
    if i == 1
        m = 1.4
    else 
        m = 1.3
    end 

    return m*x[i]
end 

function eqns!(dx, x, p, t)
    for i in 1:2
        dx[i] = indep_growth(i) + dep_growth(i) - dep_death(i) - dep_death(i)
    end 
end

And then creating the eqns!() with something like a loop:

function eqns!(dx, x, p, t)
    for i in 1:2
        dx[i] = indep_growth(i) + dep_growth(i) - dep_death(i) - dep_death(i)
    end 
end

But however, when I then try to use solve(prob), I get

ERROR: MethodError: objects of type Float64 are not callable

With a huge stacktrace.

So my question is: Is there a way to assign the actual equations in the function via a series of functions as opposed to writing them all out manually? I have re-generate my 50-eqn system like 1000 times so I can’t write it out manually each time.

Any thoughts/help is much appreciated!!

2

Instead of return (r*x[i](1-(x[i]/k))) you want return (r*x[i]*(1-(x[i]/k))). Otherwise, Julia thinks you are trying to call x[i] as a function.

1 Like
3

Thank you! Somehow incredibly frustrating but also exciting that the problem is just a typo…

5

Yeah. That one in particular is really annoying, because as a user, it’s obvious that you can’t call a number as a function, but from Julia’s perspecitve, you totally can. If you define (x::Float64)(y::Number)=x+y, you would be making it so [4.0][1](2) = 6.0

1 Like
6

I have recently faced a similar problem.

I would suggest that you consider two options:

  1. use NetworkDynamics.jl. This package makes it easy to define an ODE problem in a network.

  2. do as you suggested, but consider adding x to the signatures of the functions called from eqns!, e.g. indep_growth(i, x) as x will likely change when looking for a solution to the ODE. Also, you could define some of the parameters like a, r and k as arrays and then index the parameter accordingly instead of branching in the function. If the problem does not have many non-itineraries, you could even consider describing it using linear algebra formulas.

I hope that helps.

1 Like
7

That’s a good point re: the x in the signatures of the functions, I actually incorporated that for the same reason you mentioned after thinking about it a little bit more. So yes, totally agree, thanks!