`BifDetectEvent`

This event implements the detection of bifurcations points along a continuation curve. The detection is based on monitoring the number of unstable eigenvalues. More details are given at Detection of bifurcation points.

`FoldDetectEvent`

This event implements the detection of Fold points based on the p-component of the tangent vector to the continuation curve. It is designed to work with the predictor BorderedPred() that you pass to continuation with the keyword argument tangentAlgo.

Singleton type to trigger the computation of the jacobian using ForwardDiff.jl. It can be used for example in newton or in deflated newton.

mutable struct BTProblemMinimallyAugmented{Tprob<:BifurcationKit.AbstractBifurcationProblem, vectype, S<:BifurcationKit.AbstractLinearSolver, Sa<:BifurcationKit.AbstractLinearSolver, Sbd<:BifurcationKit.AbstractBorderedLinearSolver, Sbda<:BifurcationKit.AbstractBorderedLinearSolver, Tlens<:Setfield.Lens} <: BifurcationKit.AbstractProblemMinimallyAugmented

Structure to encode Fold / Hopf functional based on a Minimally Augmented formulation.

Fields

• prob_vf

  Functional F(x, p) - vector field - with all derivatives

• a

  close to null vector of Jᵗ

• b

  close to null vector of J

• zero

  vector zero, to avoid allocating it many times

• linsolver

  linear solver. Used to invert the jacobian of MA functional

• linsolverAdjoint

  linear solver for the jacobian adjoint

• linbdsolver

  bordered linear solver

• linbdsolverAdjoint

  linear bordered solver for the jacobian adjoint

• lens2

  second parameter axis

• usehessian

  wether to use the hessian of prob_vf

mutable struct Bautin{Tv, Tpar, Tlens, Tevr, Tevl, Tnf} <: BifurcationKit.AbstractBifurcationPoint

• x0::Any

  Bifurcation point

• params::Any

  Parameters used by the vector field.

• lens::Any

  Parameter axis used to compute the branch on which this cusp point was detected.

• ζ::Any

  Right eigenvector

• ζstar::Any

  Left eigenvector

• nf::Any

  Normal form coefficients

• type::Symbol

  Type of bifurcation

struct 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-place result = f(x, p) or inplace f(result, x, p). For type stability, the types of x and result should match

• dF::Any

  Differential of F with respect to x, signature dF(x,p,dx)

• dFad::Any

  Adjoint of the Differential of F with respect to x, signature dFad(x,p,dx)

• J::Any

  Jacobian of F at (x, p). It can assume three forms. 1. Either J is a function and J(x, p) returns a ::AbstractMatrix. In this case, the default arguments of contParams::ContinuationPar will make continuation work. 2. Or J is a function and J(x, p) returns a function taking one argument dx and returning dr of the same type as dx. In our notation, dr = J * dx. In this case, the default parameters of contParams::ContinuationPar will not work and you have to use a Matrix Free linear solver, for example GMRESIterativeSolvers, 3. Or J is a function and J(x, p) returns a variable j which can assume any type. Then, you must implement a linear solver ls as a composite type, subtype of AbstractLinearSolver which is called like ls(j, rhs) and which returns the solution of the jacobian linear system. See for example examples/SH2d-fronts-cuda.jl. This linear solver is passed to NewtonPar(linsolver = ls) which itself passed to ContinuationPar. Similarly, you have to implement an eigensolver eig as a composite type, subtype of AbstractEigenSolver.

• 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 pass Jᵗ = (x, p) -> transpose(dF(x, p)).

• d2F::Any

  Second Differential of F with respect to x, signature d2F(x,p,dx1,dx2)

• d3F::Any

  Third Differential of F with respect to x, signature d3F(x,p,dx1,dx2,dx3)

• d2Fc::Any

  [internal] Second Differential of F with respect to x which accept complex vectors dxi

• d3Fc::Any

  [internal] Third Differential of F with respect to x which accept complex vectors dxi

• isSymmetric::Any

  Whether the jacobian is auto-adjoint.

• δ::Any

  used internally to compute derivatives (with finite differences) w.r.t the parameter p.

• inplace::Bool

  optionally sets whether the function is inplace or not

Methods

• residual(pb::BifFunction, x, p) calls pb.F(x,p)
• jacobian(pb::BifFunction, x, p) calls pb.J(x, p)
• dF(pb::BifFunction, x, p, dx) calls pb.dF(x,p,dx)
• etc

struct BifurcationProblem{Tvf, Tu, Tp, Tl<:Setfield.Lens, Tplot, Trec} <: BifurcationKit.AbstractAllJetBifProblem

Structure to hold the bifurcation problem. If don't have parameters, you can pass nothing.

Fields

• VF::Any

  Vector field, typically a BifFunction

• u0::Any

  Initial guess

• params::Any

  parameters

• lens::Setfield.Lens

  Typically a Setfield.Lens. It specifies which parameter axis among params is used for continuation. For example, if par = (α = 1.0, β = 1), we can perform continuation w.r.t. α by using lens = (@lens _.α). If you have an array par = [ 1.0, 2.0] and want to perform continuation w.r.t. the first variable, you can use lens = (@lens _[1]). For more information, we refer to SetField.jl.

• plotSolution::Any

  user function to plot solutions during continuation. Signature: plotSolution(x, p; kwargs...)

• recordFromSolution::Any

  recordFromSolution = (x, p) -> norm(x) function used record a few indicators about the solution. It could be norm 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 in contres.branch (see below). Finally, the first component is used to plot in the continuation curve.

Methods

• getu0(pb) calls pb.u0
• getParams(pb) calls pb.params
• getLens(pb) calls pb.lens
• getParam(pb) calls get(pb.params, pb.lens)
• setParam(pb, p0) calls set(pb.params, pb.lens, p0)
• recordFromSolution(pb) calls pb.recordFromSolution
• plotSolution(pb) calls pb.plotSolution
• isSymmetric(pb) calls isSymmetric(pb.prob)

Constructors

• BifurcationProblem(F, u0, params, lens; J, Jᵗ, d2F, d3F, kwargs...) and kwargs are the fields above. You can pass your own jacobian with J (see [BifFunction]@ref for description of the jacobian function) and jacobian adjoint with Jᵗ. For example, this can be used to provide finite differences based jacobian.

struct BilinearMap{Tm}

This structure wraps a bilinear map to allow evaluation on Complex arguments. This is especially useful when these maps are produced by ForwardDiff.jl.

mutable struct BogdanovTakens{Tv, Tpar, Tlens, Tevr, Tevl, Tnf, Tnf2} <: BifurcationKit.AbstractBifurcationPoint

• x0::Any

  Bogdanov-Takens point

• params::Any

  Parameters used by the vector field.

• lens::Any

  Parameter axis used to compute the branch on which this BT point was detected.

• ζ::Any

  Right eigenvectors

• ζstar::Any

  Left eigenvectors

• nf::Any

  Normal form coefficients (basic)

• nfsupp::Any

  Normal form coefficients (detailed)

• type::Symbol

  Type of bifurcation

Bordered Tangent predictor

x = BorderedArray(vec1, vec2)

This defines an array (although not <: AbstractArray) to hold two arrays or an array and a scalar. This is useful when one wants to add constraints (phase, ...) to a functional for example. It is used throughout the package for the Pseudo Arc Length Continuation, for the continuation of Fold / Hopf points, for periodic orbits... It is also used to define periodic orbits as (orbit, period). As such, it is a convenient alternative to cat, vcat and friends. We chose not to make it a subtype of AbstractArray as we wish to apply the current package to general "arrays", see Requested methods for Custom State. Finally, it proves useful for the GPU where the operation x[end] can be slow.

struct BorderingBLS{S<:Union{Nothing, BifurcationKit.AbstractLinearSolver}, Ttol} <: BifurcationKit.AbstractBorderedLinearSolver

This struct is used to provide the bordered linear solver based on the Bordering Method. user the options, you can trigger a sequence of Bordering reductions to meet a precision.

Reference

This is the solver BEC + k in Govaerts, W. “Stable Solvers and Block Elimination for Bordered Systems.” SIAM Journal on Matrix Analysis and Applications 12, no. 3 (July 1, 1991): 469–83. https://doi.org/10.1137/0612034.

• solver::Union{Nothing, BifurcationKit.AbstractLinearSolver}

  Linear solver used for the Bordering method. Default: nothing

• tol::Any

  Tolerance for checking precision Default: 1.0e-12

• checkPrecision::Bool

  Check precision of the linear solve? Default: true

• k::Int64

  Number of recursions to achieve tolerance Default: 1

Constructors

• there is a simple constructor BorderingBLS(ls) where ls is a linear solver, for example ls = DefaultLS()
• you can use keyword argument to create such solver, for example BorderingBLS(solver = DefaultLS(), tol = 1e-4)

struct Branch{Tkind, Tprob, T<:Union{ContResult, Vector{ContResult}}, Tbp} <: BifurcationKit.AbstractResult{Tkind, Tprob}

A Branch is a structure which encapsulates the result of the computation of a branch bifurcating from a bifurcation point.

• γ::Union{ContResult, Vector{ContResult}}

  Set of branches branching off the bifurcation point bp

• bp::Any

  Bifurcation point. It is thought as the root of the branches in γ

struct ContIterable{Tkind<:BifurcationKit.AbstractContinuationKind, Tprob, Talg, T, S, E, TnormC, Tdot, Tfinalisesolution, TcallbackN, Tevent} <: BifurcationKit.AbstractContinuationIterable{Tkind<:BifurcationKit.AbstractContinuationKind}

Useful functions

• setParam(iter, p) set parameter with lens iter.prob.lens to p
• isEventActive(iter) whether the event detection is active
• computeEigenElements(iter) whether to compute eigen elements
• saveEigenvectors(iter) whether to save eigen vectors
• getParams(iter) get full list of params
• length(iter)
• isInDomain(iter, p) whether p in is domain [pMin, pMax]. (See ContinuationPar)
• isOnBoundary(iter, p) whether p in is {pMin, pMax}

struct 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 (recordFromSolution(u, param), param, itnewton, itlinear, ds, theta, n_unstable, n_imag, stable, step) for each continuation step i.

  • itnewton number of Newton iterations
  • itlinear total number of linear iterations during 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 for each continuation step (to detect Hopf bifurcation).
  • stable stability of the computed solution for each continuation step. Hence, stable should match eig[step] which corresponds to branch[k] for a given k.
  • step continuation step (here equal i)

