I'm using DifferentialEquations.jl to solve an ODE system as shown below. The result is not really relevant since p only contains test parameters for the purpose of producing a MWE, but the key is that I am seeing a lot of memory allocation despite using an in-place ODE function.
using DifferentialEquations
function ode_fun!(du,u,p,t)
a,b,c,d,e = p
X = @. u[1] * a * ((b-c)/b)
Y = @. u[2] * d * ((b-e)/b)
du[1] = -sum(X) + sum(Y) - u[1]*u[2]
du[2] = sum(X) - sum(Y) - u[1]*u[2]
end
#exemplary parameters
a = collect(10:-0.1:0.1)
b = a.^2
c = b*0.7
d = collect(0.01:0.01:1)
e = b*0.3
u0 = [1.0, 0.5]
p = [a,b,c,d,e]
tspan = [0.0, 100.0]
t = collect(0:0.01:100)
prob = ODEProblem(ode_fun!,u0,tspan,p,saveat=t)
@time sol = solve(prob)
1.837609 seconds (5.17 M allocations: 240.331 MiB, 2.31% gc time) #Julia 1.5.2
Since I need to solve this ODE system repeatedly I would like to reduce allocations as much as possible and am wondering if there is anything that can be done about them. I have been wondering if the issue lies with X and Y and have tried to preallocate these outside the ODE function, but have unfortunately not succeeded in reducing allocations that way.