# Using vector of vectors for initial guess in TwoPointBVProblem

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):

``````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)

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.

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.)

``````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)
``````

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

``````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 @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?

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

``````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.

1 Like

Thanks for the suggestion.

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

``````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)

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

``````#=
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

``````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.

1 Like

The following works:

``````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)
``````
1 Like

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