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