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.
Results
BifurcationKit.NonLinearSolution
— TypeStructure which holds the solution from application of Newton-Krylov algorithm to a nonlinear problem.
For example
sol = newton(prob, NewtonPar())
Fields
u::Any
: solutionprob::Any
: nonlinear problem, typically aBifurcationProblem
residuals::Any
: sequence of residualsconverged::Bool
: has algorithm converged?itnewton::Int64
: number of newton stepsitlineartot::Any
: total number of linear iterations
methods
converged(sol)
return whether the solution has converged.
BifurcationKit.ContResult
— Typestruct ContResult{Tkind<:BifurcationKit.AbstractContinuationKind, Tbr, Teigvals, Teigvec, Biftype, Tsol, Tparc, Tprob, Talg} <: BifurcationKit.AbstractResult{Tkind<:BifurcationKit.AbstractContinuationKind, Tprob}
Structure which holds the results after a call to continuation
.
You can see the propertynames of a result br
by using propertynames(br)
or propertynames(br.branch)
.
Fields
branch::StructArrays.StructArray
: holds the low-dimensional information about the branch. More precisely,branch[i+1]
contains the following information(record_from_solution(u, param), param, itnewton, itlinear, ds, θ, n_unstable, n_imag, stable, step)
for each continuation stepi
.itnewton
number of Newton iterationsitlinear
total number of linear iterations during newton (corrector)n_unstable
number of eigenvalues with positive real part for each continuation step (to detect stationary bifurcation)n_imag
number of eigenvalues with positive real part and non zero imaginary part at current continuation step (useful to detect Hopf bifurcation).stable
stability of the computed solution for each continuation step. Hence,stable
should matcheig[step]
which corresponds tobranch[k]
for a givenk
.step
continuation step (here equali
)
eig::Array{@NamedTuple{eigenvals::Teigvals, eigenvecs::Teigvec, converged::Bool, step::Int64}, 1} where {Teigvals, Teigvec}
: A vector with eigen-elements at each continuation step.sol::Any
: Vector of solutions sampled along the branch. This is set by the argumentsave_sol_every_step::Int64
(default 0) inContinuationPar
.contparams::Any
: The parameters used for the call tocontinuation
which produced this branch. Must be aContinuationPar
kind::BifurcationKit.AbstractContinuationKind
: Type of solutions computed in this branch. Default: EquilibriumCont()prob::Any
: Bifurcation problem used to compute the branch, useful for branch switching. For example, when computing periodic orbits, the functionalPeriodicOrbitTrapProblem
,ShootingProblem
... will be saved here. Default: nothingspecialpoint::Vector
: A vector holding the set of detected bifurcation points. SeeSpecialPoint
for a list of special points.alg::Any
: Continuation algorithm used for the computation of the branch
Associated methods
length(br)
number of the continuation stepsshow(br)
display information about the brancheigenvals(br, ind)
returns the eigenvalues for the ind-th continuation stepeigenvec(br, ind, indev)
returns the indev-th eigenvector for the ind-th continuation stepget_normal_form(br, ind)
compute the normal form of the ind-th points inbr.specialpoint
getlens(br)
return the parameter axis used for the branchgetlenses(br)
return the parameter two axis used for the branch when 2 parameters continuation is used (Fold, Hopf, NS, PD)br[k+1]
gives information about the k-th step. A typical run yields the following
julia> br[1]
(x = 0.0, param = 0.1, itnewton = 0, itlinear = 0, ds = -0.01, θ = 0.5, n_unstable = 2, n_imag = 2, stable = false, step = 0, eigenvals = ComplexF64[0.1 - 1.0im, 0.1 + 1.0im], eigenvecs = ComplexF64[0.7071067811865475 - 0.0im 0.7071067811865475 + 0.0im; 0.0 + 0.7071067811865475im 0.0 - 0.7071067811865475im])
which provides the value param
of the parameter of the current point, its stability, information on the newton iterations, etc. The fields can be retrieved using propertynames(br.branch)
. This information is stored in br.branch
which is a StructArray
. You can thus extract the vector of parameters along the branch as
julia> br.param
10-element Vector{Float64}:
0.1
0.08585786437626905
0.06464466094067263
0.03282485578727799
-1.2623798512809007e-5
-0.07160718539365075
-0.17899902778635765
-0.3204203840236672
-0.4618417402609767
-0.5
get_solx(br, k)
returns the k-th solution on the branchget_solp(br, k)
returns the parameter value associated with k-th solution on the branchgetparams(br)
Parameters passed to continuation and used in the equationF(x, par) = 0
.setparam(br, p0)
set the parameter valuep0
according to::Lens
for the parameters of the problembr.prob
getlens(br)
get the lens used for the computation of the branchcontinuation(br, ind)
performs automatic branch switching (aBS) from ind-th bifurcation point. Typically branching from equilibrium to equilibrium, or periodic orbit to periodic orbit.continuation(br, ind, lens2)
performs two parameters(getLens(br), lens2)
continuation of the ind-th bifurcation point.continuation(br, ind, probPO::AbstractPeriodicOrbitProblem)
performs aBS from ind-th bifurcation point (which must be a Hopf bifurcation point) to branch of periodic orbits.
Problems
BifurcationKit.BifFunction
— Typestruct BifFunction{Tf, Tdf, Tdfad, Tj, Tjad, Td2f, Td2fc, Td3f, Td3fc, Tsym, Tδ} <: BifurcationKit.AbstractBifurcationFunction
Structure to hold the vector field and its derivatives. It should rarely be called directly. Also, in essence, it is very close to SciMLBase.ODEFunction
.
Fields
F::Any
: Vector field. Function of type out-of-placeresult = f(x, p)
or inplacef(result, x, p)
. For type stability, the types ofx
andresult
should matchdF::Any
: Differential ofF
with respect tox
, signaturedF(x,p,dx)
dFad::Any
: Adjoint of the Differential ofF
with respect tox
, signaturedFad(x,p,dx)
J::Any
: Jacobian ofF
at(x, p)
. It can assume three forms. 1. EitherJ
is a function andJ(x, p)
returns a::AbstractMatrix
. In this case, the default arguments ofcontparams::ContinuationPar
will makecontinuation
work. 2. OrJ
is a function andJ(x, p)
returns a function taking one argumentdx
and returningdr
of the same type asdx
. In our notation,dr = J * dx
. In this case, the default parameters ofcontparams::ContinuationPar
will not work and you have to use a Matrix Free linear solver, for exampleGMRESIterativeSolvers
, 3. OrJ
is a function andJ(x, p)
returns a variablej
which can assume any type. Then, you must implement a linear solverls
as a composite type, subtype ofAbstractLinearSolver
which is called likels(j, rhs)
and which returns the solution of the jacobian linear system. See for exampleexamples/SH2d-fronts-cuda.jl
. This linear solver is passed toNewtonPar(linsolver = ls)
which itself passed toContinuationPar
. Similarly, you have to implement an eigensolvereig
as a composite type, subtype ofAbstractEigenSolver
.Jᵗ::Any
: jacobian adjoint, it should be implemented in an efficient manner. For matrix-free methods,transpose
is not readily available and the user must provide a dedicated method. In the case of sparse based jacobian,Jᵗ
should not be passed as it is computed internally more efficiently, i.e. it avoids recomputing the jacobian as it would be if you passJᵗ = (x, p) -> transpose(dF(x, p))
.d2F::Any
: Second Differential ofF
with respect tox
, signatured2F(x,p,dx1,dx2)
d3F::Any
: Third Differential ofF
with respect tox
, signatured3F(x,p,dx1,dx2,dx3)
d2Fc::Any
: [internal] Second Differential ofF
with respect tox
which accept complex vectors dxid3Fc::Any
: [internal] Third Differential ofF
with respect tox
which accept complex vectors dxiisSymmetric::Any
: Whether the jacobian is auto-adjoint.δ::Any
: used internally to compute derivatives (with finite differences), for example for normal form computation and codim 2 continuation.inplace::Bool
: optionally sets whether the function is inplace or not
Methods
residual(pb::BifFunction, x, p)
callspb.F(x,p)
jacobian(pb::BifFunction, x, p)
callspb.J(x, p)
dF(pb::BifFunction, x, p, dx)
callspb.dF(x,p,dx)
- etc
BifurcationKit.BifurcationProblem
— Typestruct BifurcationProblem{Tvf, Tu, Tp, Tl<:Lens, Tplot, Trec, Tgets} <: BifurcationKit.AbstractAllJetBifProblem
Structure to hold the bifurcation problem.
Fields
VF::Any
: Vector field, typically aBifFunction
u0::Any
: Initial guessparams::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:plot_solution(x, p; kwargs...)
for Plot.jl andplot_solution(ax, x, p; kwargs...)
for the Makie package(s).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
wherecontres::ContResult
is the continuation curve of the bifurcation problem. Finally, the first component is used for plotting in the continuation curve.save_solution::Any
: function to save the full solution on the branch. Some problem are mutable (like periodic orbit functional with adaptive mesh) and this function allows to save the state of the problem along with the solution itself. Signaturesave_solution(x, p)
Methods
re_make(pb; kwargs...)
modify a bifurcation problemgetu0(pb)
callspb.u0
getparams(pb)
callspb.params
getlens(pb)
callspb.lens
getparam(pb)
callsget(pb.params, pb.lens)
setparam(pb, p0)
calls_set_param(pb.params, pb.lens, p0)
record_from_solution(pb)
callspb.recordFromSolution
plot_solution(pb)
callspb.plotSolution
is_symmetric(pb)
callsis_symmetric(pb.prob)
Constructors
BifurcationProblem(F, u0, params, lens)
all derivatives are computed using ForwardDiff.BifurcationProblem(F, u0, params, lens; J, Jᵗ, d2F, d3F, kwargs...)
andkwargs
are the fields above. You can pass your own jacobian withJ
(seeBifFunction
for description of the jacobian function) and jacobian adjoint withJᵗ
. For example, this can be used to provide finite differences based jacobian usingBifurcationKit.finiteDifferences
.
BifurcationKit.DeflationOperator
— Typestruct DeflationOperator{Tp<:Real, Tdot, T<:Real, vectype} <: BifurcationKit.AbstractDeflationFactor
Structure for defining a custom distance.
This operator allows to handle the following situation. Assume you want to solve F(x)=0
with a Newton algorithm but you want to avoid the process to return some already known solutions $roots_i$. The deflation operator penalizes these roots. You can create a DeflationOperator
to define a scalar function M(u)
used to find, with Newton iterations, the zeros of the following function $F(u) \cdot Π_i(\|u - root_i\|^{-2p} + \alpha) := F(u) \cdot M(u)$ where $\|u\|^2 = dot(u, u)$. The fields of the struct DeflationOperator
are as follows:
power::Real
: powerp
. You can use anInt
for exampledot::Any
: function, this function has to be bilinear and symmetric for the linear solver to work wellα::Real
: shiftroots::Vector
: rootstmp::Any
autodiff::Bool
δ::Real
Given defOp::DeflationOperator
, one can access its roots via defOp[n]
as a shortcut for defOp.roots[n]
. Note that you can also use defOp[end]
.
Also, one can add (resp. remove) a new root by using push!(defOp, newroot)
(resp. pop!(defOp)
). Finally length(defOp)
is a shortcut for length(defOp.roots)
Constructors
DeflationOperator(p::Real, α::Real, roots::Vector{vectype}; autodiff = false)
DeflationOperator(p::Real, dt, α::Real, roots::Vector{vectype}; autodiff = false)
DeflationOperator(p::Real, α::Real, roots::Vector{vectype}, v::vectype; autodiff = false)
The option autodiff
triggers the use of automatic differentiation for the computation of the gradient of the scalar function M
. This works only on AbstractVector
for now.
Custom distance
You are asked to pass a scalar product like dot
to build a DeflationOperator
. However, in some cases, you may want to pass a custom distance dist(u, v)
. You can do this using
`DeflationOperator(p, CustomDist(dist), α, roots)`
Note that passing CustomDist(dist, true)
will trigger the use of automatic differentiation for the gradient of M
.
Linear solvers
When used with newton, you have access to the following linear solvers
- custom solver
DeflatedProblemCustomLS()
which requires solving two linear systemsJ⋅x = rhs
. - For other linear solvers
<: AbstractLinearSolver
, a matrix free method is used for the deflated functional. - if passed
Val(:autodiff)
, thenForwardDiff.jl
is used to compute the jacobian Matrix of the deflated problem - if passed
Val(:fullIterative)
, then a full matrix free method is used for the deflated problem.
BifurcationKit.DeflatedProblem
— Typepb = DeflatedProblem(prob, M::DeflationOperator, jactype)
Create a DeflatedProblem
.
This creates a deflated functional (problem) $M(u) \cdot F(u) = 0$ where M
is a DeflationOperator
which encodes the penalization term. prob
is an AbstractBifurcationProblem
which encodes the functional. It is not meant not be used directly albeit by advanced users.
Periodic orbits
BifurcationKit.PeriodicOrbitTrapProblem
— TypeThis composite type implements Finite Differences based on a Trapezoidal rule (Order 2 in time) to locate periodic orbits. More details (maths, notations, linear systems) can be found here.
Fields
prob
a bifurcation problemM::Int
number of time slicesϕ
used to set a section for the phase constraint equation, of size N*Mxπ
used in the section for the phase constraint equation, of size N*Mlinsolver: = DefaultLS()
linear solver for each time slice, i.e. to solveJ⋅sol = rhs
. This is only needed for the computation of the Floquet multipliers in a full matrix-free setting.ongpu::Bool
whether the computation takes place on the gpu (Experimental)massmatrix
a mass matrix. You can pass for example a sparse matrix. Default: identity matrix.update_section_every_step
updates the section everyupdate_section_every_step
step during continuationjacobian::Symbol
symbol which describes the type of jacobian used in Newton iterations (see below).
The scheme is as follows. We first consider a partition of $[0,1]$ given by $0<s_0<\cdots<s_m=1$ and one looks for T = x[end]
such that
$M_a\cdot\left(x_{i} - x_{i-1}\right) - \frac{T\cdot h_i}{2} \left(F(x_{i}) + F(x_{i-1})\right) = 0,\ i=1,\cdots,m-1$
with $u_{0} := u_{m-1}$ and the periodicity condition $u_{m} - u_{1} = 0$ and
where $h_1 = s_i-s_{i-1}$. $M_a$ is a mass matrix. Finally, the phase of the periodic orbit is constrained by using a section (but you could use your own)
$\sum_i\langle x_{i} - x_{\pi,i}, \phi_{i}\rangle=0.$
Constructors
The structure can be created by calling PeriodicOrbitTrapProblem(;kwargs...)
. For example, you can declare such a problem without vector field by doing
PeriodicOrbitTrapProblem(M = 100)
Orbit guess
You will see below that you can evaluate the residual of the functional (and other things) by calling pb(orbitguess, p)
on an orbit guess orbitguess
. Note that orbitguess
must be a vector of size M * N + 1 where N is the number of unknowns in the state space and orbitguess[M*N+1]
is an estimate of the period $T$ of the limit cycle. More precisely, using the above notations, orbitguess
must be $orbitguess = [x_{1},x_{2},\cdots,x_{M}, T]$.
Note that you can generate this guess from a function solution using generateSolution
.
Functional
A functional, hereby called G
, encodes this problem. The following methods are available
pb(orbitguess, p)
evaluates the functional G onorbitguess
pb(orbitguess, p, du)
evaluates the jacobiandG(orbitguess).du
functional atorbitguess
ondu
pb(Val(:JacFullSparse), orbitguess, p)
return the sparse matrix of the jacobiandG(orbitguess)
atorbitguess
without the constraints. It is calledA_γ
in the docs.pb(Val(:JacFullSparseInplace), J, orbitguess, p)
. Same aspb(Val(:JacFullSparse), orbitguess, p)
but overwritesJ
inplace. Note that the sparsity pattern must be the same independently of the values of the parameters or oforbitguess
. In this case, this is significantly faster thanpb(Val(:JacFullSparse), orbitguess, p)
.pb(Val(:JacCyclicSparse), orbitguess, p)
return the sparse cyclic matrix Jc (see the docs) of the jacobiandG(orbitguess)
atorbitguess
pb(Val(:BlockDiagSparse), orbitguess, p)
return the diagonal of the sparse matrix of the jacobiandG(orbitguess)
atorbitguess
. This allows to design Jacobi preconditioner. Useblockdiag
.
Jacobian
Specify the choice of the jacobian (and linear algorithm), jacobian
must belong to [:FullLU, :FullSparseInplace, :Dense, :DenseAD, :BorderedLU, :BorderedSparseInplace, :FullMatrixFree, :BorderedMatrixFree, :FullMatrixFreeAD]
. This is used to select a way of inverting the jacobian dG
of the functional G.
- For
jacobian = :FullLU
, we use the default linear solver based on a sparse matrix representation ofdG
. This matrix is assembled at each newton iteration. This is the default algorithm. - For
jacobian = :FullSparseInplace
, this is the same as for:FullLU
but the sparse matrixdG
is updated inplace. This method allocates much less. In some cases, this is significantly faster than using:FullLU
. Note that this method can only be used if the sparsity pattern of the jacobian is always the same. - For
jacobian = :Dense
, same as above but the matrixdG
is dense. It is also updated inplace. This option is useful to study ODE of small dimension. - For
jacobian = :DenseAD
, evaluate the jacobian using ForwardDiff - For
jacobian = :BorderedLU
, we take advantage of the bordered shape of the linear solver and use a LU decomposition to invertdG
using a bordered linear solver. - For
jacobian = :BorderedSparseInplace
, this is the same as for:BorderedLU
but the cyclic matrixdG
is updated inplace. This method allocates much less. In some cases, this is significantly faster than using:BorderedLU
. Note that this method can only be used if the sparsity pattern of the jacobian is always the same. - For
jacobian = :FullMatrixFree
, a matrix free linear solver is used fordG
: note that a preconditioner is very likely required here because of the cyclic shape ofdG
which affects negatively the convergence properties of GMRES. - For
jacobian = :BorderedMatrixFree
, a matrix free linear solver is used but forJc
only (see docs): it means thatoptions.linsolver
is used to invertJc
. These two Matrix-Free options thus expose different part of the jacobiandG
in order to use specific preconditioners. For example, an ILU preconditioner onJc
could remove the constraints indG
and lead to poor convergence. Of course, for these last two methods, a preconditioner is likely to be required. - For
jacobian = :FullMatrixFreeAD
, the evaluation map of the differential is derived using automatic differentiation. Thus, unlike the previous two cases, the user does not need to pass a Matrix-Free differential.
For these methods to work on the GPU, for example with CuArrays
in mode allowscalar(false)
, we face the issue that the function extract_period_fdtrap
won't be well defined because it is a scalar operation. Note that you must pass the option ongpu = true
for the functional to be evaluated efficiently on the gpu.
BifurcationKit.PeriodicOrbitOCollProblem
— Typepb = PeriodicOrbitOCollProblem(kwargs...)
This composite type implements an orthogonal collocation (at Gauss points) method of piecewise polynomials to locate periodic orbits. More details (maths, notations, linear systems) can be found here.
Arguments
prob
a bifurcation problemϕ::AbstractVector
used to set a section for the phase constraint equationxπ::AbstractVector
used in the section for the phase constraint equationN::Int
dimension of the state spacemesh_cache::MeshCollocationCache
cache for collocation. See docs ofMeshCollocationCache
update_section_every_step
updates the section everyupdate_section_every_step
step during continuationjacobian = DenseAnalytical()
describes the type of jacobian used in Newton iterations. Can only beAutoDiffDense(), DenseAnalytical(), FullSparse(), FullSparseInplace()
.meshadapt::Bool = false
whether to use mesh adaptationverbose_mesh_adapt::Bool = true
verbose mesh adaptation informationK::Float64 = 500
parameter for mesh adaptation, control new mesh step size. More precisely, we set max(hᵢ) / min(hᵢ) ≤ K if hᵢ denotes the time steps.
Methods
Here are some useful methods you can apply to pb
length(pb)
gives the total number of unknownssize(pb)
returns the triplet(N, m, Ntst)
getmesh(pb)
returns the mesh0 = τ0 < ... < τNtst+1 = 1
. This is useful because this mesh is born to vary during automatic mesh adaptationget_mesh_coll(pb)
returns the (static) mesh0 = σ0 < ... < σm+1 = 1
get_times(pb)
returns the vector of times (length1 + m * Ntst
) at the which the collocation is applied.generate_solution(pb, orbit, period)
generate a guess from a functiont -> orbit(t)
which approximates the periodic orbit.POSolution(pb, x)
return a function interpolating the solutionx
using a piecewise polynomials function
Orbit guess
You can evaluate the residual of the functional (and other things) by calling pb(orbitguess, p)
on an orbit guess orbitguess
. Note that orbitguess
must be of size 1 + N * (1 + m * Ntst) where N is the number of unknowns in the state space and orbitguess[end]
is an estimate of the period $T$ of the limit cycle.
Constructors
PeriodicOrbitOCollProblem(Ntst::Int, m::Int; kwargs)
creates an empty functional withNtst
andm
.
Note that you can generate this guess from a function using generate_solution
.
Functional
A functional, hereby called G
, encodes this problem. The following methods are available
pb(orbitguess, p)
evaluates the functional G onorbitguess
BifurcationKit.ShootingProblem
— Typepb = ShootingProblem(flow::Flow, ds, section; parallel = false)
Create a problem to implement the Simple / Parallel Multiple Standard Shooting method to locate periodic orbits. More details (maths, notations, linear systems) can be found here. The arguments are as follows
flow::Flow
: implements the flow of the Cauchy problem though the structureFlow
.ds
: vector of time differences for each shooting. Its length is writtenM
. IfM == 1
, then the simple shooting is implemented and the multiple one otherwise.section
: implements a phase condition. The evaluationsection(x, T)
must return a scalar number wherex
is a guess for one point on the periodic orbit andT
is the period of the guess. Also, the methodsection(x, T, dx, dT)
must be available and which returns the differential ofsection
. The type ofx
depends on what is passed to the newton solver. SeeSectionSS
for a type of section defined as a hyperplane.parallel
whether the shooting is computed in parallel (threading). Available through the use of Flows defined byEnsembleProblem
(this is automatically set up for you).par
parameters of the modellens
parameter axisupdate_section_every_step
updates the section everyupdate_section_every_step
step during continuationjacobian::Symbol
symbol which describes the type of jacobian used in Newton iterations (see below).
A functional, hereby called G
, encodes the shooting problem. For example, the following methods are available:
pb(orbitguess, par)
evaluates the functional G onorbitguess
pb(orbitguess, par, du)
evaluates the jacobiandG(orbitguess)⋅du
functional atorbitguess
ondu
.pb
(Val(:JacobianMatrixInplace), J, x, par)` compute the jacobian of the functional analytically. This is based on ForwardDiff.jl. Useful mainly for ODEs.pb(Val(:JacobianMatrix), x, par)
same as above but out-of-place.
You can then call pb(orbitguess, par)
to apply the functional to a guess. Note that orbitguess::AbstractVector
must be of size M * N + 1
where N is the number of unknowns of the state space and orbitguess[M * N + 1]
is an estimate of the period T
of the limit cycle. This form of guess is convenient for the use of the linear solvers in IterativeSolvers.jl
(for example) which only accept AbstractVector
s. Another accepted guess is of the form BorderedArray(guess, T)
where guess[i]
is the state of the orbit at the i
th time slice. This last form allows for non-vector state space which can be convenient for 2d problems for example, use GMRESKrylovKit
for the linear solver in this case.
Note that you can generate this guess from a function solution using generate_solution
.
Jacobian
jacobian
Specify the choice of the linear algorithm, which must belong to[AutoDiffMF(), MatrixFree(), AutodiffDense(), AutoDiffDenseAnalytical(), FiniteDifferences(), FiniteDifferencesMF()]
. This is used to select a way of inverting the jacobian dG- For
MatrixFree()
, matrix free jacobian, the jacobian is specified by the user inprob
. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system - For
AutoDiffMF()
, we use Automatic Differentiation (AD) to compute the (matrix-free) derivative ofx -> prob(x, p)
using a directional derivative. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system - For
AutodiffDense()
. Same as forAutoDiffMF
but the jacobian is formed as a dense Matrix. You can use a direct solver or an iterative one. - For
FiniteDifferences()
, same as forAutoDiffDense
but we use Finite Differences to compute the jacobian ofx -> prob(x, p)
using theδ = 1e-8
which can be passed as an argument. - For
AutoDiffDenseAnalytical()
. Same as forAutoDiffDense
but the jacobian is formed using a mix of AD and analytical formula. - For
FiniteDifferencesMF()
, use Finite Differences to compute the matrix-free jacobian ofx -> prob(x, p)
using theδ = 1e-8
which can be passed as an argument.
- For
Simplified constructors
- The first important constructor is the following which is used for branching to periodic orbits from Hopf bifurcation points:
pb = ShootingProblem(M::Int, prob::Union{ODEProblem, EnsembleProblem}, alg; kwargs...)
- A convenient way to build the functional is to use:
pb = ShootingProblem(prob::Union{ODEProblem, EnsembleProblem}, alg, centers::AbstractVector; kwargs...)
where prob
is an ODEProblem
(resp. EnsembleProblem
) which is used to create a flow using the ODE solver alg
(for example Tsit5()
). centers
is list of M
points close to the periodic orbit, they will be used to build a constraint for the phase. parallel = false
is an option to use Parallel simulations (Threading) to simulate the multiple trajectories in the case of multiple shooting. This is efficient when the trajectories are relatively long to compute. Finally, the arguments kwargs
are passed to the ODE solver defining the flow. Look at DifferentialEquations.jl
for more information. Note that, in this case, the derivative of the flow is computed internally using Finite Differences.
- Another way to create a Shooting problem with more options is the following where in particular, one can provide its own scalar constraint
section(x)::Number
for the phase:
pb = ShootingProblem(prob::Union{ODEProblem, EnsembleProblem}, alg, M::Int, section; parallel = false, kwargs...)
or
pb = ShootingProblem(prob::Union{ODEProblem, EnsembleProblem}, alg, ds, section; parallel = false, kwargs...)
- The next way is an elaboration of the previous one
pb = ShootingProblem(prob1::Union{ODEProblem, EnsembleProblem}, alg1, prob2::Union{ODEProblem, EnsembleProblem}, alg2, M::Int, section; parallel = false, kwargs...)
or
pb = ShootingProblem(prob1::Union{ODEProblem, EnsembleProblem}, alg1, prob2::Union{ODEProblem, EnsembleProblem}, alg2, ds, section; parallel = false, kwargs...)
where we supply now two ODEProblem
s. The first one prob1
, is used to define the flow associated to F
while the second one is a problem associated to the derivative of the flow. Hence, prob2
must implement the following vector field $\tilde F(x,y,p) = (F(x,p), dF(x,p)\cdot y)$.
BifurcationKit.PoincareShootingProblem
— Typepb = PoincareShootingProblem(flow::Flow, M, sections; δ = 1e-8, interp_points = 50, parallel = false)
This composite type implements the Poincaré Shooting method to locate periodic orbits by relying on Poincaré return maps. More details (maths, notations, linear systems) can be found here. The arguments are as follows
flow::Flow
: implements the flow of the Cauchy problem though the structureFlow
.M
: the number of Poincaré sections. IfM == 1
, then the simple shooting is implemented and the multiple one otherwise.sections
: function or callable struct which implements a Poincaré section condition. The evaluationsections(x)
must return a scalar number whenM == 1
. Otherwise, one must implement a functionsection(out, x)
which populatesout
with theM
sections. SeeSectionPS
for type of section defined as a hyperplane.δ = 1e-8
used to compute the jacobian of the functional by finite differences. If set to0
, an analytical expression of the jacobian is used instead.interp_points = 50
number of interpolation point used to define the callback (to compute the hitting of the hyperplane section)parallel = false
whether the shooting are computed in parallel (threading). Only available through the use of Flows defined byEnsembleProblem
.par
parameters of the modellens
parameter axisupdate_section_every_step
updates the section everyupdate_section_every_step
step during continuationjacobian::Symbol
symbol which describes the type of jacobian used in Newton iterations (see below).
Jacobian
jacobian
Specify the choice of the linear algorithm, which must belong to[AutoDiffMF(), MatrixFree(), AutodiffDense(), AutoDiffDenseAnalytical(), FiniteDifferences(), FiniteDifferencesMF()]
. This is used to select a way of inverting the jacobian dG- For
MatrixFree()
, matrix free jacobian, the jacobian is specified by the user inprob
. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system - For
AutoDiffMF()
, we use Automatic Differentiation (AD) to compute the (matrix-free) derivative ofx -> prob(x, p)
using a directional derivative. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system - For
AutodiffDense()
. Same as forAutoDiffMF
but the jacobian is formed as a dense Matrix. You can use a direct solver or an iterative one. - For
FiniteDifferences()
, same as forAutoDiffDense
but we use Finite Differences to compute the jacobian ofx -> prob(x, p)
using theδ = 1e-8
which can be passed as an argument. - For
AutoDiffDenseAnalytical()
. Same as forAutoDiffDense
but the jacobian is formed using a mix of AD and analytical formula. - For
FiniteDifferencesMF()
, use Finite Differences to compute the matrix-free jacobian ofx -> prob(x, p)
using theδ = 1e-8
which can be passed as an argument.
- For
Simplified constructors
The first important constructor is the following which is used for branching to periodic orbits from Hopf bifurcation points pb = PoincareShootingProblem(M::Int, prob::Union{ODEProblem, EnsembleProblem}, alg; kwargs...)
A convenient way is to create a functional is
pb = PoincareShootingProblem(prob::ODEProblem, alg, section; kwargs...)
for simple shooting or
pb = PoincareShootingProblem(prob::Union{ODEProblem, EnsembleProblem}, alg, M::Int, section; kwargs...)
for multiple shooting . Here prob
is an Union{ODEProblem, EnsembleProblem}
which is used to create a flow using the ODE solver alg
(for example Tsit5()
). Finally, the arguments kwargs
are passed to the ODE solver defining the flow. We refer to DifferentialEquations.jl
for more information.
- Another convenient call is
pb = PoincareShootingProblem(prob::Union{ODEProblem, EnsembleProblem}, alg, normals::AbstractVector, centers::AbstractVector; δ = 1e-8, kwargs...)
where normals
(resp. centers
) is a list of normals (resp. centers) which defines a list of hyperplanes $\Sigma_i$. These hyperplanes are used to define partial Poincaré return maps.
Computing the functionals
A functional, hereby called G
encodes this shooting problem. You can then call pb(orbitguess, par)
to apply the functional to a guess. Note that orbitguess::AbstractVector
must be of size M * N where N is the number of unknowns in the state space and M
is the number of Poincaré maps. Another accepted guess
is such that guess[i]
is the state of the orbit on the i
th section. This last form allows for non-vector state space which can be convenient for 2d problems for example.
Note that you can generate this guess from a function solution using generate_solution
.
pb(orbitguess, par)
evaluates the functional G onorbitguess
pb(orbitguess, par, du)
evaluates the jacobiandG(orbitguess).du
functional atorbitguess
ondu
pb
(Val(:JacobianMatrixInplace), J, x, par)` compute the jacobian of the functional analytically. This is based on ForwardDiff.jl. Useful mainly for ODEs.pb(Val(:JacobianMatrix), x, par)
same as above but out-of-place.
You can use the function getperiod(pb, sol, par)
to get the period of the solution sol
for the problem with parameters par
.
Waves
BifurcationKit.TWProblem
— TypeTWProblem(prob, ∂::Tuple, u₀; DAE = 0, jacobian::Symbol = :AutoDiff)
This composite type implements a functional for freezing symmetries in order, for example, to compute travelling waves (TW). Note that you can freeze many symmetries, not just one, by passing many Lie generators. When you call pb(x, par)
, it computes:
┌ ┐
│ f(x, par) - s⋅∂⋅x │
│ <x - u₀, ∂⋅u₀> │
└ ┘
Arguments
prob
bifurcation problem with continuous symmetries∂::Tuple = (T1, T2, ⋯)
tuple of Lie generators. In effect, each of these is an (differential) operator which can be specified as a (sparse) matrix or as an operator implementingLinearAlgebra.mul!
.u₀
reference solution
Additional Constructor(s)
pb = TWProblem(prob, ∂, u₀; kw...)
This simplified call handles the case where a single symmetry needs to be frozen.
Useful function
updatesection!(pb::TWProblem, u0)
updates the reference solution of the problem usingu0
.nb_constraints(::TWProblem)
number of constraints (or Lie generators)
Newton
BifurcationKit.newton
— Function newton(prob::AbstractBifurcationProblem, options::NewtonPar; normN = norm, callback = (;x, fx, J, residual, step, itlinear, options, x0, residuals; kwargs...) -> true, kwargs...)
This is the Newton-Krylov Solver for F(x, p0) = 0
with Jacobian w.r.t. x
written J(x, p0)
and initial guess x0
. It is important to set the linear solver options.linsolver
properly depending on your problem. This linear solver is used to solve $J(x, p_0)u = -F(x, p_0)$ in the Newton step. You can for example use linsolver = DefaultLS()
which is the operator backslash: it works well for Sparse / Dense matrices. See Linear solvers (LS) for more informations.
Arguments
prob
a::AbstractBifurcationProblem
, typically aBifurcationProblem
which holds the vector field and its jacobian. We also refer toBifFunction
for more details.options::NewtonPar
variable holding the internal parameters used by thenewton
method
Optional Arguments
normN = norm
specifies a norm for the convergence criteriacallback
function passed by the user which is called at the end of each iteration. The default one is the followingcb_default((x, fx, J, residual, step, itlinear, options, x0, residuals); k...) = true
. Can be used to update a preconditionner for example. You can use for examplecbMaxNorm
to limit the residuals norms. If yo want to specify your own, the arguments passed to the callback are as followsx
current solutionfx
current residualJ
current jacobianresidual
current norm of the residualstep
current newton stepitlinear
number of iterations to solve the linear systemoptions
a copy of the argumentoptions
passed tonewton
residuals
the history of residualskwargs
kwargs arguments, contain your initial guessx0
kwargs
arguments passed to the callback. Useful whennewton
is called fromcontinuation
Output:
solution::NonLinearSolution
, we refer toNonLinearSolution
for more information.
Make sure that the linear solver (Matrix-Free...) corresponds to your jacobian (Matrix-Free vs. Matrix based).
newton(prob, defOp, options; ...)
newton(prob, defOp, options, _linsolver; kwargs...)
This is the deflated version of the Krylov-Newton Solver for F(x, p0) = 0
.
We refer to the regular newton
for more information. It penalises the roots saved in defOp.roots
. The other arguments are as for newton
. See DeflationOperator
for more information on defOp
.
Arguments
Compared to newton
, the only different arguments are
defOp::DeflationOperator
deflation operatorlinsolver
linear solver used to invert the Jacobian of the deflated functional.- custom solver
DeflatedProblemCustomLS()
which requires solving two linear systemsJ⋅x = rhs
. - For other linear solvers
<: AbstractLinearSolver
, a matrix free method is used for the deflated functional. - if passed
Val(:autodiff)
, thenForwardDiff.jl
is used to compute the jacobian Matrix of the deflated problem - if passed
Val(:fullIterative)
, then a full matrix free method is used for the deflated problem.
- custom solver
This specific Newton-Krylov method first tries to converge to a solution sol0
close the guess x0
. It then attempts to converge from the guess x1
while avoiding the previous converged solution close to sol0
. This is very handy for branch switching. The method is based on a deflated Newton-Krylov solver.
Arguments
Compared to newton
, the only different arguments are
defOp::DeflationOperator
deflation operatorlinsolver
linear solver used to invert the Jacobian of the deflated functional.- custom solver
DeflatedProblemCustomLS()
which requires solving two linear systemsJ⋅x = rhs
. - For other linear solvers
<: AbstractLinearSolver
, a matrix free method is used for the deflated functional. - if passed
Val(:autodiff)
, thenForwardDiff.jl
is used to compute the jacobian Matrix of the deflated problem - if passed
Val(:fullIterative)
, then a full matrix free method is used for the deflated problem.
- custom solver
newton(
br,
ind_bif;
normN,
options,
start_with_eigen,
lens2,
kwargs...
)
This function turns an initial guess for a Fold / Hopf point into a solution to the Fold / Hopf problem based on a Minimally Augmented formulation.
Arguments
br
results returned after a call to continuationind_bif
bifurcation index inbr
Optional arguments:
options::NewtonPar
, default valuebr.contparams.newton_options
normN = norm
options
You can pass newton parameters different from the ones stored inbr
by using this argumentoptions
.bdlinsolver
bordered linear solver for the constraint equationstart_with_eigen = false
whether to start the Minimally Augmented problem with information from eigen elements.kwargs
keywords arguments to be passed to the regular Newton-Krylov solver
For ODE problems, it is more efficient to use the Matrix based Bordered Linear Solver passing the option bdlinsolver = MatrixBLS()
It is recommended that you use the option start_with_eigen=true
newton(prob, orbitguess, options; lens, δ, kwargs...)
This is the Newton-Krylov Solver for computing a periodic orbit using the (Standard / Poincaré) Shooting method. Note that the linear solver has to be appropriately set up in options
.
Arguments
Similar to newton
except that prob
is either a ShootingProblem
or a PoincareShootingProblem
. These two problems have specific options to be tuned, we refer to their link for more information and to the tutorials.
prob
a problem of type<: AbstractShootingProblem
encoding the shooting functional G.orbitguess
a guess for the periodic orbit. SeeShootingProblem
and SeePoincareShootingProblem
for information regarding the shape oforbitguess
.par
parameters to be passed to the functionaloptions
same as for the regularnewton
method.
Optional argument
jacobian
Specify the choice of the linear algorithm, which must belong to[AutoDiffMF(), MatrixFree(), AutodiffDense(), AutoDiffDenseAnalytical(), FiniteDifferences(), FiniteDifferencesMF()]
. This is used to select a way of inverting the jacobian dG- For
MatrixFree()
, matrix free jacobian, the jacobian is specified by the user inprob
. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system - For
AutoDiffMF()
, we use Automatic Differentiation (AD) to compute the (matrix-free) derivative ofx -> prob(x, p)
using a directional derivative. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system - For
AutodiffDense()
. Same as forAutoDiffMF
but the jacobian is formed as a dense Matrix. You can use a direct solver or an iterative one. - For
FiniteDifferences()
, same as forAutoDiffDense
but we use Finite Differences to compute the jacobian ofx -> prob(x, p)
using theδ = 1e-8
which can be passed as an argument. - For
AutoDiffDenseAnalytical()
. Same as forAutoDiffDense
but the jacobian is formed using a mix of AD and analytical formula. - For
FiniteDifferencesMF()
, use Finite Differences to compute the matrix-free jacobian ofx -> prob(x, p)
using theδ = 1e-8
which can be passed as an argument.
- For
newton(prob, orbitguess, defOp, options; lens, kwargs...)
This is the deflated Newton-Krylov Solver for computing a periodic orbit using a (Standard / Poincaré) Shooting method.
Arguments
Similar to newton
except that prob
is either a ShootingProblem
or a PoincareShootingProblem
.
Optional argument
jacobian
Specify the choice of the linear algorithm, which must belong to[AutoDiffMF(), MatrixFree(), AutodiffDense(), AutoDiffDenseAnalytical(), FiniteDifferences(), FiniteDifferencesMF()]
. This is used to select a way of inverting the jacobian dG- For
MatrixFree()
, matrix free jacobian, the jacobian is specified by the user inprob
. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system - For
AutoDiffMF()
, we use Automatic Differentiation (AD) to compute the (matrix-free) derivative ofx -> prob(x, p)
using a directional derivative. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system - For
AutodiffDense()
. Same as forAutoDiffMF
but the jacobian is formed as a dense Matrix. You can use a direct solver or an iterative one. - For
FiniteDifferences()
, same as forAutoDiffDense
but we use Finite Differences to compute the jacobian ofx -> prob(x, p)
using theδ = 1e-8
which can be passed as an argument. - For
AutoDiffDenseAnalytical()
. Same as forAutoDiffDense
but the jacobian is formed using a mix of AD and analytical formula. - For
FiniteDifferencesMF()
, use Finite Differences to compute the matrix-free jacobian ofx -> prob(x, p)
using theδ = 1e-8
which can be passed as an argument.
- For
Output:
- solution::NonLinearSolution, see
NonLinearSolution
newton(probPO, orbitguess, options; kwargs...)
This is the Krylov-Newton Solver for computing a periodic orbit using a functional G based on Finite Differences and a Trapezoidal rule.
Arguments:
prob
a problem of typePeriodicOrbitTrapProblem
encoding the functional Gorbitguess
a guess for the periodic orbit whereorbitguess[end]
is an estimate of the period of the orbit. It should be a vector of sizeN * M + 1
whereM
is the number of time slices,N
is the dimension of the phase space. This must be compatible with the numbersN, M
inprob
.par
parameters to be passed to the functionaloptions
same as for the regularnewton
method
Specify the choice of the jacobian (and linear algorithm), jacobian
must belong to [:FullLU, :FullSparseInplace, :Dense, :DenseAD, :BorderedLU, :BorderedSparseInplace, :FullMatrixFree, :BorderedMatrixFree, :FullMatrixFreeAD]
. This is used to select a way of inverting the jacobian dG
of the functional G.
- For
jacobian = :FullLU
, we use the default linear solver based on a sparse matrix representation ofdG
. This matrix is assembled at each newton iteration. This is the default algorithm. - For
jacobian = :FullSparseInplace
, this is the same as for:FullLU
but the sparse matrixdG
is updated inplace. This method allocates much less. In some cases, this is significantly faster than using:FullLU
. Note that this method can only be used if the sparsity pattern of the jacobian is always the same. - For
jacobian = :Dense
, same as above but the matrixdG
is dense. It is also updated inplace. This option is useful to study ODE of small dimension. - For
jacobian = :DenseAD
, evaluate the jacobian using ForwardDiff - For
jacobian = :BorderedLU
, we take advantage of the bordered shape of the linear solver and use a LU decomposition to invertdG
using a bordered linear solver. - For
jacobian = :BorderedSparseInplace
, this is the same as for:BorderedLU
but the cyclic matrixdG
is updated inplace. This method allocates much less. In some cases, this is significantly faster than using:BorderedLU
. Note that this method can only be used if the sparsity pattern of the jacobian is always the same. - For
jacobian = :FullMatrixFree
, a matrix free linear solver is used fordG
: note that a preconditioner is very likely required here because of the cyclic shape ofdG
which affects negatively the convergence properties of GMRES. - For
jacobian = :BorderedMatrixFree
, a matrix free linear solver is used but forJc
only (see docs): it means thatoptions.linsolver
is used to invertJc
. These two Matrix-Free options thus expose different part of the jacobiandG
in order to use specific preconditioners. For example, an ILU preconditioner onJc
could remove the constraints indG
and lead to poor convergence. Of course, for these last two methods, a preconditioner is likely to be required. - For
jacobian = :FullMatrixFreeAD
, the evaluation map of the differential is derived using automatic differentiation. Thus, unlike the previous two cases, the user does not need to pass a Matrix-Free differential.
newton(probPO, orbitguess, defOp, options; kwargs...)
This function is similar to newton(probPO, orbitguess, options, jacobianPO; kwargs...)
except that it uses deflation in order to find periodic orbits different from the ones stored in defOp
. We refer to the mentioned method for a full description of the arguments. The current method can be used in the vicinity of a Hopf bifurcation to prevent the Newton-Krylov algorithm from converging to the equilibrium point.
newton(probPO, orbitguess, options; kwargs...)
This is the Newton Solver for computing a periodic orbit using orthogonal collocation method. Note that the linear solver has to be apropriately set up in options
.
Arguments
Similar to newton
except that prob
is a PeriodicOrbitOCollProblem
.
prob
a problem of type<: PeriodicOrbitOCollProblem
encoding the shooting functional G.orbitguess
a guess for the periodic orbit.options
same as for the regularnewton
method.
Optional argument
jacobian
Specify the choice of the linear algorithm, which must belong to(AutoDiffDense(), )
. This is used to select a way of inverting the jacobian dG- For
AutoDiffDense()
. The jacobian is formed as a dense Matrix. You can use a direct solver or an iterative one usingoptions
. The jacobian is formed inplace. - For
DenseAnalytical()
Same as forAutoDiffDense
but the jacobian is formed using a mix of AD and analytical formula.
- For
newton(probPO, orbitguess, defOp, options; kwargs...)
This function is similar to newton(probPO, orbitguess, options, jacobianPO; kwargs...)
except that it uses deflation in order to find periodic orbits different from the ones stored in defOp
. We refer to the mentioned method for a full description of the arguments. The current method can be used in the vicinity of a Hopf bifurcation to prevent the Newton-Krylov algorithm from converging to the equilibrium point.
Continuation
BifurcationKit.DotTheta
— Typestruct DotTheta{Tdot, Ta}
dot::Any
: dot product used in pseudo-arclength constraintapply!::Any
: Linear operator associated with dot product, i.e. dot(x, y) = <x, Ay>, where <,> is the standard dot product on R^N. You must provide an inplace function which evaluates A. For examplex -> rmul!(x, 1/length(x))
.
This parametric type allows to define a new dot product from the one saved in dt::dot
. More precisely:
dt(u1, u2, p1::T, p2::T, theta::T) where {T <: Real}
computes, the weighted dot product $\langle (u_1,p_1), (u_2,p_2)\rangle_\theta = \theta \Re \langle u_1,u_2\rangle +(1-\theta)p_1p_2$ where $u_i\in\mathbb R^N$. The $\Re$ factor is put to ensure a real valued result despite possible complex valued arguments.
This is used in the pseudo-arclength constraint with the dot product $\frac{1}{N} \langle u_1, u_2\rangle,\quad u_i\in\mathbb R^N$
BifurcationKit.continuation
— Functioncontinuation(
prob,
alg,
contparams;
linear_algo,
bothside,
kwargs...
)
Compute the continuation curve associated to the functional F
which is stored in the bifurcation problem prob
. General information is available in Continuation methods: introduction.
Arguments:
prob::AbstractBifurcationFunction
a::AbstractBifurcationProblem
, typically aBifurcationProblem
which holds the vector field and its jacobian. We also refer toBifFunction
for more details.alg
continuation algorithm, for exampleNatural(), PALC(), Multiple(),...
. See algoscontparams::ContinuationPar
parameters for continuation. SeeContinuationPar
Optional Arguments:
plot = false
whether to plot the solution/branch/spectrum while computing the branchbothside = true
compute the branches on the two sides of the initial parameter valuep0
, merge them and return it.normC = norm
norm used in the nonlinear solvesfilename
to save the computed branch during continuation. The identifier .jld2 will be appended to this filename. This requiresusing JLD2
.callback_newton
callback for newton iterations. See docs ofnewton
. For example, it can be used to change the preconditioners.finalise_solution = (z, tau, step, contResult; kwargs...) -> true
Function called at the end of each continuation step. Can be used to alter the continuation procedure (stop it by returningfalse
), save personal data, plot... The notations arez = BorderedArray(x, p)
wherex
(resp.p
) is the current solution (resp. parameter value),tau::BorderedArray
is the tangent atz
,step::Int
is the index of the current continuation step andcontResult
is the current branch. For advanced use:- the state
state::ContState
of the continuation iterator is passed inkwargs
. This can be used for testing whether this is called from bisection for locating bifurcation points / events:in_bisection(state)
for example. This allows to escape some personal code in this case.
- the iterator
iter::ContIterable
of the continuation is passed inkwargs
.
- the state
verbosity::Int = 0
controls the amount of information printed during the continuation process. Must belong to{0,1,2,3}
. In casecontparams.newton_options.verbose = false
, the following is valid (otherwise the newton iterations are shown). Each case prints more information than the previous one:- case 0: print nothing
- case 1: print basic information about the continuation: used predictor, step size and parameter values
- case 2: print newton iterations number, stability of solution, detected bifurcations / events
- case 3: print information during bisection to locate bifurcations / events
linear_algo
set the linear solver for the continuation algorithmalg.
For example,PALC
needs a linear solver for an enlarged problem (sizen+1
instead ofn
) and one thus needs to tune the one passed incontparams.newton_options.linsolver
. This is a convenient argument to thus change thealg
linear solver and is used mostly internally. The proper way is to pass directly toalg
the correct linear solver.kind::AbstractContinuationKind
[Internal] flag to describe continuation kind (equilibrium, codim 2, ...). Default =EquilibriumCont()
Output:
contres::ContResult
composite type which contains the computed branch. SeeContResult
for more information.
Just change the sign of ds
in ContinuationPar
.
Use debug mode to access more irformation about the progression of the continuation run, like iterative solvers convergence, problem update, ...
continuation(
prob,
algdc,
contParams;
verbosity,
plot,
linear_algo,
dot_palc,
callback_newton,
filename,
normC,
kwcont...
)
This function computes the set of curves of solutions γ(s) = (x(s), p(s))
to the equation F(x,p) = 0
based on the algorithm of deflated continuation as described in Farrell, Patrick E., Casper H. L. Beentjes, and Ásgeir Birkisson. “The Computation of Disconnected Bifurcation Diagrams.” ArXiv:1603.00809 [Math], March 2, 2016. http://arxiv.org/abs/1603.00809.
Depending on the options in contParams
, it can locate the bifurcation points on each branch. Note that you can specify different predictors using alg
.
Arguments:
prob::AbstractBifurcationProblem
bifurcation problemalg::DefCont
, deflated continuation algorithm, seeDefCont
contParams
parameters for continuation. SeeContinuationPar
for more information about the options
Optional Arguments:
plot = false
whether to plot the solution while computing,callback_newton
callback for newton iterations. see docs fornewton
. Can be used to change preconditioners or affect the newton iterations. In the deflation part of the algorithm, when seeking for new branches, the callback is passed the keyword argumentfromDeflatedNewton = true
to tell the user can it is not in the continuation part (regular newton) of the algorithm,verbosity::Int
controls the amount of information printed during the continuation process. Must belong to{0,⋯,5}
,normC = norm
norm used in the Newton solves,dot_palc = (x, y) -> dot(x, y) / length(x)
, dot product used to define the weighted dot product (resp. norm) $\|(x, p)\|^2_\theta$ in the constraint $N(x, p)$ (see online docs on PALC). This argument can be used to remove the factor1/length(x)
for example in problems where the dimension of the state space changes (mesh adaptation, ...),
Outputs:
contres::DCResult
composite type which contains the computed branches. SeeContResult
for more information,
continuation(br, ind_bif, lens2; ...)
continuation(
br,
ind_bif,
lens2,
options_cont;
prob,
start_with_eigen,
detect_codim2_bifurcation,
update_minaug_every_step,
kwargs...
)
Codimension 2 continuation of Fold / Hopf points. This function turns an initial guess for a Fold / Hopf point into a curve of Fold / Hopf points based on a Minimally Augmented formulation. The arguments are as follows
br
results returned after a call to continuationind_bif
bifurcation index inbr
lens2
second parameter used for the continuation, the first one is the one used to computebr
, e.g.getlens(br)
options_cont = br.contparams
arguments to be passed to the regular continuation
Optional arguments:
bdlinsolver
bordered linear solver for the constraint equationupdate_minaug_every_step
update vectorsa, b
in Minimally Formulation everyupdate_minaug_every_step
stepsstart_with_eigen = false
whether to start the Minimally Augmented problem with information from eigen elementsdetect_codim2_bifurcation ∈ {0,1,2}
whether to detect Bogdanov-Takens, Bautin and Cusp. If equals1
non precise detection is used. If equals2
, a bisection method is used to locate the bifurcations.kwargs
keywords arguments to be passed to the regular continuation
where the parameters are as above except that you have to pass the branch br
from the result of a call to continuation
with detection of bifurcations enabled and index
is the index of Hopf point in br
you want to refine.
For ODE problems, it is more efficient to use the Matrix based Bordered Linear Solver passing the option bdlinsolver = MatrixBLS()
It is recommended that you use the option start_with_eigen = true
continuation(
prob,
x0,
par0,
x1,
p1,
alg,
lens,
contParams;
bothside,
kwargs...
)
[Internal] This function is not meant to be called directly.
This function is the analog of continuation
when the first two points on the branch are passed (instead of a single one). Hence x0
is the first point on the branch (with palc s=0
) with parameter par0
and x1
is the second point with parameter set(par0, lens, p1)
.
continuation(br, ind_bif; ...)
continuation(
br,
ind_bif,
options_cont;
alg,
δp,
ampfactor,
nev,
usedeflation,
verbosedeflation,
max_iter_deflation,
perturb,
plot_solution,
Teigvec,
scaleζ,
tol_fold,
kwargs_deflated_newton,
kwargs...
)
Automatic branch switching at branch points based on a computation of the normal form. More information is provided in Branch switching. An example of use is provided in 2d generalized Bratu–Gelfand problem.
Arguments
br
branch result from a call tocontinuation
ind_bif
index of the bifurcation point inbr
from which you want to branch fromoptions_cont
options for the call tocontinuation
Optional arguments
alg = br.alg
continuation algorithm to be used, default value:br.alg
δp
used to specify a specific value for the parameter on the bifurcated branch which is otherwise determined byoptions_cont.ds
. This allows to use a step larger thanoptions_cont.dsmax
.ampfactor = 1
factor to alter the amplitude of the bifurcated solution. Useful to magnify the bifurcated solution when the bifurcated branch is very steep.nev
number of eigenvalues to be computed to get the right eigenvectorusedeflation = false
whether to use nonlinear deflation (see Deflated problems) to help finding the guess on the bifurcatedverbosedeflation
print deflated newton iterationsmax_iter_deflation
number of newton steps in deflated newtonperturb = identity
which perturbation function to use during deflated newtonTeigvec = _getvectortype(br)
type of the eigenvector. Useful whenbr
was loaded from a file and this information was lostscaleζ = norm
pass a norm to normalize vectors during normal form computationplot_solution
change plot solution method in the problembr.prob
kwargs
optional arguments to be passed tocontinuation
, the regularcontinuation
one and toget_normal_form
.
In the case of a very large model and use of special hardware (GPU, cluster), we suggest to discouple the computation of the reduced equation, the predictor and the bifurcated branches. Have a look at methods(BifurcationKit.multicontinuation)
to see how to call these versions. These methods has been tested on GPU with very high memory pressure.
continuation(
probPO,
orbitguess,
alg,
contParams,
linear_algo;
δ,
eigsolver,
record_from_solution,
plot_solution,
kwargs...
)
This is the continuation method for computing a periodic orbit using a (Standard / Poincaré) Shooting method.
Arguments
Similar to continuation
except that probPO
is either a ShootingProblem
or a PoincareShootingProblem
. By default, it prints the period of the periodic orbit.
Optional arguments
eigsolver
specify an eigen solver for the computation of the Floquet exponents, defaults toFloquetQaD
jacobian
Specify the choice of the linear algorithm, which must belong to[AutoDiffMF(), MatrixFree(), AutodiffDense(), AutoDiffDenseAnalytical(), FiniteDifferences(), FiniteDifferencesMF()]
. This is used to select a way of inverting the jacobian dG- For
MatrixFree()
, matrix free jacobian, the jacobian is specified by the user inprob
. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system - For
AutoDiffMF()
, we use Automatic Differentiation (AD) to compute the (matrix-free) derivative ofx -> prob(x, p)
using a directional derivative. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system - For
AutodiffDense()
. Same as forAutoDiffMF
but the jacobian is formed as a dense Matrix. You can use a direct solver or an iterative one. - For
FiniteDifferences()
, same as forAutoDiffDense
but we use Finite Differences to compute the jacobian ofx -> prob(x, p)
using theδ = 1e-8
which can be passed as an argument. - For
AutoDiffDenseAnalytical()
. Same as forAutoDiffDense
but the jacobian is formed using a mix of AD and analytical formula. - For
FiniteDifferencesMF()
, use Finite Differences to compute the matrix-free jacobian ofx -> prob(x, p)
using theδ = 1e-8
which can be passed as an argument.
- For
continuation(
prob,
orbitguess,
alg,
_contParams;
linear_algo,
kwargs...
)
This is the continuation routine for computing a periodic orbit using a (Standard / Poincaré) Shooting method.
Arguments
Similar to continuation
except that prob
is either a ShootingProblem
or a PoincareShootingProblem
. By default, it prints the period of the periodic orbit.
Optional argument
linear_algo::AbstractBorderedLinearSolver
jacobian
Specify the choice of the linear algorithm, which must belong to[AutoDiffMF(), MatrixFree(), AutodiffDense(), AutoDiffDenseAnalytical(), FiniteDifferences(), FiniteDifferencesMF()]
. This is used to select a way of inverting the jacobian dG- For
MatrixFree()
, matrix free jacobian, the jacobian is specified by the user inprob
. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system - For
AutoDiffMF()
, we use Automatic Differentiation (AD) to compute the (matrix-free) derivative ofx -> prob(x, p)
using a directional derivative. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system - For
AutodiffDense()
. Same as forAutoDiffMF
but the jacobian is formed as a dense Matrix. You can use a direct solver or an iterative one. - For
FiniteDifferences()
, same as forAutoDiffDense
but we use Finite Differences to compute the jacobian ofx -> prob(x, p)
using theδ = 1e-8
which can be passed as an argument. - For
AutoDiffDenseAnalytical()
. Same as forAutoDiffDense
but the jacobian is formed using a mix of AD and analytical formula. - For
FiniteDifferencesMF()
, use Finite Differences to compute the matrix-free jacobian ofx -> prob(x, p)
using theδ = 1e-8
which can be passed as an argument.
- For
continuation(
br,
ind_bif,
_contParams,
pbPO;
prob_vf,
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 to 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.
continuation(
br,
ind_bif,
_contParams;
alg,
δp,
ampfactor,
usedeflation,
linear_algo,
detailed,
prm,
override,
kwargs...
)
Branch switching at a bifurcation point on a branch of periodic orbits (PO) specified by a br::AbstractBranchResult
. The functional used to compute the PO is br.prob
. A deflated Newton-Krylov solver can be used to improve the branch switching capabilities.
Arguments
br
branch of periodic orbits computed with aPeriodicOrbitTrapProblem
ind_bif
index of the branch point_contParams
parameters to be used by a regularcontinuation
Optional arguments
δp = 0.1
used to specify a particular guess for the parameter in the branch which is otherwise determined bycontParams.ds
. This allows to use a step larger thancontParams.dsmax
.ampfactor = 1
factor which alters the amplitude of the bifurcated solution. Useful to magnify the bifurcated solution when the bifurcated branch is very steep.detailed = false
whether to fully compute the normal form.usedeflation = true
whether to use nonlinear deflation (see Deflated problems) to help finding the guess on the bifurcated branchrecord_from_solution = (u, p) -> u[end]
, record method used in the bifurcation diagram, by default this records the period of the periodic orbit.linear_algo = BorderingBLS()
, same as forcontinuation
kwargs
keywords arguments used for a call to the regularcontinuation
and the ones specific to periodic orbits (POs).
continuation(
prob,
orbitguess,
alg,
_contParams;
record_from_solution,
linear_algo,
kwargs...
)
This is the continuation routine for computing a periodic orbit using a functional G based on Finite Differences and a Trapezoidal rule.
Arguments
prob::PeriodicOrbitTrapProblem
encodes the functional Gorbitguess
a guess for the periodic orbit whereorbitguess[end]
is an estimate of the period of the orbit. It could be a vector of sizeN * M + 1
whereM
is the number of time slices,N
is the dimension of the phase space. This must be compatible with the numbersN, M
inprob
.alg
continuation algorithmcontParams
same as for the regularcontinuation
method
Keyword arguments
linear_algo
same as incontinuation
Specify the choice of the jacobian (and linear algorithm), jacobian
must belong to [:FullLU, :FullSparseInplace, :Dense, :DenseAD, :BorderedLU, :BorderedSparseInplace, :FullMatrixFree, :BorderedMatrixFree, :FullMatrixFreeAD]
. This is used to select a way of inverting the jacobian dG
of the functional G.
- For
jacobian = :FullLU
, we use the default linear solver based on a sparse matrix representation ofdG
. This matrix is assembled at each newton iteration. This is the default algorithm. - For
jacobian = :FullSparseInplace
, this is the same as for:FullLU
but the sparse matrixdG
is updated inplace. This method allocates much less. In some cases, this is significantly faster than using:FullLU
. Note that this method can only be used if the sparsity pattern of the jacobian is always the same. - For
jacobian = :Dense
, same as above but the matrixdG
is dense. It is also updated inplace. This option is useful to study ODE of small dimension. - For
jacobian = :DenseAD
, evaluate the jacobian using ForwardDiff - For
jacobian = :BorderedLU
, we take advantage of the bordered shape of the linear solver and use a LU decomposition to invertdG
using a bordered linear solver. - For
jacobian = :BorderedSparseInplace
, this is the same as for:BorderedLU
but the cyclic matrixdG
is updated inplace. This method allocates much less. In some cases, this is significantly faster than using:BorderedLU
. Note that this method can only be used if the sparsity pattern of the jacobian is always the same. - For
jacobian = :FullMatrixFree
, a matrix free linear solver is used fordG
: note that a preconditioner is very likely required here because of the cyclic shape ofdG
which affects negatively the convergence properties of GMRES. - For
jacobian = :BorderedMatrixFree
, a matrix free linear solver is used but forJc
only (see docs): it means thatoptions.linsolver
is used to invertJc
. These two Matrix-Free options thus expose different part of the jacobiandG
in order to use specific preconditioners. For example, an ILU preconditioner onJc
could remove the constraints indG
and lead to poor convergence. Of course, for these last two methods, a preconditioner is likely to be required. - For
jacobian = :FullMatrixFreeAD
, the evaluation map of the differential is derived using automatic differentiation. Thus, unlike the previous two cases, the user does not need to pass a Matrix-Free differential.
Note that by default, the method prints the period of the periodic orbit as function of the parameter. This can be changed by providing your record_from_solution
argument.
continuation(
probPO,
orbitguess,
alg,
_contParams,
linear_algo;
δ,
eigsolver,
record_from_solution,
plot_solution,
kwargs...
)
This is the continuation method for computing a periodic orbit using an orthogonal collocation method.
Arguments
Similar to continuation
except that prob
is a PeriodicOrbitOCollProblem
. By default, it prints the period of the periodic orbit.
Keywords arguments
eigsolver
specify an eigen solver for the computation of the Floquet exponents, defaults toFloquetQaD
Continuation algorithms
BifurcationKit.Natural
— TypeNatural continuation algorithm. The predictor is the constant predictor and the parameter is incremented by `ContinuationPar().ds` at each continuation step.
BifurcationKit.Secant
— TypeSecant Tangent predictor
BifurcationKit.Bordered
— TypeBordered Tangent predictor
BifurcationKit.Polynomial
— TypePolynomial Tangent predictor
n::Int64
: Order of the polynomialk::Int64
: Length of the last solutions vector used for the polynomial fitA::Matrix{T} where T<:Real
: Matrix for the interpolationtangent::BifurcationKit.AbstractTangentComputation
: Algo for tangent when polynomial predictor is not possiblesolutions::DataStructures.CircularBuffer
: Vector of solutionsparameters::DataStructures.CircularBuffer{T} where T<:Real
: Vector of parametersarclengths::DataStructures.CircularBuffer{T} where T<:Real
: Vector of arclengthscoeffsSol::Vector
: Coefficients for the polynomials for the solutioncoeffsPar::Vector{T} where T<:Real
: Coefficients for the polynomials for the parameterupdate::Bool
: Update the predictor by adding the last point (x, p)? This can be disabled in order to just use the polynomial prediction. It is useful when the predictor is called mutiple times during bifurcation detection using bisection.
Constructor(s)
Polynomial(pred, n, k, v0)
Polynomial(n, k, v0)
n
order of the polynomialk
length of the last solutions vector used for the polynomial fitv0
example of solution to be stored. It is only used to get theeltype
of the tangent.
Can be used like
PALC(tangent = Polynomial(Bordered(), 2, 6, rand(1)))
BifurcationKit.Multiple
— TypeMultiple Tangent continuation algorithm.
alg::PALC
: Tangent predictor used Default: PALC()τ::Any
: Save the current tangentα::Real
: Damping in Newton iterations, 0 < α < 1nb::Int64
: Number of predictorscurrentind::Int64
: Index of the largest converged predictor Default: 0pmimax::Int64
: Index for lookup in residual history Default: 1imax::Int64
: Maximum index for lookup in residual history Default: 4dsfact::Real
: Factor to increase ds upon successful step Default: 1.5
Constructor(s)
Multiple(alg, x0, α, n)
Multiple(pred, x0, α, n)
Multiple(x0, α, n)
Missing docstring for PALC
. Check Documenter's build log for details.
BifurcationKit.MoorePenrose
— TypeMoore-Penrose predictor / corrector
Moore-Penrose continuation algorithm.
Additional information is available on the website.
Constructors
alg = MoorePenrose()
alg = MoorePenrose(tangent = PALC())
Fields
tangent::Any
: Tangent predictor, for examplePALC()
method::MoorePenroseLS
: Moore Penrose linear solver. Can be BifurcationKit.direct, BifurcationKit.pInv or BifurcationKit.iterativels::BifurcationKit.AbstractLinearSolver
: (Bordered) linear solver
BifurcationKit.DefCont
— Typestruct DefCont{Tdo, Talg, Tps, Tas, Tud, Tk} <: BifurcationKit.AbstractContinuationAlgorithm
Structure which holds the parameters specific to Deflated continuation.
Fields
deflation_operator::Any
: Deflation operator,::DeflationOperator
Default: nothingalg::Any
: Used as a predictor,::AbstractContinuationAlgorithm
. For examplePALC()
,Natural()
,... Default: PALC()max_branches::Int64
: maximum number of (active) branches to be computed Default: 100seek_every_step::Int64
: whether to seek new (deflated) solution at every step Default: 1max_iter_defop::Int64
: maximum number of deflated Newton iterations Default: 5perturb_solution::Any
: perturb function Default: _perturbSolutionaccept_solution::Any
: accept (solution) function Default: _acceptSolutionupdate_deflation_op::Any
: function to update the deflation operator, ie pushing new solutions Default: _updateDeflationOpjacobian::Any
: jacobian for deflated newton. Can beDeflatedProblemCustomLS()
, orVal(:autodiff)
,Val(:fullIterative)
Default: DeflatedProblemCustomLS()
Missing docstring for AsymptoticNumericalMethod.ANM
. Check Documenter's build log for details.
Events
BifurcationKit.DiscreteEvent
— Typestruct DiscreteEvent{Tcb, Tl} <: BifurcationKit.AbstractDiscreteEvent
Structure to pass a DiscreteEvent function to the continuation algorithm. A discrete call back returns a discrete value and we seek when it changes.
nb::Int64
: number of events, ie the length of the result returned by the callback functioncondition::Any
: =(iter, state) -> NTuple{nb, Int64}
callback function which at each continuation state, returns a tuple. For example, to detect a value change.computeEigenElements::Bool
: whether the event requires to compute eigen elementslabels::Any
: Labels used to display information. For examplelabels[1]
is used to qualify an event occurring in the first component. You can uselabels = ("hopf",)
orlabels = ("hopf", "fold")
. You must havelabels::Union{Nothing, NTuple{N, String}}
.
BifurcationKit.ContinuousEvent
— Typestruct ContinuousEvent{Tcb, Tl, T} <: BifurcationKit.AbstractContinuousEvent
Structure to pass a ContinuousEvent function to the continuation algorithm. A continuous call back returns a tuple/scalar value and we seek its zeros.
nb::Int64
: number of events, ie the length of the result returned by the callback functioncondition::Any
: ,(iter, state) -> NTuple{nb, T}
callback function which, at each continuation state, returns a tuple. For example, to detect crossing 1.0 and -2.0, you can pass(iter, state) -> (getp(state)+2, getx(state)[1]-1)),
. Note that the typeT
should match the one of the parameter specified by the::Lens
incontinuation
.computeEigenElements::Bool
: whether the event requires to compute eigen elementslabels::Any
: Labels used to display information. For examplelabels[1]
is used to qualify an event of the type(0, 1.3213, 3.434)
. You can uselabels = ("hopf",)
orlabels = ("hopf", "fold")
. You must havelabels::Union{Nothing, NTuple{N, String}}
.tol::Any
: Tolerance on event value to declare it as true event.
BifurcationKit.SetOfEvents
— Typestruct SetOfEvents{Tc<:Tuple, Td<:Tuple} <: BifurcationKit.AbstractEvent
Multiple events can be chained together to form a SetOfEvents
. A SetOfEvents
is constructed by passing to the constructor ContinuousEvent
, DiscreteEvent
or other SetOfEvents
instances:
SetOfEvents(cb1, cb2, cb3)
Example
BifurcationKit.SetOfEvents(BK.FoldDetectCB, BK.BifDetectCB)
You can pass as many events as you like.
eventC::Tuple
: Continuous eventeventD::Tuple
: Discrete event
BifurcationKit.PairOfEvents
— Typestruct PairOfEvents{Tc<:BifurcationKit.AbstractContinuousEvent, Td<:BifurcationKit.AbstractDiscreteEvent} <: BifurcationKit.AbstractEvent
Structure to pass a PairOfEvents function to the continuation algorithm. It is composed of a pair ContinuousEvent / DiscreteEvent. A PairOfEvents
is constructed by passing to the constructor a ContinuousEvent
and a DiscreteEvent
:
PairOfEvents(contEvent, discreteEvent)
Fields
eventC::BifurcationKit.AbstractContinuousEvent
: Continuous eventeventD::BifurcationKit.AbstractDiscreteEvent
: Discrete event
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,
pbPO;
prob_vf,
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 to 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.
Bifurcation diagram
BifurcationKit.bifurcationdiagram
— Functionbifurcationdiagram(
prob,
alg,
level,
options;
linear_algo,
kwargs...
)
Compute the bifurcation diagram associated with the problem F(x, p) = 0
recursively.
Arguments
prob::AbstractBifurcationProblem
bifurcation problemalg
continuation algorithmlevel
maximum branching (or recursion) level for computing the bifurcation diagramoptions = (x, p, level) -> contparams
this function allows to change thecontinuation
options depending on the branchinglevel
. The argumentx, p
denotes the current solution toF(x, p)=0
.kwargs
optional arguments. Look atbifurcationdiagram!
for more details.
Simplified call:
We also provide the method
bifurcationdiagram(prob, br::ContResult, level::Int, options; kwargs...)
where br
is a branch computed after a call to continuation
from which we want to compute the bifurcating branches recursively.
BifurcationKit.bifurcationdiagram!
— Functionbifurcationdiagram!(
prob,
node,
maxlevel,
options;
code,
halfbranch,
verbosediagram,
kwargs...
)
Similar to bifurcationdiagram
but you pass a previously computed node
from which you want to further compute the bifurcated branches. It is usually used with node = get_branch(diagram, code)
from a previously computed bifurcation diagram
.
Arguments
node::BifDiagNode
a node in the bifurcation diagrammaxlevel = 1
required maximal level of recursion.options = (x, p, level) -> contparams
this function allows to change thecontinuation
options depending on the branchinglevel
. The argumentx, p
denotes the current solution toF(x, p)=0
.
Optional arguments
code = "0"
code used to display iterationsusedeflation = false
halfbranch = false
for Pitchfork / Transcritical bifurcations, compute only half of the branch. Can be useful when there are symmetries.verbosediagram
verbose specific to bifurcation diagram. Print information about the branches as they are being computed.kwargs
optional arguments as forcontinuation
but also for the different versions listed in Continuation.
BifurcationKit.get_branch
— Functionget_branch(diagram, code)
Return the part of the diagram (bifurcation diagram) by recursively descending down the diagram using the Int
valued tuple code
. For example get_branch(diagram, (1,2,3,))
returns diagram.child[1].child[2].child[3]
.
BifurcationKit.get_branches_from_BP
— Functionget_branches_from_BP(diagram, indbif)
Return the part of the diagram corresponding to the indbif-th bifurcation point on the root branch.
BifurcationKit.SpecialPoint
— Typestruct SpecialPoint{T, Tp, Tv, Tvτ} <: BifurcationKit.AbstractBifurcationPoint
Structure to record special points on a curve. There are two types of special points that are recorded in this structure: bifurcation points and events (see https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/EventCallback/).
type::Symbol
: Description of the special points. In case of Events, this field records the user passed named to the event, or the default:userD
,:userC
. In case of bifurcation points, it can be one of the following:- :bp Bifurcation point, simple eigenvalue crossing the imaginary axis - :fold Fold point - :hopf Hopf point - :nd not documented bifurcation point. Detected by multiple eigenvalues crossing. Generally occurs in problems with symmetries or in cases where the continuation step size is too large and merge two different bifurcation points. - :cusp Cusp point - :gh Generalized Hopf point (also called Bautin point) - :bt Bogdanov-Takens point - :zh Zero-Hopf point - :hh Hopf-Hopf point - :ns Neimark-Sacker point - :pd Period-doubling point - :R1 Strong resonance 1:1 of periodic orbits - :R2 Strong resonance 1:2 of periodic orbits - :R3 Strong resonance 1:3 of periodic orbits - :R4 Strong resonance 1:4 of periodic orbits - :foldFlip Fold / Flip of periodic orbits - :foldNS Fold / Neimark-Sacker of periodic orbits - :pdNS Period-Doubling / Neimark-Sacker of periodic orbits - :gpd Generalized Period-Doubling of periodic orbits - :nsns Double Neimark-Sacker of periodic orbits - :ch Chenciner bifurcation of periodic orbits Default: :none
idx::Int64
: Index inbr.branch
orbr.eig
(seeContResult
) for which the bifurcation occurs. Default: 0param::Any
: Parameter value at the special point (this is an estimate). Default: 0.0norm::Any
: Norm of the equilibrium at the special point Default: 0.0printsol::Any
:printsol = record_from_solution(x, param)
whererecord_from_solution
is one of the arguments tocontinuation
Default: 0.0x::Any
: Equilibrium at the special point Default: Vector{T}(undef, 0)τ::BorderedArray{Tvτ, T} where {T, Tvτ}
: Tangent along the branch at the special point Default: BorderedArray(x, T(0))ind_ev::Int64
: Eigenvalue index responsible for detecting the special point (if applicable) Default: 0step::Int64
: Continuation step at which the special occurs Default: 0status::Symbol
:status ∈ {:converged, :guess, :guessL}
indicates whether the bisection algorithm was successful in detecting the special (bifurcation) point. Ifstatus == :guess
, the bissection algorithm failed to meet the requirements given in::ContinuationPar
. Same forstatus == :guessL
but the bissection algorithm stopped on the left of the bifurcation point. Default: :guessδ::Tuple{Int64, Int64}
:δ = (δr, δi)
where δr indicates the change in the number of unstable eigenvalues and δi indicates the change in the number of unstable eigenvalues with nonzero imaginary part.abs(δr)
is thus an estimate of the dimension of the kernel of the Jacobian at the special (bifurcation) point. Default: (0, 0)precision::Any
: Precision in the location of the special point Default: -1interval::Tuple{T, T} where T
: Interval parameter containing the special point Default: (0, 0)
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]
.
BifurcationKit.SectionSS
— Typestruct SectionSS{Tn, Tc} <: BifurcationKit.AbstractSection
This composite type (named for Section Standard Shooting) encodes a type of section implemented by a single hyperplane. It can be used in conjunction with ShootingProblem
. The hyperplane is defined by a point center
and a normal
.
normal::Any
: Normal to define hyperplanecenter::Any
: Representative point on hyperplane
BifurcationKit.SectionPS
— Typestruct SectionPS{Tn, Tc, Tnb, Tcb, Tr} <: BifurcationKit.AbstractSection
This composite type (named for SectionPoincaréShooting) encodes a type of Poincaré sections implemented by hyperplanes. It can be used in conjunction with PoincareShootingProblem
. Each hyperplane is defined par a point (one example in centers
) and a normal (one example in normals
).
M::Int64
normals::Any
centers::Any
indices::Vector{Int64}
normals_bar::Any
centers_bar::Any
radius::Any
Constructor(s)
SectionPS(normals::Vector{Tv}, centers::Vector{Tv})
Misc.
BifurcationKit.PrecPartialSchurKrylovKit
— FunctionPrecPartialSchurKrylovKit(J, x0, nev, which = :LM; krylovdim = max(2nev, 20), verbosity = 0)
Builds a preconditioner based on deflation of nev
eigenvalues chosen according to which
. A partial Schur decomposition is computed (Matrix-Free), using the package KrylovKit.jl
, from which a projection is built. The options are similar to the ones of EigKrylovKit()
.
BifurcationKit.PrecPartialSchurArnoldiMethod
— FunctionPrecPartialSchurArnoldiMethod(J, N, nev, which = LM(); tol = 1e-9, kwargs...)
Builds a preconditioner based on deflation of nev
eigenvalues chosen according to which
. A partial Schur decomposition is computed (Matrix-Free), using the package ArnoldiMethod.jl
, from which a projection is built. See the package ArnoldiMethod.jl
for how to pass the proper options.
BifurcationKit.Flow
— Typestruct Flow{TF, Tf, Tts, Tff, Td, Tad, Tse, Tprob, TprobMono, Tfs, Tcb, Tδ} <: BifurcationKit.AbstractFlow
F::Any
: The vector field(x, p) -> F(x, p)
associated to a Cauchy problem. Used for the differential of the shooting problem. Default: nothingflow::Any
: The flow (or semigroup)(x, p, t) -> flow(x, p, t)
associated to the Cauchy problem. Only the last time point must be returned in the form (u = ...) Default: nothingflowTimeSol::Any
: Flow which returns the tuple (t, u(t)). Optional, mainly used for plotting on the user side. Default: nothingflowFull::Any
: [Optional] The flow (or semigroup) associated to the Cauchy problem(x, p, t) -> flow(x, p, t)
. The whole solution on the time interval [0,t] must be returned. It is not strictly necessary to provide this, it is mainly used for plotting on the user side. Please usenothing
as default. Default: nothingjvp::Any
: The differentialdflow
of the flow w.r.t.x
,(x, p, dx, t) -> dflow(x, p, dx, t)
. One important thing is that we requiredflow(x, dx, t)
to return a Named Tuple:(t = t, u = flow(x, p, t), du = dflow(x, p, dx, t))
, the last component being the value of the derivative of the flow. Default: nothingvjp::Any
: The adjoint differentialvjpflow
of the flow w.r.t.x
,(x, p, dx, t) -> vjpflow(x, p, dx, t)
. One important thing is that we requirevjpflow(x, p, dx, t)
to return a Named Tuple:(t = t, u = flow(x, p, t), du = vjpflow(x, p, dx, t))
, the last component being the value of the derivative of the flow. Default: nothingjvpSerial::Any
: [Optional] Serial version of dflow. Used internally when using parallel multiple shooting. Please usenothing
as default. Default: nothingprob::Any
: [Internal] store the ODEProblem associated to the flow of the Cauchy problem Default: nothingprobMono::Any
: [Internal] store the ODEProblem associated to the flow of the variational problem Default: nothingflowSerial::Any
: [Internal] Serial version of the flow Default: nothingcallback::Any
: [Internal] Store possible callback Default: nothingdelta::Any
: [Internal] Default: 1.0e-8
Simplified constructor(s)
We provide a simple constructor where you only pass the vector field F
, the flow ϕ
and its differential dϕ
:
fl = Flow(F, ϕ, dϕ)
Simplified constructors for DifferentialEquations.jl
These are some simple constructors for which you only have to pass a prob::ODEProblem
or prob::EnsembleProblem
(for parallel computation) from DifferentialEquations.jl
and an ODE time stepper like Tsit5()
. Hence, you can do for example
fl = Flow(prob, Tsit5(); kwargs...)
where kwargs
is passed to SciMLBase::solve
. If your vector field depends on parameters p
, you can define a Flow
using
fl = Flow(prob, Tsit5(); kwargs...)
Finally, you can pass two ODEProblem
where the second one is used to compute the variational equation:
fl = Flow(prob1::ODEProblem, alg1, prob2::ODEProblem, alg2; kwargs...)
BifurcationKit.FloquetQaD
— Typefloquet = FloquetQaD(eigsolver::AbstractEigenSolver)
This composite type implements the computation of the eigenvalues of the monodromy matrix in the case of periodic orbits problems (based on the Shooting method or Finite Differences (Trapeze method)), also called the Floquet multipliers. The method, dubbed Quick and Dirty (QaD), is not numerically very precise for large / small Floquet exponents when the number of time sections is large because of many matrix products. It allows, nevertheless, to detect bifurcations. The arguments are as follows:
eigsolver::AbstractEigenSolver
solver used to compute the eigenvalues.
If eigsolver == DefaultEig()
, then the monodromy matrix is formed and its eigenvalues are computed. Otherwise, a Matrix-Free version of the monodromy is used.
The computation of Floquet multipliers is necessary for the detection of bifurcations of periodic orbits (which is done by analyzing the Floquet exponents obtained from the Floquet multipliers). Hence, the eigensolver eigsolver
needs to compute the eigenvalues with largest modulus (and not with largest real part which is their default behavior). This can be done by changing the option which = :LM
of eigsolver
. Nevertheless, note that for most implemented eigensolvers in the current Package, the proper option is set.
Missing docstring for guess_from_hopf(br, ind_hopf, eigsolver::AbstractEigenSolver, M, amplitude; phase = 0)
. Check Documenter's build log for details.
BifurcationKit.get_normal_form
— Functionget_normal_form(
prob,
br,
id_bif;
nev,
verbose,
ζs,
ζs_ad,
lens,
Teigvec,
scaleζ,
detailed,
autodiff,
bls,
bls_adjoint,
bls_block
)
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.bls_adjoint = bls
specify Bordered linear solver for the adjoint problem.bls_block = bls
specify Bordered linear solver when the border has dimension 2 (1 forbls
).
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,
autodiff,
detailed,
δ
)
Compute the normal form of periodic orbits. We detail the additional keyword arguments specific to periodic orbits
Optional arguments
prm = true
compute the normal form using Poincaré return map (PRM). If false, use the Iooss normal form.