Bifurcation analysis with history in ODEs

nkoch · July 11, 2022, 7:05pm

Hello, I’m relatively new to using Julia and I am trying to do bifurcation analysis.

I have a system of ODEs in which one DE depends on the history of another variable (i.e. du[2] depends on history of u[1]). From what I understand there are a few options how to implement this:

Have the history as a parameter and update it using callbacks
Have the history as a parameter, use the integrator interface and update it explicitly with each step

I haven’t been able to understand if I can (or how to) use either of these methods with BifurcationKit.jl. Are either of these methods compatible with bifurcation analysis using BifurcationKit.jl? Or is there a better way to implement this?

Thanks!
Nils

josuagrw · July 11, 2022, 10:10pm

This sounds like a typical DDE problem

https://diffeq.sciml.ai/stable/types/dde_types/

Can you say more about the way u_2 depends on the history of u_1?

rveltz · July 11, 2022, 10:57pm

Hi,

It is a DDE problem like said above. For the bifurcation analysis, it depends what you have in mind. If you want to compute the stability of stationary points and detect bifurcations, it would not take me very long. If you want to compute periodic orbits, more work is needed.

nkoch · July 11, 2022, 11:28pm

Thank you both for your responses. The DDE problem looks to be the thing I need. For the history I need a vector of u1 values from t-tau to t. If I understand the documentation correctly I could do this by setting constant_lags to a vector of lags?

And yes I at the moment I am only interested in the stability of stationary points and the detection of bifurcations.

rveltz · July 12, 2022, 5:52am

Can you post the DDE? some are much easier to handle than others

josuagrw · July 12, 2022, 7:28am

I have a strong suspicion what problem you are trying to solve but it would really help if you could explicitely say how u_1 depends on u_2 in your system.

This is your attempted solution (which intuitely makes sense). However, instead of telling us your solution it would be more helpful to know more about your problem.

See also https://xyproblem.info/

nkoch · July 12, 2022, 6:11pm

Hi,
Thank you for your suggestions. I have managed to create a working example of the kind of DDE that I would like to do bifurcation analysis on:

using Parameters, StatsBase, DifferentialEquations, Plots, DirectionalStatistics, DSP

function FHN_coupled!(du, u, p, t)
    @unpack a, b, c, l, gc = p
    du[1] = u[1] * (a - u[1]) * (u[1] - 1) - u[2] .+ gc.*(u[1] .- u[3])
    du[2] = b * u[1] - c * u[2]
    if t >= l
        du[3] = u[3] * (a - u[3]) * (u[3] - 1) - u[4] .+ gc.*(u[3] .- u[1])
        du[4] = b * u[3] - c * u[4]
        ph1 = hilbert(h[1,:] .- mean(h[1,:]))
        ph2= hilbert(h[2,:] .- mean(h[2,:]))
        u[5] =  1 .- abs.(Circular.var(angle.(ph1) .- angle.(ph2)))
    end
end

u0 = [0.5, 0.0, 0.5, 0.0, 0]
p = (a = -0.1, b =  0.01, c = 0.02, l = 50, gc = 0.025)
tspan = (0.0, 500.0)
datasize = Int(tspan[2])
tsteps = range(tspan[1], tspan[2], length=datasize)
lags = 0:1:200
h = zeros(4, length(lags))
prob = DDEProblem(FHN_coupled!, u0, h, tspan, p; constant_lags=lags)

# run DDE and plot 
sol = solve(prob, Tsit5(), saveat=tsteps)
plot(sol[1,:])
plot!(sol[3,:])
plot!(sol[5,:])

Would this DDE be easy to handle and to implement?

rveltz · July 12, 2022, 9:42pm

I dont think this has constant lags

josuagrw · July 13, 2022, 12:04pm

Please note that your function signature for FHN_coupled! does not match the function signature required for a DDEProblem.

https://diffeq.sciml.ai/stable/types/dde_types/#SciMLBase.DDEProblem

nkoch · July 13, 2022, 3:04pm

Hi,

What makes you think that this has dependent lags? I don’t think my lags depend on the state of the system.

I have attempted to fix my mistake in my DDEProblem function however I have run into some difficulties:
I cannot figure out how to pass a vector of lags in to get and array of the history.

using Parameters, StatsBase, DifferentialEquations, Plots, DirectionalStatistics, DSP

function FHN_coupled(du, u, h, p, t)
    @unpack a, b, c, l, gc, lags = p
    du[1] = u[1] * (a - u[1]) * (u[1] - 1) - u[2] .+ gc.*(u[1] .- u[3])
    du[2] = b * u[1] - c * u[2]
    if t >= l
        du[3] = u[3] * (a - u[3]) * (u[3] - 1) - u[4] .+ gc.*(u[3] .- u[1])
        du[4] = b * u[3] - c * u[4]
        hi = h(p, t .- lags)
        ph1 = hilbert(hi[1,:] .- mean(hi[1,:]))
        ph2= hilbert(hi[3,:] .- mean(hi[3,:]))
        u[5] =  1 .- abs.(Circular.var(angle.(ph1) .- angle.(ph2)))
    end
