# Interpolating solution to constrained BVP

**URL:** <https://discourse.julialang.org/t/interpolating-solution-to-constrained-bvp/137989>\
**Category:** Numerics\
**Tags:** question, differentialequation, bvp\
**Created:** [July 7, 2026, 7:01pm UTC](https://discourse.julialang.org/t/interpolating-solution-to-constrained-bvp/137989 "2026-07-07T19:01:30Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![kar15](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kar15/32/210862_2.png) [@kar15](https://discourse.julialang.org/u/kar15)\
**Post date:** [July 7, 2026, 7:01pm UTC](https://discourse.julialang.org/t/interpolating-solution-to-constrained-bvp/137989/1 "2026-07-07T19:01:30Z")

</div>

Hi all!

When trying to find interpolated values of a solution to a constrained BVP using SciMLs collocation method, I get an error. Here is a MWE:

```julia
using BoundaryValueDiffEq, OptimizationIpopt
function f!(du, u, p, t)
    du[1] = u[2]
    du[2] = u[1]
end
function bc!(res, u, p, t)
    res[1] = u(0.0)[1] - 1
    res[2] = u(1.0)[1]
end
tspan = (0.0, 1.0)
u0 = [0.0, 0.0]
bvp1 = BVProblem(BVPFunction(f!, bc!;f_prototype = zeros(2)),u0,tspan; lb = [-10.0,-10.0],ub = [10.0,10.0])
sol1 = solve(bvp1, MIRK4(; optimize = IpoptOptimizer()); dt = 0.01)
bvp2 = BVProblem(f!, bc!, u0, tspan)
sol2 = solve(bvp2, MIRK4(), dt = 0.01)
println(sol2(0.5))
println(sol1(0.5))

```

This is the output:

```julia-auto

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

This is Ipopt version 3.14.19, running with linear solver MUMPS 5.9.0.

Number of nonzeros in equality constraint Jacobian...: 802
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 0

Total number of variables............................: 202
                     variables with only lower bounds: 0
                variables with lower and upper bounds: 202
                     variables with only upper bounds: 0
Total number of equality constraints.................: 202
Total number of inequality constraints...............: 0
        inequality constraints with only lower bounds: 0
   inequality constraints with lower and upper bounds: 0
        inequality constraints with only upper bounds: 0

iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
   0 0.0000000e+00 1.00e+00 0.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
   1 0.0000000e+00 3.24e-16 0.00e+00 -1.0 1.31e+00 - 8.83e-01 1.00e+00h 1

Number of Iterations....: 1

                                   (scaled) (unscaled)
Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
Constraint violation....: 3.2439329000766293e-16 3.2439329000766293e-16
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 3.2439329000766293e-16 3.2439329000766293e-16

Number of objective function evaluations = 2
Number of objective gradient evaluations = 2
Number of equality constraint evaluations = 2
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 2
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 1
Total seconds in IPOPT = 9.040

EXIT: Optimal Solution Found.
[0.4434094419864666, -0.9595173756819017]

ArgumentError: range must be non-empty

Stacktrace:
  [1] minimum
    @ ./range.jl:863 [inlined]
  [2] copyto!(dest::Vector{Float64}, rdest::UnitRange{Int64}, src::Vector{Float64}, rsrc::UnitRange{Int64})
    @ LinearAlgebra ~/.julia/juliaup/julia-1.12.6+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/blas.jl:2241
  [3] sum_stages!(z::Vector{…}, id::BoundaryValueDiffEqMIRK.MIRKInterpolation{…}, cache::BoundaryValueDiffEqMIRK.MIRKCache{…}, w::Vector{…}, i::Int64, τ::Float64, ::Type{…})
    @ BoundaryValueDiffEqMIRK ~/.julia/packages/BoundaryValueDiffEqMIRK/z0Efo/src/interpolation.jl:116
  [4] interpolant!
    @ ~/.julia/packages/BoundaryValueDiffEqMIRK/z0Efo/src/interpolation.jl:79 [inlined]
  [5] interpolation
    @ ~/.julia/packages/BoundaryValueDiffEqMIRK/z0Efo/src/interpolation.jl:68 [inlined]
  [6] (::BoundaryValueDiffEqMIRK.MIRKInterpolation{…})(tvals::Float64, idxs::Nothing, deriv::Type, p::SciMLBase.NullParameters, continuity::Symbol)
    @ BoundaryValueDiffEqMIRK ~/.julia/packages/BoundaryValueDiffEqMIRK/z0Efo/src/interpolation.jl:13
  [7] AbstractODESolution
    @ ~/.julia/packages/SciMLBase/dpRhw/src/solutions/ode_solutions.jl:241 [inlined]
  [8] #_#625
    @ ~/.julia/packages/SciMLBase/dpRhw/src/solutions/ode_solutions.jl:225 [inlined]
  [9] AbstractODESolution
    @ ~/.julia/packages/SciMLBase/dpRhw/src/solutions/ode_solutions.jl:218 [inlined]
 [10] (::ODESolution{…})(t::Float64)
    @ SciMLBase ~/.julia/packages/SciMLBase/dpRhw/src/solutions/ode_solutions.jl:218
 [11] top-level scope
    @ In[2]:17
 [12] eval(m::Module, e::Any)
    @ Core ./boot.jl:489
Some type information was truncated. Use `show(err)` to see complete types.

```

The constrained solution has a success retcode and I can access the saved values which match those of the unconstrained version but I can’t interpolate the solution to other points. Is it a bug on my end? If not, is it possible to construct an interpolation from the saved points?

Thanks!

---

<div class="post-metadata">

**Author:** ![technocrat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/technocrat/32/220947_2.png) [@technocrat](https://discourse.julialang.org/u/technocrat)\
**Post date:** [July 7, 2026, 9:09pm UTC](https://discourse.julialang.org/t/interpolating-solution-to-constrained-bvp/137989/2 "2026-07-07T21:09:44Z")

</div>

What’s different between `bvp1` and `bvp2` is that `bvp1` is built with `BVPFunction(f!, bc!; f_prototype = zeros(2))`, while `bvp2` uses `f!`/`bc!` directly with no `f_prototype`. That’s the only structural difference that matters here; the `optimize = IpoptOptimizer()` solve algorithm itself isn’t the trigger.

Per SciMLBase’s docstring, `f_prototype` exists for optimal‑control BVPs problems where the unknown vector `u` is split into (a) true ODE/DAE state variables governed by `f`, and (b) extra “control” variables that Ipopt is free to choose, not governed by any dynamics. `f_prototype`’s length tells the solver how many of the `M` total unknowns fall into bucket (a).

In `BoundaryValueDiffEqMIRK/src/interpolation.jl`, `sum_stages!` checks:

```julia
has_control = !isnothing(cache.prob.f.f_prototype)
length_z = has_control ? length(cache.prob.f.f_prototype) : length(z)
...
if has_control
    inc = τ / dt .* (id.u[i + 1] .- id.u[i])
    copyto!(z, (length_z + 1):M, inc, (length_z + 1):M) # interpolation.jl:116
end

```

You passed `f_prototype = zeros(2)`, whose length (2) equals the full state length `M` (2) — i.e., zero actual control variables. So `(length_z+1):M` is `3:2`, an empty range. `Base.copyto!` tolerates that fine, but the specialized `copyto!(dest, rdest, src, rsrc)` method in `LinearAlgebra` (the one actually dispatched to here) calls `minimum(rdest)` unconditionally — and `minimum` of an empty `UnitRange` throws `ArgumentError: range must be non-empty`.

Note this only fires on interpolation _after_ the solve. The solve itself, the retcode, and the saved mesh values are all computed on a different code path (`__mirk_loss!`/collocation) that never touches this branch, which is why `sol1` “looks fine” until you call `sol1(0.5)`.

You don’t really need `f_prototype` because your problem has no separate control variables, only two coupled ODE states with box constraints. Just drop it:

```julia
bvp1 = BVProblem(BVPFunction(f!, bc!), u0, tspan; lb = [-10.0, -10.0], ub = [10.0, 10.0])

```

`lb`/`ub` are independent `BVProblem` fields (bounds on solution variables) and don’t require `f_prototype` — that should let interpolation fall through the normal (non-`has_control`) branch, same as `bvp2`.

---

<div class="post-metadata">

**Author:** ![kar15](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kar15/32/210862_2.png) [@kar15](https://discourse.julialang.org/u/kar15)\
**Post date:** [July 8, 2026, 5:33am UTC](https://discourse.julialang.org/t/interpolating-solution-to-constrained-bvp/137989/3 "2026-07-08T05:33:43Z")

</div>

Thanks for explaining the issue!

I had tried removing `f_prototype = zeros(2)` but then the solver does not work. I get this error:

```julia-auto
The function `copy` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  copy(::Core.PhiNode)
   @ Base expr.jl:42
  copy(::LinearAlgebra.Transpose{Bool, BitMatrix})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.12.6+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/bitarray.jl:240
  copy(::Base.JuliaSyntax.SyntaxData)
   @ Base /cache/build/builder-amdci5-4/julialang/julia-release-1-dot-12/base/JuliaSyntax/src/syntax_tree.jl:257
  ...

Stacktrace:
  [1] _broadcast_getindex_evalf
    @ ./broadcast.jl:699 [inlined]
  [2] _broadcast_getindex
    @ ./broadcast.jl:672 [inlined]
  [3] _getindex
    @ ./broadcast.jl:696 [inlined]
  [4] _broadcast_getindex
    @ ./broadcast.jl:671 [inlined]
  [5] _getindex
    @ ./broadcast.jl:620 [inlined]
  [6] getindex
    @ ./broadcast.jl:616 [inlined]
  [7] copy
    @ ./broadcast.jl:933 [inlined]
  [8] materialize
    @ ./broadcast.jl:894 [inlined]
  [9] __init(prob::BVProblem{…}, alg::MIRK4{…}; dt::Float64, abstol::Float64, adaptive::Bool, controller::BoundaryValueDiffEqCore.DefectControl{…}, nlsolve_kwargs::@NamedTuple{…}, optimize_kwargs::@NamedTuple{…}, verbose::DiffEqBase.DEVerbosity{…}, kwargs::@Kwargs{})
    @ BoundaryValueDiffEqMIRK ~/.julia/packages/BoundaryValueDiffEqMIRK/z0Efo/src/mirk.jl:129
 [10] __init
    @ ~/.julia/packages/BoundaryValueDiffEqMIRK/z0Efo/src/mirk.jl:44 [inlined]
 [11] __solve(::BVProblem{…}, ::MIRK4{…}; kwargs::@Kwargs{…})
    @ BoundaryValueDiffEqCore ~/.julia/packages/BoundaryValueDiffEqCore/CwTMw/src/BoundaryValueDiffEqCore.jl:47
 [12] solve_call(_prob::BVProblem{…}, args::MIRK4{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/F2QKS/src/solve.jl:172
 [13] solve_up(prob::BVProblem{…}, sensealg::Nothing, u0::Vector{…}, p::SciMLBase.NullParameters, args::MIRK4{…}; originator::SciMLBase.ChainRulesOriginator, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/F2QKS/src/solve.jl:646
 [14] solve_up
    @ ~/.julia/packages/DiffEqBase/F2QKS/src/solve.jl:619 [inlined]
 [15] #solve#31
    @ ~/.julia/packages/DiffEqBase/F2QKS/src/solve.jl:603 [inlined]
 [16] top-level scope
    @ In[20]:13
 [17] eval(m::Module, e::Any)
    @ Core ./boot.jl:489
Some type information was truncated. Use `show(err)` to see complete types.

```

I think when the residual is being initialized, `f_prototype` is expected when there is any constraints including box bounds : `BoundaryValueDiffEqMIRK/z0Efo/src/mirk.jl:129` has

```julia
vcat(
    [__alloc(bcresid_prototype)],__alloc.(
        copy.(
            [
                f_prototype
                    for _ in 1:Nig
            ]
        )
    )
)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 9, 2026, 7:41am UTC](https://discourse.julialang.org/t/interpolating-solution-to-constrained-bvp/137989/4 "2026-07-09T07:41:02Z")

</div>

Yeah it’s a bug.