• eig::Array{NamedTuple{(:eigenvals, :eigenvecs, :step), Tuple{Teigvals, Teigvec, 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 argument saveSolEveryStep::Int64 (default 0) in ContinuationPar.

• contparams::Any

  The parameters used for the call to continuation which produced this branch. Must be a ContinatioPar

• kind::BifurcationKit.AbstractContinuationKind

  Type of solutions computed in this branch. Default: EquilibriumCont()

• prob::Any

  Structure associated to the functional, useful for branch switching. For example, when computing periodic orbits, the functional PeriodicOrbitTrapProblem, ShootingProblem... will be saved here. Default: nothing

• specialpoint::Vector

  A vector holding the set of detected bifurcation points. See SpecialPoint for a description of the fields.

• alg::Any

  Continuation algorithm used for the computation of the branch

Associated methods

• length(br) number of the continuation steps
• eigenvals(br, ind) returns the eigenvalues for the ind-th continuation step
• eigenvec(br, ind, indev) returns the indev-th eigenvector for the ind-th continuation step
• br[k+1] gives information about the k-th step
• getSolx(br, k) returns the k-th solution on the branch
• getSolp(br, k) returns the parameter value associated with k-th solution on the branch
• getParams(br) Parameters passed to continuation and used in the equation F(x, par) = 0.
• setParam(br, p0) set the parameter value p0 according to ::Lens for the parameters of the problem br.prob
• getLens(br) get the lens used for the computation of the branch

state = ContState(ds = 1e-4,...)

Returns a variable containing the state of the continuation procedure. The fields are meant to change during the continuation procedure.

Arguments

• z_pred current solution on the branch
• converged Boolean for newton correction
• τ tangent predictor
• z previous solution
• itnewton Number of newton iteration (in corrector)
• step current continuation step
• ds step size
• theta theta parameter for constraint equation in PALC
• stopcontinuation Boolean to stop continuation

Useful functions

• copy(state) returns a copy of state
• copyto!(dest, state) returns a copy of state
• getSolution(state) returns the current solution (x, p)
• getx(state) returns the x component of the current solution
• getp(state) returns the p component of the current solution
• getpreviousp(state) returns the p component of the previous solution
• isStable(state) whether the current state is stable

options = 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.
• θ is a parameter in the arclength constraint. It is very important to tune it. See the docs of continuation. We quote them here "It should be tuned for the continuation to work properly especially in the case of large problems where the < x - x0, dx0 > component in the constraint equation might be favoured too much. Also, large thetas favour p as the corresponding term in N involves the term 1-theta."
• pMin, pMax allowed parameter range for p
• maxSteps = 100 maximum number of continuation steps
• newtonOptions::NewtonPar: options for the Newton algorithm
• saveToFile = false: save to file. A name is automatically generated or can be defined in continuation. This requires using JLD2.
• saveSolEveryStep::Int64 = 0 at which continuation steps do we save the current solution
• plotEveryStep = 10

Handling eigen elements, their computation is triggered by the argument detectBifurcation (see below)

• nev = 3 number of eigenvalues to be computed. It is automatically increased to have at least nev unstable eigenvalues. To be set for proper bifurcation detection. See Detection of bifurcation points for more informations.
• saveEigEveryStep = 1 record eigen vectors every specified steps. Important for memory limited resource, e.g. GPU.
• saveEigenvectors = true Important for memory limited resource, e.g. GPU.

Handling bifurcation detection

• tolStability = 1e-10 lower bound on the real part of the eigenvalues to test for stability of equilibria and periodic orbits
• detectFold = 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)
• detectBifurcation::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)
• dsminBisection = 1e-16 dsmin for the bisection algorithm for locating bifurcation points
• nInversion = 2 number of sign inversions in bisection algorithm
• maxBisectionSteps = 15 maximum number of bisection steps
• tolBisectionEigenvalue = 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 adapt ds in order to have a number of newton iterations per continuation step roughly constant. The higher a is, the larger the step size ds is changed at each continuation step.

Handling event detection

• detectEvent::Int ∈ {0, 1, 2} If set to 0, nothing is done. If set to 1, the event locations are seek during the continuation run, but not located precisely. If set to 2, a bisection algorithm is used to locate the event (slower).
• tolParamBisectionEvent = 1e-16 tolerance on parameter to locate event

Misc

• η = 150. parameter to estimate tangent at first point with parameter p₀ + ds / η
• finDiffEps = 1e-9 ε used in finite differences computations
• detectLoop [WORK IN PROGRESS] detect loops in the branch and stop the continuation

struct 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 function

• condition::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 type T should match the one of the parameter specified by the ::Lens in continuation.

• computeEigenElements::Bool

  whether the event requires to compute eigen elements

• labels::Any

  Labels used to display information. For example labels[1] is used to qualify an event of the type (0,1.3213,3.434). You can use labels = ("hopf",) or labels = ("hopf", "fold"). You must have labels::Union{Nothing, NTuple{N, String}}.

• tol::Any

  Tolerance on event value to declare it as true event.

mutable struct Cusp{Tv, Tpar, Tlens, Tevr, Tevl, Tnf} <: BifurcationKit.AbstractBifurcationPoint

• x0::Any

  Bifurcation point

• params::Any

  Parameters used by the vector field.

• lens::Any

  Parameter axis used to compute the branch on which this cusp point was detected.

• ζ::Any

  Right eigenvector

• ζstar::Any

  Left eigenvector

• nf::Any

  Normal form coefficients

• type::Symbol

  Type of bifurcation

Wrapper for a distance. You need to pass a function d(u, v).

struct DCResult{Tprob, Tbr, Tit, Tsol, Talg} <: BifurcationKit.AbstractBranchResult

Structure holding the result from deflated continuation.

• prob::Any

  Bifurcation problem

• branches::Any

  Branches of solution

• iter::Any

  Continuation iterator

• sol::Any

  Solutions

• alg::Any

  Algorithm

struct DefCont{Tdo, Talg, Tps, Tas, Tud, Tk} <: BifurcationKit.AbstractContinuationAlgorithm

Structure which holds the parameters specific to Deflated continuation.

Fields

• deflationOperator::Any

  Deflation operator, ::DeflationOperator Default: nothing

• alg::Any

  Used as a predictor, ::AbstractContinuationAlgorithm. For example PALC(), Natural(),... Default: PALC()

• maxBranches::Int64

  maximum number of (active) branches to be computed Default: 100

• seekEveryStep::Int64

  whether to seek new (deflated) solution at every step Default: 1

• maxIterDefOp::Int64

  maximum number of deflated Newton iterations Default: 5

• perturbSolution::Any

  perturb function Default: _perturbSolution

• acceptSolution::Any

  accept (solution) function Default: _acceptSolution

• updateDeflationOp::Any

  function to update the deflation operator Default: _updateDeflationOp

• jacobian::Any

  jacobian for deflated newton. Can be DefProbCustomLinearSolver(), or Val(:autodiff), Val(:fullIterative) Default: DefProbCustomLinearSolver()

struct DefProbCustomLinearSolver{T} <: BifurcationKit.AbstractLinearSolverForDeflation

Custom linear solver for deflated problem, very close to the Sherman-Morrison formula.

Implement the custom linear solver for the deflated problem.

struct DefProbFullIterativeLinearSolver{T} <: BifurcationKit.AbstractLinearSolverForDeflation

Full iterative linear solver for deflated problem.

The struct Default is used to provide the backslash operator to our Package

The struct Default is used to provide the backslash operator to our Package

struct DefaultLS <: BifurcationKit.AbstractDirectLinearSolver

This struct is used to provide the backslash operator. Can be used to solve (a₀ * I + a₁ * J) * x = rhs.

• useFactorization::Bool

  Whether to catch a factorization for multiple solves. Some operators may not support LU (like ApproxFun.jl) or QR factorization so it is best to let the user decides. Some matrices do not have factorize like StaticArrays.MMatrix. Default: true

struct DefaultPILS <: BifurcationKit.AbstractIterativeLinearSolver

Mainly for debugging!! This is defined as an iterative pseudo-inverse linear solver This struct is used to test Moore-Penrose continuation. Used to solve J * x = rhs.

• useFactorization::Bool

  Whether to catch a factorization for multiple solves. Some operators may not support LU (like ApproxFun.jl) or QR factorization so it is best to let the user decides. Some matrices do not have factorize like StaticArrays.MMatrix. Default: true

pb = 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.

struct DeflationOperator{Tp<:Real, Tdot, T<:Real, vectype} <: BifurcationKit.abstractDeflationFactor

Strcuture 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

  power p. You can use an Int for example

• dot::Any

  function, this function has to be bilinear and symmetric for the linear solver to work well

• α::Real

  shift

• roots::Vector

  roots

• tmp::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.

struct 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 function

• condition::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 elements

• labels::Any

  Labels used to display information. For example labels[1] is used to qualify an event occuring in the first component. You can use labels = ("hopf",) or labels = ("hopf", "fold"). You must have labels::Union{Nothing, NTuple{N, String}}.

struct DotTheta{Tdot, Ta}

• dot::Any

  dot product used in pseudo-arclength constraint

• apply!::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 example x -> 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 weigthed 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.

struct EigArnoldiMethod{T, Tby, Tw, Tkw, vectype} <: AbstractEigenSolver

• sigma::Any

  Shift for Shift-Invert method

• which::Any

  Which eigen-element to extract LR(), LM(), ...

• by::Any

  how do we sort the computed eigenvalues, defaults to real

• kwargs::Any

  Key words arguments passed to EigArpack

• x₀::Any

  Example of vector used for Krylov iterations

More information is available at ArnoldiMethod.jl. For example, you can pass the parameters tol, mindim, maxdim, restarts.

Constructor

EigArnoldiMethod(;sigma = nothing, which = ArnoldiMethod.LR(), x₀ = nothing, kwargs...)

struct EigArpack{T, Tby, Tw} <: AbstractEigenSolver

• sigma::Any

  Shift for Shift-Invert method with `(J - sigma⋅I)

• which::Symbol

  Which eigen-element to extract :LR, :LM, ...

• by::Any

  Sorting function, default to real

• kwargs::Any

  Keyword arguments passed to EigArpack

More information is available at Arpack.jl. You can pass the following parameters tol=0.0, maxiter=300, ritzvec=true, v0=zeros((0,)).

Constructor

EigArpack(sigma = nothing, which = :LR; kwargs...)

struct EigKrylovKit{T, vectype} <: BifurcationKit.AbstractMFEigenSolver

• dim::Int64

  Krylov Dimension Default: KrylovDefaults.krylovdim

• tol::Any

  Tolerance Default: 0.0001

• restart::Int64

  Number of restarts Default: 200

• maxiter::Int64

  Maximum number of iterations Default: KrylovDefaults.maxiter

• verbose::Int64

  Verbosity ∈ {0,1,2} Default: 0

• which::Symbol

  Which eigenvalues are looked for :LR (largest real), :LM, ... Default: :LR

• issymmetric::Bool

  If the linear map is symmetric, only meaningful if T<:Real Default: false

• ishermitian::Bool

  If the linear map is hermitian Default: false

• x₀::Any

  Example of vector to usen for Krylov iterations Default: nothing

Singleton type to trigger the computation of the jacobian using finite differences. It can be used for example in newton or in deflated newton.

eigfloquet = BifurcationKit.FloquetColl()

Computation of Floquet coefficients for the orthogonal collocation method. The method is based on the condensation of parameters described in [1] and used in Auto07p with a twist from [2] in which we form the monodromy matrix with a product of Ntst matrices.

This is much faster than FloquetCollGEV but less precise. The better version will use a Periodic Schur decomposition instead of the product of Ntst matrices. This will be provided in the near future.

References

[1] Doedel, Eusebius, Herbert B. Keller, et Jean Pierre Kernevez. « NUMERICAL ANALYSIS AND CONTROL OF BIFURCATION PROBLEMS (II): BIFURCATION IN INFINITE DIMENSIONS ». International Journal of Bifurcation and Chaos 01, nᵒ 04 (décembre 1991): 745‑72. https://doi.org/10.1142/S0218127491000555.

[2] Lust, Kurt. « Improved Numerical Floquet Multipliers ». International Journal of Bifurcation and Chaos 11, nᵒ 09 (septembre 2001): 2389‑2410. https://doi.org/10.1142/S0218127401003486.

Computation of Floquet coefficients for the orthogonal collocation method. The method is based on a formulation through a generalised eigenvalue problem (GEV). Relatively slow but quite precise.

This is a simplified version of [1].

Arguments

• eigls an eigensolver
• ntot total number of unknowns (without counting the period), ie length(::PeriodicOrbitOCollProblem)
• n space dimension

Example

You can create such solver like this (here n=2):

eigfloquet = BifurcationKit.FloquetCollGEV(DefaultEig(), (30*4+1)*2, 2))

References

[1] Fairgrieve, Thomas F., and Allan D. Jepson. “O. K. Floquet Multipliers.” SIAM Journal on Numerical Analysis 28, no. 5 (October 1991): 1446–62. https://doi.org/10.1137/0728075.

floquet = 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.

struct FloquetWrapper{Tpb, Tjacpb, Torbitguess, Tp}

Define a structure to interface the Jacobian of the periodic orbits Problem with the Floquet computation methods. If we use the same code as for newton (see below) but in continuation, it is difficult to tell the eigensolver that it should not use the jacobian but the monodromy matrix instead.

• pb::Any

• jacpb::Any

• x::Any

• par::Any

struct Flow{TF, Tf, Tts, Tff, Td, Tse, Tprob, TprobMono, Tfs, Tcb} <: 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: nothing

• flow::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: nothing

• flowTimeSol::Any

  Flow which returns the tuple (t, u(t)). Optional, mainly used for plotting on the user side. Default: nothing

• flowFull::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 use nothing as default. Default: nothing

• dflow::Any

  The differential dflow of the flow w.r.t. x, (x, p, dx, t) -> dflow(x, p, dx, t). One important thing is that we require dflow(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: nothing

• dfSerial::Any

  [Optional] Serial version of dflow. Used internally when using parallel multiple shooting. Please use nothing as default. Default: nothing

• prob::Any

  [Internal] store the ODEProblem associated to the flow of the Cauchy problem Default: nothing

• probMono::Any

  [Internal] store the ODEProblem associated to the flow of the variational problem Default: nothing

• flowSerial::Any

  [Internal] Serial version of the flow Default: nothing

• callback::Any

  [Internal] Store possible callback Default: nothing

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...)

mutable struct Fold{Tv, T, Tpar, Tlens<:Setfield.Lens, Tevl, Tevr, Tnf} <: BifurcationKit.AbstractSimpleBranchPoint

• x0::Any

  bifurcation point.

• p::Any

  Parameter value at the bifurcation point.

• params::Any

  Parameters used by the vector field.

• lens::Setfield.Lens

  Parameter axis used to compute the branch on which this bifurcation point was detected.

• ζ::Any

  Right eigenvector(s).

• ζstar::Any

  Left eigenvector(s).

• nf::Any

  Normal form coefficients.

• type::Symbol

  Type of bifurcation point

Predictor

You can call predictor(bp, ds; kwargs...) on such bifurcation point bp to find the zeros of the normal form polynomials.

mutable struct FoldProblemMinimallyAugmented{Tprob<:BifurcationKit.AbstractBifurcationProblem, vectype, T<:Real, S<:BifurcationKit.AbstractLinearSolver, Sa<:BifurcationKit.AbstractLinearSolver, Sbd<:BifurcationKit.AbstractBorderedLinearSolver, Sbda<:BifurcationKit.AbstractBorderedLinearSolver} <: BifurcationKit.AbstractProblemMinimallyAugmented

Structure to encode Fold / Hopf functional based on a Minimally Augmented formulation.

Fields

• prob_vf

  Functional F(x, p) - vector field - with all derivatives

• a

  close to null vector of Jᵗ

• b

  close to null vector of J

• zero

  vector zero, to avoid allocating it many times

• l1

  Lyapunov coefficient

• CP

  Cusp test value

• BT

  Bogdanov-Takens test value

• GH

  Bautin test values

• ZH

  Zero-Hopf test values

• linsolver

  linear solver. Used to invert the jacobian of MA functional

• linsolverAdjoint

  linear solver for the jacobian adjoint

• linbdsolver

  bordered linear solver

• linbdsolverAdjoint

  linear bordered solver for the jacobian adjoint

• usehessian

  wether to use the hessian of prob_vf

struct GEigArpack{T, Tby, Tw, Tb} <: BifurcationKit.AbstractGEigenSolver

• sigma::Any

  Shift for Shift-Invert method with `(J - sigma⋅I)

• which::Symbol

  Which eigen-element to extract :LR, :LM, ...

• by::Any

  Sorting function, default to real

• kwargs::Any

  Keyword arguments passed to EigArpack

• B::Any

  Mass matrix

More information is available at Arpack.jl. You can pass the following parameters tol=0.0, maxiter=300, ritzvec=true, v0=zeros((0,)).

Constructor

EigArpack(sigma = nothing, which = :LR; kwargs...)

mutable struct GMRESIterativeSolvers{T, Tl, Tr} <: BifurcationKit.AbstractIterativeLinearSolver

Linear solver based on gmres from IterativeSolvers.jl. Can be used to solve (a₀ * I + a₁ * J) * x = rhs.

• abstol::Any

  Absolute tolerance for solver Default: 0.0

• reltol::Any

  Relative tolerance for solver Default: 1.0e-8

• restart::Int64

  Number of restarts Default: 200

• maxiter::Int64

  Maximum number of iterations Default: 100

• N::Int64

  Dimension of the problem Default: 0

• verbose::Bool

  Display information during iterations Default: false

• log::Bool

  Record information Default: true

• initially_zero::Bool

  Start with zero guess Default: true

• Pl::Any

  Left preconditioner Default: IterativeSolvers.Identity()

• Pr::Any

  Right preconditioner Default: IterativeSolvers.Identity()

• ismutating::Bool

  Whether the linear operator is written inplace Default: false

mutable struct GMRESKrylovKit{T, Tl} <: BifurcationKit.AbstractIterativeLinearSolver

Create a linear solver based on GMRES from KrylovKit.jl. Can be used to solve (a₀ * I + a₁ * J) * x = rhs.

• dim::Int64

  Krylov Dimension Default: KrylovDefaults.krylovdim

• atol::Any

  Absolute tolerance for solver Default: KrylovDefaults.tol

• rtol::Any

  Relative tolerance for solver Default: KrylovDefaults.tol

• maxiter::Int64

  Maximum number of iterations Default: KrylovDefaults.maxiter

• verbose::Int64

  Verbosity ∈ {0,1,2} Default: 0

• issymmetric::Bool

  If the linear map is symmetric, only meaningful if T<:Real Default: false

• ishermitian::Bool

  If the linear map is hermitian Default: false

• isposdef::Bool

  If the linear map is positive definite Default: false

• Pl::Any

  Left preconditioner Default: nothing

mutable struct Hopf{Tv, T, Tω, Tpar, Tlens<:Setfield.Lens, Tevr, Tevl, Tnf} <: BifurcationKit.AbstractSimpleBranchPoint

• x0::Any

  Hopf point

• p::Any

  Parameter value at the Hopf point

• ω::Any

  Frequency of the Hopf point

• params::Any

  Parameters used by the vector field.

• lens::Setfield.Lens

  Parameter axis used to compute the branch on which this Hopf point was detected.

• ζ::Any

  Right eigenvector

• ζstar::Any

  Left eigenvector

• nf::Any

  Normal form coefficient ex: (a = 0., b = 1 + 1im)

• type::Symbol

  Type of Hopf bifurcation

Associated methods

Predictor

You can call predictor(bp::Hopf, ds) on such bifurcation point bp to get the guess for the periodic orbit.

mutable struct HopfHopf{Tv, Tpar, Tlens, Tevr, Tevl, Tnf} <: BifurcationKit.AbstractBifurcationPoint

• x0::Any

  Bifurcation point

• params::Any

  Parameters used by the vector field.

• lens::Any

  Parameter axis used to compute the branch on which this cusp point was detected.

• ζ::Any

  Right eigenvector

• ζstar::Any

  Left eigenvector

• nf::Any

  Normal form coefficients

• type::Symbol

  Type of bifurcation

mutable struct HopfProblemMinimallyAugmented{Tprob<:BifurcationKit.AbstractBifurcationProblem, vectype, T<:Real, S<:BifurcationKit.AbstractLinearSolver, Sa<:BifurcationKit.AbstractLinearSolver, Sbd<:BifurcationKit.AbstractBorderedLinearSolver, Sbda<:BifurcationKit.AbstractBorderedLinearSolver} <: BifurcationKit.AbstractProblemMinimallyAugmented

Structure to encode Fold / Hopf functional based on a Minimally Augmented formulation.

Fields

• prob_vf

  Functional F(x, p) - vector field - with all derivatives

• a

  close to null vector of Jᵗ

• b

  close to null vector of J

• zero

  vector zero, to avoid allocating it many times

• l1

  Lyapunov coefficient

• CP

  Cusp test value

• BT

  Bogdanov-Takens test value

• GH

  Bautin test values

• ZH

  Zero-Hopf test values

• linsolver

  linear solver. Used to invert the jacobian of MA functional

• linsolverAdjoint

  linear solver for the jacobian adjoint

• linbdsolver

  bordered linear solver

• linbdsolverAdjoint

  linear bordered solver for the jacobian adjoint

• usehessian

  wether to use the hessian of prob_vf

struct LSFromBLS{Ts} <: BifurcationKit.AbstractLinearSolver

This structure is used to provide the following linear solver. To solve (1) J⋅x = rhs, one decomposes J using Matrix by blocks and then use a bordering strategy to solve (1).

• solver::Any

  Linear solver used to solve the smaller linear systems.

struct MatrixBLS{S<:Union{Nothing, BifurcationKit.AbstractLinearSolver}} <: BifurcationKit.AbstractBorderedLinearSolver

This struct is used to provide the bordered linear solver based on inverting the full matrix.

• solver::Union{Nothing, BifurcationKit.AbstractLinearSolver}

  Linear solver used to invert the full matrix.

struct MatrixFreeBLS{S<:Union{Nothing, BifurcationKit.AbstractLinearSolver}} <: BifurcationKit.AbstractBorderedLinearSolver

This struct is used to provide the bordered linear solver based a matrix free operator for the full system in (x, p).

• solver::Union{Nothing, BifurcationKit.AbstractLinearSolver}

  Linear solver used to solve the extended linear system

• useBorderedArray::Bool

  What is the structure used to hold (x, p). If true, this is achieved using BorderedArray. If false, a Vector is used.

cache = MeshCollocationCache(Ntst::Int, m::Int, Ty = Float64)

Structure to hold the cache for the collocation method.

• Ntst::Int64

  Coarse mesh size

• degree::Int64

  Collocationn degree, usually called m

• lagrange_vals::Matrix

  Lagrange matrix

• lagrange_∂::Matrix

  Lagrange matrix for derivative

• gauss_nodes::Vector

  Gauss nodes

• gauss_weight::Vector

  Gauss weights

• mesh::Vector

  Values for the coarse mesh, call τj. This can be adapted.

• mesh_coll::LinRange

  Values for collocation poinnts, call σj. These are fixed.

• full_mesh::Vector

  Full mesh containing both the coarse mesh and the collocation points.

Constructor

 MeshCollocationCache(Ntst::Int, m::Int, Ty = Float64)

• Ntst number of time steps
• m degree of the collocation polynomials
• Ty type of the time variable

Moore-Penrose predictor / corrector

Moore-Penrose continuation algorithm.

Additional information is available on the website.

Constructors

alg = MoorePenrose()

alg = MoorePenrose(pred::AbstractPredictor)

• tangent::Any

  Tangent predictor, example PALC()

• method::MoorePenroseLS

  Use a direct linear solver

• ls::BifurcationKit.AbstractLinearSolver

  (Bordered) linear solver

MoorePenrose(; tangent, method, ls)

Multiple Tangent continiation algorithm.

• alg::PALC

  Tangent predictor used Default: PALC()

• τ::Any

  Save the current tangent

• α::Real

  Damping in Newton iterations, 0 < α < 1

• nb::Int64

  Number of predictors

• currentind::Int64

  Index of the largest converged predictor Default: 0

• pmimax::Int64

  Index for lookup in residual history Default: 1

• imax::Int64

  Maximum index for lookup in residual history Default: 4

• dsfact::Real

  Factor to increase ds upon successful step Default: 1.5

Constructor(s)

Multiple(alg, x0, α, n)

Multiple(pred, x0, α, n)

Multiple(x0, α, n)

Natural continuation algorithm.

This is a type which holds information for the bifurcation points of equilibria with dim(Ker)>1.

mutable struct NdBranchPoint{Tv, T, Tpar, Tlens<:Setfield.Lens, Tevl, Tevr, Tnf} <: BifurcationKit.AbstractBranchPoint

• x0::Any

  bifurcation point

• p::Any

  Parameter value at the bifurcation point

• params::Any

  Parameters used by the vector field.

• lens::Setfield.Lens

  Parameter axis used to compute the branch on which this bifurcation point was detected.

• ζ::Any

  Right eigenvectors

• ζstar::Any

  Left eigenvectors

• nf::Any

  Normal form coefficients

• type::Symbol

  Type of bifurcation point

Associated methods

You can call type(bp::NdBranchPoint), length(bp::NdBranchPoint).

Predictor

You can call predictor(bp, ds) on such bifurcation point bp to find the zeros of the normal form polynomials.

Manipulating the normal form

• You can use bp(Val(:reducedForm), x, p) to evaluate the normal form polynomials on the vector x for (scalar) parameter p.

• You can use bp(x, δp::Real) to get the (large dimensional guess) associated to the low dimensional vector x. Note that we must have length(x) == length(bp).

• You can use BifurcationKit.nf(bp; kwargs...) to pretty print the normal form with a string.

struct 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 for F(x) Default: 1.0e-12

• maxIter::Int64

  number of Newton iterations Default: 25

• verbose::Bool

  display Newton iterations? Default: false

• linsolver::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 damping alpha

Structure to hold the solution from application of Newton-Krylov algorithm to a nonlinear problem.

• u::Any

  solution

• prob::Any

  nonlinear problem

• residuals::Any

  sequence of residuals

• converged::Bool

  has algorithm converged?

• itnewton::Int64

  number of newton iterations

• itlineartot::Any

  total number of linear iterations

struct ODEBifProblem{Tvf, Tu, Tp, Tl<:Setfield.Lens, Tplot, Trec} <: BifurcationKit.AbstractAllJetBifProblem

Structure to hold the bifurcation problem. If don't have parameters, you can pass nothing.

Fields

• VF::Any

  Vector field, typically a BifFunction

• u0::Any

  Initial guess

• params::Any

  parameters

• lens::Setfield.Lens

  Typically a Setfield.Lens. It specifies which parameter axis among params is used for continuation. For example, if par = (α = 1.0, β = 1), we can perform continuation w.r.t. α by using lens = (@lens _.α). If you have an array par = [ 1.0, 2.0] and want to perform continuation w.r.t. the first variable, you can use lens = (@lens _[1]). For more information, we refer to SetField.jl.

• plotSolution::Any

  user function to plot solutions during continuation. Signature: plotSolution(x, p; kwargs...)

• recordFromSolution::Any

  recordFromSolution = (x, p) -> norm(x) function used record a few indicators about the solution. It could be norm 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 in contres.branch (see below). Finally, the first component is used to plot in the continuation curve.

Methods

• getu0(pb) calls pb.u0
• getParams(pb) calls pb.params
• getLens(pb) calls pb.lens
• getParam(pb) calls get(pb.params, pb.lens)
• setParam(pb, p0) calls set(pb.params, pb.lens, p0)
• recordFromSolution(pb) calls pb.recordFromSolution
• plotSolution(pb) calls pb.plotSolution
• isSymmetric(pb) calls isSymmetric(pb.prob)

Constructors

• BifurcationProblem(F, u0, params, lens; J, Jᵗ, d2F, d3F, kwargs...) and kwargs are the fields above. You can pass your own jacobian with J (see [BifFunction]@ref for description of the jacobian function) and jacobian adjoint with Jᵗ. For example, this can be used to provide finite differences based jacobian.

struct PALC{Ttang<:BifurcationKit.AbstractTangentComputation, Tbls<:BifurcationKit.AbstractLinearSolver, T} <: BifurcationKit.AbstractContinuationAlgorithm

Pseudo-arclength continuation algorithm.

Additional information is available on the website.

Fields

• tangent::BifurcationKit.AbstractTangentComputation

  Tangent predictor, must be a subtype of AbstractTangentComputation. For example Secant() or Bordered(), Default: Secant()

• bothside::Bool

  [internal], Default: false

• bls::BifurcationKit.AbstractLinearSolver

  Bordered linear solver used to invert the jacobian of the newton bordered problem. It is also used to compute the tangent for the predictor Bordered(), Default: MatrixBLS()

• doArcLengthScaling::Bool

  Unused for now, Default: false

• gGoal::Any

  Unused for now, Default: 0.5

• gMax::Any

  Unused for now, Default: 0.8

• θMin::Any

  Unused for now, Default: 0.001

Associated methods

struct PDEBifProblem{Tvf, Tu, Tp, Tl<:Setfield.Lens, Tplot, Trec} <: BifurcationKit.AbstractAllJetBifProblem

Structure to hold the bifurcation problem. If don't have parameters, you can pass nothing.

Fields

• VF::Any

  Vector field, typically a BifFunction

• u0::Any

  Initial guess

• params::Any

  parameters

• lens::Setfield.Lens

  Typically a Setfield.Lens. It specifies which parameter axis among params is used for continuation. For example, if par = (α = 1.0, β = 1), we can perform continuation w.r.t. α by using lens = (@lens _.α). If you have an array par = [ 1.0, 2.0] and want to perform continuation w.r.t. the first variable, you can use lens = (@lens _[1]). For more information, we refer to SetField.jl.

• plotSolution::Any

  user function to plot solutions during continuation. Signature: plotSolution(x, p; kwargs...)

• recordFromSolution::Any

  recordFromSolution = (x, p) -> norm(x) function used record a few indicators about the solution. It could be norm 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 in contres.branch (see below). Finally, the first component is used to plot in the continuation curve.

Methods

• getu0(pb) calls pb.u0
• getParams(pb) calls pb.params
• getLens(pb) calls pb.lens
• getParam(pb) calls get(pb.params, pb.lens)
• setParam(pb, p0) calls set(pb.params, pb.lens, p0)
• recordFromSolution(pb) calls pb.recordFromSolution
• plotSolution(pb) calls pb.plotSolution
• isSymmetric(pb) calls isSymmetric(pb.prob)

Constructors

• BifurcationProblem(F, u0, params, lens; J, Jᵗ, d2F, d3F, kwargs...) and kwargs are the fields above. You can pass your own jacobian with J (see [BifFunction]@ref for description of the jacobian function) and jacobian adjoint with Jᵗ. For example, this can be used to provide finite differences based jacobian.

cache to remove allocations from PeriodicOrbitOCollProblem

Structure to encode the solution associated to a functional ::PeriodicOrbitOCollProblem. In particular, this allows to use the collocation polynomials to interpolate the solution. Hence, if sol::POOCollSolution, one can call

sol = BifurcationKit.POOCollSolution(prob_coll, x)
sol(t)

on any time t.

struct 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.

• eventC::BifurcationKit.AbstractContinuousEvent

  Continuous event

• eventD::BifurcationKit.AbstractDiscreteEvent

  Discrete event

pb = 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 equation
• xπ::AbstractVector used in the section for the phase constraint equation
• N::Int dimension of the state space
• mesh_cache::MeshCollocationCache cache for collocation. See docs of MeshCollocationCache
• updateSectionEveryStep updates the section every updateSectionEveryStep step during continuation
• jacobian::Symbol symbol which describes the type of jacobian used in Newton iterations. Can only be :autodiffDense.
• meshadapt::Bool = false whether to use mesh adaptation
• verboseMeshAdapt::Bool = true verbose mesh adaptation information
• K::Float64 = 500 parameter for mesh adaptation, control new mesh step size

Methods

Here are some useful methods you can apply to pb

• length(pb) gives the total number of unknowns
• size(pb) returns the triplet (N, m, Ntst)
• getMesh(pb) returns the mesh 0 = τ0 < ... < τNtst+1 = 1. This is useful because this mesh is born to vary by automatic mesh adaptation
• getMeshColl(pb) returns the (static) mesh 0 = σ0 < ... < σm+1 = 1
• getTimes(pb) returns the vector of times (length 1 + m * Ntst) at the which the collocation is applied.
• generateSolution(pb, orbit, period) generate a guess from a function t -> orbit(t) which approximates the periodic orbit.
• POOCollSolution(pb, x) return a function interpolating the solution x using a piecewise polynomials function

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 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 with Ntstand m.

Note that you can generate this guess from a function using generateSolution.

Functional

A functional, hereby called G, encodes this problem. The following methods are available

• pb(orbitguess, p) evaluates the functional G on orbitguess

This 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.

Arguments

• prob a bifurcation problem
• M::Int number of time slices
• ϕ used to set a section for the phase constraint equation, of size N*M
• xπ used in the section for the phase constraint equation, of size N*M
• linsolver: = DefaultLS() linear solver for each time slice, i.e. to solve J⋅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.
• updateSectionEveryStep updates the section every updateSectionEveryStep step during continuation
• jacobian::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.$

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 on orbitguess
• pb(orbitguess, p, du) evaluates the jacobian dG(orbitguess).du functional at orbitguess on du
• pb(Val(:JacFullSparse), orbitguess, p) return the sparse matrix of the jacobian dG(orbitguess) at orbitguess without the constraints. It is called A_γ in the docs.
• pb(Val(:JacFullSparseInplace), J, orbitguess, p). Same as pb(Val(:JacFullSparse), orbitguess, p) but overwrites J inplace. Note that the sparsity pattern must be the same independantly of the values of the parameters or of orbitguess. In this case, this is significantly faster than pb(Val(:JacFullSparse), orbitguess, p).
• pb(Val(:JacCyclicSparse), orbitguess, p) return the sparse cyclic matrix Jc (see the docs) of the jacobian dG(orbitguess) at orbitguess
• pb(Val(:BlockDiagSparse), orbitguess, p) return the diagonal of the sparse matrix of the jacobian dG(orbitguess) at orbitguess. This allows to design Jacobi preconditioner. Use blockdiag.

Jacobian

• jacobian = :FullLU. Specify the choice of the jacobian (and linear algorithm), jacobian must belong to [:FullLU, :FullSparseInplace, :BorderedLU, :FullMatrixFree, :BorderedMatrixFree, :FullSparseInplace]. 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 of dG. 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 matrix dG 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 matrix dG is dense. It is also updated inplace. This option is useful to study ODE of small dimension.
• For jacobian = :BorderedLU, we take advantage of the bordered shape of the linear solver and use a LU decomposition to invert dG using a bordered linear solver.
• For jacobian = :BorderedSparseInplace, this is the same as for :BorderedLU but the cyclic matrix dG 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 for dG: note that a preconditioner is very likely required here because of the cyclic shape of dG which affects negatively the convergence properties of GMRES.
• For jacobian = :BorderedMatrixFree, a matrix free linear solver is used but for Jc only (see docs): it means that options.linsolver is used to invert Jc. These two Matrix-Free options thus expose different part of the jacobian dG in order to use specific preconditioners. For example, an ILU preconditioner on Jc could remove the constraints in dG and lead to poor convergence. Of course, for these last two methods, a preconditioner is likely to be required.
• For jacobian = :FullMatrixFreeAD, the evalution 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.

This method returns the jacobian of the functional G encoded in PeriodicOrbitTrapProblem using a Sparse representation.

mutable struct Pitchfork{Tv, T, Tpar, Tlens<:Setfield.Lens, Tevl, Tevr, Tnf} <: BifurcationKit.AbstractSimpleBranchPoint

• x0::Any

  bifurcation point.

• p::Any

  Parameter value at the bifurcation point.

• params::Any

  Parameters used by the vector field.

• lens::Setfield.Lens

  Parameter axis used to compute the branch on which this bifurcation point was detected.

• ζ::Any

  Right eigenvector(s).

• ζstar::Any

  Left eigenvector(s).

• nf::Any

  Normal form coefficients.

• type::Symbol

  Type of bifurcation point

Predictor

You can call predictor(bp, ds; kwargs...) on such bifurcation point bp to find the zeros of the normal form polynomials.

pb = 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 structure Flow.
• M: the number of Poincaré sections. If M == 1, then the simple shooting is implemented and the multiple one otherwise.
• sections: function or callable struct which implements a Poincaré section condition. The evaluation sections(x) must return a scalar number when M == 1. Otherwise, one must implement a function section(out, x) which populates out with the M sections. See SectionPS for type of section defined as a hyperplane.
• δ = 1e-8 used to compute the jacobian of the functional by finite differences. If set to 0, 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 by EnsembleProblem.
• par parameters of the model
• lens parameter axis
• updateSectionEveryStep updates the section every updateSectionEveryStep step during continuation
• jacobian::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]. This is used to select a way of inverting the jacobian dG
  • For :MatrixFree, we use an iterative solver (e.g. GMRES) to solve the linear system. The jacobian was specified by the user in prob.
  • For :autodiffMF, we use iterative solver (e.g. GMRES) to solve the linear system. We use Automatic Differentiation to compute the (matrix-free) derivative of x -> prob(x, p).
  • For :autodiffDense. Same as for :autodiffMF but the jacobian is formed as a dense Matrix. You can use a direct solver or an iterative one using options.
  • For :FiniteDifferencesDense, same as for :autodiffDense but we use Finite Differences to compute the jacobian of x -> prob(x, p) using the δ = 1e-8 which can be passed as an argument.
  • For :autodiffDenseAnalytical. Same as for :autodiffDense but the jacobian is using a mix of AD and analytical formula.
  • For :FiniteDifferences, use Finite Differences to compute the jacobian of x -> prob(x, p) using the δ = 1e-8 which can be passed as an argument.