end


u0 = [0.5, 0.0, 0.5, 0.0, 0]
tspan = (0.0, 500.0)
datasize = Int(tspan[2])
tsteps = range(tspan[1], tspan[2], length=datasize*100)
lags =collect(1.:200.)
p = (a = -0.1, b =  0.01, c = 0.02, l = 50, gc = 0.025, lags=lags)
h(p,t) = zeros(5, length(lags))
tau = [1,2]
prob = DDEProblem(FHN_coupled,u0,h,tspan,p; constant_lags=lags) 

ERROR: LoadError: MethodError: no method matching isless(::Vector{Float64}, ::Float64)
Closest candidates are:
  isless(!Matched::T, ::T) where T<:Union{Float16, Float32, Float64} at C:\Users\nilsk\AppData\Local\Programs\Julia-1.7.2\share\julia\base\float.jl:460
  isless(!Matched::Union{StatsBase.PValue, StatsBase.TestStat}, ::AbstractFloat) at C:\Users\nilsk\.julia\packages\StatsBase\pnfAI\src\statmodels.jl:99
  isless(!Matched::Union{StatsBase.PValue, StatsBase.TestStat}, ::Real) at C:\Users\nilsk\.julia\packages\StatsBase\pnfAI\src\statmodels.jl:90
  ...
Stacktrace:
  [1] <(x::Vector{Float64}, y::Float64)
    @ Base .\operators.jl:352
  [2] (::DelayDiffEq.HistoryFunction{typeof(h), DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}})(p::NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, t::Vector{Float64}, ::Type{Val{0}}; idxs::Nothing)
    @ DelayDiffEq C:\Users\nilsk\.julia\packages\DelayDiffEq\5rXFy\src\history_function.jl:26
  [3] (::DelayDiffEq.HistoryFunction{typeof(h), DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}})(p::NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, t::Vector{Float64}, ::Type{Val{0}}) (repeats 2 times)
    @ DelayDiffEq C:\Users\nilsk\.julia\packages\DelayDiffEq\5rXFy\src\history_function.jl:21
  [4] FHN_coupled(du::Vector{Float64}, u::Vector{Float64}, h::DelayDiffEq.HistoryFunction{typeof(h), DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}}, p::NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, t::Float64)
    @ Main C:\Users\nilsk\AppData\Roaming\JetBrains\PyCharm2021.2\scratches\scratch_52.jl:17
  [5] ODEFunctionWrapper
    @ C:\Users\nilsk\.julia\packages\DelayDiffEq\5rXFy\src\functionwrapper.jl:59 [inlined]
  [6] perform_step!(integrator::DelayDiffEq.DDEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, DDEProblem{Vector{Float64}, Tuple{Float64, Float64}, Vector{Float64}, Tuple{}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DDEFunction{true, typeof(FHN_coupled), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, typeof(h), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), DelayDiffEq.HistoryFunction{typeof(h), DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DelayDiffEq.DDEStats}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), DelayDiffEq.HistoryFunction{typeof(h), DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DelayDiffEq.FPSolver{NLFunctional{Rational{Int64}, Rational{Int64}}, true, Float64, DelayDiffEq.FPFunctionalCache{Vector{Float64}, Vector{Float64}}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Discontinuity{Float64, Int64}, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Tuple{}}, Float64, Int64, DelayDiffEq.HistoryFunction{typeof(h), DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}}, DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Discontinuity{Float64, Int64}}, Vector{Float64}, Float64, Nothing}, cache::OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, repeat_step::Bool)
    @ OrdinaryDiffEq C:\Users\nilsk\.julia\packages\OrdinaryDiffEq\ZgJ9s\src\perform_step\low_order_rk_perform_step.jl:704
  [7] perform_step!
    @ C:\Users\nilsk\.julia\packages\OrdinaryDiffEq\ZgJ9s\src\perform_step\low_order_rk_perform_step.jl:675 [inlined]
  [8] perform_step!
    @ C:\Users\nilsk\.julia\packages\DelayDiffEq\5rXFy\src\integrators\interface.jl:125 [inlined]
  [9] solve!(integrator::DelayDiffEq.DDEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, DDEProblem{Vector{Float64}, Tuple{Float64, Float64}, Vector{Float64}, Tuple{}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DDEFunction{true, typeof(FHN_coupled), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, typeof(h), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardDDEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), DelayDiffEq.HistoryFunction{typeof(h), DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DelayDiffEq.DDEStats}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), DelayDiffEq.HistoryFunction{typeof(h), DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DelayDiffEq.FPSolver{NLFunctional{Rational{Int64}, Rational{Int64}}, true, Float64, DelayDiffEq.FPFunctionalCache{Vector{Float64}, Vector{Float64}}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Float64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Discontinuity{Float64, Int64}, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Tuple{}}, Float64, Int64, DelayDiffEq.HistoryFunction{typeof(h), DelayDiffEq.HistoryODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, NamedTuple{(:a, :b, :c, :l, :gc, :lags), Tuple{Float64, Float64, Float64, Int64, Float64, Vector{Float64}}}, DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{DelayDiffEq.ODEFunctionWrapper{true, typeof(FHN_coupled), typeof(h), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}}, DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Discontinuity{Float64, Int64}}, Vector{Float64}, Float64, Nothing})
    @ DelayDiffEq C:\Users\nilsk\.julia\packages\DelayDiffEq\5rXFy\src\solve.jl:358
 [10] #__solve#32
    @ C:\Users\nilsk\.julia\packages\DelayDiffEq\5rXFy\src\solve.jl:5 [inlined]
 [11] #solve_call#28
    @ C:\Users\nilsk\.julia\packages\DiffEqBase\RHAWf\src\solve.jl:429 [inlined]
 [12] #solve_up#34
    @ C:\Users\nilsk\.julia\packages\DiffEqBase\RHAWf\src\solve.jl:767 [inlined]
 [13] #solve#33
    @ C:\Users\nilsk\.julia\packages\DiffEqBase\RHAWf\src\solve.jl:752 [inlined]
 [14] top-level scope
    @ C:\Users\nilsk\AppData\Roaming\JetBrains\PyCharm2021.2\scratches\scratch_52.jl:36
 [15] include(fname::String)
    @ Base.MainInclude .\client.jl:451
 [16] top-level scope
    @ none:1
