Skip to content

Commit

Permalink
correct docs
Browse files Browse the repository at this point in the history
  • Loading branch information
@rveltz
rveltz committed Sep 24, 2023
1 parent a1df562 commit 923d05e
Show file tree
Hide file tree
Showing 5 changed files with 74 additions and 80 deletions.
2 changes: 0 additions & 2 deletions docs/Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
[deps]
BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
DDEBifurcationKit = "954e4062-bdb8-4e3f-9eee-d47105dd3e65"
DiffEqOperators = "9fdde737-9c7f-55bf-ade8-46b3f136cc48"
DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
Expand All @@ -16,7 +15,6 @@ SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"

[compat]
BifurcationKit = "0.3"
DiffEqOperators = "4.34.0"
DifferentialEquations = "6.19.0, 7.1.0"
Documenter = "0.27"
Setfield = "0.5.0, 0.7.0, 0.8.0, 1.0.0"
Expand Down
51 changes: 25 additions & 26 deletions docs/src/tutorials/dde/Humphries.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ $$x^{\prime}(t)=-\gamma x(t)-\kappa_1 x\left(t-a_1-c x(t)\right)-\kappa_2 x\left
We first instantiate the model

```@example TUTHumphries
using Revise, DDEBifurcationKit, Parameters, Setfield, Plots
using Revise, DDEBifurcationKit, Parameters, Plots
using BifurcationKit
const BK = BifurcationKit
Expand Down Expand Up @@ -43,7 +43,7 @@ We then compute the branch

```@example TUTHumphries
optn = NewtonPar(verbose = true, eigsolver = DDE_DefaultEig())
opts = ContinuationPar(p_max = 13., p_min = 0., newtonOptions = optn, ds = -0.01, detectBifurcation = 3, nev = 3, )
opts = ContinuationPar(p_max = 13., p_min = 0., newton_options = optn, ds = -0.01, detect_bifurcation = 3, nev = 3, )
br = continuation(prob, PALC(), opts; verbosity = 0, bothside = true)
```

Expand All @@ -60,13 +60,13 @@ We tell the solver to consider br.specialpoint[2] and continue it.

```@example TUTHumphries
brhopf = continuation(br, 2, (@lens _.κ2),
setproperties(br.contparams, detectBifurcation = 2, dsmax = 0.04, max_steps = 230, p_max = 5., p_min = -1.,ds = -0.02);
setproperties(br.contparams, detect_bifurcation = 2, dsmax = 0.04, max_steps = 230, p_max = 5., p_min = -1.,ds = -0.02);
verbosity = 0, plot = false,
# we disable detection of Bautin bifurcation as the
# Hopf normal form is not implemented for SD-DDE
detectCodim2Bifurcation = 0,
detect_codim2_bifurcation = 0,
bothside = true,
startWithEigen = true)
start_with_eigen = true)
scene = plot(brhopf, vars = (:κ1, :κ2))
```
Expand All @@ -78,26 +78,26 @@ We compute the branch of periodic orbits from the Hopf bifurcation points using
```julia
# 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, detectBifurcation = 0, plot_every_step = 20, save_sol_every_step=1)
@set! opts_po_cont.newtonOptions.tol = 1e-9
@set! opts_po_cont.newtonOptions.verbose = true

# arguments for periodic orbits
args_po = ( recordFromSolution = (x, p) -> begin
xtt = BK.get_periodic_orbit(p.prob, x, nothing)
_max = maximum(xtt[1,:])
_min = minimum(xtt[1,:])
return (amp = _max - _min,
period = getperiod(p.prob, x, nothing))
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

# arguments for periodic orbits
args_po = ( record_from_solution = (x, p) -> begin
xtt = BK.get_periodic_orbit(p.prob, x, nothing)
_max = maximum(xtt[1,:])
_min = minimum(xtt[1,:])
return (amp = _max - _min,
period = getperiod(p.prob, x, nothing))
end,
plot_solution = (x, p; k...) -> begin
xtt = BK.get_periodic_orbit(p.prob, x, nothing)
plot!(xtt.t, xtt[1,:]; label = "x", k...)
plot!(br; subplot = 1, putspecialptlegend = false)
end,
plot_solution = (x, p; k...) -> begin
xtt = BK.get_periodic_orbit(p.prob, x, nothing)
plot!(xtt.t, xtt[1,:]; label = "x", k...)
plot!(br; subplot = 1, putspecialptlegend = false)
end,
normC = norminf)

probpo = PeriodicOrbitOCollProblem(200, 2; N = 1)
normC = norminf)

