Greetings,
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(θ)
end
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
end
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
end
#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!