Issue with StaticArrays in Boundary Value Problems


Recently I’ve been using the amazing DifferentialEquations package. However, few days ago I noticed that the package used for Boundary Value Problems BoundaryValueDiffEq has some issues when using Static Arrays. For instance, using the example from the documentation will return a solution with retcode: Failure.
Here’s a MWE:

using StaticArrays, BoundaryValueDiffEq

const g = 9.81
L = 1.0
tspan = (0.0, pi / 2)
function simplependulum!(du, u, p, t)
    θ = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -(g / L) * sin(θ)

function bc2a!(resid_a, u_a, p) # u_a is at the beginning of the time span
    resid_a[1] = u_a[1] + pi / 2 # the solution at the beginning of the time span should be -pi/2
function bc2b!(resid_b, u_b, p) # u_b is at the ending of the time span
    resid_b[1] = u_b[1] - pi / 2 # the solution at the end of the time span should be pi/2

#This is the same problem that appears in the docs, but using 'MVector' for the initial condition.
bvp_SA = TwoPointBVProblem(simplependulum!, (bc2a!, bc2b!), MVector{2}([pi / 2, pi / 2]), tspan;  
                   bcresid_prototype = (zeros(1), zeros(1)))
sol_SA = solve(bvp_SA, MIRK4(), dt = 0.05)
println("Retcode using MVector: ", sol_SA.retcode) #This will printout 'Failure'

Besides, the solution is very different form the correct one (i.e. the one reported in the documentation, obtained using ‘Vector’ types). I also tried using bcresid_prototype = (MVector{1}(zeros(1)), MVector{1}(zeros(1))) when creating the BVProblem, but that didn’t solve the issue.

Finally, trying to spot the problem from the relevant lines (here and here), if I understood correctly, it seems that the solver first returns ‘Failure’ because the defect norm is too large. Then, after several calls to the nonlinear solver, it also happens that the mesh is incompatible with the number of maximum subintervals. I think that’s what triggers the last ‘Failure’, but I’m not completely sure.

I don’t know if I’m doing something wrong when creating or solving the problem. Thanks in advance for any help!

Open an issue on BoundaryValueDiffEq.jl. I think this might just be a limitation in the current MIRK4 implementation that we should track and get tested.