Using vector of vectors for initial guess in TwoPointBVProblem

New to Julia
,
1

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

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

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?

3

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
4

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)
@ Base .\loading.jl:1428

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
5

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)

with the patch Fix vector of vector initial conditions with GeneralMIRK4 by ChrisRackauckas · Pull Request #82 · SciML/BoundaryValueDiffEq.jl · GitHub

2 Likes
6

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