Hi there!

I’m working on a problem where I must solve a system of ODEs repeatedly with different parameters and initial conditions. In particular, I’m not really interested in the solution as a whole, I only need to know the final state of the system once a certain condition is met.

Here’s a working example: Let’s suppose that I want to solve the Lorenz system 10 times, each time with different initial positions and system parameters. I want to evolve the system until `t = 100.0`

or the particle reaches a sphere of a certain radius `r`

, that is, it’s `x^2 + y^2 + z^2 - r^2 = 0`

. Furthermore I don’t want the solver to save all intermediate steps and interpolation information, I’m only interested in the particle’s position at the final integration time. For that end, one could use the following code:

```
using DifferentialEquations
function parameterized_lorenz!(du,u,p,t)
du[1] = p[1]*(u[2]-u[1])
du[2] = u[1]*(p[2]-u[3]) - u[2]
du[3] = u[1]*u[2] - p[3]*u[3]
end
function main()
initialData = [[rand(-10.0:10.0), rand(-10.0:10.0), rand(-10.0:10.0)] for i in 1:10]
parameters = [[rand(0.0:10.0), rand(0.0:28.0), rand(0.0:10.0/3.0)] for i in 1:10]
tspan = (0.0,100.0)
maxRadius = 5.0
condition(u, t, integrator) = (u[1]^2 + u[2]^2 + u[3]^2) - maxRadius^2
affect!(integrator) = terminate!(integrator)
cb = ContinuousCallback(condition, affect!)
for i in 1:10
u0 = initialData[i]
p = parameters[i]
prob = ODEProblem(parameterized_lorenz!,u0,tspan,p)
sol = solve(prob,
Tsit5(),
adaptive = true,
maxiters = 1e5,
abstol = 1.0e-10,
retol = 1.0e-10,
dtmin = 1.0e-12,
callback = cb,
dense = false,
save_on = false,
save_start = false,
save_end = true,
alias_u0 = true,
calck = false
)
println("Endpoint ", i, ": ", sol[1][1], ", ", sol[1][2], ", ", sol[1][3])
end
end
main()
```

Running this code with `--track-allocation=user`

I see that `solve`

and `println`

are the only two instructions allocating memory inside the loop. Since `println`

is only there for testing purposes, I would like to know what else can be done in order to minimize allocations (again, given that I’m only interested in the position values at the final time step).

More specifically, is it possible to pre-allocate memory for the solver’s output so that it does not need to create new solution objects at each iteration?

I think this use case of repeatedly solving a system is fairly common, so if I’m missing any existing discussion on the topic please let me know.