sciml / boundaryvaluediffeq.jl Goto Github PK

Boundary value problem (BVP) solvers for scientific machine learning (SciML)

License: Other

Julia 100.00%

differentialequations differential-equations bvp neural-bvp neural-ode neural-differential-equations gpu sciml scientific-machine-learning

boundaryvaluediffeq.jl's People

Contributors

Stargazers

Watchers

boundaryvaluediffeq.jl's Issues

DimensionMismatch("dimensions must match") with Shooting()

I get a solution when I solve a problem with

sol1 = solve(bvp1, p=(Cm, Cp));

but when I run

sol1 = solve(bvp1, Shooting(Tsit5()), p=(Cm, Cp));

I get

DimensionMismatch("dimensions must match")

Stacktrace:
 [1] promote_shape at ./indices.jl:154 [inlined]
 [2] promote_shape at ./indices.jl:145 [inlined]
 [3] -(::ODESolution{Float64,2,Core.Array{Core.Array{Float64,1},1},Nothing,Nothing,Core.Array{Float64,1},Core.Array{Core.Array{Core.Array{Float64,1},1},1},BVProblem{Core.Array{Float64,1},Tuple{Float64,Float64},false,Tuple{Core.Array{Float64,1},Core.Array{Float64,1}},ODEFunction{false,typeof(rod_rhs),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(bc1!),DiffEqBase.StandardBVProblem,Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{false,typeof(rod_rhs),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Core.Array{Core.Array{Float64,1},1},Core.Array{Float64,1},Core.Array{Core.Array{Core.Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},DiffEqBase.DEStats}, ::Core.Array{Float64,1}) at ./arraymath.jl:38
 [4] bc1!(::Core.Array{Float64,1}, ::ODESolution{Float64,2,Core.Array{Core.Array{Float64,1},1},Nothing,Nothing,Core.Array{Float64,1},Core.Array{Core.Array{Core.Array{Float64,1},1},1},BVProblem{Core.Array{Float64,1},Tuple{Float64,Float64},false,Tuple{Core.Array{Float64,1},Core.Array{Float64,1}},ODEFunction{false,typeof(rod_rhs),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(bc1!),DiffEqBase.StandardBVProblem,Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{false,typeof(rod_rhs),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Core.Array{Core.Array{Float64,1},1},Core.Array{Float64,1},Core.Array{Core.Array{Core.Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}},DiffEqBase.DEStats}, ::Tuple{Core.Array{Float64,1},Core.Array{Float64,1}}, ::Core.Array{Float64,1}) at ./In[15]:11
 [5] (::getfield(BoundaryValueDiffEq, Symbol("##18#19")){Base.Iterators.Pairs{Symbol,Tuple{Core.Array{Float64,1},Core.Array{Float64,1}},Tuple{Symbol},NamedTuple{(:p,),Tuple{Tuple{Core.Array{Float64,1},Core.Array{Float64,1}}}}},BVProblem{Core.Array{Float64,1},Tuple{Float64,Float64},false,Tuple{Core.Array{Float64,1},Core.Array{Float64,1}},ODEFunction{false,typeof(rod_rhs),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(bc1!),DiffEqBase.StandardBVProblem,Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}},Shooting{Tsit5,getfield(BoundaryValueDiffEq, Symbol("##1#2"))},typeof(bc1!)})(::Core.Array{Float64,1}, ::Core.Array{Float64,1}) at /home/user/.julia/packages/BoundaryValueDiffEq/XvS8h/src/solve.jl:12
 [6] (::getfield(NLSolversBase, Symbol("#fj_finitediff!#21")){getfield(BoundaryValueDiffEq, Symbol("##18#19")){Base.Iterators.Pairs{Symbol,Tuple{Core.Array{Float64,1},Core.Array{Float64,1}},Tuple{Symbol},NamedTuple{(:p,),Tuple{Tuple{Core.Array{Float64,1},Core.Array{Float64,1}}}}},BVProblem{Core.Array{Float64,1},Tuple{Float64,Float64},false,Tuple{Core.Array{Float64,1},Core.Array{Float64,1}},ODEFunction{false,typeof(rod_rhs),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(bc1!),DiffEqBase.StandardBVProblem,Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}},Shooting{Tsit5,getfield(BoundaryValueDiffEq, Symbol("##1#2"))},typeof(bc1!)},DiffEqDiffTools.JacobianCache{Core.Array{Float64,1},Core.Array{Float64,1},Core.Array{Float64,1},UnitRange{Int64},Nothing,Val{:central},Float64}})(::Core.Array{Float64,1}, ::Core.Array{Float64,2}, ::Core.Array{Float64,1}) at /home/user/.julia/packages/NLSolversBase/NsXIC/src/objective_types/oncedifferentiable.jl:137
 [7] value_jacobian!!(::NLSolversBase.OnceDifferentiable{Core.Array{Float64,1},Core.Array{Float64,2},Core.Array{Float64,1}}, ::Core.Array{Float64,1}, ::Core.Array{Float64,2}, ::Core.Array{Float64,1}) at /home/user/.julia/packages/NLSolversBase/NsXIC/src/interface.jl:124
 [8] value_jacobian!! at /home/user/.julia/packages/NLSolversBase/NsXIC/src/interface.jl:122 [inlined]
 [9] trust_region_(::NLSolversBase.OnceDifferentiable{Core.Array{Float64,1},Core.Array{Float64,2},Core.Array{Float64,1}}, ::Core.Array{Float64,1}, ::Float64, ::Float64, ::Int64, ::Bool, ::Bool, ::Bool, ::Float64, ::Bool, ::NLsolve.NewtonTrustRegionCache{Core.Array{Float64,1}}) at /home/user/.julia/packages/NLsolve/W6Ptc/src/solvers/trust_region.jl:119
 [10] trust_region at /home/user/.julia/packages/NLsolve/W6Ptc/src/solvers/trust_region.jl:229 [inlined] (repeats 2 times)
 [11] #nlsolve#25(::Symbol, ::Float64, ::Float64, ::Int64, ::Bool, ::Bool, ::Bool, ::LineSearches.Static, ::getfield(NLsolve, Symbol("##31#33")), ::Float64, ::Bool, ::Int64, ::Int64, ::Int64, ::Float64, ::typeof(NLsolve.nlsolve), ::NLSolversBase.OnceDifferentiable{Core.Array{Float64,1},Core.Array{Float64,2},Core.Array{Float64,1}}, ::Core.Array{Float64,1}) at /home/user/.julia/packages/NLsolve/W6Ptc/src/nlsolve/nlsolve.jl:26
 [12] (::getfield(NLsolve, Symbol("#kw##nlsolve")))(::NamedTuple{(:method, :xtol, :ftol, :iterations, :store_trace, :show_trace, :extended_trace, :linesearch, :factor, :autoscale, :m, :beta, :aa_start, :droptol, :linsolve),Tuple{Symbol,Float64,Float64,Int64,Bool,Bool,Bool,LineSearches.Static,Float64,Bool,Int64,Int64,Int64,Float64,getfield(NLsolve, Symbol("##31#33"))}}, ::typeof(NLsolve.nlsolve), ::NLSolversBase.OnceDifferentiable{Core.Array{Float64,1},Core.Array{Float64,2},Core.Array{Float64,1}}, ::Core.Array{Float64,1}) at ./none:0
 [13] #nlsolve#30(::Symbol, ::Float64, ::Float64, ::Int64, ::Bool, ::Bool, ::Bool, ::LineSearches.Static, ::Float64, ::Bool, ::Int64, ::Int64, ::Int64, ::Float64, ::Symbol, ::getfield(NLsolve, Symbol("##31#33")), ::Bool, ::typeof(NLsolve.nlsolve), ::getfield(BoundaryValueDiffEq, Symbol("##18#19")){Base.Iterators.Pairs{Symbol,Tuple{Core.Array{Float64,1},Core.Array{Float64,1}},Tuple{Symbol},NamedTuple{(:p,),Tuple{Tuple{Core.Array{Float64,1},Core.Array{Float64,1}}}}},BVProblem{Core.Array{Float64,1},Tuple{Float64,Float64},false,Tuple{Core.Array{Float64,1},Core.Array{Float64,1}},ODEFunction{false,typeof(rod_rhs),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(bc1!),DiffEqBase.StandardBVProblem,Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}},Shooting{Tsit5,getfield(BoundaryValueDiffEq, Symbol("##1#2"))},typeof(bc1!)}, ::Core.Array{Float64,1}) at /home/user/.julia/packages/NLsolve/W6Ptc/src/nlsolve/nlsolve.jl:73
 [14] nlsolve(::Function, ::Core.Array{Float64,1}) at /home/user/.julia/packages/NLsolve/W6Ptc/src/nlsolve/nlsolve.jl:59
 [15] #1 at /home/user/.julia/packages/BoundaryValueDiffEq/XvS8h/src/algorithms.jl:11 [inlined]
 [16] #__solve#17(::Base.Iterators.Pairs{Symbol,Tuple{Core.Array{Float64,1},Core.Array{Float64,1}},Tuple{Symbol},NamedTuple{(:p,),Tuple{Tuple{Core.Array{Float64,1},Core.Array{Float64,1}}}}}, ::typeof(DiffEqBase.__solve), ::BVProblem{Core.Array{Float64,1},Tuple{Float64,Float64},false,Tuple{Core.Array{Float64,1},Core.Array{Float64,1}},ODEFunction{false,typeof(rod_rhs),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(bc1!),DiffEqBase.StandardBVProblem,Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}, ::Shooting{Tsit5,getfield(BoundaryValueDiffEq, Symbol("##1#2"))}) at /home/user/.julia/packages/BoundaryValueDiffEq/XvS8h/src/solve.jl:15
 [17] #__solve at ./none:0 [inlined]
 [18] #solve_call#433(::Base.Iterators.Pairs{Symbol,Tuple{Core.Array{Float64,1},Core.Array{Float64,1}},Tuple{Symbol},NamedTuple{(:p,),Tuple{Tuple{Core.Array{Float64,1},Core.Array{Float64,1}}}}}, ::typeof(DiffEqBase.solve_call), ::BVProblem{Core.Array{Float64,1},Tuple{Float64,Float64},false,Tuple{Core.Array{Float64,1},Core.Array{Float64,1}},ODEFunction{false,typeof(rod_rhs),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(bc1!),DiffEqBase.StandardBVProblem,Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}, ::Shooting{Tsit5,getfield(BoundaryValueDiffEq, Symbol("##1#2"))}) at /home/user/.julia/packages/DiffEqBase/IDJcQ/src/solve.jl:38
 [19] #solve_call at ./none:0 [inlined]
 [20] #solve#434(::Base.Iterators.Pairs{Symbol,Tuple{Core.Array{Float64,1},Core.Array{Float64,1}},Tuple{Symbol},NamedTuple{(:p,),Tuple{Tuple{Core.Array{Float64,1},Core.Array{Float64,1}}}}}, ::typeof(solve), ::BVProblem{Core.Array{Float64,1},Tuple{Float64,Float64},false,DiffEqBase.NullParameters,ODEFunction{false,typeof(rod_rhs),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(bc1!),DiffEqBase.StandardBVProblem,Base.Iterators.Pairs{Symbol,Tuple{Core.Array{Float64,1},Core.Array{Float64,1}},Tuple{Symbol},NamedTuple{(:p,),Tuple{Tuple{Core.Array{Float64,1},Core.Array{Float64,1}}}}}}, ::Shooting{Tsit5,getfield(BoundaryValueDiffEq, Symbol("##1#2"))}) at /home/user/.julia/packages/DiffEqBase/IDJcQ/src/solve.jl:63
 [21] (::getfield(DiffEqBase, Symbol("#kw##solve")))(::NamedTuple{(:p,),Tuple{Tuple{Core.Array{Float64,1},Core.Array{Float64,1}}}}, ::typeof(solve), ::BVProblem{Core.Array{Float64,1},Tuple{Float64,Float64},false,DiffEqBase.NullParameters,ODEFunction{false,typeof(rod_rhs),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(bc1!),DiffEqBase.StandardBVProblem,Base.Iterators.Pairs{Symbol,Tuple{Core.Array{Float64,1},Core.Array{Float64,1}},Tuple{Symbol},NamedTuple{(:p,),Tuple{Tuple{Core.Array{Float64,1},Core.Array{Float64,1}}}}}}, ::Shooting{Tsit5,getfield(BoundaryValueDiffEq, Symbol("##1#2"))}) at ./none:0
 [22] top-level scope at In[19]:1

Which dimensions do not match if they match in the 1st case? Any workaround for this?

Add nonstandard array support

It seems that BVPSystem requires plain old Vectors. Any way to lift this restriction for ComponentArrays and other nonstandard arrays?

Adaptivity in MIRK solver

(I'm not sure to what extent this overlaps with #12 which mentions defects)

The possibility of having this was mentioned in this blog post:

Adaptivity in the MIRK BVP solvers

Lack of adaptivity strongly limits the usability in any real use. (For example, I currently rely on Mathematica for BVPs.)

Any hopes we will have adaptivity in a near future?

Solution interpolation in BC with MIRK solver not working?

I follow the simple pendulum example 3

using BoundaryValueDiffEq
using Plots
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

using OrdinaryDiffEq
u₀_2 = [-1.6, -1.7] # the initial guess
function bc3!(residual, sol, p, t)
    residual[1] = sol(pi / 4)[1] + pi / 2 # use the interpolation here, since indexing will be wrong for adaptive methods
    residual[2] = sol(pi / 2)[1] - pi / 2
end

bvp3 = BVProblem(simplependulum!, bc3!, u₀_2, tspan)
sol3 = solve(bvp3, Shooting(Vern7()))

which works as expected.

However, when I try to solve with MIRK4(), the solution interpolation in the BC does not seem to work. Is this expected?

sol3 = solve(bvp3, MIRK4(), dt=0.05)

MethodError: objects of type Vector{Vector{Float64}} are not callable
Use square brackets [] for indexing an Array.

Issue #107 looks to me as if it should be working?

Pkg.status()

  [764a87c0] BoundaryValueDiffEq v5.4.0
  [1dea7af3] OrdinaryDiffEq v6.59.2
  [91a5bcdd] Plots v1.39.0
  [aa65fe97] SnoopCompile v2.10.8

How to input a tolerance in BVP

Is there a way to pass the ftol argument to the BVProblem? As an example, say that I want the residuals of this problem to be zero up to 25 decimals:

using BoundaryValueDiffEq, OrdinaryDiffEq
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/p[1])*sin(θ)
end
u₀_2 = BigFloat.(["-0.47", "-4.79"]) # the initial guess
tspan = BigFloat.(tspan)
function bc3!(residual, sol, p, t)
    residual[1] = sol(pi/4)[1] + pi/2 # use the interpolation here, since indexing will be wrong for adaptive methods
    residual[2] = sol(pi/2)[1] - pi/2
end
bvp3 = BVProblem(simplependulum!, bc3!, u₀_2, tspan, [L])
sol3 = solve(bvp3, Shooting(Feagin14()), abstol = 1e-40, reltol = 1e-35)

@show sol3(pi/2)[1] - pi/2

How can I do that?

Continuous Benchmarking

Add more serious workloads
Add MIRK Methods

Free-Final-Time TPBVP/generalized BVP

A common problem in optimal control is when you want to achieve some objective state but do not have a specific time by which that objective must be attained. One example of such a problem would be rocket landing guidance ("we want to land at a point without running out of fuel, but don't care about how long it takes"). A similar problem for a generalized BVP with intermediate state constraints is when the system needs to attain an intermediate point at some time, such as to go through a waypoin.

Right now, as far as I'm aware, this needs to be ``baked into'' the problem; you can build in a time dilation parameter into your model that maps a virtual time as proposed to BVP to real-world time which can then be optimized over by the solver. This is not general, though, and doesn't take advantage of solver adaptation.

It would be useful then to have an interface by which free-final-time could be specified in a BVProblem/TPBVProblem and then be solved taking advantage of adaptation.

Cleanup test files

https://github.com/JuliaDiffEq/BoundaryValueDiffEq.jl/blob/master/test/runtests.jl#L7

These shooting method parts should probably be their own file, which is also included.

Also,

https://github.com/JuliaDiffEq/BoundaryValueDiffEq.jl/blob/master/test/runtests.jl#L5

I prefer to wrap each included file into a @testset. This will make it run the entire testset, instead of erroring at the first failed test in the file. It also gives a very nice named printout to make it easier to follow what's going on. Example:

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/test/runtests.jl

(though don't follow the ode folder thing: that's a relic of the past that needs to be dealt with... you can group related tests in folders, but you don't have enough tests for that to matter quite yet).

Utilize reinit fast path in shooting

Regression with (General)MIRK4 on BoundaryValueDiffEq.jl v4.0

The following code works on BoundaryValueDiffEq.jl v3.0:

using Manifolds, BoundaryValueDiffEq
using Manifolds: affine_connection

M = Manifolds.EmbeddedTorus(3, 2)
A = Manifolds.DefaultTorusAtlas()

a1 = [0.5, -1.2]
a2 = [-0.5, 0.3]
i = (0, 0)
solver=GeneralMIRK4()
dt=0.05
tspan = (0.0, 1.0)
function bc1!(residual, u, p, t)
    mid = div(length(u[1]), 2)
    residual[1:mid] = u[1][1:mid] - a1
    return residual[(mid + 1):end] = u[end][1:mid] - a2
end
function chart_log_problem!(du, u, params, t)
    M, A, i = params
    mid = div(length(u), 2)
    a = u[1:mid]
    dx = u[(mid + 1):end]
    ddx = -affine_connection(M, A, i, a, dx, dx)
    du[1:mid] .= dx
    du[(mid + 1):end] .= ddx
    return du
end
bvp1 = BVProblem(
    chart_log_problem!,
    bc1!,
    (p, t) -> vcat(t * a1 + (1 - t) * a2, zero(a1)),
    tspan,
    (M, A, i),
)
sol1 = solve(bvp1, solver, dt=dt)

I tried upgrading it to the recently released version v4.0 by replacing solver=GeneralMIRK4() with solver=MIRK4() but I get the following error:

ERROR: DimensionMismatch: new dimensions (1, 40) must be consistent with array size 20
Stacktrace:
  [1] (::Base.var"#throw_dmrsa#318")(dims::Tuple{Int64, Int64}, len::Int64)
    @ Base ./reshapedarray.jl:41
  [2] reshape
    @ ./reshapedarray.jl:45 [inlined]
  [3] reshape
    @ ./reshapedarray.jl:117 [inlined]
  [4] interp_setup(S::BoundaryValueDiffEq.BVPSystem{ODEFunction{true, SciMLBase.FullSpecialize, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ManifoldsBoundaryValueDiffEqExt.var"#bc1!#3"{Vector{Float64}, Vector{Float64}}, Vector{Float64}, Vector{Float64}, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}}, mesh::Vector{Float64}, y::Matrix{Float64}, ITU::BoundaryValueDiffEq.MIRKInterpTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}, k_discrete::Array{Float64, 3}, p::Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, dt_::Vector{Float64})
    @ BoundaryValueDiffEq ~/.julia/dev/BoundaryValueDiffEq/src/adaptivity.jl:209
  [5] defect_estimate(S::BoundaryValueDiffEq.BVPSystem{ODEFunction{true, SciMLBase.FullSpecialize, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ManifoldsBoundaryValueDiffEqExt.var"#bc1!#3"{Vector{Float64}, Vector{Float64}}, Vector{Float64}, Vector{Float64}, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}}, cache::BoundaryValueDiffEq.MIRK4Cache{Array{Float64, 3}}, alg::MIRK4{NewtonRaphson{0, true, Val{:forward}, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), true, nothing}}, ITU::BoundaryValueDiffEq.MIRKInterpTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}, y::Matrix{Float64}, dt::Vector{Float64})
    @ BoundaryValueDiffEq ~/.julia/dev/BoundaryValueDiffEq/src/adaptivity.jl:166
  [6] __solve(prob::BVProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ManifoldsBoundaryValueDiffEqExt.var"#bc1!#3"{Vector{Float64}, Vector{Float64}}, SciMLBase.StandardBVProblem, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, alg::MIRK4{NewtonRaphson{0, true, Val{:forward}, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), true, nothing}}; dt::Float64, abstol::Float64, adaptive::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ BoundaryValueDiffEq ~/.julia/dev/BoundaryValueDiffEq/src/solve.jl:59
  [7] __solve
    @ ~/.julia/dev/BoundaryValueDiffEq/src/solve.jl:33 [inlined]
  [8] solve_call(_prob::BVProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ManifoldsBoundaryValueDiffEqExt.var"#bc1!#3"{Vector{Float64}, Vector{Float64}}, SciMLBase.StandardBVProblem, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, args::MIRK4{NewtonRaphson{0, true, Val{:forward}, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), true, nothing}}; merge_callbacks::Bool, kwargshandle::KeywordArgError, kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:dt,), Tuple{Float64}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/rVTlI/src/solve.jl:511
  [9] solve_up(prob::BVProblem{ManifoldsBoundaryValueDiffEqExt.var"#2#4"{Vector{Float64}, Vector{Float64}}, Tuple{Float64, Float64}, true, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, ODEFunction{true, SciMLBase.FullSpecialize, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, ManifoldsBoundaryValueDiffEqExt.var"#bc1!#3"{Vector{Float64}, Vector{Float64}}, SciMLBase.StandardBVProblem, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, sensealg::Nothing, u0::Function, p::Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, args::MIRK4{NewtonRaphson{0, true, Val{:forward}, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), true, nothing}}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:dt,), Tuple{Float64}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/rVTlI/src/solve.jl:972
 [10] solve_up
    @ ~/.julia/packages/DiffEqBase/rVTlI/src/solve.jl:945 [inlined]
 [11] #solve#40
    @ ~/.julia/packages/DiffEqBase/rVTlI/src/solve.jl:882 [inlined]

I tried debugging it but I couldn't see what may be wrong.

Invalid Bandind error / Manifolds.jl

BoundaryValueDiffEq.jl is again causing test failures in Manifolds.jl:

  Invalid Bandind
  Stacktrace:
    [1] _bandsize
      @ ~/.julia/packages/ArrayInterface/Ic7Rb/ext/ArrayInterfaceBandedMatricesExt.jl:46 [inlined]
    [2] #1
      @ ./none:0 [inlined]
    [3] iterate(g::Base.Generator, s::Vararg{Any})
      @ Base ./generator.jl:47 [inlined]
    [4] collect(itr::Base.Generator{StepRange{Int64, Int64}, ArrayInterfaceBandedMatricesExt.var"#1#2"{Int64, Int64}})
      @ Base ./array.jl:834
    [5] ArrayInterface.BandedMatrixIndex(rowsize::Int64, colsize::Int64, lowerbandwidth::Int64, upperbandwidth::Int64, isrow::Bool)
      @ ArrayInterfaceBandedMatricesExt ~/.julia/packages/ArrayInterface/Ic7Rb/ext/ArrayInterfaceBandedMatricesExt.jl:59
    [6] findstructralnz
      @ ~/.julia/packages/ArrayInterface/Ic7Rb/ext/ArrayInterfaceBandedMatricesExt.jl:66 [inlined]
    [7] #_#94
      @ ~/.julia/packages/SparseDiffTools/PsgQF/src/highlevel/coloring.jl:29 [inlined]
    [8] PrecomputedJacobianColorvec
      @ ~/.julia/packages/SparseDiffTools/PsgQF/src/highlevel/coloring.jl:26 [inlined]
    [9] sparse_jacobian_cache(ad::AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, sd::PrecomputedJacobianColorvec{BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, Vector{Int64}, Vector{Int64}}, f!::BoundaryValueDiffEq.var"#135#141"{BoundaryValueDiffEq.MIRKCache{true, Float64, Tuple{Int64}, BVPFunction{true, SciMLBase.FullSpecialize, false, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, BVProblem{Vector{Vector{Float64}}, Tuple{Float64, Float64}, true, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, BVPFunction{true, SciMLBase.FullSpecialize, false, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, SciMLBase.StandardBVProblem, @Kwargs{}}, SciMLBase.StandardBVProblem, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, MIRK4{Nothing, BVPJacobianAlgorithm{AutoForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}}}, BoundaryValueDiffEq.MIRKTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Matrix{Float64}}, BoundaryValueDiffEq.MIRKInterpTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{PreallocationTools.DiffCache{Matrix{Float64}, Vector{Float64}}}, Vector{Matrix{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, Vector{Vector{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}, Vector{Float64}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Tuple{Int64}, @NamedTuple{defect_threshold::Float64, MxNsub::Int64, abstol::Float64, dt::Float64, adaptive::Bool}}, BoundaryValueDiffEq.var"#117#123"{BoundaryValueDiffEq.MIRKCache{true, Float64, Tuple{Int64}, BVPFunction{true, SciMLBase.FullSpecialize, false, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, BVProblem{Vector{Vector{Float64}}, Tuple{Float64, Float64}, true, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, BVPFunction{true, SciMLBase.FullSpecialize, false, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, SciMLBase.StandardBVProblem, @Kwargs{}}, SciMLBase.StandardBVProblem, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, MIRK4{Nothing, BVPJacobianAlgorithm{AutoForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}}}, BoundaryValueDiffEq.MIRKTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Matrix{Float64}}, BoundaryValueDiffEq.MIRKInterpTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{PreallocationTools.DiffCache{Matrix{Float64}, Vector{Float64}}}, Vector{Matrix{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, Vector{Vector{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}, Vector{Float64}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Tuple{Int64}, @NamedTuple{defect_threshold::Float64, MxNsub::Int64, abstol::Float64, dt::Float64, adaptive::Bool}}}}, fx::Vector{Float64}, x::Vector{Float64})
      @ SparseDiffTools ~/.julia/packages/SparseDiffTools/PsgQF/src/highlevel/forward_mode.jl:41
   [10] __sparse_jacobian_cache
      @ ~/.julia/packages/BoundaryValueDiffEq/ComOY/src/sparse_jacobians.jl:17 [inlined]
   [11] __construct_nlproblem(cache::BoundaryValueDiffEq.MIRKCache{true, Float64, Tuple{Int64}, BVPFunction{true, SciMLBase.FullSpecialize, false, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, BVProblem{Vector{Vector{Float64}}, Tuple{Float64, Float64}, true, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, BVPFunction{true, SciMLBase.FullSpecialize, false, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, SciMLBase.StandardBVProblem, @Kwargs{}}, SciMLBase.StandardBVProblem, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, MIRK4{Nothing, BVPJacobianAlgorithm{AutoForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}}}, BoundaryValueDiffEq.MIRKTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Matrix{Float64}}, BoundaryValueDiffEq.MIRKInterpTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{PreallocationTools.DiffCache{Matrix{Float64}, Vector{Float64}}}, Vector{Matrix{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, Vector{Vector{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}, Vector{Float64}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Tuple{Int64}, @NamedTuple{defect_threshold::Float64, MxNsub::Int64, abstol::Float64, dt::Float64, adaptive::Bool}}, y::Vector{Float64}, loss_bc::BoundaryValueDiffEq.var"#115#121"{BoundaryValueDiffEq.MIRKCache{true, Float64, Tuple{Int64}, BVPFunction{true, SciMLBase.FullSpecialize, false, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, BVProblem{Vector{Vector{Float64}}, Tuple{Float64, Float64}, true, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, BVPFunction{true, SciMLBase.FullSpecialize, false, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, SciMLBase.StandardBVProblem, @Kwargs{}}, SciMLBase.StandardBVProblem, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, MIRK4{Nothing, BVPJacobianAlgorithm{AutoForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}}}, BoundaryValueDiffEq.MIRKTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Matrix{Float64}}, BoundaryValueDiffEq.MIRKInterpTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{PreallocationTools.DiffCache{Matrix{Float64}, Vector{Float64}}}, Vector{Matrix{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, Vector{Vector{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}, Vector{Float64}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Tuple{Int64}, @NamedTuple{defect_threshold::Float64, MxNsub::Int64, abstol::Float64, dt::Float64, adaptive::Bool}}, SciMLBase.StandardBVProblem}, loss_collocation::BoundaryValueDiffEq.var"#117#123"{BoundaryValueDiffEq.MIRKCache{true, Float64, Tuple{Int64}, BVPFunction{true, SciMLBase.FullSpecialize, false, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, BVProblem{Vector{Vector{Float64}}, Tuple{Float64, Float64}, true, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, BVPFunction{true, SciMLBase.FullSpecialize, false, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, SciMLBase.StandardBVProblem, @Kwargs{}}, SciMLBase.StandardBVProblem, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, MIRK4{Nothing, BVPJacobianAlgorithm{AutoForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}}}, BoundaryValueDiffEq.MIRKTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Matrix{Float64}}, BoundaryValueDiffEq.MIRKInterpTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{PreallocationTools.DiffCache{Matrix{Float64}, Vector{Float64}}}, Vector{Matrix{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, Vector{Vector{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}, Vector{Float64}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Tuple{Int64}, @NamedTuple{defect_threshold::Float64, MxNsub::Int64, abstol::Float64, dt::Float64, adaptive::Bool}}}, loss::BoundaryValueDiffEq.var"#119#125"{BoundaryValueDiffEq.MIRKCache{true, Float64, Tuple{Int64}, BVPFunction{true, SciMLBase.FullSpecialize, false, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, BVProblem{Vector{Vector{Float64}}, Tuple{Float64, Float64}, true, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, BVPFunction{true, SciMLBase.FullSpecialize, false, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, SciMLBase.StandardBVProblem, @Kwargs{}}, SciMLBase.StandardBVProblem, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, MIRK4{Nothing, BVPJacobianAlgorithm{AutoForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}}}, BoundaryValueDiffEq.MIRKTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Matrix{Float64}}, BoundaryValueDiffEq.MIRKInterpTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{PreallocationTools.DiffCache{Matrix{Float64}, Vector{Float64}}}, Vector{Matrix{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, Vector{Vector{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}, Vector{Float64}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Tuple{Int64}, @NamedTuple{defect_threshold::Float64, MxNsub::Int64, abstol::Float64, dt::Float64, adaptive::Bool}}, SciMLBase.StandardBVProblem}, ::SciMLBase.StandardBVProblem)
      @ BoundaryValueDiffEq ~/.julia/packages/BoundaryValueDiffEq/ComOY/src/solve/mirk.jl:315
   [12] __construct_nlproblem
      @ BoundaryValueDiffEq ~/.julia/packages/BoundaryValueDiffEq/ComOY/src/solve/mirk.jl:214 [inlined]
   [13] solve!(cache::BoundaryValueDiffEq.MIRKCache{true, Float64, Tuple{Int64}, BVPFunction{true, SciMLBase.FullSpecialize, false, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, BVProblem{Vector{Vector{Float64}}, Tuple{Float64, Float64}, true, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, BVPFunction{true, SciMLBase.FullSpecialize, false, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, SciMLBase.StandardBVProblem, @Kwargs{}}, SciMLBase.StandardBVProblem, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, MIRK4{Nothing, BVPJacobianAlgorithm{AutoForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}}}, BoundaryValueDiffEq.MIRKTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Matrix{Float64}}, BoundaryValueDiffEq.MIRKInterpTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{PreallocationTools.DiffCache{Matrix{Float64}, Vector{Float64}}}, Vector{Matrix{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, Vector{Vector{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}, Vector{Float64}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Tuple{Int64}, @NamedTuple{defect_threshold::Float64, MxNsub::Int64, abstol::Float64, dt::Float64, adaptive::Bool}})
      @ BoundaryValueDiffEq ~/.julia/packages/BoundaryValueDiffEq/ComOY/src/solve/mirk.jl:144
   [14] #__solve#157
      @ ~/.julia/packages/BoundaryValueDiffEq/ComOY/src/BoundaryValueDiffEq.jl:47 [inlined]
   [15] __solve
      @ ~/.julia/packages/BoundaryValueDiffEq/ComOY/src/BoundaryValueDiffEq.jl:45 [inlined]
   [16] #solve_call#34
      @ ~/.julia/packages/DiffEqBase/a6p43/src/solve.jl:561 [inlined]
   [17] solve_call
      @ ~/.julia/packages/DiffEqBase/a6p43/src/solve.jl:527 [inlined]
   [18] solve_up(prob::BVProblem{Vector{Vector{Float64}}, Tuple{Float64, Float64}, true, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, BVPFunction{true, SciMLBase.FullSpecialize, false, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), ManifoldsBoundaryValueDiffEqExt.var"#bc1!#2"{Vector{Float64}, Vector{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, SciMLBase.StandardBVProblem, @Kwargs{}}, sensealg::Nothing, u0::Vector{Vector{Float64}}, p::Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, args::MIRK4{Nothing, BVPJacobianAlgorithm{Nothing, Nothing, Nothing}}; kwargs::@Kwargs{dt::Float64})
      @ DiffEqBase ~/.julia/packages/DiffEqBase/a6p43/src/solve.jl:1010
   [19] solve_up
      @ DiffEqBase ~/.julia/packages/DiffEqBase/a6p43/src/solve.jl:996 [inlined]
   [20] #solve#40
      @ DiffEqBase ~/.julia/packages/DiffEqBase/a6p43/src/solve.jl:933 [inlined]
   [21] solve_chart_log_bvp(M::Manifolds.EmbeddedTorus{Int64}, a1::Vector{Float64}, a2::Vector{Float64}, A::Manifolds.DefaultTorusAtlas, i::Tuple{Int64, Int64}; solver::MIRK4{Nothing, BVPJacobianAlgorithm{Nothing, Nothing, Nothing}}, dt::Float64, kwargs::@Kwargs{})
      @ ManifoldsBoundaryValueDiffEqExt ~/work/Manifolds.jl/Manifolds.jl/ext/ManifoldsBoundaryValueDiffEqExt.jl:70
   [22] solve_chart_log_bvp(M::Manifolds.EmbeddedTorus{Int64}, a1::Vector{Float64}, a2::Vector{Float64}, A::Manifolds.DefaultTorusAtlas, i::Tuple{Int64, Int64})
      @ ManifoldsBoundaryValueDiffEqExt ~/work/Manifolds.jl/Manifolds.jl/ext/ManifoldsBoundaryValueDiffEqExt.jl:51
   [23] macro expansion
      @ ~/work/Manifolds.jl/Manifolds.jl/test/manifolds/embedded_torus.jl:92 [inlined]
   [24] macro expansion
      @ /opt/hostedtoolcache/julia/1.10.0-rc2/x64/share/julia/stdlib/v1.10/Test/src/Test.jl:1577 [inlined]
   [25] top-level scope
      @ ~/work/Manifolds.jl/Manifolds.jl/test/manifolds/embedded_torus.jl:11
   [26] include
      @ ./client.jl:489 [inlined]
   [27] macro expansion
      @ ./timing.jl:279 [inlined]
   [28] include_test(path::String)
      @ Main ~/work/Manifolds.jl/Manifolds.jl/test/utils.jl:27
   [29] macro expansion
      @ ~/work/Manifolds.jl/Manifolds.jl/test/runtests.jl:191 [inlined]
   [30] macro expansion
      @ /opt/hostedtoolcache/julia/1.10.0-rc2/x64/share/julia/stdlib/v1.10/Test/src/Test.jl:1577 [inlined]
   [31] top-level scope
      @ ~/work/Manifolds.jl/Manifolds.jl/test/runtests.jl:11
   [32] include(fname::String)
      @ Base.MainInclude ./client.jl:489
   [33] top-level scope
      @ none:6
   [34] eval
      @ Core ./boot.jl:385 [inlined]
   [35] exec_options(opts::Base.JLOptions)
      @ Base ./client.jl:291
   [36] _start()
      @ Base ./client.jl:552

See here: https://github.com/JuliaManifolds/Manifolds.jl/actions/runs/7211150273/job/19646434007?pr=692 .
I could write an MWE but at this point I think our use case should be just added to your test suite? It keeps failing so often and it's not too complex.

Defect control adaptivity

http://epubs.siam.org/doi/pdf/10.1137/S1064827593251496
http://cs.stmarys.ca/~muir/ShampineMuirXu2006.pdf
http://epubs.siam.org/doi/pdf/10.1137/0912052

The next step should be adaptivity. A really nice way to handle the interpolations in a way that can be used internally is to simply add them to the solution object, and use that solution object on the inside. You can actually see this done in OrdinaryDiffEq.jl here:

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/solve.jl#L250

So the steps for a defect control are:

Add a way to call an interpolation scheme on the sol type (which is algorithm-dependent). You might want to match OrdinaryDiffEq.jl for this, though it's quite complicated. We should discuss this.
Change the timing that the solution is built so it's built before the rest of the routine, and your solver routine just updates the sol.u and sol.t parts (this should actually require no change if the pointers are setup correctly!)
Then, the parts of the adaptive algorithm which require using the continuous extension for the defect, you can just call sol(t) (just like a user would).

The nice thing about this is that if you get the interpolation on the sol type right, then it'll work both internally and for users just the same. Also, this means you can pass the "current sol" to the bc!, making it exactly the same as the Shooting setup and allowing conditions on the interpolation.

saveat doesn't work when specifying a negative part of tspan

When solving a BVP over [-1, 1] (ie tspan = (-1.0, 1.0) ) and specifying a saveat=0.01, only the values of x (ie t) over 0.0 are stored. This can be seen in the following script drop.jl. This is with the development version of BoundaryValueDiffEq and the following version of Julia:

Version 0.7.0-DEV.66 (2017-05-08 17:16 UTC)
Commit 25f241c* (2 days old master)

drop.jl

MATLAB Guide, Higham and Higham example on p. 214 (Third Edition)

using BoundaryValueDiffEq
using DiffEqBase, OrdinaryDiffEq
using PyPlot

function drop_d!(x, y, dy)
dy[1] = y[2]
dy[2] = (y[1] - 1) * ((1 + y[2]^2)^(3/2))
end

function bc!(resid, sol)
resid[1] = sol[1][1]
resid[2] = sol[end][1]
end

this will show the problem

tspan = (-1.0, 1.0)

use the following and it works

#tspan = (0.0, 2.0)

u0 = [0., 1.]
bvp = BVProblem(drop_d!, bc!, u0, tspan)
resid_f = Array(Float64, 2)

sol = solve(bvp, Shooting(Tsit5()))
bc!(resid_f, sol)
println("Norm of residual: ", norm(resid_f))

sol2 = solve(bvp, Shooting(Tsit5()), saveat=0.01)
bc!(resid_f, sol2)
println("Norm of residual: ", norm(resid_f))

plot(sol.t, sol[1, :], label="No saveat")
plot(sol2.t, sol2[1, :], "*", label="Using saveat")
legend()
title("Water droplet BVP solved by BVProblem")
axis("equal")

Passing parameters to the residuals method

It would be convenient to have access to the parameters of the differential equations within the residuals method ( as well as the timespan ).
In the optimal case even for the initial value vector ( 3rd parameter of TwoPointBVProblem ).

True banded matrix support

Open this issue for keeping track of upstream issue/PR. JuliaNLSolvers/NLsolve.jl#113 & JuliaNLSolvers/NLsolve.jl#125.

Functions in Types

Every function has its own type. But callable types are not <:Function. So you might want to let

https://github.com/JuliaDiffEq/BoundaryValueDiffEq.jl/blob/master/src/BoundaryValueDiffEq.jl#L38

those roam free without a restriction.

Optimization-based shooting solutions

Using the parameter estimation functionality:

http://docs.juliadiffeq.org/latest/analysis/parameter_estimation.html#Parameter-Estimation-1

we can solve BVPs with more equations than solutions by changing problem_new_parameters

http://docs.juliadiffeq.org/latest/analysis/parameter_estimation.html#The-Problem-Generator-Function-1

to change the initial condition (and maybe some parameters) as part of the optimization, and then

http://docs.juliadiffeq.org/latest/analysis/parameter_estimation.html#Note-About-Loss-Functions-1

defining a loss function which is a sum of squares of the BVP condition residuals.

This would allow one to solve BVPs with more conditions than equations in a way that simply minimizes sum squared error.

Autoselect Correct Solver for Nonlinear Least Squares

Currently, we fix the solver at the time of algorithm creation. This means for a NLS problem our defaults are incorrect.

Easy solution:

Keep nlsolve as nothing
Run a concretize step
- If not nothing then return the solver
- If nothing choose the solver based on the generated problem

More TODOs

With #89 finally merged, let's push BoundaryValueDiffEq.jl to the next level!🎉🎉🎉🎉

Here are some issues we need to address in the near future:

The code style is a little weird, we need to use a better format style for the package
Currently, the solving speed is very slow and didn't use any special tricks for performance, we need to speed up the defect estimation process
The code is not optimized, we need to do a lot of refactors to make it match the DiffEq interface
Add benchmarks for the existing BVP solvers

GE control

http://cs.smu.ca/tech_reports/txt2012_001.pdf

Very long compilation time for solve()

I've upgraded to v5.1 today and found very long compilation times when trying the simple pendulum example.

My base environment:
[6e4b80f9] BenchmarkTools v1.3.2
[7073ff75] IJulia v1.24.2
[295af30f] Revise v3.5.6
[0f1e0344] WebIO v0.8.21
[de0858da] Printf

I created a minimal environment for testing:
[764a87c0] BoundaryValueDiffEq v5.1.0
[91a5bcdd] Plots v1.39.0

I'm using Julia 1.9.3.

I run "Example 1: Simple pendulum" from the docs and find

@time sol1 = solve(bvp1, MIRK4(), dt = 0.05)
225.577779 seconds (289.71 M allocations: 60.160 GiB, 6.70% gc time, 99.98% compilation time)

Second execution is fast:
@time sol1 = solve(bvp1, MIRK4(), dt = 0.05)
0.001343 seconds (6.28 k allocations: 879.165 KiB)

I'm quite confident this did not happend when I played around 1-2 weeks ago with an earlier version. I tried reverting to v5.0.0, but I'm not experienced in doing stuff like this and obtained an error message upon solve(), so I can't tell whether this also happens for me with v5.0.0 now.

I managed to revert to v4.0.1
⌃ [764a87c0] BoundaryValueDiffEq v4.0.1

And I get more reasonable compilation times:
@time sol1 = solve(bvp1, MIRK4(), dt = 0.05)
28.803127 seconds (42.95 M allocations: 6.039 GiB, 6.34% gc time, 99.76% compilation time)

DiffEq solution and testing tools

building the DiffEq solution isn't documented anywhere, but you can see it in action here:

https://github.com/JuliaDiffEq/Sundials.jl/blob/master/src/common.jl#L471

and it's defined in DiffEqBase:

https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/solutions/ode_solutions.jl#L17

Notice that it has a special check for the dispatch f(Val{:analytic},t,u0) which calculates the error at time t given an initial condition u0. You'd probably need to update u0 at the end for that to work, or set that to have a separate path which uses the first timeseries value.

But once that's setup, all of the testing tools will work. The testing tools are poorly documented here:

http://devdocs.juliadiffeq.org/latest/index.html

or just look at the code here:

https://github.com/JuliaDiffEq/DiffEqDevTools.jl

Sensitivity Analysis of BVProblems

I am trying to implement the "shooting method" to solve a simple differential equation (\ddot{x} = -x*(1-x^2) ) as a boundary value problem. I also need to optimize the initial guess when I solve the differential equation.

When I used DifferentialEquations.jl and Optimization.jl for this problem, I encountered the following error. (Please see optimization_diff in my code attached below.)

ERROR: Incompatible problem+solver pairing.
For example, this can occur if an ODE solver is passed with an SDEProblem.
Solvers are only capable of handling specific problem types. Please double
check that the chosen pairing is capable for handling the given problems.

Problem type: ODEProblem
Solver type: Shooting
Problem types compatible with the chosen solver: nothing

On the other hand, if I solved the differential equation without the optimization (the function diff_eq), I did not enter any problem. Also, when I printed the properties of my problem, I found it is Problem type of differential equation: BVProblem.... So, I wonder if the error comes from the incompatibility between Optimization.jl and DifferentialEquations.jl.

Below is my environment.

MacOS (M2Max)
Julia: 1.10.0
Zygote v0.6.68
DifferentialEquations v7.12.0
Optimization v3.20.2
OptimizationOptimJL v0.1.14
SciMLSensitivity v7.51.0

using DifferentialEquations, OptimizationOptimJL, Plots

function diff_eq(u0, initial, final, xspan=(0.0, 1.0))
    function f(du, u, p, x)
        du[1] = u[2]
        du[2] = -u[1] * (1 - u[1]^2)
    end

    function bc!(residual, u, p, x)
        residual[1] = u[1][1] - initial
        residual[2] = u[end][1] - final
    end

    bvp = BVProblem(f, bc!, u0, xspan)
    println("Problem type of differential equation: $bvp")
    sol = solve(bvp, Shooting(Vern7()))
    return sol
end


function optimization_diff(initial_guess, initial, final, xspan=(0.0, 1.0))
    x0 = initial_guess
    function objective(v, p)
        u0 = [initial, v[1]]
        sol = diff_eq(u0, initial, final, xspan)
        maximum(sol)
    end
    optprob = OptimizationFunction(objective, Optimization.AutoZygote())
    prob = OptimizationProblem(optprob, x0, [1.0])
    sol = solve(prob, BFGS())
end

diffeq_test = diff_eq([0.0, 0.0], 0.0, 1.0, (0.0, 1.0))
opt_test = optimization_diff([0.0], 0.0, 1.0, (0.0, 1.0))

Move problem types

The problem types should be defined in DiffEqBase so others can build off of it without importing the solvers.

MIRK4() BVP solver gives wrong results

Here is a MWE. Solving an ODE very similar to the Blasius problem (boundary layers in fluid mechanics), as written below:

using DifferentialEquations
using Plots
# ODE
function bleq!(du, u, p, t)
    g, h, k = u
    dg = h
    dh = k
    dk = -g*k  + h^2  - 1
    du[1] = dg
    du[2] = dh
    du[3] = dk
end

# IC: let g = 0, h = 0. 
# BC: at large η, g' = 1
function bc!(resid, u, p, t)
#     k, h, g = 
    resid[1] = u[1][1] # g=0 at t=0
    resid[2] = u[1][2] # g'=0 at t=0
    resid[3] = u[end][2] - 1.0 #g'=1 at t=end
end
# Guessed IC?
icg = [0, 0, 1.2]
ηsp = (0.0, 10.0)
bvp = BVProblem(bleq!, bc!, icg, ηsp)

@time sol = solve(bvp, Shooting(Tsit5()), dt=0.5)

plot(sol, labels=permutedims(["g", "g'", "g''"]), legend=:topleft)

As written, with Shooting(Tsit5()), the problem solves correctly and my plot results look like

When I replace the solver with MIRK4(), as recommended by https://diffeq.sciml.ai/stable/solvers/bvp_solve/ , I get

Singularity handling in MIRK methods

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Residual vector dimension

In BVProblem, can the residual vector have dimension larger that the dimension of u?

Alternatively, if I want to solve a TwoPointBVProblem for a state vector of dimension 10, and I have BCs on 6 variables at t[0] and on 6 variables at t[end], how do I code the residual vector of dimension 12?

For example, if I add an additional condition to the pendulum from the docs

using DifferentialEquations
using Plots

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 bc1!(residual, u, p, t) # u[1] is the beginning of the time span, and u[end] is the ending
    residual[1] = u[1][1] + 1 
    residual[2] = u[end][1] - 1 
    residual[3] = u[end][2] + 2 # from the plot you can see that this is realistic
end

bvp1 = BVProblem(simplependulum!, bc1!, [pi/2,pi/2], tspan)
sol1 = solve(bvp1, Shooting(Tsit5()))
plot(sol1)

the solver casts BoundsError: attempt to access 2-element Core.Array{Float64,1} at index [3]

Add Mass Matrices to the collocation method

and also test DAEs

BV Problem : No possibility to provide non constant initial guess

According to the tutorial on BVP "The third argument of BVProblem is the initial guess of the solution, which is constant in this example.". For the pendulum example (which is a 2nd order dgl) you have to provide two initial guesses, e.g. u0=[0,5]. But what to do if one would like to provide a non constant initial guess such as u0[sin(t),cos(t)]? Then one would provide u0 as u0=[[0,0],[1,1],[n,n]] where n is the desired length of the solution. But this approach throws an error "DimensionMismatch("tried to assign n elements to 2 destinations")". Here is a working example that shows this error

nonConstantIntialGuessError.txt

Throw error on supplying callback for MIRK4

As MIRK4 is fully implicit, callbacks are not supported and therefore supplying one should yield an error.

BigFloat support for TwoPointBVProblem using MIRK4

Using BigFloat with a TwoPointBVProblem and MIRK4 solver currently fails with the following error:

ArgumentError: matrix type SparseArrays.SparseMatrixCSC{BigFloat,Int64}not supported. Try lu(convert(SparseMatrixCSC{Float64/ComplexF64,Int}, A)) for sparse floating point LU using UMFPACK or lu(Array(A)) for generic dense LU.

Multiple Shooting Methods

Multiple shooting methods can be faster and more robust than shooting methods, and just use a simple extension to shooting methods. This requires constrained optimization, but we have all of JuliaOpt to rely on so we can make this pretty good. It would be useful to setup.

Abstract typing

In the DiffEq-ecosystem, the solvers internally match the input type given by the user. In the BVP case, this is a little more complicated but for now let's keep the idea of u0 for this purpose. If that's the case, you don't necessarily have a Vector{Vector}, but a Vector{AbstractArray} for some kind of AbstractArray determined by u0.

It looks like you don't have much required to expand you typing, though these function signatures are overtyped:

https://github.com/JuliaDiffEq/BoundaryValueDiffEq.jl/blob/master/src/collocation.jl#L8

(and the ::Function is unnecessary too)

https://github.com/JuliaDiffEq/BoundaryValueDiffEq.jl/blob/master/src/collocation.jl#L15

I don't see the reason to force type matching here. There can be very good reasons to use a separate type for independent and dependent variables. One case which comes to mind is using rationals for the independent variables. dx=1//10 allows you to make sure you don't have the floating point issues with dx=1/10 having dx = 1/10 + 1/10 + ... + 1/10 = 0.9999994 != 1.0, and using rationals here really doesn't have a computational cost. However, using rationals for the dependent variables can have a massive computational cost since those have to update, while the x are constant. That's a long explanation of: if you use the types the user gives you, it can be very useful to the user.

https://github.com/JuliaDiffEq/BoundaryValueDiffEq.jl/blob/master/src/algorithms.jl#L19

That's another spot which assumes Vector{Vector}.

Improve initial guess interface

BVProblem and TwoPointBVProblem should accept a function initialGuess( t ) and/or a solution structure in order to reuse solutions as initial guess for continuation.
The documentation on this should be improved/ updated accordingly.

Feature request: BVP-DAE support

It would be nice, if solving of boundary value problems for implicit or semi-explicit DAEs was supported. Such solvers exist in FORTRAN (COLDAE).

Dimension mismatch in solve() for MIRK methods providing an initial solution

BVPSolve is allocating more memory than neccessary when provided a initial guess. It is basically allocating a n times n Vector of Vectors, where n is the number of timepoints.

Additionally i get a DimensionMismatch in https://github.com/JuliaDiffEq/BoundaryValueDiffEq.jl/blob/d6b2bcace6c743250366bfea67d2366c4d79d018/src/solve.jl#L55.

It would be good if the dimensions of the initial solution where checked to have the right dimensions and throw a corresponding error if not.

That could be done when the BVProblem gets defined. My testing code can be found here.

Support Initial Guess Function

As pointed out in #117 this is not currently supported.

MIRK
Multiple Shooting
Shooting (should just default to picking the initial condition)

TODO list of GSoC

This is a list of TODOs for my GSoC project.

High priority:

Low priority:

Register

@JuliaRegistrator register()

Interoplation for MIRK methods

We need to implement interpolations for MIRK methods so that we can approximate the solution between mesh points.

Feature request: Derivative boundary condition for BVP ODEs

At the moment BoundaryValueDiffEq doesn't seem to support derivative boundary conditions(neumann, robin) at the boundaries in ODEs of BVP form. Requesting maintainers for same !

Error afer updating to BoundaryValueDiffEq v5: BoundsError in `interpolation`

After updating to BoundaryValueDiffEq v5 CI shows the following error:

Torus in ℝ³: Error During Test at /home/runner/work/Manifolds.jl/Manifolds.jl/test/manifolds/embedded_torus.jl:93
  Test threw exception
  Expression: (sol_log(0.0))[1:2] ≈ p0x
  BoundsError: attempt to access 81-element Vector{Float64} at index [0]
  Stacktrace:
    [1] getindex
      @ ./essentials.jl:13 [inlined]
    [2] interpolation(tval::Float64, id::BoundaryValueDiffEq.MIRKInterpolation{Vector{Float64}, Vector{Vector{Float64}}}, idxs::Nothing, deriv::Type{Val{0}}, p::Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, continuity::Symbol)
      @ BoundaryValueDiffEq ~/.julia/packages/BoundaryValueDiffEq/fSVnX/src/interpolation.jl:81
    [3] (::BoundaryValueDiffEq.MIRKInterpolation{Vector{Float64}, Vector{Vector{Float64}}})(tvals::Float64, idxs::Nothing, deriv::Type, p::Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, continuity::Symbol)
      @ BoundaryValueDiffEq ~/.julia/packages/BoundaryValueDiffEq/fSVnX/src/interpolation.jl:8
    [4] AbstractODESolution
      @ ~/.julia/packages/SciMLBase/McEqc/src/solutions/ode_solutions.jl:75 [inlined]
    [5] #_#433
      @ ~/.julia/packages/SciMLBase/McEqc/src/solutions/ode_solutions.jl:66 [inlined]
    [6] AbstractODESolution
      @ ~/.julia/packages/SciMLBase/McEqc/src/solutions/ode_solutions.jl:64 [inlined]
    [7] (::ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Nothing, BVProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Tuple{Manifolds.EmbeddedTorus{Int64}, Manifolds.DefaultTorusAtlas, Tuple{Int64, Int64}}, BVPFunction{true, SciMLBase.FullSpecialize, false, typeof(ManifoldsBoundaryValueDiffEqExt.chart_log_problem!), ManifoldsBoundaryValueDiffEqExt.var"#bc1!#3"{Vector{Float64}, Vector{Float64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, SciMLBase.StandardBVProblem, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, MIRK4{NewtonRaphson{nothing, AutoForwardDiff{0, Bool}, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), LineSearch{Static, AutoFiniteDiff{Val{:forward}, Val{:forward}, Val{:hcentral}}, Bool}}, MIRKJacobianComputationAlgorithm{AutoForwardDiff{nothing, Nothing}, AutoSparseForwardDiff{nothing, Nothing}, AutoSparseForwardDiff{nothing, Nothing}}}, BoundaryValueDiffEq.MIRKInterpolation{Vector{Float64}, Vector{Vector{Float64}}}, Nothing, Nothing})(t::Float64)
      @ SciMLBase ~/.julia/packages/SciMLBase/McEqc/src/solutions/ode_solutions.jl:64
    [8] macro expansion
      @ /opt/hostedtoolcache/julia/1.9.3/x64/share/julia/stdlib/v1.9/Test/src/Test.jl:478 [inlined]
    [9] macro expansion
      @ ~/work/Manifolds.jl/Manifolds.jl/test/manifolds/embedded_torus.jl:93 [inlined]
   [10] macro expansion
      @ /opt/hostedtoolcache/julia/1.9.3/x64/share/julia/stdlib/v1.9/Test/src/Test.jl:1498 [inlined]
   [11] top-level scope
      @ ~/work/Manifolds.jl/Manifolds.jl/test/manifolds/embedded_torus.jl:11

See here: JuliaManifolds/Manifolds.jl#657 . The v4 version works fine.

Clean Files

https://github.com/JuliaDiffEq/BoundaryValueDiffEq.jl/blob/master/src/algorithms.jl#L17

These functions definitely don't belong in there.

Unexpected kwargs passed to new warning system

This commit to DiffEqBase causes warnings to appear due to unexpected kwargs being passed to the solver: :default_set and :second_time. Either this package or that one ought to be updated so that the default behavior does not produce warnings. (I wasn't sure which repo to create the issue in; forgive me if it ought to be in DiffEqBase.)

DifferentialEquations v7.1.0
BoundaryValueDiffEq v2.7.1

Simple example:

julia> using DifferentialEquations, BoundaryValueDiffEq

julia> dudt(u,p,t) = [u[2] ; -t]  # y'' = -t
dudt (generic function with 1 method)

julia> function boundary!(residual,u,p,t)
           residual[1] = u[1][1]  # u(0) = 0
           residual[2] = u[end][2]  # u'(1) = 0
       end
boundary! (generic function with 1 method)

julia> solve(BVProblem(dudt, boundary!, zeros(2), (0.0, 1.0)))
┌ Warning: Warning: Unrecognized keyword arguments found. Future versions will error.
│ The only allowed keyword arguments to `solve` are: 
│ (:dense, :saveat, :save_idxs, :tstops, :d_discontinuities, :save_everystep, :save_on, :save_start, :save_end, :initialize_save, :adaptive, :abstol, :reltol, :dt, :dtmax, :dtmin, :force_dtmin, :internalnorm, :controller, :gamma, :beta1, :beta2, :qmax, :qmin, :qsteady_min, :qsteady_max, :qoldinit, :failfactor, :calck, :alias_u0, :maxiters, :callback, :isoutofdomain, :unstable_check, :verbose, :merge_callbacks, :progress, :progress_steps, :progress_name, :progress_message, :timeseries_errors, :dense_errors, :calculate_errors, :initializealg, :alg, :save_noise, :delta, :seed, :alg_hints, :kwargshandle, :trajectories, :batch_size, :sensealg, :advance_to_tstop, :stop_at_next_tstop)
│ See https://diffeq.sciml.ai/stable/basics/common_solver_opts/ for more details.
│ 
│ Set kwargshandle=KeywordArgError for an error message and more details.
│ Set kwargshandle=KeywordArgSilent to ignore this message.
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/ziNGu/src/solve.jl:210

The warning is repeated many times, presumably for each Shooting iteration.

dt

dt is defined as a keyword argument for the common interface:

http://docs.juliadiffeq.org/latest/basics/common_solver_opts.html#Basic-Stepsize-Control-1

To match the rest of the DiffEq-universe, you should have it be a keyword argument. That means you just need to change:

https://github.com/JuliaDiffEq/BoundaryValueDiffEq.jl/blob/master/src/solve.jl#L20
https://github.com/JuliaDiffEq/BoundaryValueDiffEq.jl/blob/master/src/algorithms.jl#L13

You should look out for standards like this when developing the API, since even if you're working in a new area, matching the API of similar tools makes it much easier to transfer knowledge.

Note that in that same spot:

http://docs.juliadiffeq.org/latest/basics/common_solver_opts.html#Basic-Stepsize-Control-1

There's an adaptive keyword argument for turning on/off adaptivity, along with tolerances. While those have slightly different meanings in the ODE/SDE/DDE/DAE/PDE/BVP senses, it's clear what it means in each case, so when you go to adaptivity you should follow this same setup.

Solution has no du component

When using TwoPointBVProblem with MIRK4, the output solution doesn't contain a du field.

Doesn't support MVector

I was hoping to use an in-place method on a statically sized array, but BVProblem seemed unhappy about this. MWE:

using BoundaryValueDiffEq, StaticArrays
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 bc1!(residual, u, p, t)
    residual[1] = u[end÷2][1] + pi/2 # the solution at the middle of the time span should be -pi/2
    residual[2] = u[end][1] - pi/2 # the solution at the end of the time span should be pi/2
end

u0 = MVector{2}([pi/2,pi/2]),

bvp1 = BVProblem(simplependulum!, bc1!, u0, tspan)
sol1 = solve(bvp1, GeneralMIRK4(), dt=0.05)

#+RESULTS:
u0 must be a Vector or Vector of Arrays

Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:33
 [2] __solve(prob::BVProblem{MVector{2, Float64}, Tuple{Float64, Float64}, true, DiffEqBase.NullParameters, ODEFunction{true, typeof(simplependulum!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(bc1!), DiffEqBase.StandardBVProblem, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, alg::GeneralMIRK4; dt::Float64, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ BoundaryValueDiffEq ~/.julia/packages/BoundaryValueDiffEq/SKDiI/src/solve.jl:38
 [3] solve_call(_prob::BVProblem{MVector{2, Float64}, Tuple{Float64, Float64}, true, DiffEqBase.NullParameters, ODEFunction{true, typeof(simplependulum!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(bc1!), DiffEqBase.StandardBVProblem, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, args::GeneralMIRK4; merge_callbacks::Bool, kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:dt,), Tuple{Float64}}})
   @ DiffEqBase ~/.julia/packages/DiffEqBase/znIav/src/solve.jl:65
 [4] solve_up(prob::BVProblem{MVector{2, Float64}, Tuple{Float64, Float64}, true, DiffEqBase.NullParameters, ODEFunction{true, typeof(simplependulum!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(bc1!), DiffEqBase.StandardBVProblem, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, sensealg::Nothing, u0::MVector{2, Float64}, p::DiffEqBase.NullParameters, args::GeneralMIRK4; kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:dt,), Tuple{Float64}}})
   @ DiffEqBase ~/.julia/packages/DiffEqBase/znIav/src/solve.jl:92
 [5] solve(prob::BVProblem{MVector{2, Float64}, Tuple{Float64, Float64}, true, DiffEqBase.NullParameters, ODEFunction{true, typeof(simplependulum!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, typeof(bc1!), DiffEqBase.StandardBVProblem, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, args::GeneralMIRK4; sensealg::Nothing, u0::Nothing, p::Nothing, kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:dt,), Tuple{Float64}}})
   @ DiffEqBase ~/.julia/packages/DiffEqBase/znIav/src/solve.jl:74
 [6] top-level scope
   @ In[6]:17
 [7] eval
   @ ./boot.jl:360 [inlined]
 [8] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1090

Switching u0 to a Vector makes the above code work fine.

MIRK methods failed on swirling flow III when viscosity parameter is small

Describe the bug 🐞

$\epsilon$ is the viscosity parameter, when it is set to 0.1, we can successfully get the solution, but when we set $\epsilon$ as 0.01 or smaller, MIRK methods throw error.

Minimal Reproducible Example 👇

using BoundaryValueDiffEq
eps = 0.01
function test!(du, u, p, t)
    eps = p
    du[1] = u[2]
    du[2] = (u[1]*u[4] - u[3]*u[2])/eps
    du[3] = u[4]
    du[4] = u[5]
    du[5] = u[6]
    du[6] = (-u[3]*u[6] - u[1]*u[2])/eps
end

function bc!(res, u, p, t)
    res[1] = u[1][1] + 1.0
    res[2] = u[1][3]
    res[3] = u[1][4]
    res[4] = u[end][1] - 1.0
    res[5] = u[end][3]
    res[6] = u[end][4]
end
tspan = (0.0, 1.0)
u0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
prob = BVProblem(test!, bc!, u0, tspan, eps)
sol=solve(prob, MIRK4(), dt=0.01)

Error & Stacktrace ⚠️

julia> sol = solve(prob, MIRK4(), dt=0.01)
ERROR: MethodError: no method matching (::BoundaryValueDiffEq.var"#119#125"{BoundaryValueDiffEq.MIRKCache{true, Float64, Tuple{Int64}, BVPFunction{true, SciMLBase.FullSpecialize, true, typeof(test!), Tuple{typeof(bca!), typeof(bcb!)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, Tuple{typeof(bca!), typeof(bcb!)}, BVProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, BVPFunction{true, SciMLBase.FullSpecialize, true, typeof(test!), Tuple{typeof(bca!), typeof(bcb!)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, TwoPointBVProblem{true}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, TwoPointBVProblem{true}, SciMLBase.NullParameters, MIRK4{Nothing, BVPJacobianAlgorithm{AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}}}, BoundaryValueDiffEq.MIRKTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Matrix{Float64}}, BoundaryValueDiffEq.MIRKInterpTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{PreallocationTools.DiffCache{Matrix{Float64}, Vector{Float64}}}, Vector{Matrix{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, Vector{Vector{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}, Vector{Float64}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Tuple{Tuple{Int64}, Tuple{Int64}}, NamedTuple{(:defect_threshold, :MxNsub, :abstol, :dt, :adaptive), Tuple{Float64, Int64, Float64, Float64, Bool}}}, TwoPointBVProblem{true}})(::Vector{Float64}, ::Vector{Float64})

Closest candidates are:
  (::BoundaryValueDiffEq.var"#119#125")(::Any, ::Any, ::Any)
   @ BoundaryValueDiffEq ~/.julia/packages/BoundaryValueDiffEq/SrG4p/src/solve/mirk.jl:208

Stacktrace:
  [1] (::NonlinearFunction{true, SciMLBase.FullSpecialize, BoundaryValueDiffEq.var"#119#125"{BoundaryValueDiffEq.MIRKCache{true, Float64, Tuple{Int64}, BVPFunction{true, SciMLBase.FullSpecialize, true, typeof(test!), Tuple{typeof(bca!), typeof(bcb!)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, Tuple{typeof(bca!), typeof(bcb!)}, BVProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, BVPFunction{true, SciMLBase.FullSpecialize, true, typeof(test!), Tuple{typeof(bca!), typeof(bcb!)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, TwoPointBVProblem{true}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, TwoPointBVProblem{true}, SciMLBase.NullParameters, MIRK4{Nothing, BVPJacobianAlgorithm{AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}}}, BoundaryValueDiffEq.MIRKTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Matrix{Float64}}, BoundaryValueDiffEq.MIRKInterpTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{PreallocationTools.DiffCache{Matrix{Float64}, Vector{Float64}}}, Vector{Matrix{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, Vector{Vector{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}, Vector{Float64}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Tuple{Tuple{Int64}, Tuple{Int64}}, NamedTuple{(:defect_threshold, :MxNsub, :abstol, :dt, :adaptive), Tuple{Float64, Int64, Float64, Float64, Bool}}}, TwoPointBVProblem{true}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, BoundaryValueDiffEq.var"#147#151"{SparseDiffTools.ForwardDiffJacobianCache{SparseDiffTools.MatrixColoringResult{Vector{Int64}, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, ArrayInterface.BandedMatrixIndex, ArrayInterface.BandedMatrixIndex}, ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{BoundaryValueDiffEq.BoundaryValueDiffEqTag, Float64}, Float64, 8}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{BoundaryValueDiffEq.BoundaryValueDiffEqTag, Float64}, Float64, 8}}, Vector{Float64}, Vector{Vector{NTuple{8, Float64}}}, Vector{Int64}, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, Vector{Float64}, Vector{Float64}}, Vector{Float64}, BoundaryValueDiffEq.var"#145#149"{BoundaryValueDiffEq.MIRKCache{true, Float64, Tuple{Int64}, BVPFunction{true, SciMLBase.FullSpecialize, true, typeof(test!), Tuple{typeof(bca!), typeof(bcb!)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, Tuple{typeof(bca!), typeof(bcb!)}, BVProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, BVPFunction{true, SciMLBase.FullSpecialize, true, typeof(test!), Tuple{typeof(bca!), typeof(bcb!)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, TwoPointBVProblem{true}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, TwoPointBVProblem{true}, SciMLBase.NullParameters, MIRK4{Nothing, BVPJacobianAlgorithm{AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}}}, BoundaryValueDiffEq.MIRKTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Matrix{Float64}}, BoundaryValueDiffEq.MIRKInterpTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{PreallocationTools.DiffCache{Matrix{Float64}, Vector{Float64}}}, Vector{Matrix{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, Vector{Vector{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}, Vector{Float64}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Tuple{Tuple{Int64}, Tuple{Int64}}, NamedTuple{(:defect_threshold, :MxNsub, :abstol, :dt, :adaptive), Tuple{Float64, Int64, Float64, Float64, Bool}}}, BoundaryValueDiffEq.var"#119#125"{BoundaryValueDiffEq.MIRKCache{true, Float64, Tuple{Int64}, BVPFunction{true, SciMLBase.FullSpecialize, true, typeof(test!), Tuple{typeof(bca!), typeof(bcb!)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, Tuple{typeof(bca!), typeof(bcb!)}, BVProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, BVPFunction{true, SciMLBase.FullSpecialize, true, typeof(test!), Tuple{typeof(bca!), typeof(bcb!)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, RecursiveArrayTools.ArrayPartition{Float64, Tuple{Vector{Float64}, Vector{Float64}}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, TwoPointBVProblem{true}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, TwoPointBVProblem{true}, SciMLBase.NullParameters, MIRK4{Nothing, BVPJacobianAlgorithm{AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}}}, BoundaryValueDiffEq.MIRKTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Matrix{Float64}}, BoundaryValueDiffEq.MIRKInterpTableau{Int64, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{PreallocationTools.DiffCache{Matrix{Float64}, Vector{Float64}}}, Vector{Matrix{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, Vector{Vector{Float64}}, Vector{PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}}, PreallocationTools.DiffCache{Vector{Float64}, Vector{Float64}}, Vector{Float64}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Tuple{Tuple{Int64}, Tuple{Int64}}, NamedTuple{(:defect_threshold, :MxNsub, :abstol, :dt, :adaptive), Tuple{Float64, Int64, Float64, Float64, Bool}}}, TwoPointBVProblem{true}}}, BVPJacobianAlgorithm{AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}, AutoSparseForwardDiff{nothing, BoundaryValueDiffEq.BoundaryValueDiffEqTag}}}, Nothing, Nothing, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, BandedMatrices.BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Vector{Int64}, Nothing, Vector{Float64}})(::Vector{Float64}, ::Vararg{Vector{Float64}})

Environment (please complete the following information):

Output of using Pkg; Pkg.status()

julia> Pkg.status()
Status `~/.julia/environments/v1.9/Project.toml`
  [ded0fc24] BVProblemLibrary v0.1.4
⌃ [6e4b80f9] BenchmarkTools v1.3.2
  [764a87c0] BoundaryValueDiffEq v5.6.0
  [54ca160b] ODEInterface v0.5.0

Output of versioninfo()

julia> versioninfo()
Julia Version 1.9.2
Commit e4ee485e909 (2023-07-05 09:39 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × Intel(R) Core(TM) i5-10210U CPU @ 1.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-14.0.6 (ORCJIT, skylake)
  Threads: 1 on 8 virtual cores

README example should use interpolation instead of indexing

While trying MIRK solvers and Shooting solvers with example problem in README, the two kinds of solvers gave two different solutions:

using BoundaryValueDiffEq, OrdinaryDiffEq
function simplependulum!(du, u, p, t)
    θ = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -9.81 * sin(θ)
end
function bc!(residual, u, p, t)
    residual[1] = u[end ÷ 2][1] + pi / 2
    residual[2] = u[end][1] - pi / 2
end
tspan = (0.0, pi / 2)
prob = BVProblem(simplependulum!, bc!, [pi / 2, pi / 2], tspan)
mirk_sol = solve(prob, MIRK4(), dt = 0.05)
shooting_sol = solve(prob, Shooting(Tsit5()), dt=0.05)

Solution from MIRK(Correct):

Solution from Shooting(Wrong):