# Using vector of vectors for initial guess in TwoPointBVProblem

**URL:** <https://discourse.julialang.org/t/using-vector-of-vectors-for-initial-guess-in-twopointbvproblem/88362>\
**Category:** New to Julia\
**Tags:** diffeq, bvp\
**Created:** [October 6, 2022, 4:49pm UTC](https://discourse.julialang.org/t/using-vector-of-vectors-for-initial-guess-in-twopointbvproblem/88362 "2022-10-06T16:49:46Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![morningglory05](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/morningglory05/32/43282_2.png) [@morningglory05](https://discourse.julialang.org/u/morningglory05)\
**Post date:** [October 6, 2022, 4:49pm UTC](https://discourse.julialang.org/t/using-vector-of-vectors-for-initial-guess-in-twopointbvproblem/88362/1 "2022-10-06T16:49:46Z")

</div>

Hello! I am trying to run a bvp to obtain a periodic orbit. A forward integration IVP works with the current parameters, but there is an instability that causes problems with different parameters. This can be solved by using a bvp where the initial guesses at the collocation nodes are the points found using forward integration. Basically, I need to use a collocation bvp method. The output is incorrect when the initial guess is a constant. My understanding of `TwoPointBVProblem` is that I should be able to use a vector of arrays to give different inital guesses at each collocation node.

When I run the following cell block (the code leading up to this cell block is included further down):

```julia
bvp1 = TwoPointBVProblem(TC!, bc_po!, sol.u, tspan)
sol6 = solve(bvp1, GeneralMIRK4(), dt = .1) 

```

I get the following error:

> MethodError: no method matching zero(::Type{Vector{Float64}})  
> Closest candidates are:  
> zero(::Union{Type{P}, P}) where P\<:Dates.Period at C:\Users\johns.julia\juliaup\julia-1.8.2+0.x64\share\julia\stdlib\v1.8\Dates\src\periods.jl:53  
> zero(::SparseArrays.AbstractSparseArray) at C:\Users\johns.julia\juliaup\julia-1.8.2+0.x64\share\julia\stdlib\v1.8\SparseArrays\src\SparseArrays.jl:42  
> zero(::VectorizationBase.Vec{W, T}) where {W, T} at C:\Users\johns.julia\packages\VectorizationBase\cTMlH\src\llvm\_intrin\vbroadcast.jl:181  
> …
> 
> Stacktrace:  
> [1] zeros(#unused#::Type{Vector{Float64}}, dims::Tuple{Int64})  
> @ Base .\array.jl:589  
> [2] zeros(#unused#::Type{Vector{Float64}}, dims::Int64)  
> @ Base .\array.jl:584  
> [3] vector\_alloc  
> @ C:\Users\johns.julia\packages\BoundaryValueDiffEq\g0ius\src\vector\_auxiliary.jl:5 [inlined]  
> [4] BoundaryValueDiffEq.BVPSystem(prob::BVProblem{Vector{Vector{Float64}}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, typeof(TC!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT\_OBSERVED), Nothing, Nothing}, SciMLBase.TwoPointBVPFunction{typeof(bc\_po!)}, SciMLBase.StandardBVProblem, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, x::Vector{Float64}, order::Int64)  
> @ BoundaryValueDiffEq C:\Users\johns.julia\packages\BoundaryValueDiffEq\g0ius\src\collocation.jl:23  
> [5] \_\_solve(prob::BVProblem{Vector{Vector{Float64}}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, typeof(TC!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT\_OBSERVED), Nothing, Nothing}, SciMLBase.TwoPointBVPFunction{typeof(bc\_po!)}, SciMLBase.StandardBVProblem, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, alg::GeneralMIRK4; dt::Float64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})  
> @ BoundaryValueDiffEq C:\Users\johns.julia\packages\BoundaryValueDiffEq\g0ius\src\solve.jl:33  
> [6] solve\_call(\_prob::BVProblem{Vector{Vector{Float64}}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, typeof(TC!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT\_OBSERVED), Nothing, Nothing}, SciMLBase.TwoPointBVPFunction{typeof(bc\_po!)}, SciMLBase.StandardBVProblem, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, args::GeneralMIRK4; merge\_callbacks::Bool, kwargshandle::KeywordArgError, kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:dt,), Tuple{Float64}}})  
> @ DiffEqBase C:\Users\johns.julia\packages\DiffEqBase\5rKYk\src\solve.jl:472  
> [7] solve\_up(prob::BVProblem{Vector{Vector{Float64}}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, typeof(TC!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT\_OBSERVED), Nothing, Nothing}, SciMLBase.TwoPointBVPFunction{typeof(bc\_po!)}, SciMLBase.StandardBVProblem, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, sensealg::Nothing, u0::Vector{Vector{Float64}}, p::SciMLBase.NullParameters, args::GeneralMIRK4; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:dt,), Tuple{Float64}}})  
> @ DiffEqBase C:\Users\johns.julia\packages\DiffEqBase\5rKYk\src\solve.jl:834  
> [8] solve(prob::BVProblem{Vector{Vector{Float64}}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, typeof(TC!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT\_OBSERVED), Nothing, Nothing}, SciMLBase.TwoPointBVPFunction{typeof(bc\_po!)}, SciMLBase.StandardBVProblem, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, args::GeneralMIRK4; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{true}, kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:dt,), Tuple{Float64}}})  
> @ DiffEqBase C:\Users\johns.julia\packages\DiffEqBase\5rKYk\src\solve.jl:801  
> [9] top-level scope  
> @ In[10]:2  
> [10] eval  
> @ .\boot.jl:368 [inlined]  
> [11] include\_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)  
> @ Base .\loading.jl:1428

