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!

[ChrisRackauckas]

I think this is just a missing dispatch, we should support AbstractVectorOfArray and AbstractDiffEqArray u0s. For now just do the sol3.u

[avik-pal]

This is a SciMLBase fix that is needed here. AbstractDiffEq array is already supported inside MIRK (at least on #155), however, here a solution object is being passed to u0 which should be converted to a DiffEqArray during problem construction (I guess there could be some value in letting the solution type propagate so that we could use interpolation to initialize on an arbitrary mesh)

[avik-pal]

To add to Chris' solution, note that just passing sol3.u would make an implicit assumption that the initial guess mesh is uniform (which happens to be true in this case). But in general you would want to pass a DiffEqArray in most cases.