probpo = PeriodicOrbitOCollProblem(200, 2; N = 1, jacobian = BK.AutoDiffDense())
br_pocoll = continuation(
br, 2, opts_po_cont,
probpo;
Expand All @@ -107,8 +107,7 @@ br_pocoll = continuation(
args_po...,
ampfactor = 1/0.467829783456199 * 0.1,
δp = 0.01,
callback_newton = BK.cbMaxNorm(10.0)
end
callback_newton = BK.cbMaxNorm(10.0),
)
```

Expand Down
23 changes: 12 additions & 11 deletions docs/src/tutorials/dde/HutchinsonDiff.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,15 +18,10 @@ where $a>0,d>0$.
We start by discretizing the above PDE based on finite differences.

```@example TUTHut
using Revise, DDEBifurcationKit, Parameters, Setfield, LinearAlgebra, Plots, SparseArrays
using Revise, DDEBifurcationKit, Parameters, LinearAlgebra, Plots, SparseArrays
using BifurcationKit
const BK = BifurcationKit
# sup norm
norminf(x) = norm(x, Inf)
using DiffEqOperators
function Hutchinson(u, ud, p)
@unpack a,d,Δ = p
d .* (Δ*u) .- a .* ud[1] .* (1 .+ u)
Expand All @@ -44,9 +39,15 @@ We can now instantiate the model
Nx = 200; Lx = pi/2;
X = -Lx .+ 2Lx/Nx*(0:Nx-1) |> collect
# boundary condition
Q = Neumann0BC(X[2]-X[1])
Δ = sparse(CenteredDifference(2, 2, X[2]-X[1], Nx) * Q)[1]
function DiffOp(N, lx; order = 2)
hx = 2lx/N
Δ = spdiagm(0 => -2ones(N), 1 => ones(N-1), -1 => ones(N-1) )
Δ[1,1]=-1; Δ[end,end]=-1
Δ = Δ / hx^2
return Δ
end
Δ = DiffOp(Nx, Lx)
nothing #hide
```

Expand All @@ -60,7 +61,7 @@ x0 = zeros(Nx)
prob = ConstantDDEBifProblem(Hutchinson, delaysF, x0, pars, (@lens _.a))
optn = NewtonPar(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)
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)
br = continuation(prob, PALC(), opts; verbosity = 0, plot = false, normC = norminf)
br
```
Expand Down Expand Up @@ -91,7 +92,7 @@ end
prob2 = ConstantDDEBifProblem(Hutchinson, delaysF, x0, pars, (@lens _.a); J = JacHutchinson)
optn = NewtonPar(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)
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)
br = continuation(prob2, PALC(), opts; verbosity = 1, plot = true, normC = norminf)
br
```
Expand Down
51 changes: 23 additions & 28 deletions docs/src/tutorials/dde/neuron.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,13 +17,10 @@ $$\left\{\begin{array}{l}
We first instantiate the model

```@example TUTneuron
using Revise, DDEBifurcationKit, Parameters, Setfield, LinearAlgebra, Plots
using Revise, DDEBifurcationKit, Parameters, LinearAlgebra, Plots
using BifurcationKit
const BK = BifurcationKit
# sup norm
norminf(x) = norm(x, Inf)
function neuronVF(x, xd, p)
@unpack κ, β, a12, a21, τs, τ1, τ2 = p
[
Expand All @@ -40,7 +37,7 @@ x0 = [0.01, 0.001]
prob = ConstantDDEBifProblem(neuronVF, delaysF, x0, pars, (@lens _.τs))
optn = NewtonPar(verbose = true, eigsolver = DDE_DefaultEig())
opts = ContinuationPar(p_max = 13., p_min = 0., newtonOptions = optn, ds = -0.01, detectBifurcation = 3, nev = 5, dsmax = 0.2, n_inversion = 4)
opts = ContinuationPar(p_max = 13., 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 = 1, plot = true, bothside = true, normC = norminf)
```

Expand All @@ -55,7 +52,7 @@ scene = plot(br)
As in [BifurcationKit.jl](https://github.com/rveltz/BifurcationKit.jl), it is straightforward to compute the normal forms.

```@example TUTneuron
hopfpt = BK.getNormalForm(br, 2)
hopfpt = BK.get_normal_form(br, 2)
```

## Continuation of Hopf points
Expand All @@ -64,18 +61,16 @@ We follow the Hopf points in the parameter plane $(a_{21},\tau_s)$. We tell the
```@example TUTneuron
# continuation of the first Hopf point
brhopf = continuation(br, 3, (@lens _.a21),
setproperties(br.contparams, detectBifurcation = 1, dsmax = 0.04, max_steps = 230, p_max = 15., p_min = -1.,ds = -0.02);
verbosity = 2, plot = true,
detectCodim2Bifurcation = 2,
setproperties(br.contparams, detect_bifurcation = 1, dsmax = 0.04, max_steps = 230, p_max = 15., p_min = -1.,ds = -0.02);
detect_codim2_bifurcation = 2,
# bothside = true,
startWithEigen = true)
start_with_eigen = true)
# continuation of the second Hopf point
brhopf2 = continuation(br, 2, (@lens _.a21),
setproperties(br.contparams, detectBifurcation = 1, dsmax = 0.1, max_steps = 56, p_max = 15., p_min = -1.,ds = -0.01, n_inversion = 4);
verbosity = 2, plot = true,
detectCodim2Bifurcation = 2,
startWithEigen = true,
setproperties(br.contparams, detect_bifurcation = 1, dsmax = 0.1, max_steps = 56, p_max = 15., p_min = -1.,ds = -0.01, n_inversion = 4);
detect_codim2_bifurcation = 2,
start_with_eigen = true,
bothside=true)
scene = plot(brhopf, brhopf2, legend = :top)
Expand All @@ -94,12 +89,12 @@ We then compute the branch of periodic orbits from the Hopf bifurcation points u

