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
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] end end function boundary_condition!(residual, u, p, t) _, b = p @. residual = u - b end 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?