Skip to content

Commit

Permalink
update for [email protected]
Browse files Browse the repository at this point in the history
  • Loading branch information
@rveltz
rveltz committed Oct 6, 2024
1 parent 6636498 commit b0a48cd
Show file tree
Hide file tree
Showing 20 changed files with 130 additions and 142 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"


[compat]
BifurcationKit = "^0.3.5"
BifurcationKit = "^0.4.0"
DocStringExtensions = "^0.9"
ForwardDiff = "^0.10"
NonlinearEigenproblems = "^1.0.1"
Expand Down
4 changes: 2 additions & 2 deletions examples/HutchinsonDiffusion.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ Q = Neumann0BC(X[2]-X[1])
pars = (a = 0.5, d = 1, τ = 1.0, Δ = Δ, N = Nx)
x0 = zeros(Nx)

prob = ConstantDDEBifProblem(Hutchinson, delaysF, x0, pars, (@lens _.a))
prob = ConstantDDEBifProblem(Hutchinson, delaysF, x0, pars, (@optic _.a))

optn = NewtonPar(verbose = true, eigsolver = DDE_DefaultEig())
opts = ContinuationPar(p_max = 10., p_min = 0., newtonOptions = optn, ds = 0.01, detectBifurcation = 3, nev = 5, dsmax = 0.2, n_inversion = 4)
Expand All @@ -46,7 +46,7 @@ function JacHutchinson(u, p)
return J0, [J1]
end

prob2 = ConstantDDEBifProblem(Hutchinson, delaysF, x0, pars, (@lens _.a); J = JacHutchinson)
prob2 = ConstantDDEBifProblem(Hutchinson, delaysF, x0, pars, (@optic _.a); J = JacHutchinson)

optn = NewtonPar(verbose = true, eigsolver = DDE_DefaultEig())
opts = ContinuationPar(p_max = 10., p_min = 0., newtonOptions = optn, ds = 0.01, detectBifurcation = 3, nev = 5, dsmax = 0.2, n_inversion = 4)
Expand Down
22 changes: 8 additions & 14 deletions examples/TravisWebb.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,16 +3,13 @@ cd(@__DIR__)
cd("..")
# using Pkg, LinearAlgebra, Test
# pkg"activate ."
using Revise, DDEBifurcationKit, Parameters, Setfield, LinearAlgebra, Plots, SparseArrays
using Revise, DDEBifurcationKit, LinearAlgebra, Plots, SparseArrays
using BifurcationKit
const BK = BifurcationKit
const DDEBK = DDEBifurcationKit


using DiffEqOperators

function TW(u, ud, p)
@unpack a,Δ = p
(;a,Δ) = p
(Δ*u) .+ u .- a .* ud[1] .* (1 .+ u)
end

Expand All @@ -23,13 +20,13 @@ Nx = 500; Lx = pi/2;
X = -Lx .+ 2Lx/Nx*(0:Nx-1) |> collect

# boundary condition
Q = Dirichlet0BC(X[2]-X[1]|>typeof)
Δ = sparse(CenteredDifference(2, 2, X[2]-X[1], Nx) * Q)[1]
h = 2Lx/Nx
Δ = spdiagm(0 => -2ones(Nx), 1 => ones(Nx-1), -1 => ones(Nx-1) ) / h^2; Δ[1,1]=Δ

pars = (a = 0.5, τ = 1.0, Δ = Δ, N = Nx)
x0 = zeros(Nx)

prob = ConstantDDEBifProblem(TW, delaysF, x0, pars, (@lens _.a))
prob = ConstantDDEBifProblem(TW, delaysF, x0, pars, (@optic _.a))

optn = NewtonPar(verbose = true, eigsolver = DDE_DefaultEig())
opts = ContinuationPar(p_max = 10., p_min = 0., newton_options = optn, ds = 0.01, detect_bifurcation = 3, nev = 5, dsmax = 0.2, n_inversion = 4)
Expand All @@ -42,14 +39,14 @@ plot(hopfpt.ζ |> real)
################################################################################
# case where we specify the jacobian
function JacTW(u, p)
@unpack a,Δ = p
(;a,Δ) = p
# we compute the jacobian at the steady state
J0 = Δ + I .- a .* Diagonal(u)
J1 = -a .* Diagonal(1 .+ u)
return J0, [J1]
end

