Solving via Continuation: Interpolation error when using previous solution as initial guess

[sebastiantk]

Hi,

I am trying to solve a two point boundary value problem via continuation but I am encountering an error. I replicate the problem using the example problem.

I first solve the problem from the end point to the middle of the time span:

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


function bc2a!(resid_a, u_a, p) # u_a is at the beginning of the time span
    x0 = p
    resid_a[1] = u_a[1] - x0 # 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
    x0 = p
    resid_b[1] = u_b[1] - pi / 2 # the solution at the end of the time span should be pi/2
end


bvp3 = TwoPointBVProblem(simplependulum!, (bc2a!, bc2b!), [pi / 2, pi / 2], (pi / 4, pi / 2), -pi / 2;
    bcresid_prototype=(zeros(1), zeros(1)))
sol3 = solve(bvp3, MIRK4(), dt=0.05)

which yields the correct solution

I then try to solve the problem over the original time span using the previous solution as the initial guess

bvp4 = TwoPointBVProblem(simplependulum!, (bc2a!, bc2b!), sol3, (0, pi / 2), pi / 2;
    bcresid_prototype=(zeros(1), zeros(1)))
sol4 = solve(bvp4, MIRK4(), dt=0.05)

This however yields an error

ERROR: MethodError: no method matching sortperm(::Float64; rev::Bool)

Closest candidates are:
  sortperm(::StaticArraysCore.StaticArray{Tuple{N}, T, 1} where {N, T}; alg, lt, by, rev, order)
   @ StaticArrays ~/.julia/packages/StaticArrays/oOCPP/src/sort.jl:26
  sortperm(::AbstractUnitRange) got unsupported keyword argument "rev"
   @ Base range.jl:1415
  sortperm(::AbstractRange) got unsupported keyword argument "rev"
   @ Base range.jl:1416
  ...

Stacktrace:
  [1] interpolation!
    @ ~/.julia/packages/BoundaryValueDiffEq/spgaR/src/interpolation.jl:47 [inlined]
  [2] (::BoundaryValueDiffEq.MIRKInterpolation{…})(val::Float64, tvals::Float64, idxs::Nothing, deriv::Type, p::Float64, continuity::Symbol)
    @ BoundaryValueDiffEq ~/.julia/packages/BoundaryValueDiffEq/spgaR/src/interpolation.jl:16
  [3] (::ODESolution{…})(v::Float64, t::Float64, ::Type{…}; idxs::Nothing, continuity::Symbol)
    @ SciMLBase ~/.julia/packages/SciMLBase/aft1j/src/solutions/ode_solutions.jl:145
  [4] (::ODESolution{…})(v::Float64, t::Float64, ::Type{…})
    @ SciMLBase ~/.julia/packages/SciMLBase/aft1j/src/solutions/ode_solutions.jl:143
  [5] __initial_guess(f::ODESolution{…}, p::Float64, t::Float64)
    @ BoundaryValueDiffEq ~/.julia/packages/BoundaryValueDiffEq/spgaR/src/utils.jl:153
  [6] concrete_jacobian_algorithm(jac_alg::BVPJacobianAlgorithm{…}, prob_type::TwoPointBVProblem{…}, prob::BVProblem{…}, alg::MIRK4{…})
    @ BoundaryValueDiffEq ~/.julia/packages/BoundaryValueDiffEq/spgaR/src/types.jl:79
  [7] concrete_jacobian_algorithm(jac_alg::BVPJacobianAlgorithm{…}, prob::BVProblem{…}, alg::MIRK4{…})
    @ BoundaryValueDiffEq ~/.julia/packages/BoundaryValueDiffEq/spgaR/src/types.jl:74
  [8] macro expansion
    @ ~/.julia/packages/Setfield/PdKfV/src/sugar.jl:197 [inlined]
  [9] __init(prob::BVProblem{…}, alg::MIRK4{…}; dt::Float64, abstol::Float64, adaptive::Bool, kwargs::@Kwargs{})
    @ BoundaryValueDiffEq ~/.julia/packages/BoundaryValueDiffEq/spgaR/src/solve/mirk.jl:36
 [10] init_call(_prob::BVProblem{…}, args::MIRK4{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/a6p43/src/solve.jl:489
 [11] init_up(prob::BVProblem{…}, sensealg::Nothing, u0::ODESolution{…}, p::Float64, args::MIRK4{…}; kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/a6p43/src/solve.jl:522
 [12] init(prob::BVProblem{…}, args::MIRK4{…}; sensealg::Nothing, u0::Nothing, p::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/a6p43/src/solve.jl:502
 [13] __solve(::BVProblem{…}, ::MIRK4{…}; kwargs::@Kwargs{…})
    @ BoundaryValueDiffEq ~/.julia/packages/BoundaryValueDiffEq/spgaR/src/BoundaryValueDiffEq.jl:46
 [14] solve_call(_prob::BVProblem{…}, args::MIRK4{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/a6p43/src/solve.jl:561
 [15] solve_up(prob::BVProblem{…}, sensealg::Nothing, u0::ODESolution{…}, p::Float64, args::MIRK4{…}; kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/a6p43/src/solve.jl:1010
 [16] solve(prob::BVProblem{…}, args::MIRK4{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/a6p43/src/solve.jl:933

I am using Julia v1.10.0 and BoundaryValueDiffEq v5.6.1

Should the initial guess be provided differently?

Thank you!