It is of note that another user tried something similar and got the following error:

> “DimensionMismatch(“tried to assign n elements to 2 destinations”)”

When I run their code I get my error, not theirs. Their question is posted on the [github issues](https://github.com/SciML/BoundaryValueDiffEq.jl/issues/73).

Here is the complete code necessary to see my error: (Note, I’m new to coding, so I’m still learning how to make things easier to read.)

```julia
using ForwardDiff
using BoundaryValueDiffEq
using OrdinaryDiffEq

#System Constants
ss = 1 #excitatory parameter
sj = 0 #inhibitory parameter
glb = 0.05
el = -70
gnab = 3
ena = 50
gkb = 5
ek = -90
gtb = 2 
et = 90 
gex = 0.1 
vex = 0
gsyn = 0.13
vsyn = -85
iext = 0.41
eps = 1
qht = 2.5

#System functions
function f(v,h,r)
    -(glb*(v-el)+gnab*(1/(1+exp(-(v+37)/7)))^3*h*(v-ena)+gkb*(0.75*(1-h))^4*(v-ek)+gtb*(1/(1+exp(-(v+60)/6.2)))^2*r*(v-et))-gex*ss*(v-vex)-gsyn*sj*(v-vsyn)+iext
end

function g(v,h)
    eps*((1/(1+exp((v+41)/4)))-h)/(1/((0.128*exp(-(v+46)/18))+(4/(1+exp(-(v+23)/5)))))
end

function k(v,r)
    qht*((1/(1+exp((v+84)/4)))-r)/((28+exp(-(v+25)/10.5)))
end

#Dynamical System
function TC!(du,u,p,t)
    v,h,r = u
    
    du[1] = dv = f(v,h,r)
    du[2] = dh = g(v,h)
    du[3] = dr = k(v,r)
end

#Finding initial guesses by forward integration
T = 7.588145762136627 #orbit length
u0 = [-40.296570996984855, 0.7298857398191566, 0.0011680534089275774]
tspan = (0.0,T)
prob = ODEProblem(TC!,u0,tspan,dt=0.01)
sol = solve(prob, Rodas4P(), reltol = 1e-30, abstol = 1e-15)

#The BVP set up
function bc_po!(residual, u, p, t)
    residual[1] = u[1][1] - u[end][1]
    residual[2] = u[1][2] - u[end][2]
    residual[3] = u[1][3] - u[end][3]
end

#This is the part of the code that has problems
bvp1 = TwoPointBVProblem(TC!, bc_po!, sol.u, tspan)
sol6 = solve(bvp1, GeneralMIRK4(), dt = .1) 

```

---

<div class="post-metadata">

**Author:** ![liamfdoherty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/liamfdoherty/32/24516_2.png) [@liamfdoherty](https://discourse.julialang.org/u/liamfdoherty)\
**Post date:** [October 11, 2022, 8:02pm UTC](https://discourse.julialang.org/t/using-vector-of-vectors-for-initial-guess-in-twopointbvproblem/88362/2 "2022-10-11T20:02:30Z")

</div>

Following the stacktrace, the error can be reproduced in the line

```julia
zeros(eltype(sol.u), length(sol.u))

```

and the error appears to be coming from the fact that Julia doesn’t know what to do about making a vector of zeros with a `Vector` element type (maybe because of issues of not having information of size to allocate?). My initial thought would be to use StaticArrays, but in [this post](https://discourse.julialang.org/t/problem-with-differentialequations-staticarrays/39129) @ChrisRackauckas indicates that they aren’t appropriate for ODE solving. Maybe somebody has ideas if there is a way to package the initial guess `sol.u`, maybe as a Matrix instead of as `Vector{Vector{Float64}}` or something similar?

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petersimon/32/25193_2.png) [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)\
**Post date:** [October 11, 2022, 8:50pm UTC](https://discourse.julialang.org/t/using-vector-of-vectors-for-initial-guess-in-twopointbvproblem/88362/3 "2022-10-11T20:50:53Z")

</div>

You could try using an `Marray` or `MVector` from `StaticArrays`. They are mutable and `zeros` has a method defined to work with them:

```julia
julia> using StaticArrays

julia> t = [@MVector [1.0,2.0] for _ in 1:4]
4-element Vector{MVector{2, Float64}}:
 [1.0, 2.0]
 [1.0, 2.0]
 [1.0, 2.0]
 [1.0, 2.0]

julia> zeros(eltype(t), 3)
3-element Vector{MVector{2, Float64}}:
 [0.0, 0.0]
 [0.0, 0.0]
 [0.0, 0.0]

```

Edit: Disclaimer: I know nothing about the differential equations packages. This suggestion only addresses the surface problem you reported.

---

<div class="post-metadata">

**Author:** ![morningglory05](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/morningglory05/32/43282_2.png) [@morningglory05](https://discourse.julialang.org/u/morningglory05)\
**Post date:** [October 12, 2022, 5:55pm UTC](https://discourse.julialang.org/t/using-vector-of-vectors-for-initial-guess-in-twopointbvproblem/88362/4 "2022-10-12T17:55:28Z")

</div>

Thanks for the suggestion.

When I use an MVector instead (as in the following changes to my original code):

```julia
using StaticArrays

init_1 = [@MVector [sol.u[i]] for i in 1:length(sol.u)]
bvp1 = TwoPointBVProblem(TC!, bc_po!, init_1, tspan)
sol6 = solve(bvp1, GeneralMIRK4(), dt = .1) 

```

I get the following error:

> MethodError: no method matching one(::Type{Vector{Float64}})  
> Closest candidates are:  
> one(::Union{Type{T}, T}) where T\<:AbstractString at strings/basic.jl:262  
> one(::Union{Type{P}, P}) where P\<:Dates.Period at C:\Users\johns.julia\juliaup\julia-1.8.2+0.x64\share\julia\stdlib\v1.8\Dates\src\periods.jl:54  
> one(::VectorizationBase.Vec{W, T}) where {W, T} at C:\Users\johns.julia\packages\VectorizationBase\cTMlH\src\llvm\_intrin\vbroadcast.jl:182  
> …
> 
> Stacktrace:  
> [1] recursive\_unitless\_eltype(a::Type{MVector{1, Vector{Float64}}})  
> @ RecursiveArrayTools C:\Users\johns.julia\packages\RecursiveArrayTools\IgoeN\src\utils.jl:217  
> [2] recursive\_unitless\_eltype(a::Vector{MVector{1, Vector{Float64}}})  
> @ RecursiveArrayTools C:\Users\johns.julia\packages\RecursiveArrayTools\IgoeN\src\utils.jl:210  
> [3] solve\_call(\_prob::BVProblem{Vector{MVector{1, Vector{Float64}}}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, typeof(TC!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT\_OBSERVED), Nothing, Nothing}, SciMLBase.TwoPointBVPFunction{typeof(bc\_po!)}, SciMLBase.StandardBVProblem, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, args::GeneralMIRK4; merge\_callbacks::Bool, kwargshandle::KeywordArgError, kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:dt,), Tuple{Float64}}})  
> @ DiffEqBase C:\Users\johns.julia\packages\DiffEqBase\5rKYk\src\solve.jl:458  
> [4] solve\_up(prob::BVProblem{Vector{MVector{1, Vector{Float64}}}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, typeof(TC!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT\_OBSERVED), Nothing, Nothing}, SciMLBase.TwoPointBVPFunction{typeof(bc\_po!)}, SciMLBase.StandardBVProblem, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, sensealg::Nothing, u0::Vector{MVector{1, Vector{Float64}}}, p::SciMLBase.NullParameters, args::GeneralMIRK4; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:dt,), Tuple{Float64}}})  
> @ DiffEqBase C:\Users\johns.julia\packages\DiffEqBase\5rKYk\src\solve.jl:834  
> [5] solve(prob::BVProblem{Vector{MVector{1, Vector{Float64}}}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.FullSpecialize, typeof(TC!), UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT\_OBSERVED), Nothing, Nothing}, SciMLBase.TwoPointBVPFunction{typeof(bc\_po!)}, SciMLBase.StandardBVProblem, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, args::GeneralMIRK4; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{true}, kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol}, NamedTuple{(:dt,), Tuple{Float64}}})  
> @ DiffEqBase C:\Users\johns.julia\packages\DiffEqBase\5rKYk\src\solve.jl:801  
> [6] top-level scope  
> @ In[21]:2  
> [7] eval  
> @ .\boot.jl:368 [inlined]  
> [8] include\_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)  
> @ Base .\loading.jl:1428