prob2 = ConstantDDEBifProblem(TW, delaysF, x0, pars, (@lens _.a); J = JacTW)
prob2 = ConstantDDEBifProblem(TW, delaysF, x0, pars, (@optic _.a); J = JacTW)

optn = NewtonPar(verbose = false, eigsolver = DDE_DefaultEig())
opts = ContinuationPar(p_max = 2., p_min = 0., newton_options = optn, ds = 0.01, detect_bifurcation = 3, nev = 5, dsmax = 0.2, n_inversion = 6)
Expand All @@ -63,9 +60,6 @@ function TW_DE(du,u,h,p,t)
end

h(p, t) = real(hopfpt.ζ) .*(1+ 0.01cos(t/4))
prob_de = DDEProblem(TW_DE,h(pars,0),h,(0.,120.),setproperties(pars, a = br.specialpoint[1].param + 0.1); constant_lags=delaysF(pars))
alg = MethodOfSteps(Rosenbrock23())
sol = solve(prob_de,alg, progress=true)
plot(sol.t, sol[1,:])


heatmap(sol.t,1:Nx,reduce(hcat,sol.u))
15 changes: 8 additions & 7 deletions examples/delayed-log.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ using Plots

function delayedlogVF(x, xd, p)
(;λ) = p
y = xd[1][1]
y = xd[1,1]
[
(λ - y) * x[1]
]
Expand All @@ -26,24 +26,25 @@ end
pars = (λ=1.1,b=0.)
x0 = [1.]

prob = ConstantDDEBifProblem(delayedlogVF, delaysF, x0, pars, (@lens _.λ), record_from_solution=(x,p)-> (x=x[1], _x=1))
prob = ConstantDDEBifProblem(delayedlogVF, delaysF, x0, pars, (@optic _.λ), record_from_solution=(x,p;k...)-> (x=x[1], _x=1))

optn = NewtonPar(verbose = false, eigsolver = DDE_DefaultEig())
opts = ContinuationPar(p_max = 2., p_min = 0., newton_options = optn, ds = 0.01, detect_bifurcation = 3, nev = 6, n_inversion = 6 )
br = BK.continuation(prob, PALC(), opts; verbosity = 1, plot = true, bothside = false)
plot(br)

get_normal_form(br, 1)
################################################################################
# computation periodic orbit

# continuation parameters
opts_po_cont = ContinuationPar(dsmax = 0.1, ds= 0.001, dsmin = 1e-4, p_max = 10., p_min=-5., max_steps = 130,
nev = 3, tol_stability = 1e-8, detect_bifurcation = 0, plot_every_step = 2, save_sol_every_step=1)
@set! opts_po_cont.newton_options.tol = 1e-8
@set! opts_po_cont.newton_options.verbose = true
@reset opts_po_cont.newton_options.tol = 1e-8
@reset opts_po_cont.newton_options.verbose = true