Simplified constructors

• A simpler 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 ith 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 generateSolution.

• pb(orbitguess, par) evaluates the functional G on orbitguess
• pb(orbitguess, par, du) evaluates the jacobian dG(orbitguess).du functional at orbitguess on du
• 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.

Polynomial Tangent predictor

• n::Int64

  Order of the polynomial

• k::Int64

  Length of the last solutions vector used for the polynomial fit

• A::Matrix{T} where T<:Real

  Matrix for the interpolation

• tangent::BifurcationKit.AbstractTangentComputation

  Algo for tangent when polynomial predictor is not possible

• solutions::DataStructures.CircularBuffer

  Vector of solutions

• parameters::DataStructures.CircularBuffer{T} where T<:Real

  Vector of parameters

• arclengths::DataStructures.CircularBuffer{T} where T<:Real

  Vector of arclengths

• coeffsSol::Vector

  Coefficients for the polynomials for the solution

• coeffsPar::Vector{T} where T<:Real

  Coefficients for the polynomials for the parameter

• update::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 polynomial
• k length of the last solutions vector used for the polynomial fit
• v0 example of solution to be stored. It is only used to get the eltype of the tangent!!

Secant Tangent predictor

struct SectionPS{Tn, Tc, Tnb, Tcb} <: 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

