High allocations count for TwoPointBVProblem in DifferentialEquations.jl


I’m currently trying to port some simulation code from Python to Julia. The Python version uses scipy’s solve_bvp function to solve a two point boundary value problem.

I am using DifferentialEquations.jl to accomplish this task in Julia. It works and the result matches the solve_bvp output just fine. However, the Julia version is much, much slower compared to solve_bvp. I have benchmarked my rhs function and boundary condition to make them as fast as possible, but it seems DifferentialEquations.jl allocates a lot.

Additionally, it seems that TwoPointBVProblems are not correctly handled, as in BoundaryValueDiffEq.jl/solve.jl at master · SciML/BoundaryValueDiffEq.jl · GitHub the isa condition will always result to false, because the constructor of TwoPointBVProblem always sets the field problem_type to StandardBVProblem.

I have boiled down my simulation to a short MWE, which showcases the problem (the MWE is not a real two point BV problem mathematically, but this should not be relevant here):

using DifferentialEquations, BenchmarkTools

a = [1, 2, 3]
b = [1, 1, 1]

function rhs!(du, u, p, t)
    a, _ = p
    for i in axes(u, 1)
        du[i] = u[i] * a[i]

function boundary_condition!(residual, u, p, t)
    _, b = p
    @. residual = u[1] - b

tspan = (0, 1)

u0 = zeros(3, 11)
u0 = [u0[:,i] for i in axes(u0, 2)]
p = [a, b]
bvp = TwoPointBVProblem(rhs!, boundary_condition!, u0, tspan, p)

@benchmark solve($bvp, MIRK4(), dt=0.1, save_everystep=false)


BenchmarkTools.Trial: 898 samples with 1 evaluation.
 Range (min … max):  4.562 ms … 8.521 ms  ┊ GC (min … max):  0.00% … 37.20%
 Time  (median):     4.968 ms             ┊ GC (median):     0.00%
 Time  (mean ± σ):   5.567 ms ± 1.299 ms  ┊ GC (mean ± σ):  11.16% ± 14.98%

  ▁    █▆                                                    
  █▄▅▆▆███▅▃▃▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▃▂▃▂▃▃▃▄▅▄▃ ▃
  4.56 ms        Histogram: frequency by time       8.34 ms <

 Memory estimate: 9.16 MiB, allocs estimate: 115818.

Is there anything I can improve in terms of speed for this kind of problem?

just out of curoisity, have you tried Shooting?

Not really, as it is not the ideal algorithm type for the real problem (aside from the MWE).

The BVP solvers will be going through a complete overhaul this spring and summer, so I’m just going to hold off here and say it’s going to completely change and be fixed by that.

1 Like

Thanks @ChrisRackauckas ! I guess then I will use SciPy via Py(thon)Call for now and wait for the overhaul.