# arguments for periodic orbits
args_po = ( record_from_solution = (x, p) -> begin
args_po = ( record_from_solution = (x, p;k...) -> begin
xtt = BK.get_periodic_orbit(p.prob, x, nothing)
return (max = maximum(xtt[1,:]),
min = minimum(xtt[1,:]),
Expand All @@ -57,8 +58,8 @@ args_po = ( record_from_solution = (x, p) -> begin
normC = norminf)

probpo = PeriodicOrbitOCollProblem(40, 4; N = 1)
# probpo = PeriodicOrbitTrapProblem(M = 2000, jacobian = :DenseAD, N = 2)
br_pocoll = @time continuation(
# probpo = PeriodicOrbitTrapProblem(M = 2000, jacobian = :DenseAD, N = 2)
br_pocoll = @time continuation(
br, 1, opts_po_cont,
# PeriodicOrbitOCollProblem(100, 4);
probpo;
Expand Down
10 changes: 5 additions & 5 deletions examples/humpries.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ end
pars = (κ1=0., κ2=2.3, a1=1.3, a2=6, γ=4.75, c=1.)
x0 = zeros(1)

prob = DDEBK.SDDDEBifProblem(humpriesVF, delaysF, x0, pars, (@lens _.κ1))
prob = DDEBK.SDDDEBifProblem(humpriesVF, delaysF, x0, pars, (@optic _.κ1))

optn = NewtonPar(verbose = true, eigsolver = DDE_DefaultEig())
opts = ContinuationPar(p_max = 13., p_min = 0., newton_options = optn, ds = -0.01, detect_bifurcation = 3, nev = 3, )
Expand All @@ -36,7 +36,7 @@ br = continuation(prob, PALC(), opts; verbosity = 1, plot = true, bothside = tru

plot(br)
################################################################################
brhopf = continuation(br, 2, (@lens _.κ2),
brhopf = continuation(br, 2, (@optic _.κ2),
ContinuationPar(br.contparams, detect_bifurcation = 2, dsmax = 0.04, max_steps = 230, p_max = 5., p_min = -1.,ds = -0.02);
verbosity = 2, plot = true,
detect_codim2_bifurcation = 0,
Expand All @@ -50,11 +50,11 @@ plot(brhopf, vars = (:κ1, :κ2))
# continuation parameters
opts_po_cont = ContinuationPar(dsmax = 0.05, ds= 0.001, dsmin = 1e-4, p_max = 12., p_min=-5., max_steps = 3000,
nev = 3, tol_stability = 1e-8, detect_bifurcation = 0, plot_every_step = 20, save_sol_every_step=1)
@set! opts_po_cont.newton_options.tol = 1e-9
@set! opts_po_cont.newton_options.verbose = true
@reset opts_po_cont.newton_options.tol = 1e-9
@reset opts_po_cont.newton_options.verbose = true

# arguments for periodic orbits
args_po = ( record_from_solution = (x, p) -> begin
args_po = ( record_from_solution = (x, p; k...) -> begin
xtt = BK.get_periodic_orbit(p.prob, x, nothing)
_max = maximum(xtt[1,:])
_min = minimum(xtt[1,:])
Expand Down
22 changes: 11 additions & 11 deletions examples/ikeda.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ delaysF(par) = [1.0]
pars = (Λ=0.1,b=0.)
x0 = [-sqrt(pi)]

prob = ConstantDDEBifProblem(ikedaVF, delaysF, x0, pars, (@lens _.Λ), record_from_solution=(x,p)-> (x=x[1], _x=1))
prob = ConstantDDEBifProblem(ikedaVF, delaysF, x0, pars, (@optic _.Λ), record_from_solution=(x,p;k...)-> (x=x[1], _x=1))

optn = NewtonPar(verbose = false, eigsolver = DDE_DefaultEig())
opts = ContinuationPar(p_max = 2., p_min = 0., newton_options = optn, ds = 0.01, detect_bifurcation = 3, nev = 4, n_inversion = 12 )
Expand All @@ -37,11 +37,11 @@ BK.get_normal_form(br, 1) # l1= -0.0591623057, b = 0.09293196762669392
# continuation parameters
opts_po_cont = ContinuationPar(dsmax = 0.2, ds= -0.001, dsmin = 1e-4, p_max = 10., p_min=-5., max_steps = 40,
nev = 3, tol_stability = 1e-8, detect_bifurcation = 0, plot_every_step = 1, save_sol_every_step=1)
@set! opts_po_cont.newton_options.tol = 1e-7
@set! opts_po_cont.newton_options.verbose = true
@reset opts_po_cont.newton_options.tol = 1e-7
@reset opts_po_cont.newton_options.verbose = true

# arguments for periodic orbits
args_po = ( record_from_solution = (x, p) -> begin
args_po = ( record_from_solution = (x, p; k...) -> begin
xtt = BK.get_periodic_orbit(p.prob, x, nothing)
return (max = maximum(xtt[1,:]),
min = minimum(xtt[1,:]),
Expand Down Expand Up @@ -79,18 +79,18 @@ plot(br, br_pocoll)
plot(br_pocoll, vars = (:param, :period))

################################################################################
using DifferentialEquations
using DifferentialEquations

function ikedaVF_DE(du,u,h,p,t)
(; Λ) = p
du[1] = -pi/2 + Λ/2 * h(p,t-1)[1]^2;
end

u0 = -2ones(1)
h(p, t) = -2ones(1) .+ 0.001cos(t/4)
h(p,t) = br_pocoll.orbit(t)
prob_de = DDEProblem(ikedaVF_DE,u0,h,(0.,1240.),setproperties(pars, Λ = br.specialpoint[1].param + 0.01); constant_lags=delaysF(pars))
alg = MethodOfSteps(Rosenbrock23())
sol = solve(prob_de,alg)
plot(plot(sol, xlims = (1000,1020)), plot(sol))
h(p, t) = -2ones(1) .+ 0.001cos(t/4)
h(p,t) = br_pocoll.orbit(t)
prob_de = DDEProblem(ikedaVF_DE,u0,h,(0.,1240.),setproperties(pars, Λ = br.specialpoint[1].param + 0.01); constant_lags=delaysF(pars))
alg = MethodOfSteps(Rosenbrock23())
sol = solve(prob_de,alg)
plot(plot(sol, xlims = (1000,1020)), plot(sol))

32 changes: 15 additions & 17 deletions examples/neuron.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,28 +18,28 @@ delaysF(par) = [par.τ1, par.τ2, par.τs]
pars = (κ = 0.5, β = -1, a12 = 1, a21 = 0.5, τ1 = 0.2, τ2 = 0.2, τs = 1.5)
x0 = [0.01, 0.001]

prob = ConstantDDEBifProblem(neuronVF, delaysF, x0, pars, (@lens _.τs); record_from_solution = (x,p)->(x1 = x[1], x2 = x[2]))
prob = ConstantDDEBifProblem(neuronVF, delaysF, x0, pars, (@optic _.τs); record_from_solution = (x,p; k...) -> (x1 = x[1], x2 = x[2]))

optn = NewtonPar(verbose = false, eigsolver = DDE_DefaultEig())
opts = ContinuationPar(p_max = 5., p_min = 0., newton_options = optn, ds = -0.01, detect_bifurcation = 3, nev = 5, dsmax = 0.2, n_inversion = 4)
br = continuation(prob, PALC(), opts; verbosity = 0, plot = true, bothside = true, normC = norminf)

plot(br)
################################################################################
prob2 = ConstantDDEBifProblem(neuronVF, delaysF, x0, pars, (@lens _.a21); record_from_solution = prob.recordFromSolution)
prob2 = ConstantDDEBifProblem(neuronVF, delaysF, x0, pars, (@optic _.a21); record_from_solution = prob.recordFromSolution)
br2 = BK.continuation(prob2, PALC(), ContinuationPar(opts, ds = 0.1, p_max = 4., n_inversion=8); verbosity = 0, plot = true, normC = norminf)

# @set! br2.contparams.newton_options.eigsolver.σ = 1e-5
# @reset br2.contparams.newton_options.eigsolver.σ = 1e-5
BK.get_normal_form(br2, 1)
#Hopf l1 ≈ −0.0601.
BK.get_normal_form(br2, 2)

opts_pitch = ContinuationPar(br2.contparams, ds = 0.005, dsmax = 0.03, n_inversion = 4)
@set! opts_pitch.newton_options.eigsolver.γ = 0.001
@reset opts_pitch.newton_options.eigsolver.γ = 0.001
br_pitchfork = continuation(br2, 2, opts_pitch)
plot(br2, br_pitchfork)
################################################################################
brhopf = continuation(br, 2, (@lens _.a21),
brhopf = continuation(br, 2, (@optic _.a21),
ContinuationPar(br.contparams, detect_bifurcation = 1, dsmax = 0.04, max_steps = 230, p_max = 15., p_min = -1.,ds = -0.02);
verbosity = 2, plot = true,
detect_codim2_bifurcation = 2,
Expand All @@ -49,34 +49,32 @@ brhopf = continuation(br, 2, (@lens _.a21),
plot(brhopf, vars = (:a21, :τs))
plot(brhopf, vars = (:τs, :ω))

brhopf2 = continuation(br, 2, (@lens _.a21),
brhopf2 = continuation(br, 2, (@optic _.a21),
ContinuationPar(br.contparams, detect_bifurcation = 1, dsmax = 0.1, max_steps = 56, p_max = 1.5, p_min = -1.,ds = -0.01, n_inversion = 4);
verbosity = 2, plot = true,
detect_codim2_bifurcation = 2,
start_with_eigen = true,
bothside=true)

plot(brhopf, brhopf2, legend = :top)


################################################################################
# computation periodic orbit
opts_po_cont = ContinuationPar(dsmax = 0.1, ds= 0.0001, dsmin = 1e-4, p_max = 10., p_min=-5., max_steps = 20,
nev = 3, tol_stability = 1e-8, detect_bifurcation = 0, plot_every_step = 2, save_sol_every_step=1, detect_fold = true)
@set! opts_po_cont.newton_options.tol = 1e-8
@set! opts_po_cont.newton_options.verbose = true
@reset opts_po_cont.newton_options.tol = 1e-8
@reset opts_po_cont.newton_options.verbose = true

# arguments for periodic orbits
args_po = ( record_from_solution = (x, p) -> begin
args_po = ( record_from_solution = (x, p; k...) -> begin
_par = BK.getparams(p.prob)
xtt = BK.get_periodic_orbit(p.prob, x, @set _par.a21=p.p)
xtt = BK.get_periodic_orbit(p.prob, x, p.p)
return (max = maximum(xtt[1,:]),
min = minimum(xtt[1,:]),
period = getperiod(p.prob, x, @set _par.a21=p.p))
period = getperiod(p.prob, x, p.p))
end,
plot_solution = (x, p; k...) -> begin
_par = BK.getparams(p.prob)
xtt = BK.get_periodic_orbit(p.prob, x, @set _par.a21=p.p)
xtt = BK.get_periodic_orbit(p.prob, x, p.p)
plot!(xtt.t, xtt[1,:]; label = "V1", k...)
plot!(xtt.t, xtt[2,:]; label = "V2", k...)
plot!(br2; subplot = 1, putspecialptlegend = false)
Expand Down Expand Up @@ -124,14 +122,14 @@ plot(br2, br_pocoll, br_pitchfork);plot!(br_pocoll, vars = (:param,:min))
using DifferentialEquations

function neuron_DE(du,u,h,p,t)
@unpack κ, β, a12, a21, τs, τ1, τ2 = p
(;κ, β, a12, a21, τs, τ1, τ2) = p
du[1] = -κ * u[1] + β * tanh(h(p, t-τs)[1]) + a12 * tanh(h(p, t-τ2)[2])
du[2] = -κ * u[2] + β * tanh(h(p, t-τs)[2]) + a21 * tanh(h(p, t-τ1)[1])
end

h(p, t) = -0*ones(2) .+ 0.25sin(t/4)
prob_de = DDEProblem(neuron_DE,h(pars, 0),h,(0.,20000.),ContinuationPar(pars, a21 = br2.specialpoint[1].param + 0.01); constant_lags=delaysF(pars), abstol = 1e-10, reltol = 1e-9)
prob_de = DDEProblem(neuron_DE,h(pars, 0),h,(0.,20000.),(pars..., a21 = br2.specialpoint[1].param + 0.01); constant_lags=delaysF(pars), abstol = 1e-10, reltol = 1e-9)
alg = MethodOfSteps(Tsit5())
sol = solve(prob_de,alg)
plot(plot(sol, xlims = (sol.t[end]-100,sol.t[end]), vars=(:t,:1),ylims=(-0.1,0.1)), plot(sol.t,sol[1,:]),title = "a21=$(sol.prob.p.a21)")
plot(plot(sol, xlims = (sol.t[end]-100,sol.t[end]), vars=(:t,:1),ylims=(-0.1,0.1)), plot(sol.t,sol[1,:]),title = "a21=$(sol.prob.p.a21)")

Loading

0 comments on commit b0a48cd

Please sign in to comment.