```@example TUTneuron
# continuation parameters
opts_po_cont = ContinuationPar(dsmax = 0.1, ds= -0.0001, dsmin = 1e-4, p_max = 10., p_min=-0., max_steps = 120, detectBifurcation = 0, save_sol_every_step=1)
@set! opts_po_cont.newtonOptions.tol = 1e-8
@set! opts_po_cont.newtonOptions.verbose = false
opts_po_cont = ContinuationPar(dsmax = 0.1, ds= -0.0001, dsmin = 1e-4, p_max = 10., p_min=-0., max_steps = 120, detect_bifurcation = 0, save_sol_every_step=1)
@set! opts_po_cont.newton_options.tol = 1e-8
@set! opts_po_cont.newton_options.verbose = false
# arguments for periodic orbits
args_po = ( recordFromSolution = (x, p) -> begin
args_po = ( record_from_solution = (x, p) -> begin
xtt = BK.get_periodic_orbit(p.prob, x, nothing)
return (max = maximum(xtt[1,:]),
min = minimum(xtt[1,:]),
Expand All @@ -113,16 +108,16 @@ args_po = ( recordFromSolution = (x, p) -> begin
end,
normC = norminf)
probpo = PeriodicOrbitOCollProblem(60, 4; N = 2)
br_pocoll = @time continuation(
br2, 1, opts_po_cont,
probpo;
verbosity = 0, plot = false,
args_po...,
δp = 0.001,
normC = norminf,
)
scene = plot(br2, br_pocoll)
probpo = PeriodicOrbitOCollProblem(60, 4; N = 2, jacobian = BK.AutoDiffDense())
br_pocoll = @time continuation(
br2, 1, opts_po_cont,
probpo;
verbosity = 0, plot = false,
args_po...,
δp = 0.001,
normC = norminf,
)
scene = plot(br2, br_pocoll)
```

We can plot the periodic orbit as they approach the homoclinic point.
Expand Down
27 changes: 14 additions & 13 deletions docs/src/tutorials/dde/neuronV2.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ where $g(z)=[\tanh (z-1)+\tanh (1)] \cosh (1)^2$.
We first instantiate the model

