# Efficient way to simulate large system of ODEs

**URL:** <https://discourse.julialang.org/t/efficient-way-to-simulate-large-system-of-odes/63057>\
**Category:** Performance\
**Tags:** question, diffeq, ode\
**Created:** [June 17, 2021, 1:15am UTC](https://discourse.julialang.org/t/efficient-way-to-simulate-large-system-of-odes/63057 "2021-06-17T01:15:47Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![colebrookson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/colebrookson/32/19567_2.png) [@colebrookson](https://discourse.julialang.org/u/colebrookson)\
**Post date:** [June 17, 2021, 1:15am UTC](https://discourse.julialang.org/t/efficient-way-to-simulate-large-system-of-odes/63057/1 "2021-06-17T01:15:47Z")

</div>

# 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:

```julia-auto
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](https://global.discourse-cdn.com/julialang/original/3X/9/9/9993986b465274571ee7c6a60d55765a48b71f85.png)

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:

```julia-auto
 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:

```julia-auto
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

```julia-auto
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!!

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [June 17, 2021, 2:03am UTC](https://discourse.julialang.org/t/efficient-way-to-simulate-large-system-of-odes/63057/2 "2021-06-17T02:03:36Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![colebrookson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/colebrookson/32/19567_2.png) [@colebrookson](https://discourse.julialang.org/u/colebrookson)\
**Post date:** [June 17, 2021, 2:06am UTC](https://discourse.julialang.org/t/efficient-way-to-simulate-large-system-of-odes/63057/3 "2021-06-17T02:06:17Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [June 17, 2021, 2:11am UTC](https://discourse.julialang.org/t/efficient-way-to-simulate-large-system-of-odes/63057/5 "2021-06-17T02:11:02Z")

</div>

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`

---

<div class="post-metadata">

**Author:** ![gzagatti](https://avatars.discourse-cdn.com/v4/letter/g/97f17d/32.png) [@gzagatti](https://discourse.julialang.org/u/gzagatti)\
**Post date:** [June 18, 2021, 11:45am UTC](https://discourse.julialang.org/t/efficient-way-to-simulate-large-system-of-odes/63057/6 "2021-06-18T11:45:07Z")

</div>

I have recently faced a similar problem.

I would suggest that you consider two options:

1. use [`NetworkDynamics.jl`](https://pik-icone.github.io/NetworkDynamics.jl/dev/getting_started_with_network_dynamics/). 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.

---

<div class="post-metadata">

**Author:** ![colebrookson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/colebrookson/32/19567_2.png) [@colebrookson](https://discourse.julialang.org/u/colebrookson)\
**Post date:** [June 18, 2021, 10:18pm UTC](https://discourse.julialang.org/t/efficient-way-to-simulate-large-system-of-odes/63057/7 "2021-06-18T22:18:49Z")

</div>

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!