Based on the first line of the stacktrace, the github documentation tells me:

```julia
#=
recursive_unitless_eltype(a)

Grabs the unitless element type. For example, if
ones has a `Array{Array{Float64,N},N}`, this will return `Array{Float64,N}`.
=#

recursive_unitless_eltype(a) = recursive_unitless_eltype(eltype(a))
recursive_unitless_eltype(a::Type{Any}) = Any

```

Which is part of `RecursiveArrayTools`. However, just running the line

```julia
recursive_unitless_eltype([@MVector [sol.u[i]] for i in 1:length(sol.u)])

```

Gives me the same error as before. This inability to mutate using the `RecursiveArrayTools` seems to be what was referenced in the post mentioned by @liamfdoherty.

I don’t understand static arrays enough to know if there are options that might work.

---

<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:** [December 25, 2022, 6:12am UTC](https://discourse.julialang.org/t/using-vector-of-vectors-for-initial-guess-in-twopointbvproblem/88362/5 "2022-12-25T06:12:08Z")

</div>

The following works:

```julia
using BoundaryValueDiffEq
using OrdinaryDiffEq

#System Constants
ss = 1 #excitatory parameter
sj = 0 #inhibitory parameter
glb = 0.05
el = -70
gnab = 3
ena = 50
gkb = 5
ek = -90
gtb = 2 
et = 90 
gex = 0.1 
vex = 0
gsyn = 0.13
vsyn = -85
iext = 0.41
eps = 1
qht = 2.5

#System functions
function f(v,h,r)
    -(glb*(v-el)+gnab*(1/(1+exp(-(v+37)/7)))^3*h*(v-ena)+gkb*(0.75*(1-h))^4*(v-ek)+gtb*(1/(1+exp(-(v+60)/6.2)))^2*r*(v-et))-gex*ss*(v-vex)-gsyn*sj*(v-vsyn)+iext
end

function g(v,h)
    eps*((1/(1+exp((v+41)/4)))-h)/(1/((0.128*exp(-(v+46)/18))+(4/(1+exp(-(v+23)/5)))))
end

function k(v,r)
    qht*((1/(1+exp((v+84)/4)))-r)/((28+exp(-(v+25)/10.5)))
end

#Dynamical System
function TC!(du,u,p,t)
    v,h,r = u
    
    du[1] = dv = f(v,h,r)
    du[2] = dh = g(v,h)
    du[3] = dr = k(v,r)
end

#Finding initial guesses by forward integration
T = 7.588145762136627 #orbit length
u0 = [-40.296570996984855, 0.7298857398191566, 0.0011680534089275774]
tspan = (0.0,T)
prob = ODEProblem(TC!,u0,tspan,dt=0.01)
sol = solve(prob, Rodas4P(), reltol = 1e-12, abstol = 1e-12, saveat = 0.5)

#The BVP set up
function bc_po!(residual, u, p, t)
    residual[1] = u[1][1] - u[end][1]
    residual[2] = u[1][2] - u[end][2]
    residual[3] = u[1][3] - u[end][3]
end

#This is the part of the code that has problems
bvp1 = TwoPointBVProblem(TC!, bc_po!, sol.u, tspan)
sol6 = solve(bvp1, GeneralMIRK4(), dt = 0.5) 

```

with the patch [Fix vector of vector initial conditions with GeneralMIRK4 by ChrisRackauckas · Pull Request #82 · SciML/BoundaryValueDiffEq.jl · GitHub](https://github.com/SciML/BoundaryValueDiffEq.jl/pull/82)

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [December 25, 2022, 7:39am UTC](https://discourse.julialang.org/t/using-vector-of-vectors-for-initial-guess-in-twopointbvproblem/88362/6 "2022-12-25T07:39:41Z")

</div>

Hi

This is a common problem in bifurcation theory. You need to add a phase condition to remove the singularity. This is treated at length in BifurcationKit.jl where 4 different methods are provided