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