in expression starting at C:\Users\nilsk\FHN_coupled_lag.jl:36

My current solution to this is to pass the lags as a parameter and to get the history for each lag and construct the array in that manner. I don’t think that this is ideal though, is there a better way to do this?

using Parameters, StatsBase, DifferentialEquations, Plots, DirectionalStatistics, DSP

function FHN_coupled(du, u, h, p, t)
    @unpack a, b, c, l, gc, lags = p
    du[1] = u[1] * (a - u[1]) * (u[1] - 1) - u[2] .+ gc.*(u[1] .- u[3])
    du[2] = b * u[1] - c * u[2]
    if t >= l
        du[3] = u[3] * (a - u[3]) * (u[3] - 1) - u[4] .+ gc.*(u[3] .- u[1])
        du[4] = b * u[3] - c * u[4]
        hist = zeros(2, length(lags))
        T = typeof(t)
        for i in 1:length(lags)
            hist[1,i] = h(p, t - T(lags[i]))[1]
            hist[2,i] = h(p, t - T(lags[i]))[3]
        end
        ph1 = hilbert(hist[1,:] .- mean(hist[1,:]))
        ph2 = hilbert(hist[2,:] .- mean(hist[2,:]))
        u[5] =  1 .- abs.(Circular.var(angle.(ph1) .- angle.(ph2)))
    end
end


u0 = [0.5, 0.0, 0.5, 0.0, 0]
tspan = (0.0, 500.0)
datasize = Int(tspan[2])
tsteps = range(tspan[1], tspan[2], length=datasize*100)
lags = collect(1:200)
p = (a = -0.1, b =  0.01, c = 0.02, l = 50, gc = 0.025, lags=lags)
h(p,t) = zeros(5)

prob = DDEProblem(FHN_coupled,u0,h,tspan,p; constant_lags=lags)

# run DDE and plot
sol = solve(prob,MethodOfSteps(Tsit5()), saveat=tsteps)
plot(sol[1,:], dpi=300 )
plot!(sol[3,:], dpi=300 )
plot!(sol[5,:], dpi=300 )

josuagrw · July 13, 2022, 4:04pm

My first impression is that

        for i in 1:length(lags)
            hist[1,i] = h(p, t - T(lags[1]))[1]
            hist[2,i] = h(p, t - T(lags[1]))[3]
        end

should probably be

        for i in 1:length(lags)
            hist[1,i] = h(p, t - T(lags[i]))[1]
            hist[2,i] = h(p, t - T(lags[i]))[3]
        end

Is that correct?

Whether this approach is correct in your case, depends on the internals of the function hilbert. Could you explain what that does?

Do you really intend to have 200 discrete lags or are you really approximating an integral over the history?

nkoch · July 13, 2022, 4:32pm

Yes, my apologies, that is correct. I have edited the code in my previous reply to reflect this.

As for the function hilbert(x) it returns the analytic signal for x, namely x_A = x +j \hat{x} where \hat{x} is the Hilbert transform of x. The exact function and its source code can be found here:
https://docs.juliadsp.org/stable/util/#DSP.Util.hilbert

github.com

JuliaDSP/DSP.jl/blob/f54471f002551331cf3dbbee3e4e130fddf08204/src/util.jl#L45-L51


      
          """
              hilbert(x)
          
          
Computes the analytic representation of x, ``x_a = x + j
          \\hat{x}``, where ``\\hat{x}`` is the Hilbert transform of x,
          along the first dimension of x.
          """

I do intend to have this many lags, as I would like to have a decent estimate of the analytic signal based on the history.

josuagrw · July 13, 2022, 5:18pm

Ok, since you want to use a discrete Hilbert transform, it makes sense to produce discrete samples of the history.

It’s up to you if you want to use quadrature (e.g. QuadGK.jl) to calculate the time average (mean) of the history more accurately.