Constructor(s)

SectionPS(normals::Vector{Tv}, centers::Vector{Tv})

struct 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 hyperplane

• center::Any

  Representative point on hyperplane

struct SetOfEvents{Tc<:Tuple, Td<:Tuple} <: BifurcationKit.AbstractEvent

Multiple events can be chained together to form a SetOfEvents. A SetOfEvents is constructed by passing 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 event

• eventD::Tuple

  Discrete event

pb = ShootingProblem(flow::Flow, ds, section; parallel = false)

Create a problem to implement the Standard 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 structure Flow.
• ds: vector of time differences for each shooting. Its length is written M. If M == 1, then the simple shooting is implemented and the multiple one otherwise.
• section: implements a phase condition. The evaluation section(x, T) must return a scalar number where x is a guess for one point on the periodic orbit and T is the period of the guess. Also, the method section(x, T, dx, dT) must be available and which returns the differential of section. The type of x depends on what is passed to the newton solver. See SectionSS 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 by EnsembleProblem (this is automatically set up for you).
• par parameters of the model
• lens parameter axis
• updateSectionEveryStep updates the section every updateSectionEveryStep step during continuation
• jacobian::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 on orbitguess
• pb(orbitguess, par, du; δ = 1e-9) evaluates the jacobian dG(orbitguess)⋅du functional at orbitguess on du. The optional argument δ is used to compute a finite difference approximation of the derivative of the section.
• 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 AbstractVectors. Another accepted guess is of the form BorderedArray(guess, T) where guess[i] is the state of the orbit at the ith 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 generateSolution.

