Solving same DiffEq repeatedly but with different initial conditions

I have a set of coupled ODEs that I’d like to solve, which is simple enough using DifferentialEquations.jl. My problem is that I need to solve the same ODE repeatedly for many different initial conditions parametrized by an array x.

Currently I just do this in a loop. I’m sure there’s a smarter way but I haven’t quite found it in the documentation yet. This is very inefficient of course since I need to keep setting up the problem again every time.

Could someone give me some pointers?

Here is a minimal non-working example of what I’m doing:

using Elliptic.Jacobi

function dn_ab!(du,u,p,t)
    du[1] .= +im*p[1].*u[1] .+ im.*u[2].*dn.(t,real(p[2]))
    du[2] .= -im*p[1].*u[2] .+ im.*u[1].*dn.(t,real(p[2]))
end

function ab_dn_t0(x)
    A = +χ .+ 0.5.*Ω*λ.*(x - xs) .- π/4
    B = -χ .+ 0.5.*Ω*λ.*(x - xs).- π/4
    return [2*im*sin(A); 2*cos(B)]
end

# How do you get rid of this whole loop?
tspan = (0.0,abs(minimum(c.box.t)))
for i in eachindex(x)
    u0 = ab_dn_t0(x[i])
    prob = ODEProblem(dn_ab!, u0, tspan, [λ, m])
    sol = solve(prob, Tsit5(), saveat=t) # t is an array
   # do some processing on sol and save the result
end
# do some more processing on the saved result for all x

You can use the integrator interface with reinit.

Or you can use the ensemble interface to auto-parallelize it.

This includes automated GPU parallelism with DiffEqGPU.jl

4 Likes

Thanks Chris, that resolves it for my use case.

That’s wonderful, I’ll definitely move to this at some point if my system gets large enough.