Library
Parameters
BifurcationKit.NewtonPar
— Typestruct NewtonPar{T, L<:BifurcationKit.AbstractLinearSolver, E<:AbstractEigenSolver}
Returns a variable containing parameters to affect the newton
algorithm when solving F(x) = 0
.
Arguments (with default values):
tol::Any
: absolute tolerance forF(x)
Default: 1.0e-12max_iterations::Int64
: number of Newton iterations Default: 25verbose::Bool
: display Newton iterations? Default: falselinsolver::BifurcationKit.AbstractLinearSolver
: linear solver, must be<: AbstractLinearSolver
Default: DefaultLS()eigsolver::AbstractEigenSolver
: eigen solver, must be<: AbstractEigenSolver
Default: DefaultEig()linesearch::Bool
: Default: falseα::Any
: Default: convert(typeof(tol), 1.0)αmin::Any
: Default: convert(typeof(tol), 0.001)
Arguments for line search (Armijo)
linesearch = false
: use line search algorithm (i.e. Newton with Armijo's rule)α = 1.0
: initial value of α (damping) parameter for line search algorithmαmin = 0.001
: minimal value of the dampingalpha
For performance reasons, we decided to use an immutable structure to hold the parameters. One can use the package Setfield.jl
to drastically simplify the mutation of different fields. See the tutorials for examples.
BifurcationKit.ContinuationPar
— Typeoptions = ContinuationPar(dsmin = 1e-4,...)
Returns a variable containing parameters to affect the continuation
algorithm used to solve F(x,p) = 0
.
Arguments
dsmin, dsmax
are the minimum, maximum arclength allowed value. It controls the density of points in the computed branch of solutions.ds = 0.01
is the initial arclength.p_min, p_max
allowed parameter range forp
max_steps = 100
maximum number of continuation stepsnewton_options::NewtonPar
: options for the Newton algorithmsave_to_file = false
: save to file. A name is automatically generated or can be defined incontinuation
. This requiresusing JLD2
.save_sol_every_step::Int64 = 0
at which continuation steps do we save the current solutionplot_every_step = 10
at which continuation steps do we plot the current solution
Handling eigen elements, their computation is triggered by the argument detect_bifurcation
(see below)
nev = 3
number of eigenvalues to be computed. It is automatically increased to have at leastnev
unstable eigenvalues. To be set for proper bifurcation detection. See Detection of bifurcation points of Equilibria for more informations.save_eig_every_step = 1
record eigen vectors every specified steps. Important for memory limited resource, e.g. GPU.save_eigenvectors = true
Important for memory limited resource, e.g. GPU.
Handling bifurcation detection
tol_stability = 1e-10
lower bound on the real part of the eigenvalues to test for stability of equilibria and periodic orbitsdetect_fold = true
detect Fold bifurcations? It is a useful option although the detection of Fold is cheap. Indeed, it may happen that there is a lot of Fold points and this can saturate the memory in memory limited devices (e.g. on GPU)detect_bifurcation::Int
∈ {0, 1, 2, 3} If set to 0, nothing is done. If set to 1, the eigen-elements are computed. If set to 2, the bifurcations points are detected during the continuation run, but not located precisely. If set to 3, a bisection algorithm is used to locate the bifurcations points (slower). The possibility to switch off detection is a useful option. Indeed, it may happen that there are a lot of bifurcation points and this can saturate the memory of memory limited devices (e.g. on GPU)dsmin_bisection = 1e-16
dsmin for the bisection algorithm for locating bifurcation pointsn_inversion = 2
number of sign inversions in bisection algorithmmax_bisection_steps = 15
maximum number of bisection stepstol_bisection_eigenvalue = 1e-16
tolerance on real part of eigenvalue to detect bifurcation points in the bisection steps
Handling ds
adaptation (see continuation
for more information)
a = 0.5
aggressiveness factor. It is used to adaptds
in order to have a number of newton iterations per continuation step roughly constant. The highera
is, the larger the step sizeds
is changed at each continuation step.
Handling event detection
detect_event::Int
∈ {0, 1, 2} If set to 0, nothing is done. If set to 1, the event locations are sought during the continuation run, but not located precisely. If set to 2, a bisection algorithm is used to locate the event (slower).tol_param_bisection_event = 1e-16
tolerance on parameter to locate event
Misc
η = 150.
parameter to estimate tangent at first point with parameter p₀ + ds / ηdetect_loop
[WORK IN PROGRESS] detect loops in the branch and stop the continuation
For performance reasons, we decided to use an immutable structure to hold the parameters. One can use the package Setfield.jl
to drastically simplify the mutation of different fields. See tutorials for more examples.
Problems
DDEBifurcationKit.ConstantDDEBifProblem
— Typestruct ConstantDDEBifProblem{Tvf, Tdf, Tu, Td, Tp, Tl<:Lens, Tplot, Trec, Tδ} <: DDEBifurcationKit.AbstractDDEBifurcationProblem
Structure to hold the bifurcation problem. If don't have parameters, you can pass nothing
.
Fields
VF::Any
: Vector field, typically aBifFunction
delays::Any
: function delays. It takes the parameters and return the non-zero delays in anAsbtractVector
form. Example:delays = par -> [1.]
u0::Any
: Initial guessdelays0::Any
: initial delays (set internally by the constructor)params::Any
: parameterslens::Lens
: Typically aSetfield.Lens
. It specifies which parameter axis amongparams
is used for continuation. For example, ifpar = (α = 1.0, β = 1)
, we can perform continuation w.r.t.α
by usinglens = (@lens _.α)
. If you have an arraypar = [ 1.0, 2.0]
and want to perform continuation w.r.t. the first variable, you can uselens = (@lens _[1])
. For more information, we refer toSetField.jl
.plotSolution::Any
: user function to plot solutions during continuation. Signature:plotSolution(x, p; kwargs...)
recordFromSolution::Any
:record_from_solution = (x, p) -> norm(x)
function used record a few indicators about the solution. It could benorm
or(x, p) -> x[1]
. This is also useful when saving several huge vectors is not possible for memory reasons (for example on GPU...). This function can return pretty much everything but you should keep it small. For example, you can do(x, p) -> (x1 = x[1], x2 = x[2], nrm = norm(x))
or simply(x, p) -> (sum(x), 1)
. This will be stored incontres.branch
(see below). Finally, the first component is used to plot in the continuation curve.δ::Any
: delta for Finite differences
Methods
getu0(pb)
callspb.u0
getparams(pb)
callspb.params
getlens(pb)
callspb.lens
getParam(pb)
callsget(pb.params, pb.lens)
setParam(pb, p0)
callsset(pb.params, pb.lens, p0)
record_from_solution(pb)
callspb.record_from_solution
plotSolution(pb)
callspb.plotSolution
isSymmetric(pb)
callsisSymmetric(pb.prob)
Constructors
ConstantDDEBifProblem(F, delays, u0, params, lens; J, Jᵗ, d2F, d3F, kwargs...)
andkwargs
are the fields above.
DDEBifurcationKit.SDDDEBifProblem
— Typestruct SDDDEBifProblem{Tvf, Tdf, Tu, Td, Tp, Tl<:Lens, Tplot, Trec, Tδ} <: DDEBifurcationKit.AbstractDDEBifurcationProblem
Structure to hold the bifurcation problem. If don't have parameters, you can pass nothing
.
Fields
VF::Any
: Vector field, typically aBifFunction
delays::Any
: function delays. It takes the parameters and the state and return the non-zero delays in anAsbtractVector
form. Example:delays = (par, u) -> [1. + u[1]^2]
u0::Any
: Initial guessdelays0::Any
: initial delays (set internally by the constructor)params::Any
: parameterslens::Lens
: Typically aSetfield.Lens
. It specifies which parameter axis amongparams
is used for continuation. For example, ifpar = (α = 1.0, β = 1)
, we can perform continuation w.r.t.α
by usinglens = (@lens _.α)
. If you have an arraypar = [ 1.0, 2.0]
and want to perform continuation w.r.t. the first variable, you can uselens = (@lens _[1])
. For more information, we refer toSetField.jl
.plotSolution::Any
: user function to plot solutions during continuation. Signature:plotSolution(x, p; kwargs...)
recordFromSolution::Any
:record_from_solution = (x, p) -> norm(x)
function used record a few indicators about the solution. It could benorm
or(x, p) -> x[1]
. This is also useful when saving several huge vectors is not possible for memory reasons (for example on GPU...). This function can return pretty much everything but you should keep it small. For example, you can do(x, p) -> (x1 = x[1], x2 = x[2], nrm = norm(x))
or simply(x, p) -> (sum(x), 1)
. This will be stored incontres.branch
(see below). Finally, the first component is used to plot in the continuation curve.δ::Any
: delta for Finite differences
Methods
getu0(pb)
callspb.u0
getparams(pb)
callspb.params
getlens(pb)
callspb.lens
getParam(pb)
callsget(pb.params, pb.lens)
setParam(pb, p0)
callsset(pb.params, pb.lens, p0)
record_from_solution(pb)
callspb.record_from_solution
plotSolution(pb)
callspb.plotSolution
isSymmetric(pb)
callsisSymmetric(pb.prob)
Constructors
SDDDEBifProblem(F, delays, u0, params, lens; J, Jᵗ, d2F, d3F, kwargs...)
andkwargs
are the fields above.
Eigen solvers
DDEBifurcationKit.DDE_DefaultEig
— Typemutable struct DDE_DefaultEig{T, Tw, Tv} <: DDEBifurcationKit.AbstractDDEEigenSolver
Default eigen solver for DDEBifurcationKit based on the julia package NonlinearEigenproblems.jl. ore precisely, we rely on NonlinearEigenproblems.iar_chebyshev
for the computation of eigenvalues.
Fields
maxit::Int64
: Default: 100which::Any
: Default: realσ::Any
: Default: 0.0γ::Any
: Default: 1.0tol::Any
: Default: 1.0e-10logger::Int64
: Default: 0check_error_every::Int64
: Default: 1v::Any
: Default: nothing
Constructors
DDE_DefaultEig(; kwargs...)
andkwargs
are the fields above.
Branch switching (branch point)
Missing docstring for continuation(br::ContResult, ind_bif::Int, optionsCont::ContinuationPar ; kwargs...)
. Check Documenter's build log for details.
Branch switching (Hopf point)
BifurcationKit.continuation
— Methodcontinuation(
br,
ind_bif,
_contParams,
probPO;
alg,
δp,
ampfactor,
usedeflation,
nev,
kwargs...
)
Perform automatic branch switching from a Hopf bifurcation point labelled ind_bif
in the list of the bifurcated points of a previously computed branch br::ContResult
. It first computes a Hopf normal form.
Arguments
br
branch result from a call tocontinuation
ind_hopf
index of the bifurcation point inbr
contParams
parameters for the call tocontinuation
probPO
problem used to specify the way the periodic orbit is computed. It can bePeriodicOrbitTrapProblem
,ShootingProblem
orPoincareShootingProblem
.
Optional arguments
alg = br.alg
continuation algorithmδp
used to specify a particular guess for the parameter on the bifurcated branch which is otherwise determined bycontParams.ds
. This allows to use a step larger thancontParams.dsmax
.ampfactor = 1
factor which alter the amplitude of the bifurcated solution. Useful to magnify the bifurcated solution when the bifurcated branch is very steep.usedeflation = true
whether to use nonlinear deflation (see Deflated problems) to help finding the guess on the bifurcated branchnev
number of eigenvalues to be computed to get the right eigenvector- all
kwargs
fromcontinuation
A modified version of prob
is passed to plot_solution
and finalise_solution
.
You have to be careful about the options contParams.newton_options.linsolver
. In the case of Matrix-Free solver, you have to pass the right number of unknowns N * M + 1
. Note that the options for the preconditioner are not accessible yet.
Utils for periodic orbits
BifurcationKit.getperiod
— Functiongetperiod(, x)
getperiod(, x, par)
Compute the period of the periodic orbit associated to x
.
getperiod(prob, x, p)
Compute the period of the periodic orbit associated to x
.
getperiod(psh, x_bar, par)
Compute the period of the periodic orbit associated to x_bar
.
BifurcationKit.getamplitude
— Functiongetamplitude(prob, x, p; ratio)
Compute the amplitude of the periodic orbit associated to x
. The keyword argument ratio = 1
is used as follows. If length(x) = 1 + ratio * n
, the call returns the amplitude over x[1:n]
.
getamplitude(prob, x, p; ratio)
Compute the amplitude of the periodic orbit associated to x
. The keyword argument ratio = 1
is used as follows. If length(x) = ratio * n
, the call returns the amplitude over x[1:n]
.
BifurcationKit.getmaximum
— Functiongetmaximum(prob, x, p; ratio)
Compute the maximum of the periodic orbit associated to x
. The keyword argument ratio = 1
is used as follows. If length(x) = 1 + ratio * n
, the call returns the amplitude over x[1:n]
.
getmaximum(prob, x, p)
Compute the maximum of the periodic orbit associated to x
.
getmaximum(prob, x, p; ratio)
Compute the maximum of the periodic orbit associated to x
. The keyword argument ratio = 1
is used as follows. If length(x) = ratio * n
, the call returns the amplitude over x[1:n]
.
Misc.
BifurcationKit.guess_from_hopf
— Methodguess_from_hopf(
br,
ind_hopf,
eigsolver,
M,
amplitude;
phase
)
This function returns several useful quantities regarding a Hopf bifurcation point. More precisely, it returns:
- the parameter value at which a Hopf bifurcation occurs
- the period of the bifurcated periodic orbit
- a guess for the bifurcated periodic orbit
- the equilibrium at the Hopf bifurcation point
- the eigenvector at the Hopf bifurcation point.
The arguments are
br
: the continuation branch which lists the Hopf bifurcation pointsind_hopf
: index of the bifurcation branch, as inbr.specialpoint
eigsolver
: the eigen solver used to find the eigenvectorsM
number of time slices in the periodic orbit guessamplitude
: amplitude of the periodic orbit guess
BifurcationKit.get_normal_form
— Functionget_normal_form(
prob,
br,
id_bif;
nev,
verbose,
ζs,
lens,
Teigvec,
scaleζ,
detailed,
autodiff,
bls
)
Compute the normal form of the bifurcation point located at br.specialpoint[ind_bif]
.
Arguments
prob::AbstractBifurcationProblem
br
result from a call tocontinuation
ind_bif
index of the bifurcation point inbr.specialpoint
Optional arguments
nev
number of eigenvalues used to compute the spectral projection. This number has to be adjusted when used with iterative methods.verbose
whether to display informationζs
list of vectors spanning the kernel ofdF
at the bifurcation point. Useful to enforce the basis for the normal form.lens::Lens
specify which parameter to take the partial derivative ∂pFscaleζ
function to normalise the kernel basis. Indeed, when used with large vectors andnorm
, it results in ζs and the normal form coefficient being super small.autodiff = true
whether to use ForwardDiff for the differentiations w.r.t the parameters that are required to compute the normal form. Used for example for Bogdanov-Takens point. You can set toautodiff = false
if you wish.detailed = true
whether to compute only a simplified normal form. Used for example for Bogdanov-Takens point.bls = MatrixBLS()
specify Bordered linear solver. Used for example for Bogdanov-Takens point.
Available method
You can directly call
get_normal_form(br, ind_bif ; kwargs...)
which is a shortcut for get_normal_form(getProb(br), br, ind_bif ; kwargs...)
.
Once the normal form nf
has been computed, you can call predictor(nf, δp)
to obtain an estimate of the bifurcating branch.
get_normal_form(
prob,
br,
id_bif;
nev,
verbose,
ζs,
lens,
Teigvec,
scaleζ,
prm,
δ,
detailed
)
Compute the normal form of periodic orbits. Same arguments as the function getNormalForm
for equilibria. We detail the additional keyword arguments specific to periodic orbits
Optional arguments
prm = true
compute the normal form using Poincaré return map. For collocation, there will be another way to compute the normal form in the future.