Jacobian

• jacobian Specify the choice of the linear algorithm, which must belong to [:autodiffMF, :MatrixFree, :autodiffDense, :autodiffDenseAnalytical, :FiniteDifferences]. This is used to select a way of inverting the jacobian dG
  • For :MatrixFree, we use an iterative solver (e.g. GMRES) to solve the linear system. The jacobian was specified by the user in prob.
  • For :autodiffMF, we use iterative solver (e.g. GMRES) to solve the linear system. We use Automatic Differentiation to compute the (matrix-free) derivative of x -> prob(x, p).
  • For :autodiffDense. Same as for :autodiffMF but the jacobian is formed as a dense Matrix. You can use a direct solver or an iterative one using options.
  • For :FiniteDifferencesDense, same as for :autodiffDense but we use Finite Differences to compute the jacobian of x -> prob(x, p) using the δ = 1e-8 which can be passed as an argument.
  • For :autodiffDenseAnalytical. Same as for :autodiffDense but the jacobian is using a mix of AD and analytical formula.
  • For :FiniteDifferences, use Finite Differences to compute the jacobian of x -> prob(x, p) using the δ = 1e-8 which can be passed as an argument.

Simplified constructors

• A simpler 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 ODEProblems. 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)$.

struct SpecialPoint{T, Tp, Tv, Tvτ} <: BifurcationKit.AbstractBifurcationPoint

Structure to record a generic special (bifurcation) point.

• type::Symbol

  Bifurcation type, :hopf, :bp.... Default: :none

• idx::Int64

  Index in br.eig (see ContResult) for which the bifurcation occurs. Default: 0

• param::Any

  Parameter value at the special (bifurcation) point, this is an estimate. Default: 0.0

• norm::Any

  Norm of the equilibrium at the special (bifurcation) point Default: 0.0

• printsol::Any

  printsol = recordFromSolution(x, param) where recordFromSolution is one of the arguments to continuation Default: 0.0

• x::Any

  Equilibrium at the special (bifurcation) point Default: Vector{T}(undef, 0)

• τ::BorderedArray{Tvτ, T} where {T, Tvτ}

  Tangent along the branch at the special (bifurcation) point Default: BorderedArray(x, T(0))

• ind_ev::Int64

  Eigenvalue index responsible for the special (bifurcation) (if applicable) Default: 0

• step::Int64

  Continuation step at which the special (bifurcation) occurs Default: 0