```@example TUTneuron2
using Revise, DDEBifurcationKit, Parameters, Setfield, Plots
using Revise, DDEBifurcationKit, Parameters, Plots
using BifurcationKit
const BK = BifurcationKit
Expand All @@ -42,7 +42,7 @@ x0 = [0.01, 0.001]
prob = ConstantDDEBifProblem(neuron2VF, delaysF, x0, pars, (@lens _.a))
optn = NewtonPar(verbose = false, eigsolver = DDE_DefaultEig(maxit=100))
opts = ContinuationPar(p_max = 1., p_min = 0., newtonOptions = optn, ds = 0.01, detectBifurcation = 3, nev = 9, dsmax = 0.2, n_inversion = 4)
opts = ContinuationPar(p_max = 1., p_min = 0., newton_options = optn, ds = 0.01, detect_bifurcation = 3, nev = 9, dsmax = 0.2, n_inversion = 4)
br = continuation(prob, PALC(tangent=Bordered()), opts)
```

Expand All @@ -57,7 +57,7 @@ scene = plot(br)
As in [BifurcationKit.jl](https://github.com/rveltz/BifurcationKit.jl), it is straightforward to compute the normal forms.

```@example TUTneuron2
hopfpt = BK.getNormalForm(br, 2)
hopfpt = BK.get_normal_form(br, 2)
```

## Continuation of Hopf points
Expand All @@ -66,18 +66,18 @@ We follow the Hopf points in the parameter plane $(a,c)$. We tell the solver to
```@example TUTneuron2
# continuation of the first Hopf point
brhopf = continuation(br, 1, (@lens _.c),
setproperties(br.contparams, detectBifurcation = 1, dsmax = 0.01, max_steps = 100, p_max = 1.1, p_min = -0.1,ds = 0.01, n_inversion = 2);
setproperties(br.contparams, detect_bifurcation = 1, dsmax = 0.01, max_steps = 100, p_max = 1.1, p_min = -0.1,ds = 0.01, n_inversion = 2);
verbosity = 0,
detectCodim2Bifurcation = 2,
detect_codim2_bifurcation = 2,
bothside = true,
startWithEigen = true)
start_with_eigen = true)
brhopf2 = continuation(br, 2, (@lens _.c),
setproperties(br.contparams, detectBifurcation = 1, dsmax = 0.01, max_steps = 100, p_max = 1.1, p_min = -0.1,ds = -0.01);
setproperties(br.contparams, detect_bifurcation = 1, dsmax = 0.01, max_steps = 100, p_max = 1.1, p_min = -0.1,ds = -0.01);
verbosity = 0,
detectCodim2Bifurcation = 2,
detect_codim2_bifurcation = 2,
bothside = true,
startWithEigen = true)
start_with_eigen = true)
scene = plot(brhopf, vars = (:a, :c), xlims = (0,0.7), ylims = (0,1))
plot!(scene, brhopf2, vars = (:a, :c), xlims = (-0,0.7), ylims = (-0.1,1))
Expand All @@ -94,14 +94,15 @@ br2 = continuation(prob2, PALC(), setproperties(opts, p_max = 1.22);)
# change tolerance for avoiding error computation of the EV
opts_fold = br.contparams
@set! opts_fold.newtonOptions.eigsolver.σ = 1e-7
@set! opts_fold.newton_options.eigsolver.σ = 1e-7
brfold = continuation(br2, 3, (@lens _.a),
setproperties(opts_fold; detectBifurcation = 1, dsmax = 0.01, max_steps = 100, p_max = 0.6, p_min = -0.6,ds = -0.01, n_inversion = 2, tol_stability = 1e-6);
setproperties(opts_fold; detect_bifurcation = 1, dsmax = 0.01, max_steps = 70, p_max = 0.6, p_min = -0.6,ds = -0.01, n_inversion = 2, tol_stability = 1e-6);
verbosity = 1, plot = true,
detectCodim2Bifurcation = 2,
detect_codim2_bifurcation = 2,
update_minaug_every_step = 1,
bothside = false,
startWithEigen = true)
start_with_eigen = true)
scene = plot(brfold, vars = (:a, :c), branchlabel = "Fold")
plot!(scene, brhopf, vars = (:a, :c), branchlabel = "Hopf")
Expand Down

0 comments on commit 923d05e

Please sign in to comment.