• status::Symbol

  status ∈ {:converged, :guess} indicates whether the bisection algorithm was successful in detecting the special (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 (bifurcation) point Default: -1

• interval::Tuple{T, T} where T

  Interval containing the special (bifurcation) point Default: (0, 0)

TWProblem(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 implementing LinearAlgebra.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 using u0.
• nbConstraints(::TWProblem) number of constraints (or Lie generators)

mutable struct Transcritical{Tv, T, Tpar, Tlens<:Setfield.Lens, Tevl, Tevr, Tnf} <: BifurcationKit.AbstractSimpleBranchPoint

• x0::Any

  bifurcation point.

• p::Any

  Parameter value at the bifurcation point.

• params::Any

  Parameters used by the vector field.

• lens::Setfield.Lens

  Parameter axis used to compute the branch on which this bifurcation point was detected.

• ζ::Any

  Right eigenvector(s).

• ζstar::Any

  Left eigenvector(s).

• nf::Any

  Normal form coefficients.

• type::Symbol

  Type of bifurcation point

Predictor

You can call predictor(bp, ds; kwargs...) on such bifurcation point bp to find the zeros of the normal form polynomials.

struct TrilinearMap{Tm}

This structure wraps a trilinear map to allow evaluation on Complex arguments. This is especially useful when these maps are produced by ForwardDiff.jl.

mutable struct ZeroHopf{Tv, Tpar, Tlens, Tevr, Tevl, Tnf} <: BifurcationKit.AbstractBifurcationPoint

• x0::Any

  Bifurcation point

• params::Any

  Parameters used by the vector field.

• lens::Any

  Parameter axis used to compute the branch on which this cusp point was detected.

• ζ::Any

  Right eigenvector

• ζstar::Any

  Left eigenvector

• nf::Any

  Normal form coefficients

• type::Symbol

  Type of bifurcation

cb = cbMaxNorm(maxres)

Create a callback used to reject residals larger than cb.maxres in the Newton iterations. See docs for newton.

size(tree)
size(tree, code)

Return the size of the bifurcation diagram. The argument code is the same as in getBranch.

The method size returns (n, m, Ntst) when applied to a PeriodicOrbitOCollProblem

Function to compute the Matrix-Free version of Aγ, see docs for its expression.

For an initial guess from the index of a Fold bifurcation point located in ContResult.specialpoint, returns a point which will be refined using newtonFold.

For an initial guess from the index of a Fold bifurcation point located in ContResult.specialpoint, returns a point which will be refined using newtonFold.

For an initial guess from the index of a Hopf bifurcation point located in ContResult.specialpoint, returns a point which will be refined using newtonHopf.

Function to compute the Matrix-Free version of the cyclic matrix Jc, see docs for its expression.

This function implements the functional for finding periodic orbits based on finite differences using the Trapezoidal rule. It works for inplace / out of place vector fields pb.F

Matrix free expression of the Jacobian of the problem for computing periodic obits when evaluated at u and applied to du.

PrecPartialSchurArnoldiMethod(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.

PrecPartialSchurKrylovKit(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().

`SaveAtEvent(positions::Tuple)`

This event implements the detection of when the parameter values, used during continuation, equals one of the values in positions. This state is then saved in the branch.

For example, you can use it like continuation(args...; event = SaveAtEvent((1., 2., -3.)))

_ContResult(
    prob,
    alg,
    printsol,
    br,
    x0,
    τ,
    eiginfo,
    contParams,
    computeEigElements,
    kind
)

Function is used to initialize the composite type ContResult according to the options contained in contParams

Arguments

• br result from getStateSummary
• par: parameters
• lens: lens to specify the continuation parameter
• eiginfo: eigen-elements (eigvals, eigvecs)

This function returns a₀ * I + a₁ * J and ensures that we don't perform unnecessary computations like 0I + 1J.

This function implements the operator a₀ * I + a₁ * J and ensures that we don't perform unnecessary computations like 0I + 1J.

_cat!(br, br2)

Merge two ContResults and put the result in br.

Same as _cat! but determine the ordering so that the branches merge properly

• ss tuple of speeds
• D tuple of Lie generators

bifurcationdiagram!(
    prob,
    node,
    maxlevel,
    options;
    code,
    usedeflation,
    halfbranch,
    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 = getBranch(diagram, code) from a previously computed bifurcation diagram.

Arguments

• node::BifDiagNode a node in the bifurcation diagram
• maxlevel = 1 required maximal level of recursion.
• options = (x, p, level) -> contparams this function allows to change the continuation options depending on the branching level. The argument x, p denotes the current solution to F(x, p)=0.

Optional arguments

• code = "0" code used to display iterations
• usedeflation = false
• halfbranch = false for Pitchfork/Transcritical bifurcations, compute only half of the branch. Can be useful when there are symmetries.
• kwargs optional arguments as for continuation but also for the different versions listed in Continuation.

bifurcationdiagram(prob, alg, level, options; kwargs...)

Compute the bifurcation diagram associated with the problem F(x, p) = 0 recursively.

Arguments

• prob::AbstractBifurcationProblem bifurcation problem
• alg continuation algorithm
• level maximum branching (or recursion) level for computing the bifurcation diagram
• options = (x, p, level) -> contparams this function allows to change the continuation options depending on the branching level. The argument x, p denotes the current solution to F(x, p)=0.
• kwargs optional arguments. Look at bifurcationdiagram! 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.

biorthogonalise(ζs, ζstars, verbose)

Bi-orthogonalise the arguments.

bogdanovTakensNormalForm(
    prob_ma,
    L,
    pt;
    δ,
    verbose,
    detailed,
    autodiff,
    bls
)

Compute the Bogdanov-Takens normal form.

Arguments

• prob_ma a FoldProblemMinimallyAugmented or HopfProblemMinimallyAugmented
• pt::BogdanovTakens BogdanovTakens bifurcation point
• ls linear solver

Optional arguments

• δ = 1e-8 used for finite differences
• verbose bool to print information
• autodiff = true only for Bogdanov-Takens point. Whether to use ForwardDiff for the many differentiations that are required to compute the normal form.
• detailed = true only for Bogdanov-Takens point. Whether to compute only a simplified normal form.

bogdanovTakensNormalForm(
    _prob,
    br,
    ind_bif;
    δ,
    nev,
    verbose,
    ζs,
    lens,
    Teigvec,
    scaleζ,
    bls,
    detailed,
    autodiff
)

Compute the Bogdanov-Takens normal form.

Arguments

• prob bifurcation problem, typically br.prob
• br branch result from a call to continuation
• ind_bif index of the bifurcation point in br
• options options for the Newton solver

Optional arguments

• δ = 1e-8 used for finite differences
• nev = 5 number of eigenvalues to compute to estimate the spectral projector
• verbose bool to print information
• autodiff = true only for Bogdanov-Takens point. Whether to use ForwardDiff for the many differentiations that are required to compute the normal form.
• detailed = true only for Bogdanov-Takens point. Whether to compute only a simplified normal form.

computeError(pb, x; normE, verbosity, K, kw...)

Ascher, Uri M., Robert M. M. Mattheij, and Robert D. Russell. Numerical Solution of Boundary Value Problems for Ordinary Differential Equations. Society for Industrial and Applied Mathematics, 1995. https://doi.org/10.1137/1.9781611971231.

p. 368

continuation(br, ind_bif, lens2)
continuation(
    br,
    ind_bif,
    lens2,
    options_cont;
    startWithEigen,
    detectCodim2Bifurcation,
    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 continuation
• ind_bif bifurcation index in br
• lens2 second parameter used for the continuation, the first one is the one used to compute br, 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 equation
• updateMinAugEveryStep update vectors a,b in Minimally Formulation every updateMinAugEveryStep steps
• startWithEigen = false whether to start the Minimally Augmented problem with information from eigen elements
• detectCodim2Bifurcation ∈ {0,1,2} whether to detect Bogdanov-Takens, Bautin and Cusp. If equals 1 non precise detection is used. If equals 2, 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.

continuation(
    probPO,
    orbitguess,
    alg,
    contParams,
    linearAlgo;
    δ,
    eigsolver,
    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 to FloquetQaD
• jacobian Specify the choice of the linear algorithm, which must belong to [:autodiffMF, :MatrixFree, :autodiffDense, :autodiffDenseAnalytical, :FiniteDifferences]. This is used to select a way of inverting the jacobian dG
  • For :MatrixFree, we use an iterative solver (e.g. GMRES) to solve the linear system. The jacobian was specified by the user in prob.
  • For :autodiffMF, we use iterative solver (e.g. GMRES) to solve the linear system. We use Automatic Differentiation to compute the (matrix-free) derivative of x -> prob(x, p).
  • For :autodiffDense. Same as for :autodiffMF but the jacobian is formed as a dense Matrix. You can use a direct solver or an iterative one using options.
  • For :FiniteDifferencesDense, same as for :autodiffDense but we use Finite Differences to compute the jacobian of x -> prob(x, p) using the δ = 1e-8 which can be passed as an argument.
  • For :autodiffDenseAnalytical. Same as for :autodiffDense but the jacobian is using a mix of AD and analytical formula.
  • For :FiniteDifferences, use Finite Differences to compute the jacobian of x -> prob(x, p) using the δ = 1e-8 which can be passed as an argument.

continuation(
    prob,
    alg,
    contParams;
    linearAlgo,
    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 a BifurcationProblem which holds the vector field and its jacobian. We also refer to BifFunction for more details.
• alg continuation algorithm, for example Natural(), PALC(), Multiple(),.... See algos
• contParams parameters for continuation. See ContinuationPar

Optional Arguments:

• plot = false whether to plot the solution/branch/spectrum while computing
• bothside = true compute the branches on the two sides of p0, merge them and return it.
• finaliseSolution = (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 returning false), saving personal data, plotting... The notations are $z=(x, p)$, tau is the tangent at z, step is the index of the current continuation step and ContResult is the current branch. For advanced use, the current state::ContState of the continuation is passed in kwargs. Note that you can have a better control over the continuation procedure by using an iterator, see Iterator Interface.
• verbosity::Int = 0 controls the amount of information printed during the continuation process. Must belong to {0,1,2,3}. In case contParams.newtonOptions.verbose = false, the following is valid (Otherwise the newton iterations are shown). Each case prints more information then 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 bifurcation / events
• normC = norm norm used in the Newton solves
• dotPALC = DotTheta(), this sets up a dot product (x, y) -> dot(x, y) / length(x) 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 factor 1/length(x) for example in problems where the dimension of the state space changes (mesh adaptation, ...)
• filename name of a file to save the computed branch during continuation. The identifier .jld2 will be appended to this filename. This requires using JLD2.
• callbackN callback for newton iterations. See docs for newton. For example, it can be used to change preconditioners.
• kind::AbstractContinuationKind [Internal] flag to describe continuation kind (equilibrium, codim 2, ...). Default = EquilibriumCont()

Output:

• contres::ContResult composite type which contains the computed branch. See ContResult for more information.

continuation(
    prob,
    algdc,
    contParams;
    verbosity,
    plot,
    linearAlgo,
    315,
    dotPALC,
    callbackN,
    filename,
    normN
)

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 problem
• alg::DefCont, deflated continuation algorithm, see DefCont
• contParams parameters for continuation. See ContinuationPar for more information about the options

Optional Arguments:

• plot = false whether to plot the solution while computing,
• callbackN callback for newton iterations. see docs for newton. 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 argument fromDeflatedNewton = 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},
• normN = norm norm used in the Newton solves,
• dotPALC = (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 factor 1/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. See ContResult for more information,

continuation(
    br,
    ind_bif,
    _contParams,
    probPO;
    alg,
    δp,
    ampfactor,
    usedeflation,
    nev,
    kwargs...
)

Perform automatic branch switching from a Hopf bifurcation point labelled ind_bif in the list of the bifurcated points of a previously computed branch br::ContResult. It first computes a Hopf normal form.

Arguments

• br branch result from a call to continuation
• ind_hopf index of the bifurcation point in br
• contParams parameters for the call to continuation
• probPO problem used to specify the way the periodic orbit is computed. It can be PeriodicOrbitTrapProblem, ShootingProblem or PoincareShootingProblem .

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 by contParams.ds. This allows to use a step larger than contParams.dsmax.
• ampfactor = 1 factor which alter the amplitude of the bifurcated solution. Useful to magnify the bifurcated solution when the bifurcated branch is very steep.
• usedeflation = true whether to use nonlinear deflation (see Deflated problems) to help finding the guess on the bifurcated branch
• nev number of eigenvalues to be computed to get the right eigenvector
• all kwargs from continuation

A modified version of prob is passed to plotSolution and finaliseSolution.

continuation(
    prob,
    orbitguess,
    alg,
    _contParams;
    linearAlgo,
    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

• linearAlgo::AbstractBorderedLinearSolver
• jacobian Specify the choice of the linear algorithm, which must belong to [:autodiffMF, :MatrixFree, :autodiffDense, :autodiffDenseAnalytical, :FiniteDifferences]. This is used to select a way of inverting the jacobian dG
  • For :MatrixFree, we use an iterative solver (e.g. GMRES) to solve the linear system. The jacobian was specified by the user in prob.
  • For :autodiffMF, we use iterative solver (e.g. GMRES) to solve the linear system. We use Automatic Differentiation to compute the (matrix-free) derivative of x -> prob(x, p).
  • For :autodiffDense. Same as for :autodiffMF but the jacobian is formed as a dense Matrix. You can use a direct solver or an iterative one using options.
  • For :FiniteDifferencesDense, same as for :autodiffDense but we use Finite Differences to compute the jacobian of x -> prob(x, p) using the δ = 1e-8 which can be passed as an argument.
  • For :autodiffDenseAnalytical. Same as for :autodiffDense but the jacobian is using a mix of AD and analytical formula.
  • For :FiniteDifferences, use Finite Differences to compute the jacobian of x -> prob(x, p) using the δ = 1e-8 which can be passed as an argument.

continuation(
    probPO,
    orbitguess,
    alg,
    _contParams,
    linearAlgo;
    eigsolver,
    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 to FloquetQaD

continuation(
    prob,
    orbitguess,
    alg,
    _contParams;
    recordFromSolution,
    linearAlgo,
    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 G
• orbitguess a guess for the periodic orbit where orbitguess[end] is an estimate of the period of the orbit. It could be a vector of size N * M + 1 where M is the number of time slices, N is the dimension of the phase space. This must be compatible with the numbers N, M in prob.
• alg continuation algorithm
• contParams same as for the regular continuation method

Keyword arguments

• linearAlgo same as in continuation
• jacobian = :FullLU. Specify the choice of the jacobian (and linear algorithm), jacobian must belong to [:FullLU, :FullSparseInplace, :BorderedLU, :FullMatrixFree, :BorderedMatrixFree, :FullSparseInplace]. 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 of dG. 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 matrix dG 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 matrix dG is dense. It is also updated inplace. This option is useful to study ODE of small dimension.
• For jacobian = :BorderedLU, we take advantage of the bordered shape of the linear solver and use a LU decomposition to invert dG using a bordered linear solver.
• For jacobian = :BorderedSparseInplace, this is the same as for :BorderedLU but the cyclic matrix dG 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 for dG: note that a preconditioner is very likely required here because of the cyclic shape of dG which affects negatively the convergence properties of GMRES.
• For jacobian = :BorderedMatrixFree, a matrix free linear solver is used but for Jc only (see docs): it means that options.linsolver is used to invert Jc. These two Matrix-Free options thus expose different part of the jacobian dG in order to use specific preconditioners. For example, an ILU preconditioner on Jc could remove the constraints in dG and lead to poor convergence. Of course, for these last two methods, a preconditioner is likely to be required.
• For jacobian = :FullMatrixFreeAD, the evalution 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 recordFromSolution argument.

continuation(br, ind_bif)
continuation(
    br,
    ind_bif,
    optionsCont;
    alg,
    δp,
    ampfactor,
    nev,
    usedeflation,
    Teigvec,
    scaleζ,
    verbosedeflation,
    maxIterDeflation,
    perturb,
    plotSolution,
    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 to continuation
• ind_bif index of the bifurcation point in br from which you want to branch from
• optionsCont options for the call to continuation

Optional arguments

• alg = br.alg continuation algorithm to be used, default value: br.alg
• δp used to specify a particular guess for the parameter on the bifurcated branch which is otherwise determined by optionsCont.ds. This allows to use a step larger than optionsCont.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.
• nev number of eigenvalues to be computed to get the right eigenvector
• usedeflation = true whether to use nonlinear deflation (see Deflated problems) to help finding the guess on the bifurcated
• plotSolution change plot solution method in the problem br.prob
• kwargs optional arguments to be passed to continuation, the regular continuation one.

continuation(
    br,
    ind_bif,
    _contParams;
    δp,
    ampfactor,
    usedeflation,
    linearAlgo,
    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 a PeriodicOrbitTrapProblem
• ind_bif index of the branch point
• _contParams parameters to be used by a regular continuation

Optional arguments

• δp = 0.1 used to specify a particular guess for the parameter in the branch which is otherwise determined by contParams.ds. This allows to use a step larger than contParams.dsmax.
• ampfactor = 1 factor which alter the amplitude of the bifurcated solution. Useful to magnify the bifurcated solution when the bifurcated branch is very steep.
• usedeflation = true whether to use nonlinear deflation (see Deflated problems) to help finding the guess on the bifurcated branch
• recordFromSolution = (u, p) -> u[end], record method used in the bifurcation diagram, by default this records the period of the periodic orbit.
• linearAlgo = BorderingBLS(), same as for continuation
• kwargs keywords arguments used for a call to the regular continuation and the ones specific to periodic orbits (POs).

continuation(
    prob,
    x0,
    par0,
    x1,
    p1,
    alg,
    lens,
    contParams;
    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).

continuationFold(
    prob,
    alg,
    foldpointguess,
    par,
    lens1,
    lens2,
    eigenvec,
    eigenvec_ad,
    options_cont;
    normC,
    updateMinAugEveryStep,
    bdlinsolver,
    jacobian_ma,
    computeEigenElements,
    usehessian,
    kwargs...
)

Codim 2 continuation of Fold points. This function turns an initial guess for a Fold point into a curve of Fold points based on a Minimally Augmented formulation. The arguments are as follows

• prob::AbstractBifurcationFunction
• foldpointguess initial guess (x0, p10) for the Fold point. It should be a BorderedArray as returned by the function FoldPoint
• par set of parameters
• lens1 parameter axis for parameter 1
• lens2 parameter axis for parameter 2
• eigenvec guess for the 0 eigenvector at p1_0
• eigenvec_ad guess for the 0 adjoint eigenvector
• options_cont arguments to be passed to the regular continuation

Optional arguments:

• bdlinsolver bordered linear solver for the constraint equation
• updateMinAugEveryStep update vectors a, b in Minimally Formulation every updateMinAugEveryStep steps
• computeEigenElements = false whether to compute eigenelements. If options_cont.detecttEvent>0, it allows the detection of ZH points.
• kwargs keywords arguments to be passed to the regular continuation

Simplified call

continuationFold(br::AbstractBranchResult, ind_fold::Int64, lens2::Lens, options_cont::ContinuationPar ; kwargs...)

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 Fold point in br that you want to continue.

codim 2 continuation of Hopf points. This function turns an initial guess for a Hopf point into a curve of Hopf points based on a Minimally Augmented formulation. The arguments are as follows

• prob::AbstractBifurcationProblem
• hopfpointguess initial guess (x0, p10) for the Hopf point. It should be a Vector or a BorderedArray
• par set of parameters
• lens1 parameter axis for parameter 1
• lens2 parameter axis for parameter 2
• eigenvec guess for the iω eigenvector at p1_0
• eigenvec_ad guess for the -iω eigenvector at p1_0
• options_cont keywords arguments to be passed to the regular continuation

Optional arguments:

• bdlinsolver bordered linear solver for the constraint equation
• updateMinAugEveryStep update vectors a,b in Minimally Formulation every updateMinAugEveryStep steps
• computeEigenElements = false whether to compute eigenelements. If options_cont.detecttEvent>0, it allows the detection of ZH, HH points.
• kwargs keywords arguments to be passed to the regular continuation

Simplified call:

continuationHopf(br::AbstractBranchResult, ind_hopf::Int, lens2::Lens, options_cont::ContinuationPar ;  kwargs...)

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 that you want to refine.

continuationPOTrap(
    prob,
    orbitguess,
    alg,
    contParams,
    linearAlgo;
    eigsolver,
    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 G
• orbitguess a guess for the periodic orbit where orbitguess[end] is an estimate of the period of the orbit. It could be a vector of size N * M + 1 where M is the number of time slices, N is the dimension of the phase space. This must be compatible with the numbers N, M in prob.
• alg conntinuation algorithm
• contParams same as for the regular continuation method
• linearAlgo same as in continuation

Keywords arguments

• eigsolver specify an eigen solver for the computation of the Floquet exponents, defaults to FloquetQaD

• jacobian = :FullLU. Specify the choice of the jacobian (and linear algorithm), jacobian must belong to [:FullLU, :FullSparseInplace, :BorderedLU, :FullMatrixFree, :BorderedMatrixFree, :FullSparseInplace]. 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 of dG. 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 matrix dG 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 matrix dG is dense. It is also updated inplace. This option is useful to study ODE of small dimension.

• For jacobian = :BorderedLU, we take advantage of the bordered shape of the linear solver and use a LU decomposition to invert dG using a bordered linear solver.

• For jacobian = :BorderedSparseInplace, this is the same as for :BorderedLU but the cyclic matrix dG 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 for dG: note that a preconditioner is very likely required here because of the cyclic shape of dG which affects negatively the convergence properties of GMRES.

• For jacobian = :BorderedMatrixFree, a matrix free linear solver is used but for Jc only (see docs): it means that options.linsolver is used to invert Jc. These two Matrix-Free options thus expose different part of the jacobian dG in order to use specific preconditioners. For example, an ILU preconditioner on Jc could remove the constraints in dG and lead to poor convergence. Of course, for these last two methods, a preconditioner is likely to be required.

• For jacobian = :FullMatrixFreeAD, the evalution 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 recordFromSolution argument.

correctBifurcation(contres)

This function uses information in the branch to detect codim 2 bifurcations like BT, ZH and Cusp.

cuspNormalForm(
    _prob,
    br,
    ind_bif;
    δ,
    nev,
    verbose,
    ζs,
    lens,
    Teigvec,
    scaleζ
)

Compute the Cusp normal form.

Arguments

• prob bifurcation problem
• pt::Cusp Cusp bifurcation point
• ls linear solver

Optional arguments

• δ = 1e-8 used for finite differences
• verbose bool to print information

This function populates Jc with the cyclic matrix using the different Jacobians

detectLoop(br, x, p; rtol, verbose)

Function to detect continuation branches which loop on themselve.

This function computes the derivative of the Poincare return map Π(x) = ϕ(t(x),x) where t(x) is the return time of x to the section.

eigenvals(br, ind)
eigenvals(br, ind, verbose)

Return the eigenvalues of the ind-th continuation step. verbose is used to tell the number of unstable eigen elements.

eigenvalsfrombif(br, ind)

Return the eigenvalues of the ind-th bifurcation point.

eigenvec(br, ind, indev)

Return the indev-th eigenvectors of the ind-th continuation step.

finiteDifferences!(F, J, x::AbstractVector; δ = 1e-9)

Compute a Jacobian by Finite Differences, update J. Use the centered formula (f(x+δ)-f(x-δ))/2δ

finiteDifferences(F, x::AbstractVector; δ = 1e-9)

Compute a Jacobian by Finite Differences

from(br)

Return the bifurcation point of a ::Branch.

generateHomoclinicSolution(pb, orbit, T)

This function generates an initial guess for the solution of the problem pb based on the orbit t -> orbit(t) for t ∈ [0,1] and half time return T.

generateSolution(pb, orbit, period)

This function generates an initial guess for the solution of the problem pb based on the orbit t -> orbit(t) for t ∈ [0,1] and the period period.

generateSolution(pb, orbit, period)

This function generates an initial guess for the solution of the problem pb based on the orbit t -> orbit(t) for t ∈ [0,1] and the period.

getAdjointBasis(Lstar, λ, eigsolver; nev, verbose)

Return a left eigenvector for an eigenvalue closest to λ. nev indicates how many eigenvalues must be computed by the eigensolver. Indeed, for iterative solvers, it may be needed to compute more eigenvalues than necessary.

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].

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) = 1 + ratio * n, the call returns the amplitude over x[1:n].

Function for coarse detection of bifurcation points.

getBlocks(A, N, M)

This function extracts the indices of the blocks composing the matrix A which is a M x M Block matrix where each block N x N has the same sparsity.

getBranch(tree, code)

Return the part of the tree (bifurcation diagram) by recursively descending down the tree using the Int valued tuple code. For example getBranch(tree, (1,2,3,)) returns tree.child[1].child[2].child[3].

getBranchesFromBP(tree, indbif)

Return the part of the tree corresponding to the indbif-th bifurcation point on the root branch.

getFirstPointsOnBranch(br, bpnf, solfromRE)
getFirstPointsOnBranch(
    br,
    bpnf,
    solfromRE,
    optionsCont;
    δp,
    Teigvec,
    verbosedeflation,
    maxIterDeflation,
    lsdefop,
    perturbGuess,
    kwargs...
)

Function to transform predictors solfromRE in the normal form coordinates of bpnf into solutions. Note that solfromRE = (before = Vector{vectype}, after = Vector{vectype}).

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].

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) = 1 + ratio * n, the call returns the amplitude over x[1:n].

Returns the vector of size m+1, 0 = τ1 < τ1 < ... < τm+1 = 1

getNormalForm(
    prob,
    br,
    id_bif;
    nev,
    verbose,
    ζs,
    lens,
    Teigvec,
    scaleζ,
    detailed,
    autodiff
)

Compute the normal form of the bifurcation point located at br.specialpoint[ind_bif].

Arguments

• prob::AbstractBifurcationProblem
• br result from a call to continuation
• ind_bif index of the bifurcation point in br.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 of dF at the bifurcation point. Useful to enforce the basis for the normal form.
• lens::Lens specify which parameter to take the partial derivative ∂pF
• scaleζ function to normalise the kernel basis. Indeed, when used with large vectors and norm, it results in ζs and the normal form coefficient being super small.
• autodiff = true only for Bogdanov-Takens point. Whether to use ForwardDiff for the many differentiations that are required to compute the normal form.
• detailed = true only for Bogdanov-Takens point. Whether to compute only a simplified normal form.

Based on Golubitsky, Martin, David G Schaeffer, and Ian Stewart. Singularities and Groups in Bifurcation Theory. New York: Springer-Verlag, 1985, VI.1.d page 295.

Available method

Once the normal form nf has been computed, you can call predictor(nf, δp) to obtain an estimate of the bifurcating periodic orbit.v

getNormalForm1d(
    prob,
    br,
    ind_bif;
    nev,
    verbose,
    lens,
    Teigvec,
    tolFold,
    scaleζ
)

Compute a normal form based on Golubitsky, Martin, David G Schaeffer, and Ian Stewart. Singularities and Groups in Bifurcation Theory. New York: Springer-Verlag, 1985, VI.1.d page 295.

getPeriod(, 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.

getPeriodicOrbit(prob, u, p)

Compute the full periodic orbit associated to x. Mainly for plotting purposes.

getPeriodicOrbit(prob, u, p)

Compute the full periodic orbit associated to x. Mainly for plotting purposes.

getPeriodicOrbit(prob, x_bar, p)

Compute the full periodic orbit associated to x. Mainly for plotting purposes.

getPeriodicOrbit(prob, x, par)

Compute the full periodic orbit associated to x. Mainly for plotting purposes.

getTimeDiff(pb, u)

Compute norm(du/dt)

getTimes(pb)

Return all the times at which the problem is evaluated.

guessFromHopf(br, ind_hopf, eigsolver, M, amplitude; phase)

This function returns several useful quantities regarding a Hopf bifurcation point. More precisely, it returns:

• the parameter value at which a Hopf bifurcation occurs
• the period of the bifurcated periodic orbit
• a guess for the bifurcated periodic orbit
• the equilibrium at the Hopf bifurcation point
• the eigenvector at the Hopf bifurcation point.

The arguments are

• br: the continuation branch which lists the Hopf bifurcation points
• ind_hopf: index of the bifurcation branch, as in br.specialpoint
• eigsolver: the eigen solver used to find the eigenvectors
• M number of time slices in the periodic orbit guess
• amplitude: amplitude of the periodic orbit guess

This function solves the linear problem associated with a linearization of the minimally augmented formulation of the Hopf bifurcation point. The keyword debugArray is used to debug the routine by returning several key quantities.

hopfNormalForm(
    prob,
    br,
    ind_hopf;
    nev,
    verbose,
    lens,
    Teigvec,
    scaleζ
)

Compute the Hopf normal form.

Arguments

• prob::AbstractBifurcationProblem bifurcation problem
• br branch result from a call to continuation
• ind_hopf index of the bifurcation point in br
• options options for the Newton solver

Optional arguments

• nev = 5 number of eigenvalues to compute to estimate the spectral projector
• verbose bool to print information

Available method

Once the normal form hopfnf has been computed, you can call predictor(hopfnf, ds) to obtain an estimate of the bifurcating periodic orbit.

hopfNormalForm(prob, pt, ls; verbose)

Compute the Hopf normal form.

Arguments

• prob::AbstractBifurcationProblem bifurcation problem
• pt::Hopf Hopf bifurcation point
• ls linear solver

Optional arguments

• verbose bool to print information

This function checks whether the solution with eigenvalues eigvalues is stable and also compute the number of unstable eigenvalues with nonzero imaginary part

Matrix by blocks expression of the Jacobian for the PO functional computed at the space-time guess: u0

kernelDim(bp)

Return the dimension of the kernel of the special point.

Function to locate precisely bifurcation points using a bisection algorithm. We make sure that at the end of the algorithm, the state is just after the bifurcation point (in the s coordinate).

modCounter(step, everyN)

This function implements a counter. If everyN == 0, it returns false. Otherwise, it returns true when step is a multiple of everyN

multicontinuation(br, ind_bif)
multicontinuation(
    br,
    ind_bif,
    optionsCont;
    δp,
    ampfactor,
    nev,
    Teigvec,
    ζs,
    verbosedeflation,
    scaleζ,
    perturbGuess,
    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 to continuation
• ind_bif index of the bifurcation point in br from which you want to branch from
• optionsCont options for the call to continuation

Optional arguments

• alg = br.alg continuation algorithm to be used, default value: br.alg
• δp used to specify a particular guess for the parameter on the bifurcated branch which is otherwise determined by optionsCont.ds. This allows to use a step larger than optionsCont.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.
• nev number of eigenvalues to be computed to get the right eigenvector
• verbosedeflation = true whether to display the nonlinear deflation iterations (see Deflated problems) to help finding the guess on the bifurcated branch
• scaleζ norm used to normalize eigenbasis when computing the reduced equation
• Teigvec type of the eigenvector. Useful when br was loaded from a file and this information was lost
• ζs basis of the kernel
• perturbGuess = identity perturb the guess from the predictor just before the deflated-newton correction
• kwargs optional arguments to be passed to continuation, the regular continuation one.

multicontinuation(br, bpnf, defOpm, defOpp)
multicontinuation(
    br,
    bpnf,
    defOpm,
    defOpp,
    optionsCont;
    alg,
    δp,
    Teigvec,
    verbosedeflation,
    maxIterDeflation,
    lsdefop,
    plotSolution,
    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 to continuation
• bpnf normal form
• defOpm::DeflationOperator, defOpp::DeflationOperator to specify converged points on nonn-trivial branches before/after the bifurcation points.

The rest is as the regular multicontinuation function.

newton(prob, orbitguess, options; lens, kwargs...)

This is the Newton-Krylov Solver for computing a periodic orbit using (Standard / Poincaré) Shooting method. Note that the linear solver has to be apropriately 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. See ShootingProblem and See PoincareShootingProblem for information regarding the shape of orbitguess.
• par parameters to be passed to the functional
• options same as for the regular newton method.

Optional argument

• jacobian Specify the choice of the linear algorithm, which must belong to [:autodiffMF, :MatrixFree, :autodiffDense, :autodiffDenseAnalytical, :FiniteDifferences]. This is used to select a way of inverting the jacobian dG
  • For :MatrixFree, we use an iterative solver (e.g. GMRES) to solve the linear system. The jacobian was specified by the user in prob.
  • For :autodiffMF, we use iterative solver (e.g. GMRES) to solve the linear system. We use Automatic Differentiation to compute the (matrix-free) derivative of x -> prob(x, p).
  • For :autodiffDense. Same as for :autodiffMF but the jacobian is formed as a dense Matrix. You can use a direct solver or an iterative one using options.
  • For :FiniteDifferencesDense, same as for :autodiffDense but we use Finite Differences to compute the jacobian of x -> prob(x, p) using the δ = 1e-8 which can be passed as an argument.
  • For :autodiffDenseAnalytical. Same as for :autodiffDense but the jacobian is using a mix of AD and analytical formula.
  • For :FiniteDifferences, use Finite Differences to compute the jacobian of x -> prob(x, p) using the δ = 1e-8 which can be passed as an argument.

	newton(prob::AbstractBifurcationProblem, options::NewtonPar; normN = norm, callback = (;x, f, J, res, iteration, itlinear, optionsN; 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. The function normN allows to specify a norm for the convergence criteria. 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 a BifurcationProblem which holds the vector field and its jacobian. We also refer to BifFunction for more details.
• options::NewtonPar variable holding the internal parameters used by the newton method
• callback function passed by the user which is called at the end of each iteration. The default one is the following cbDefault((x, f, J, res, it, itlinear, options); k...) = true. Can be used to update a preconditionner for example. You can use for example cbMaxNorm to limit the residuals norms. If yo want to specify your own, the arguments passed to the callback are as follows
  • x current solution
  • f current residual
  • J current jacobian
  • res current norm of the residual
  • iteration current newton iteration
  • itlinear number of iterations to solve the linear system
  • optionsN a copy of the argument options passed to newton
  • kwargs kwargs arguments, contain your initial guess x0
• kwargs arguments passed to the callback. Useful when newton is called from continuation

Output:

• solution::NonLinearSolution, we refer to NonLinearSolution for more information.

newton(
    br,
    ind_bif;
    normN,
    options,
    startWithEigen,
    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 continuation
• ind_bif bifurcation index in br

Optional arguments:

• options::NewtonPar, default value br.contparams.newtonOptions
• normN = norm
• options You can pass newton parameters different from the ones stored in br by using this argument options.
• bdlinsolver bordered linear solver for the constraint equation
• startWithEigen = 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

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 regular newton 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 using options. The jacobian is formed inplace.

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 type PeriodicOrbitTrapProblem encoding the functional G
• orbitguess a guess for the periodic orbit where orbitguess[end] is an estimate of the period of the orbit. It should be a vector of size N * M + 1 where M is the number of time slices, N is the dimension of the phase space. This must be compatible with the numbers N, M in prob.
• par parameters to be passed to the functional
• options same as for the regular newton method
• jacobian = :FullLU. Specify the choice of the jacobian (and linear algorithm), jacobian must belong to [:FullLU, :FullSparseInplace, :BorderedLU, :FullMatrixFree, :BorderedMatrixFree, :FullSparseInplace]. 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 of dG. 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 matrix dG 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 matrix dG is dense. It is also updated inplace. This option is useful to study ODE of small dimension.
• For jacobian = :BorderedLU, we take advantage of the bordered shape of the linear solver and use a LU decomposition to invert dG using a bordered linear solver.
• For jacobian = :BorderedSparseInplace, this is the same as for :BorderedLU but the cyclic matrix dG 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 for dG: note that a preconditioner is very likely required here because of the cyclic shape of dG which affects negatively the convergence properties of GMRES.
• For jacobian = :BorderedMatrixFree, a matrix free linear solver is used but for Jc only (see docs): it means that options.linsolver is used to invert Jc. These two Matrix-Free options thus expose different part of the jacobian dG in order to use specific preconditioners. For example, an ILU preconditioner on Jc could remove the constraints in dG and lead to poor convergence. Of course, for these last two methods, a preconditioner is likely to be required.
• For jacobian = :FullMatrixFreeAD, the evalution 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.

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 coonverged solution close to sol0. This is very handy for branch switching. The method is based on a deflated Newton-Krylov solver.

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]. This is used to select a way of inverting the jacobian dG
  • For :MatrixFree, we use an iterative solver (e.g. GMRES) to solve the linear system. The jacobian was specified by the user in prob.
  • For :autodiffMF, we use iterative solver (e.g. GMRES) to solve the linear system. We use Automatic Differentiation to compute the (matrix-free) derivative of x -> prob(x, p).
  • For :autodiffDense. Same as for :autodiffMF but the jacobian is formed as a dense Matrix. You can use a direct solver or an iterative one using options.
  • For :FiniteDifferencesDense, same as for :autodiffDense but we use Finite Differences to compute the jacobian of x -> prob(x, p) using the δ = 1e-8 which can be passed as an argument.
  • For :autodiffDenseAnalytical. Same as for :autodiffDense but the jacobian is using a mix of AD and analytical formula.
  • For :FiniteDifferences, use Finite Differences to compute the jacobian of x -> prob(x, p) using the δ = 1e-8 which can be passed as an argument.

Output:

• solution::NonLinearSolution, see NonLinearSolution

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 operator
• linsolver linear solver used to invert the Jacobian of the deflated functional.
  • custom solver DefProbCustomLinearSolver() with requires solving two linear systems J⋅x = rhs.
  • For other linear solvers <: AbstractLinearSolver, a matrix free method is used for the deflated functional.
  • if passed Val(:autodiff), then ForwardDiff.jl is used to compute the jacobian of the deflated problem
  • if passed Val(:fullIterative), then a full matrix free method is used.

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.

newtonBT(
    prob,
    btpointguess,
    par,
    lens2,
    eigenvec,
    eigenvec_ad,
    options;
    normN,
    jacobian_ma,
    usehessian,
    bdlinsolver,
    kwargs...
)

This function turns an initial guess for a BT point into a solution to the BT problem based on a Minimally Augmented formulation. The arguments are as follows

• prob::AbstractBifurcationFunction
• btpointguess initial guess (x0, p0) for the BT point. It should be a BorderedArray as returned by the function BTPoint
• par parameters used for the vector field
• eigenvec guess for the 0 eigenvector
• eigenvec_ad guess for the 0 adjoint eigenvector
• options::NewtonPar options for the Newton-Krylov algorithm, see NewtonPar.

Optional arguments:

• normN = norm
• bdlinsolver bordered linear solver for the constraint equation
• jacobian_ma::Symbol = true specify the way the (newton) linear system is solved. Can be (:autodiff, :finitedifferences, :minaug)
• kwargs keywords arguments to be passed to the regular Newton-Krylov solver

Simplified call

Simplified call to refine an initial guess for a Fold point. More precisely, the call is as follows

newton(br::AbstractBranchResult, ind_bt::Int; options = br.contparams.newtonOptions, kwargs...)

The parameters / options are as usual 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 bifurcation point in br you want to refine. You can pass newton parameters different from the ones stored in br by using the argument options.

newtonBT(
    br,
    ind_bt;
    probvf,
    normN,
    options,
    nev,
    startWithEigen,
    bdlinsolver,
    kwargs...
)

This function turns an initial guess for a Bogdanov-Takens point into a solution to the Bogdanov-Takens problem based on a Minimally Augmented formulation.

Arguments

• br results returned after a call to continuation
• ind_bif bifurcation index in br

Optional arguments:

• options::NewtonPar, default value br.contparams.newtonOptions
• normN = norm
• options You can pass newton parameters different from the ones stored in br by using this argument options.
• jacobian_ma::Symbol = true specify the way the (newton) linear system is solved. Can be (:autodiff, :finitedifferences, :minaug)
• bdlinsolver bordered linear solver for the constraint equation
• startWithEigen = 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

newtonFold(
    prob,
    foldpointguess,
    par,
    eigenvec,
    eigenvec_ad,
    options;
    normN,
    bdlinsolver,
    usehessian,
    kwargs...
)

This function turns an initial guess for a Fold point into a solution to the Fold problem based on a Minimally Augmented formulation. The arguments are as follows

• prob::AbstractBifurcationFunction
• foldpointguess initial guess (x0, p0) for the Fold point. It should be a BorderedArray as returned by the function FoldPoint
• par parameters used for the vector field
• eigenvec guess for the 0 eigenvector
• eigenvec_ad guess for the 0 adjoint eigenvector
• options::NewtonPar options for the Newton-Krylov algorithm, see NewtonPar.

Optional arguments:

• normN = norm
• bdlinsolver bordered linear solver for the constraint equation
• kwargs keywords arguments to be passed to the regular Newton-Krylov solver

Simplified call

Simplified call to refine an initial guess for a Fold point. More precisely, the call is as follows

newtonFold(br::AbstractBranchResult, ind_fold::Int; options = br.contparams.newtonOptions, kwargs...)

The parameters / options are as usual 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 bifurcation point in br you want to refine. You can pass newton parameters different from the ones stored in br by using the argument options.

newtonHopf(
    prob,
    hopfpointguess,
    par,
    eigenvec,
    eigenvec_ad,
    options;
    normN,
    bdlinsolver,
    usehessian,
    kwargs...
)

This function turns an initial guess for a Hopf point into a solution to the Hopf problem based on a Minimally Augmented formulation. The arguments are as follows

• prob::AbstractBifurcationProblem where p is a set of parameters.
• hopfpointguess initial guess (x0, p0) for the Hopf point. It should a BorderedArray as returned by the function HopfPoint.
• par parameters used for the vector field
• eigenvec guess for the iω eigenvector
• eigenvec_ad guess for the -iω eigenvector
• options::NewtonPar options for the Newton-Krylov algorithm, see NewtonPar.

Optional arguments:

• normN = norm
• bdlinsolver bordered linear solver for the constraint equation
• kwargs keywords arguments to be passed to the regular Newton-Krylov solver

Simplified call:

Simplified call to refine an initial guess for a Hopf point. More precisely, the call is as follows

newtonHopf(br::AbstractBranchResult, ind_hopf::Int; normN = norm, options = br.contparams.newtonOptions, kwargs...)

The parameters / options are as usual 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 bifurcation point in br you want to refine. You can pass newton parameters different from the ones stored in br by using the argument options.

This is the classical Newton-Krylov solver for F(x, p) = 0 together with the scalar condition n(x, p) ≡ θ ⋅ <x - x0, τx> + (1-θ) ⋅ (p - p0) * τp - n0 = 0. This makes a problem of dimension N + 1.

The initial guess for the newton method is located in state.z_pred

nf(bp; tol, digits)

Print the normal form bp with a nice string.

phaseCondition(pb, , )

[INTERNAL] Implementation of phase condition.

Arguments

• uj n x (m + 1)
• guj n x m

predictor(hp, ds; verbose, ampfactor)

This function provides prediction for the orbits of the Hopf bifurcation point.

Arguments

• bp::Hopf the bifurcation point
• ds at with distance relative to the bifurcation point do you want the prediction. Can be negative. Basically the parameter is p = bp.p + ds

Optional arguments

• verbose display information
• ampfactor = 1 factor multiplied to the amplitude of the periodic orbit.

predictor(
    bp,
    δp;
    verbose,
    ampfactor,
    nbfailures,
    maxiter,
    perturb,
    J,
    normN,
    optn
)

This function provides prediction for what the zeros of the reduced equation / normal form should be. The algorithm for finding these zeros is based on deflated newton.

predictor(bp, ds; verbose, ampfactor)

This function provides prediction for the zeros of the Pitchfork bifurcation point.

Arguments

• bp::Pitchfork the bifurcation point
• ds at with distance relative to the bifurcation point do you want the prediction. Based on the criticality of the Picthfork, its sign is enforced no matter what you pass. Basically the parameter is bp.p + abs(ds) * dsfactor where dsfactor = ±1 depending on the criticality.

Optional arguments

• verbose display information
• ampfactor = 1 factor multiplying prediction

predictor(bp, ds; verbose, ampfactor)

This function provides prediction for the zeros of the Transcritical bifurcation point.

Arguments

• bp::Transcritical the bifurcation point
• ds at with distance relative to the bifurcation point do you want the prediction. Can be negative. Basically the parameter is p = bp.p + ds

Optional arguments

• verbose display information
• ampfactor = 1 factor multiplying prediction

projection(psh, x)

Compute the projection of each vector (x[i, :] is a Vector) on the Poincaré section.

projection(psh, x)

Compute the projection of each vector (x[i] is a Vector) on the Poincaré section.

reMake(
    prob;
    u0,
    params,
    lens,
    recordFromSolution,
    plotSolution,
    J,
    d2F,
    d3F
)

This function allows to change the fields of a problem ::AbstractBifurcationProblem. For example, you can change the initial condition by doing

reMake(prob; u0 = new_u0)

split_events(cs, ds, args...)

Split comma seperated callbacks into sets of continous and discrete callbacks. Inspired by DiffEqBase.

update!(hyp::SectionPS, normals, centers)

Update the hyperplanes saved in hyp.

update!(pb::PoincareShootingProblem, centers_bar; _norm = norm)

This function updates the normals and centers of the hyperplanes defining the Poincaré sections.