Library

struct NewtonPar{T, L<:BifurcationKit.AbstractLinearSolver, E<:AbstractEigenSolver}

Returns a variable containing parameters to affect the newton algorithm when solving F(x) = 0.

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

Internal fields

• tol::Any: absolute tolerance for F(x). Default: 1.0e-12

• max_iterations::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)

struct ContinuationPar{T, S<:BifurcationKit.AbstractLinearSolver, E<:AbstractEigenSolver}

Returns a variable containing the parameters to affect the continuation algorithm used to solve F(x, p) = 0.

Arguments

• dsmin, dsmax are the minimum, maximum allowed arc-length value. It controls the density of points in the computed branch of solutions.
• ds = 0.01 is the initial arc-length.
• p_min, p_max allowed parameter range for p
• max_steps = 100 maximum number of continuation steps
• newton_options::NewtonPar: options for the Newton algorithm
• save_to_file = false: save to file. A name is automatically generated or can be defined in continuation. This requires using JLD2.
• save_sol_every_step::Int64 = 1 at which continuation steps do we save the current solution
• plot_every_step = 10 at which continuation steps do we plot the current solution

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

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

Handling bifurcation detection

• tol_stability = 1e-10 lower bound on the real part of the eigenvalues to test for stability of equilibria and periodic orbits
• detect_fold = true detect Fold bifurcations? It is a useful option although the detection of Fold is cheap. Indeed, it may happen that there is a lot of Fold points and this can saturate the memory in memory limited devices (e.g. on GPU)
• detect_bifurcation::Int ∈ {0, 1, 2, 3} If set to 0, nothing is done. If set to 1, the eigen-elements are computed. If set to 2, the bifurcations points are detected during the continuation run, but not located precisely. If set to 3, a bisection algorithm is used to locate the bifurcations points (slower). The possibility to switch off detection is a useful option. Indeed, it may happen that there are a lot of bifurcation points and this can saturate the memory of memory limited devices (e.g. on GPU)
• dsmin_bisection = 1e-16 minimal ds for the bisection algorithm for locating bifurcation points
• n_inversion = 2 number of sign inversions in bisection algorithm
• max_bisection_steps = 25 maximum number of bisection steps
• tol_bisection_eigenvalue = 1e-16 tolerance on real part of eigenvalue to detect bifurcation points in the bisection steps

Handling ds adaptation (see continuation for more information)

• a = 0.5 aggressiveness factor. It is used to 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

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

Misc

• η = 150. parameter to estimate tangent at first point with parameter p₀ + ds / η
• detect_loop [WORK IN PROGRESS] detect loops in the branch and stop the continuation

cb = cbMaxNorm(maxres)

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

cb = cbMaxNormAndΔp(maxres, δp)

Create a callback used to reject residuals larger than cb.maxres or parameter step larger than δp in the Newton iterations. See docs for newton.

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

For example

sol = newton(prob, NewtonPar())

methods

• converged(sol) return whether the solution has converged.

Internal fields

• u::Any: solution.

• prob::Any: nonlinear problem, typically a BifurcationProblem.

• residuals::Any: sequence of residuals.

• converged::Bool: has algorithm converged?

• itnewton::Int64: number of newton steps.

• itlineartot::Any: total number of linear iterations.

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

Internal fields

• branch::StructArrays.StructArray: holds the low-dimensional information about the branch. More precisely, branch[i+1] contains the following information (record_from_solution(u, param), param, itnewton, itlinear, ds, n_unstable, n_imag, stable, step) for each continuation step i.

  • itnewton number of Newton iterations.
  • itlinear total number of linear iterations during newton (corrector).
  • n_unstable number of eigenvalues with positive real part for each continuation step (to detect stationary bifurcation).
  • n_imag number of eigenvalues with positive real part and non zero imaginary part at current continuation step (useful to detect Hopf bifurcation).
  • stable stability of the computed solution for each continuation step. Hence, stable should match eig[step] which corresponds to branch[k] for a given k.
  • step continuation step (here equal i).

• eig::Array{@NamedTuple{eigenvals::Teigvals, eigenvecs::Teigvec, converged::Bool, step::Int64}, 1} where {Teigvals, Teigvec}: A vector with eigen-elements at each continuation step.

• sol::Any: Vector of solutions sampled along the branch. This is set by the argument save_sol_every_step::Int64 (default 1) in ContinuationPar.

• contparams::Any: The parameters used for the call to continuation which produced this branch. Must be a ContinuationPar.

• kind::BifurcationKit.AbstractContinuationKind: Type of solutions computed in this branch. Default: EquilibriumCont()

• prob::Any: Bifurcation problem used to compute the branch, useful for branch switching. For example, when computing periodic orbits, the functional Trapeze, Shooting... will be saved here. Default: nothing

• specialpoint::Vector: A vector holding the list of detected bifurcation points. See SpecialPoint for a list of special points.

• alg::Any: Continuation algorithm used for the computation of the branch

Associated methods

• length(br) number of the continuation steps.
• show(br) display information about the branch.
• propertynames(br) give the propertynames of a result.
• 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.
• get_normal_form(br, ind) compute the normal form of the ind-th points in br.specialpoint.
• getlens(br) return the parameter axis used for the branch.
• get_lenses(br) return the parameter two axis used for the branch when 2 parameters continuation is used (Fold, Hopf, NS, PD).
• get_solx(br, k) returns the k-th solution on the branch.
• get_solp(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.
• getparams(br, ind) Parameters passed to continuation and used in the equation F(x, par) = 0 for the ind-th continuation step.
• setparam(br, p0) set the parameter value p0 according to ::Lens for the parameters of the problem getprob(br).
• getlens(br) get the lens used for the computation of the branch.
• eigenvals(br, ind) give the eigenvalues at continuation step ind.
• eigenvalsfrombif(br, ind) give the eigenvalues at bifurcation point index ind.
• BifurcationKit._merge(br1, br2) merge twos branches into a single ContResult.
• type(br, ind) returns the type of the ind-th bifurcation point.
• br[k+1] gives information about the k-th step. A typical run yields the following.
• get_solution(br, ind) returns the ind-th solution.

julia> br[1]
(x = 0.0, param = 0.1, itnewton = 0, itlinear = 0, ds = -0.01, n_unstable = 2, n_imag = 2, stable = false, step = 0, eigenvals = ComplexF64[0.1 - 1.0im, 0.1 + 1.0im], eigenvecs = ComplexF64[0.7071067811865475 - 0.0im 0.7071067811865475 + 0.0im; 0.0 + 0.7071067811865475im 0.0 - 0.7071067811865475im])

which provides the value param of the parameter of the current point, its stability, information on the newton iterations, etc. The fields can be retrieved using propertynames(br.branch). This information is stored in br.branch which is a StructArray. You can thus extract the vector of parameters along the branch as

julia> br.param
10-element Vector{Float64}:
 0.1
 0.08585786437626905
 0.06464466094067263
 0.03282485578727799
-1.2623798512809007e-5
-0.07160718539365075
-0.17899902778635765
-0.3204203840236672
-0.4618417402609767
-0.5

• continuation(br, ind) performs automatic branch switching (aBS) from ind-th bifurcation point. Typically branching from equilibrium to equilibrium, or periodic orbit to periodic orbit.
• continuation(br, ind, lens2) performs two parameters (getlens(br), lens2) continuation of the ind-th bifurcation point.
• continuation(br, ind, probPO::AbstractBoundaryValueDiscretization) performs aBS from ind-th bifurcation point (which must be a Hopf bifurcation point) to branch of periodic orbits.

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 SpecialPoint{T, Tp, Tv, Tvτ} <: BifurcationKit.AbstractBifurcationPoint

Structure to record special points on a curve. There are two types of special points that are recorded in this structure: bifurcation points and events (see https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/EventCallback/).

Internal fields

• type::Symbol: Description of the special points. In case of Events, this field records the user passed name to the event, or the default :userD, :userC. In case of bifurcation points, it can be one of the following names:

  - :bp Bifurcation point, simple eigenvalue crossing the imaginary axis
  - :fold Fold point
  - :hopf Hopf point
  - :nd not documented bifurcation point. Detected by multiple eigenvalues crossing. Generally occurs in problems with symmetries or in cases where the continuation step size is too large and merge two different bifurcation points.
  - :cusp Cusp point
  - :gh Generalized Hopf point (also called Bautin point)
  - :bt Bogdanov-Takens point
  - :zh Zero-Hopf point
  - :hh Hopf-Hopf point
  - :ns Neimark-Sacker point
  - :pd Period-doubling point
  - :R1 Strong resonance 1:1 of periodic orbits
  - :R2 Strong resonance 1:2 of periodic orbits
  - :R3 Strong resonance 1:3 of periodic orbits
  - :R4 Strong resonance 1:4 of periodic orbits
  - :foldFlip Fold / Flip of periodic orbits
  - :foldNS Fold / Neimark-Sacker of periodic orbits
  - :pdNS  Period-Doubling / Neimark-Sacker of periodic orbits
  - :gpd Generalized Period-Doubling of periodic orbits
  - :nsns Double Neimark-Sacker of periodic orbits
  - :ch Chenciner bifurcation of periodic orbits
   Default: :none

• idx::Int64: Index in br.branch or br.eig (see ContResult) for which the bifurcation occurs. Default: 0

• param::Any: Parameter value at the special point (this is an estimate). Default: 0.0

• norm::Any: Norm of the equilibrium at the special point. Default: 0.0

• printsol::Any: printsol = record_from_solution(x, param) where record_from_solution is one of the arguments to continuation. Default: 0.0

• x::Any: Solution at the special point. Default: Vector{T}(undef, 0)

• τ::BorderedArray{Tvτ, T} where {T, Tvτ}: Tangent along the branch at the special point. Default: BorderedArray(x, zero(T))

• ind_ev::Int64: Eigenvalue index responsible for detecting the special point (if applicable). Default: 0

• step::Int64: Continuation step at which the special occurs. Default: 0

• status::Symbol: status ∈ {:converged, :guess, :guessL} indicates whether the bisection algorithm was successful in detecting the special (bifurcation) point. If status == :guess, the bisection algorithm failed to meet the requirements given in ::ContinuationPar. Same for status == :guessL but the bisection algorithm stopped on the left of the bifurcation point. Default: :guess

• δ::Tuple{Int64, Int64}: δ = (δr, δi) where δr indicates the change in the number of unstable eigenvalues and δi indicates the change in the number of unstable eigenvalues with nonzero imaginary part. abs(δr) is thus an estimate of the dimension of the kernel of the Jacobian at the special (bifurcation) point. Default: (0, 0)

• precision::Any: Precision in the location of the special point Default: -1

• interval::Tuple{T, T} where T: Interval parameter containing the special point Default: (0, 0)

Associated methods

• BifurcationKit.type(::SpecialPoint) returns the bifurcation type (::Symbol)

Structure to hold a connected component of a bifurcation diagram which is encoded as a tree of BifDiagNode(s).

Internal fields

• level::Int64: current recursion level in the tree.

• code::Int64: code for finding the current node in the tree, this is the index of the bifurcation point from which γ branches off.

• γ::Any: branch associated to the current node.

• child::Any: children of current node. These are the different branches off the bifurcation point in γ.

Methods

• hasbranch(diagram)
• get_branch(diagram) return the AbstractBranch stored inside the current node.
• from(diagram) return the parent bifurcation point.
• diagram[code] For example diagram[1,2,3] returns diagram.child[1].child[2].child[3]. This is essentially the Ulam-Harris-Neveu labeling.

struct BifFunction{Tf, TFinp, Tdf, Tdfad, Tj, Tjad, TJinp, Td2f, Td2fc, Td3f, Td3fc, Tδ, Tjet} <: 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.

Internal 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

• F!::Any: Same as F but inplace with signature F!(result, x, p)

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

• J!::Any: Inplace jacobian

• 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::Bool: Whether the jacobian is auto-adjoint.

• δ::Any: used internally to compute derivatives (with finite differences), for example for normal form computation and codim 2 continuation.

• inplace::Bool: optionally sets whether the function is inplace or not. You can use in_bisection(state) to inquire whether the current state is in bisection mode.

• jet::Any: jet of the vector field.

Methods

• residual(pb::BifFunction, x, p) calls pb.F(x,p)
• residual!(pb::BifFunction, o, x, p) calls pb.F(o, x, p)
• jacobian(pb::BifFunction, x, p) calls pb.J(x, p)
• dF(pb::BifFunction, x, p, dx) calls pb.dF(x, p, dx)
• R21(pb::BifFunction, x, p, dx1, dx2, dp1) calls pb.jet.R21(x, p, dx1, dx2, dp1). Same for the other jet functions.
• etc

struct BifurcationProblem{Tvf, Tu, Tp, Tl<:Union{typeof(identity), IndexLens, PropertyLens, ComposedFunction}, Tplot, Trec, Tgets, Tupdate} <: BifurcationKit.AbstractAllJetBifProblem

Structure to hold a bifurcation problem. Generic case, the user has to set most options.

Methods

• re_make(pb; kwargs...) modify a bifurcation problem
• 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)
• record_from_solution(pb) calls pb.recordFromSolution
• plot_solution(pb) calls pb.plotSolution
• is_symmetric(pb) calls is_symmetric(pb.VF)
• getdelta(prob)

Constructors

• BifurcationProblem(F, u0, params, lens) all derivatives are computed using ForwardDiff.
• 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 for description of the jacobian function) and jacobian adjoint with Jᵗ. For example, this can be used to provide finite differences based jacobian using BifurcationKit.finite_differences. You can also pass
  • record_from_solution see above
  • plot_solution see above
  • issymmetric[=false] whether the jacobian is symmetric, this remove the need of providing an adjoint
  • jvp jacobian-vector product, signature jvp(x, p, dx)
  • vjp vector-jacobian product (adjoint of jvp), signature vjp(x, p,dx)
  • d2F second Differential of F with respect to x, signature d2F(x, p, dx1, dx2)
  • d3F third Differential of F with respect to x, signature d3F(x, p, dx1, dx2, dx3)
  • save_solution specify a particular way to record solution which are written in br.sol. This can be useful in very particular situations and we recommend using record_from_solution instead. For example, it is used internally to record the mesh in the collocation method because this mesh can be modified.

Internal fields

• VF::Any: Vector field, typically a BifFunction.

• u0::Any: Initial guess.

• params::Any: Parameters.

• lens::Union{typeof(identity), IndexLens, PropertyLens, ComposedFunction}: Typically a Accessors.PropertyLens. It specifies which parameter axis among params is used for continuation. For example, if par = (α = 1.0, β = 1.78), we can perform continuation w.r.t. α by using lens = (@optic _.α). 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 = (@optic _[1]) or pass directly lens = 1. For more information, we refer to Accessors.jl.

• plotSolution::Any: user function to plot solutions during continuation. Signature: plot_solution(x, p; kwargs...) for Plot.jl and plot_solution(ax, x, p; ax1 = nothing, kwargs...) for the Makie.jl.

• recordFromSolution::Any: record_from_solution = (x, p; k...) -> norm(x) function used to record a few indicators about the solution. It could be norm or (x, p; k...) -> 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; k...) -> (x1 = x[1], x2 = x[2], nrm = norm(x)) or simply (x, p; k...) -> (sum(x), 1) or (x, p; k...) -> [x[1], x[2]]. This will be stored in contres.branch where contres::AbstractBranchResult is the continuation curve of the bifurcation problem. Finally, the first component is used for plotting in the continuation curve.

• save_solution::Any: Function to save the full solution on the branch. Some problem are updated during computation (like periodic orbit functional with adaptive mesh) and this function allows to save the state of the problem along with the solution itself. Note that this should allocate the output (i.e. not as a view). Signature: save_solution(x, p). Defaults to save_solution_default(x, p) = x. The saved solution can be retrieved back using saved_solution with signature saved_solution(x_saved). Finally, a function restore_problem!(prob, x_saved, pars) which defaults to restore_problem!(prob, x_saved, pars) = prob should be provided. It allows to update the problem to the state it was at x_saved

• update!::Any: Function to update the problem after each continuation step. Defaults to update_default. It has signature update!(prob, iter, state) where prob is the current bifurcation problem, iter::ContIterable the current iterable and state::ContState the current state of the continuation. It should return true if the update was successful and false otherwise. The continuation will stop if false is returned. This type of function is useful for example to update the problem internals like meshes, preconditioners, etc. Note that you can extract the type of continuation from iter (for example FoldCont) and thus modify the method update! accordingly.

struct ODEBifProblem{Tvf, Tu, Tp, Tl<:Union{typeof(identity), IndexLens, PropertyLens, ComposedFunction}, Tplot, Trec, Tgets, Tupdate} <: BifurcationKit.AbstractAllJetBifProblem

Structure to hold a bifurcation problem. Specific to Ordinary Differential Equations (ODE). The options are set accordingly. 🚧🚧 This is work in progress 🚧🚧.

Methods

• re_make(pb; kwargs...) modify a bifurcation problem
• 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)
• record_from_solution(pb) calls pb.recordFromSolution
• plot_solution(pb) calls pb.plotSolution
• is_symmetric(pb) calls is_symmetric(pb.VF)
• getdelta(prob)

Constructors

• ODEBifProblem(F, u0, params, lens) all derivatives are computed using ForwardDiff.
• ODEBifProblem(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 for description of the jacobian function) and jacobian adjoint with Jᵗ. For example, this can be used to provide finite differences based jacobian using BifurcationKit.finite_differences. You can also pass
  • record_from_solution see above
  • plot_solution see above
  • issymmetric[=false] whether the jacobian is symmetric, this remove the need of providing an adjoint
  • jvp jacobian-vector product, signature jvp(x, p, dx)
  • vjp vector-jacobian product (adjoint of jvp), signature vjp(x, p,dx)
  • d2F second Differential of F with respect to x, signature d2F(x, p, dx1, dx2)
  • d3F third Differential of F with respect to x, signature d3F(x, p, dx1, dx2, dx3)
  • save_solution specify a particular way to record solution which are written in br.sol. This can be useful in very particular situations and we recommend using record_from_solution instead. For example, it is used internally to record the mesh in the collocation method because this mesh can be modified.

Internal fields

• VF::Any: Vector field, typically a BifFunction.

• u0::Any: Initial guess.

• params::Any: Parameters.

• lens::Union{typeof(identity), IndexLens, PropertyLens, ComposedFunction}: Typically a Accessors.PropertyLens. It specifies which parameter axis among params is used for continuation. For example, if par = (α = 1.0, β = 1.78), we can perform continuation w.r.t. α by using lens = (@optic _.α). 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 = (@optic _[1]) or pass directly lens = 1. For more information, we refer to Accessors.jl.

• plotSolution::Any: user function to plot solutions during continuation. Signature: plot_solution(x, p; kwargs...) for Plot.jl and plot_solution(ax, x, p; ax1 = nothing, kwargs...) for the Makie.jl.

• recordFromSolution::Any: record_from_solution = (x, p; k...) -> norm(x) function used to record a few indicators about the solution. It could be norm or (x, p; k...) -> 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; k...) -> (x1 = x[1], x2 = x[2], nrm = norm(x)) or simply (x, p; k...) -> (sum(x), 1) or (x, p; k...) -> [x[1], x[2]]. This will be stored in contres.branch where contres::AbstractBranchResult is the continuation curve of the bifurcation problem. Finally, the first component is used for plotting in the continuation curve.

• save_solution::Any: Function to save the full solution on the branch. Some problem are updated during computation (like periodic orbit functional with adaptive mesh) and this function allows to save the state of the problem along with the solution itself. Note that this should allocate the output (i.e. not as a view). Signature: save_solution(x, p). Defaults to save_solution_default(x, p) = x. The saved solution can be retrieved back using saved_solution with signature saved_solution(x_saved). Finally, a function restore_problem!(prob, x_saved, pars) which defaults to restore_problem!(prob, x_saved, pars) = prob should be provided. It allows to update the problem to the state it was at x_saved

• update!::Any: Function to update the problem after each continuation step. Defaults to update_default. It has signature update!(prob, iter, state) where prob is the current bifurcation problem, iter::ContIterable the current iterable and state::ContState the current state of the continuation. It should return true if the update was successful and false otherwise. The continuation will stop if false is returned. This type of function is useful for example to update the problem internals like meshes, preconditioners, etc. Note that you can extract the type of continuation from iter (for example FoldCont) and thus modify the method update! accordingly.

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

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) ⋅ Πᵢ(||u - rootᵢ||₂⁻²ᵖ + α) := F(u) ⋅ M(u)$

where $||u||₂² = dot(u, u)$. The fields of the struct DeflationOperator are as follows:

Internal fields

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

• accumulator::Any: accumulator ∈(Val(:Prod), Val(:Mean))

• tmp::Any: [internal] to reduce allocations during computation.

• autodiff::Bool: [internal] to reduce allocations during computation.

• δ::Real: [internal] for finite differences.

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(power::Real, α::T, roots::Vector{vectype}; autodiff = false, δ = convert(T, 1e-8), accumulator = Val(:Prod))
• DeflationOperator(power::Real, dot, α::Real, roots::Vector{vectype}; autodiff = false, δ = convert(VI.scalartype(roots[1]), 1e-8) , accumulator = Val(:Prod))
• DeflationOperator(power::Real, α::T, roots::Vector{vectype}, tmp::vectype; autodiff = false, accumulator = Val(:Prod))

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

The option autodiff of the constructors above also applies when a custom distance is used.

Linear solvers / jacobians

When used with newton, you have access to the following linear solvers:

• custom solver DeflatedProblemCustomLS() which requires solving two linear systems J⋅x = rhs.
• if passed Val(:autodiff), then ForwardDiff.jl is used to compute the jacobian Matrix of the deflated problem.
• if passed Val(:fullIterative), then a full matrix free method is used for the deflated problem.

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.

Arguments

• jactype selects the jacobian for the newton solve. Can be Val(:autodiff), Val(:fullIterative), Val(:Custom)

struct Trapeze{Tprob, vectype, Tls<:BifurcationKit.AbstractLinearSolver, T, Tmass, Tjac} <: BifurcationKit.AbstractFiniteDifferencesDiscretization

This composite type implements a finite-difference discretization based on the Trapeze rule (aka Crank-Nicolson, order 2 in time) to locate periodic orbits / BVP. More details (maths, notations, linear systems) can be found here.

The scheme is as follows. We discretize the time on M slices $x_{1},\cdots,x_{M}$, each of dimension N, and one looks for the period T = x[end] such that the following Crank-Nicolson relations hold for $i = 1, \cdots, M-1$

$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,$

where we used the cyclic convention $x_{0} := x_{M-1}$ and where $h_i = s_i - s_{i-1}$ are the normalized steps of the mesh mesh ($s_0 = 0 < s_1 < \cdots < s_M = 1$). The orbit is finally closed by the periodicity condition $x_{M} - x_{1} = 0$. Here $M_a$ is a mass matrix (the identity by default, see field massmatrix) and $F$ stands for the residual of the vector field encoded in prob_vf.

The phase of the periodic orbit, which removes the indeterminacy due to the invariance of periodic orbits under time shifts, is constrained by using a section (but you could use your own)

$\sum_i\langle x_{i} - x_{\pi,i}, \phi_{i}\rangle=0.$

The pair $(\phi, x_{\pi})$ is stored in the fields ϕ and xπ: $x_\pi$ is a reference state on the orbit while $\phi$ are the (normalized) normal vectors of the section. It is updated automatically during continuation by updatesection! every update_section_every_step steps.

Internal fields

• prob_vf::Any: Vector field (or bifurcation problem) whose residual F and jacobian are used to assemble the functional G. nothing is allowed to build a bare discretization. Default: nothing

• ϕ::Any: Normal vectors ϕ of the phase constraint, of size N * M (see the section equation in the documentation of Trapeze). Default: nothing

• xπ::Any: Reference point xπ of the phase constraint, of size N * M (see the section equation in the documentation of Trapeze). Default: nothing

• M::Int64: Number of time slices. Default: 0

• mesh::BifurcationKit.TimeMesh: Mesh of (normalized) time steps, see TimeMesh. Default: TimeMesh(M)

• N::Int64: Dimension of the problem in case of an AbstractVector state space. Default: 0

• linsolver::BifurcationKit.AbstractLinearSolver: Linear solver used to invert the jacobian of a single time slice, i.e. to solve J⋅sol = rhs. Only needed in a matrix-free setting (e.g. BorderedMatrixFree()) or for the computation of the Floquet multipliers. Default: DefaultLS()

• ongpu::Bool: Whether the computation takes place on the gpu (Experimental). When true, the functional returns vcat(out[begin:end-1], phase_cond) which is compatible with CuArrays in the mode allowscalar(false). Default: false

• isautonomous::Bool: Whether the vector field is autonomous, i.e. does not depend explicitly on time. Default: true

• massmatrix::Any: Mass matrix $M_a$ of the time discretization. You can pass for example a sparse matrix. Default: nothing, i.e. the identity matrix. Default: nothing

• update_section_every_step::UInt64: Frequency at which the phase constraint is updated during continuation, see updatesection!. Default: 1

• jacobian::Any: Type of jacobian used in Newton iterations (see the # Jacobian section in the documentation of Trapeze). Default: Dense()

Methods

Here are some useful methods you can apply to pb::Trapeze:

• length(pb) gives the number M * N of unknowns of the time-discretized state (without the period T, which is stored as x[end]).
• get_mesh_size(pb) returns the number M of time slices.
• get_state_dim(pb) returns the dimension N of a time slice.
• get_times(pb) returns the normalized times sᵢ at which the orbit is discretized, i.e. the cumulative sum of the mesh steps.
• get_time_slices(pb, x) returns the state part of the guess x (i.e. x[1:M*N], the period is dropped) reshaped as an N x M matrix.
• get_time_step(pb, i) returns the i-th normalized mesh step hᵢ.
• get_mass_matrix(pb) returns the mass matrix, defaulting to a sparse identity matrix if none was provided. Passing Val(true) as a second argument returns instead an identity matrix of the form I(N).
• hasmassmatrix(pb) returns true if a mass matrix was provided.
• getparams(pb), getlens(pb) and setparam(pb, p) give access to the parameters of the underlying vector field.
• getperiod(pb, x) returns the period T = x[end] of the guess x.
• getdelta(pb) returns the step δ used for finite differences.
• generate_solution(pb, orbit, period) generates a guess from a function t -> orbit(t) for t ∈ [0, 2π] and a period period.
• generate_ci_problem(pb, bifprob, sol, tspan) generates a Trapeze problem together with a guess from an ODE solution sol.
• get_periodic_orbit(pb, x, pars) computes the full periodic orbit, mainly for plotting purposes.

Constructors

The structure can be created by calling Trapeze(;kwargs...). For example, you can declare such a problem without vector field by doing

Trapeze(M = 100)

A more realistic way to build the problem is to provide the bifurcation problem together with the number of time slices M (or a non-uniform mesh given by a vector of steps) and the dimension N of a time slice:

Trapeze(prob_vf, M::Int, N::Int)
Trapeze(prob_vf, ϕ, xπ, M::Int, N::Int, ls = DefaultLS(); kwargs...)

In the second form, ϕ and xπ (see above) provide the initial section; they are stored into vectors of length N * M, the extra entries (if the provided vectors are shorter) being set to 0. The keyword massmatrix allows to specify a mass matrix. When the discretization is created with a vector field, the residual F of prob_vf and its jacobian are used to assemble the functional G.

Orbit guess

An orbit guess 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 generate_solution or from an ODEProblem solution using generate_ci_problem. You can evaluate the residual of the functional G on an orbit guess using po_residual(pb, orbitguess, p) and its jacobian with the methods listed in the # Functional section below.

Functional

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

• po_residual(pb, orbitguess, p) evaluates the functional G on orbitguess
• po_residual!(pb, out, orbitguess, p) same as po_residual but writes the result into out
• po_jvp(pb, orbitguess, p, du) evaluates the jacobian dG(orbitguess)⋅du functional at orbitguess on du
• po_jacobian_sparse(pb, orbitguess, p) returns the sparse matrix of the jacobian dG(orbitguess) at orbitguess. It is called $A_γ$ in the docs.
• po_jacobian_sparse!(pb, J, orbitguess, p). Same as po_jacobian_sparse but overwrites J inplace. Note that the sparsity pattern must be the same independently of the values of the parameters or of orbitguess. In this case, this is significantly faster than po_jacobian_sparse.
• jacobian_cyclic_sparse(pb, orbitguess, p) returns the sparse cyclic matrix $J_c$ (see the docs) of the jacobian dG(orbitguess) at orbitguess
• jacobian_block_diag(pb, orbitguess, p) returns the block diagonal of the sparse matrix of the jacobian dG(orbitguess) at orbitguess, i.e. the matrices $I - (T\,h_i/2)\,J(x_i)$ associated to each slice (the last block being the identity). Its inverse is a natural block Jacobi preconditioner.

Jacobian

These methods only differ in the linear algebra used to invert the jacobian dG of the functional G (see Trapeze); the discretization is otherwise the same. The value of jacobian must belong to (BifurcationKit.Dense(), BifurcationKit.AutoDiffDense(), BifurcationKit.FullLU(), BifurcationKit.FullMatrixFree(), BifurcationKit.BorderedLU(), BifurcationKit.BorderedMatrixFree(), BifurcationKit.FullSparseInplace(), BifurcationKit.BorderedSparseInplace(), BifurcationKit.AutoDiffMF()).

• 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 right choice when the sparsity pattern can change.
• For jacobian = FullSparseInplace(), this is the same as for FullLU() but the sparse matrix dG is updated inplace. This method allocates much less and, in some cases, is significantly faster than 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, and it is also updated inplace. This option is useful to study ODEs of small dimension.
• For jacobian = AutoDiffDense(), the jacobian is evaluated using automatic differentiation (ForwardDiff).
• For jacobian = BorderedLU(), we take advantage of the bordered shape of dG and invert it with a bordered linear solver based on a LU decomposition of the cyclic matrix.
• For jacobian = BorderedSparseInplace(), this is the same as for BorderedLU() but the cyclic matrix Jc is updated inplace. This method allocates much less and, in some cases, is significantly faster than 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 (given by options.linsolver) is used to invert dG: note that a preconditioner is very likely required here because of the cyclic shape of dG which negatively affects the convergence properties of GMRES.
• For jacobian = BorderedMatrixFree(), a matrix-free linear solver is used as well but only for Jc (see the docs): options.linsolver is then used to invert Jc. These two matrix-free options thus expose different parts of the jacobian dG in order to apply 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 = AutoDiffMF(), the evaluation map of the differential is derived using automatic differentiation. Thus, unlike the previous two cases, the user does not need to pass a matrix-free differential.

struct Collocation{Tprob<:Union{Nothing, BifurcationKit.AbstractBifurcationProblem}, Tjac<:BifurcationKit.AbstractJacobianType, 𝒯, vectype, ∂vectype, Tmass} <: BifurcationKit.AbstractDifferentialDiscretization

This composite type implements an orthogonal collocation (at Gauss points) method of piecewise polynomials to locate periodic orbits / BVP. More details (maths, notations, linear systems) can be found online.

Internal fields

• prob_vf::Union{Nothing, BifurcationKit.AbstractBifurcationProblem}: Bifurcation problem. Default: nothing

• ϕ::Any: Used to set a section for the phase constraint equation. Default: nothing

• xπ::Any: Used in the section for the phase constraint equation. Default: nothing

• ∂ϕ::Any: We store the derivative of ϕ, no need to recompute it each time. Default: nothing

• N::Int64: Dimension of the state space. Default: 0

• isautonomous::Bool: Whether the vector field is autonomous, i.e. does not depend explicitly on time. Default: true

• massmatrix::Any: Mass matrix to handle differential-algebraic (DAE) problems. Defaults to nothing for standard ODE problems. Default: nothing

• update_section_every_step::UInt64: Update the section every update_section_every_step step during continuation. Default: 1

• jacobian::BifurcationKit.AbstractJacobianType: Describes the type of jacobian used in Newton/PALC/etc iterations. See below for more information. Default: DenseAnalytical()

• mesh_cache::BifurcationKit.MeshCollocationCache: Cache for collocation. See docs of MeshCollocationCache. Default: nothing

• cache::BifurcationKit.POCollCache: Cache for allocation free computations. See docs of POCollCache. Default: nothing

• meshadapt::Bool: Whether to use mesh adaptation. Default: false

• verbose_mesh_adapt::Bool: Verbose mesh adaptation information. Default: false

• K::Float64: Parameter for mesh adaptation, control new mesh step size. More precisely, we set max(hᵢ) / min(hᵢ) ≤ K if hᵢ denotes the time steps. Default: 100

Methods

Here are some useful methods you can apply to coll::Collocation:

• length(coll) gives the total number of unknowns.
• size(coll) returns the triplet (N, m, Ntst).
• getmesh(coll) returns the mesh 0 = τ₁ < ... < τₙₜₛₜ₊₁ = 1. This is useful because this mesh is bound to vary during automatic mesh adaptation.
• get_mesh_coll(coll) returns the (static) mesh -1 = σ₁ < ... < σₘ₊₁ = 1.
• get_times(coll) returns the vector of times (length 1 + m * Ntst) at which the collocation is applied.
• generate_solution(coll, orbit, period) generates a guess from a function t -> orbit(t) which approximates the periodic orbit.
• POInterpolation(coll, x) returns a function interpolating the solution x using a piecewise polynomial function.
• getperiod(coll, po, p) returns the period of the periodic orbit po.

Orbit guess

An orbit guess 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.

Collocation is a discretization, not a functional. You can evaluate the residual G and its jacobian on an orbit guess using po_residual(coll, orbitguess, p), po_residual!(coll, out, orbitguess, p) and po_analytical_jacobian(coll, orbitguess, p).

Note that you can generate this guess from a function using generate_solution or generate_ci_problem.

Jacobian

Specify the choice of the jacobian (and linear algorithm), jacobian must belong to (BifurcationKit.AutoDiffDense(), BifurcationKit.DenseAnalytical(), BifurcationKit.FullSparse(), BifurcationKit.DenseAnalyticalInplace(), BifurcationKit.FullSparseInplace(), BifurcationKit.AutoDiffMF()).

This is used to select a way of inverting the jacobian dG of the functional G. See website for more information.

Mesh adaptation

Mesh adaptation is activated by setting meshadapt = true. It then modifies the mesh, i.e. the distribution of the Ntst intervals over [0, 1], so as to equilibrate an estimate of the discretization error along the orbit. More precisely, it is performed every update_section_every_step continuation steps after a successful Newton convergence (and not during bisection).

At each adaptation, the local error is estimated from the jump of the m-th derivative of the polynomial interpolating the orbit, which stands in for the (m + 1)-th derivative of the true solution. From this, a monitor function ϕ is built and a new mesh is computed by equipartition of ∫ϕ, i.e. the intervals are redistributed so that each of them carries the same amount of error.

The new mesh is enforced to satisfy max(hᵢ) / min(hᵢ) ≤ K, where hᵢ denotes the time steps; this bound is set through the parameter K. Set verbose_mesh_adapt = true to print information about the new mesh and the monitor function at each adaptation.

Note that mesh adaptation modifies getmesh(coll) in place. The mesh is stored in the solutions saved along the branch (see POSavedSolutionAndState), so that it can be restored, e.g. when starting a new branch from a bifurcation point.

Constructors

• Collocation(Ntst::Int, m::Int; kwargs) creates an empty functional with Ntst and m.

Functional

A functional, hereby called G, encodes this problem. Collocation is a discretization of G; it is wrapped into the functional PeriodicOrbit(coll). The following methods are available on this functional

• residual(PeriodicOrbit(coll), orbitguess, p) evaluates the functional G on orbitguess
• residual!(PeriodicOrbit(coll), out, orbitguess, p) evaluates the functional G on orbitguess
• jacobian(probPO, orbitguess, p) evaluates the jacobian dG of the functional G on orbitguess, where probPO = PeriodicOrbitFunctionalColl(coll) is the wrapped problem.

struct Shooting{Tf<:BifurcationKit.AbstractFlow, Tjac<:BifurcationKit.AbstractJacobianType, Ts, Tsection, Tpar, Tlens} <: BifurcationKit.AbstractShootingDiscretization

Create a problem to implement the Simple / Parallel Multiple Standard Shooting method to locate periodic orbits / BVP. More details (maths, notations, linear systems) can be found here. The arguments are as described below.

A functional, hereby called G, encodes the shooting problem. For example, the following methods are available:

• po_residual(pb, orbitguess, par) evaluates the functional G on orbitguess
• po_jvp(pb, orbitguess, par, du) evaluates the jacobian dG(orbitguess)⋅du functional at orbitguess on du.
• po_jacobian(pb, orbitguess, par) computes the matrix of the jacobian dG(orbitguess) analytically, based on monodromy matrices. Useful mainly for ODEs.
• po_jacobian!(pb, J, orbitguess, par) same as above but overwrites J inplace.

You can then call po_residual(pb, orbitguess, par) to apply the functional to a guess. Note that you can generate this guess from a function solution using generate_solution or generate_ci_problem.

Allowed types

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

Internal fields

• M::Int64: 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. Default: 0

• flow::BifurcationKit.AbstractFlow: flow::Flow: implements the flow of the Cauchy problem though the structure Flow. Default: Flow()

• ds::Any: ds: vector of the time differences for each shooting, of length M. Defaults to a uniform partition of [0, 1]. Default: diff(LinRange(0, 1, M + 1))

• section::Any: 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 must return 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. Default: nothing

• parallel::Bool: 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). Default: false

• par::Any: par parameters of the model Default: nothing

• lens::Any: lens parameter axis. Default: nothing

• update_section_every_step::UInt64: updates the section every update_section_every_step step during continuation Default: 1

• jacobian::BifurcationKit.AbstractJacobianType: Describes the type of jacobian used in Newton iterations (see below). Default: AutoDiffDense()

Methods

Here are some useful methods you can apply to pb::Shooting:

• get_mesh_size(pb) returns the number M of time slices.
• get_time_slices(pb, x) returns the time slices, either as an N x M matrix or as a BorderedArray.
• getparams(pb), getlens(pb) and setparam(pb, p) give access to the parameters of the model.
• isparallel(pb) returns true if the multiple trajectories are simulated in parallel (threading).
• generate_solution(pb, orbit, period) generates a guess from a function t -> orbit(t).
• POInterpolation(pb, x) returns a function interpolating the periodic orbit x.
• get_periodic_orbit(pb, x, pars) computes the full periodic orbit, mainly for plotting purposes.

Jacobian

jacobian Specify the choice of the linear algorithm, which must belong to (BifurcationKit.AutoDiffMF(), BifurcationKit.MatrixFree(), BifurcationKit.AutoDiffDense(), BifurcationKit.AutoDiffDenseAnalytical(), BifurcationKit.FiniteDifferences(), BifurcationKit.FiniteDifferencesMF()). This is used to select a way of inverting the jacobian dG

1. For `MatrixFree()`, matrix free jacobian, the jacobian is specified by the user in `prob`. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system
2. For `AutoDiffMF()`, we use Automatic Differentiation (AD) to compute the (matrix-free) derivative of `x -> prob(x, p)` using a directional derivative, also called JVP product. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system
3. 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.
4. For `FiniteDifferences()`, 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.
5. For `AutoDiffDenseAnalytical()`. Same as for `AutoDiffDense` but the jacobian is formed using a mix of AD and analytical formula.
6. For `FiniteDifferencesMF()`, use Finite Differences to compute the matrix-free jacobian of `x -> prob(x, p)` using the `δ = 1e-8` which can be passed as an argument.

Simplified constructors

• The first important constructor is the following which is used for branching to periodic orbits from Hopf bifurcation points:

pb = Shooting(M::Int, prob::Union{ODEProblem, EnsembleProblem}, alg; kwargs...)

• A convenient way to build the functional is to use:

pb = Shooting(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 = Shooting(prob::Union{ODEProblem, EnsembleProblem}, alg, M::Int, section; parallel = false, kwargs...)

or

pb = Shooting(prob::Union{ODEProblem, EnsembleProblem}, alg, ds, section; parallel = false, kwargs...)

• The next way is an elaboration of the previous one

pb = Shooting(prob1::Union{ODEProblem, EnsembleProblem}, alg1, prob2::Union{ODEProblem, EnsembleProblem}, alg2, M::Int, section; parallel = false, kwargs...)

or

pb = Shooting(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 PoincareShooting{Tf, Tjac<:BifurcationKit.AbstractJacobianType, Tsection<:SectionPS, Tpar, Tlens} <: BifurcationKit.AbstractPoincareShootingDiscretization

This composite type implements the Poincaré Shooting method to locate periodic orbits / BVP by relying on Poincaré return maps. More details (maths, notations, linear systems) can be found here. The arguments are as described below.

Internal fields

• M::Int64: M: the number of Poincaré sections. If M == 1, then the simple shooting is implemented and the multiple one otherwise. Default: 0

• flow::Any: flow::Flow: implements the flow of the Cauchy problem though the structure Flow. Default: Flow()

• section::SectionPS: 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. Default: SectionPS(M)

• δ::Float64: δ = 1e-8 used to compute the jacobian of the functional by finite differences. Default: 1.0e-8

• parallel::Bool: parallel = false whether the shooting are computed in parallel (threading). Only available through the use of Flows defined by EnsembleProblem. Default: false

• par::Any: par parameters of the model Default: nothing

• lens::Any: lens parameter axis Default: nothing

• update_section_every_step::UInt64: update_section_every_step updates the section every update_section_every_step step during continuation Default: 1

• jacobian::BifurcationKit.AbstractJacobianType: Describes the type of jacobian used in Newton iterations (see below). Default: AutoDiffDenseAnalytical()

Jacobian

jacobian Specify the choice of the linear algorithm, which must belong to (BifurcationKit.AutoDiffMF(), BifurcationKit.MatrixFree(), BifurcationKit.AutoDiffDense(), BifurcationKit.AutoDiffDenseAnalytical(), BifurcationKit.FiniteDifferences(), BifurcationKit.FiniteDifferencesMF()). This is used to select a way of inverting the jacobian dG

1. For `MatrixFree()`, matrix free jacobian, the jacobian is specified by the user in `prob`. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system
2. For `AutoDiffMF()`, we use Automatic Differentiation (AD) to compute the (matrix-free) derivative of `x -> prob(x, p)` using a directional derivative, also called JVP product. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system
3. 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.
4. For `FiniteDifferences()`, 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.
5. For `AutoDiffDenseAnalytical()`. Same as for `AutoDiffDense` but the jacobian is formed using a mix of AD and analytical formula.
6. For `FiniteDifferencesMF()`, use Finite Differences to compute the matrix-free jacobian of `x -> prob(x, p)` using the `δ = 1e-8` which can be passed as an argument.

Simplified constructors

• The first important constructor is the following which is used for branching to periodic orbits from Hopf bifurcation points pb = PoincareShooting(M::Int, prob::Union{ODEProblem, EnsembleProblem}, alg; kwargs...)

• A convenient way is to create a functional is

pb = PoincareShooting(prob::ODEProblem, alg, section; kwargs...)

for simple shooting or

pb = PoincareShooting(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 = PoincareShooting(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 po_residual(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 in the state space and M is the number of Poincaré maps (the coordinates are the projections on the Poincaré sections, hence the - 1). 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 generate_solution.

• po_residual(pb, orbitguess, par) evaluates the functional G on orbitguess
• po_jvp(pb, orbitguess, par, du) evaluates the jacobian dG(orbitguess).du functional at orbitguess on du
• po_jacobian(pb, orbitguess, par) computes the matrix of the jacobian dG(orbitguess) analytically, based on monodromy matrices. Useful mainly for ODEs.
• po_jacobian!(pb, J, orbitguess, par) same as above but overwrites J inplace.

struct TWModel{Tprob, Tu0, TDu0, TD, Tj} <: BifurcationKit.AbstractTravelingWaveDiscretization

This composite type implements a functional for freezing symmetries in order, for example, to compute traveling waves (TW). Note that you can freeze many symmetries, not just one, by passing many Lie generators. TWModel is a discretization: the residual of the frozen system is obtained by wrapping pb into a TravellingWave functional, residual(TravellingWave(pb), x, par), which computes:

                ┌                   ┐
                │ f(x, par) - s⋅∂⋅x │
                │   <x - u₀, ∂⋅u₀>  │
                └                   ┘

Arguments

• prob bifurcation problem with continuous symmetries
• ∂::Tuple 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 = TWModel(prob, ∂, u₀; kw...)

This simplified call handles the case where a single symmetry needs to be frozen.

Useful function

• updatesection!(pb::TWModel, u0) updates the reference solution of the problem using u0.
• nb_constraints(::TWModel) number of constraints (or Lie generators)

Internal fields

• prob_vf::Any: vector field, must be AbstractBifurcationProblem.

• ∂::Any: Infinitesimal generator of symmetries, differential operator.

• u₀::Any: reference solution, we only need one!

• ∂u₀::Any: Default: (∂ * u₀,)

• DAE::Int64: Default: 0

• nc::Int64: [Internal] number of constraints Default: 1

• jacobian::Any: Default: AutoDiff()

Basic eigen solver to compute the stability of the wave based on the eigenvalues of `J + η * ∂`.

Eigen solver to compute the stability of the wave based on the eigenvalues of the GEV, see [documentation](https://bifurcationkit.github.io/BifurcationKitDocs.jl/dev/intro_wave/#Wave-stability).

Specific continuation routine for wave problems.

Arguments

struct DefaultLS <: BifurcationKit.AbstractDirectLinearSolver

This struct is used to provide the backslash operator `. Can be used to solveJ * x = rhs`.

J can be an AbstractMatrix, a linear operator or a ShiftedOperator (possibly wrapping a MassAndJacobian) to solve shifted systems.

Internal fields

• 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 solver is used to test Moore-Penrose continuation. It is a direct solver based on the backslash operator, used to solve J * x = rhs.

J can be an AbstractMatrix, a linear operator or a ShiftedOperator (possibly wrapping a MassAndJacobian) to solve shifted systems.

Internal fields

• 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

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

Linear solver based on gmres from IterativeSolvers.jl. Can be used to solve J * x = rhs.

J can be an AbstractMatrix, a linear operator or a ShiftedOperator (possibly wrapping a MassAndJacobian) to solve shifted systems.

The struct is mutable so that you can modify the preconditioners.

Internal fields

• 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{𝒯, 𝒯l} <: BifurcationKit.AbstractIterativeLinearSolver

Create a linear solver based on linsolve from KrylovKit.jl. Can be used to solve J * x = rhs.

J can be an AbstractMatrix, a linear operator or a ShiftedOperator (possibly wrapping a MassAndJacobian) to solve shifted systems.

The struct is mutable so that you can modify the preconditioners.

Internal fields

• 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 KrylovLS{K, 𝒯l, 𝒯r} <: BifurcationKit.AbstractIterativeLinearSolver

Create a linear solver based on Krylov.jl. Can be used to solve J * x = rhs. You have access to cg, cr, gmres, symmlq, cg_lanczos, cg_lanczos_shift_seq...

J can be an AbstractMatrix, a linear operator or a ShiftedOperator (possibly wrapping a MassAndJacobian) to solve shifted systems.

The struct is mutable so that you can modify the preconditioners.

Example

You can create a Krylov solver with the following code:

KrylovLS(atol = 1e-11, rtol = 1e-8)

Other methods

Look at KrylovLSInplace for a method where the Krylov space is kept in memory

Internal fields

• KrylovAlg::Symbol: Krylov method

• kwargs::Any: Arguments passed to the linear solver

• Pl::Any: Left preconditioner

• Pr::Any: Right preconditioner

mutable struct KrylovLSInplace{F, K, 𝒯l, 𝒯r} <: BifurcationKit.AbstractIterativeLinearSolver

Create an inplace linear solver based on Krylov.jl. Can be used to solve J * x = rhs.

J can be an AbstractMatrix, a linear operator or a ShiftedOperator (possibly wrapping a MassAndJacobian) to solve shifted systems.

The Krylov space is pre-allocated. This is really great for GPU but also for CPU.

The struct is mutable so that you can modify the preconditioners.

Internal fields

• workspace::Any: Can be Krylov.GmresWorkspace for example.

• KrylovAlg::Symbol: Krylov method.

• kwargs::Any: Arguments passed to the linear solver.

• Pl::Any: Left preconditioner.

• Pr::Any: Right preconditioner.

• is_inplace::Bool: Is the linear mapping inplace.

struct DefaultEig{T} <: BifurcationKit.AbstractDirectEigenSolver

The struct DefaultEig is used to provide the eigen method to BifurcationKit.

Internal fields

• which::Any: How do we sort the computed eigenvalues. Default: real

Constructors

Just pass the above fields like DefaultEig(; which = abs)

struct EigArpack{T, Tby, Tw} <: BifurcationKit.AbstractIterativeEigenSolver

Create an eigen solver based on Arpack.jl.

Internal fields

• 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

Constructor

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

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

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

Create an eigen solver based on KrylovKit.jl.

Internal fields

• 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

Constructors

Just pass the above fields like EigKrylovKit(;dim=2)

struct EigArnoldiMethod{T, Tby, Tw, Tkw, vectype} <: BifurcationKit.AbstractIterativeEigenSolver

Internal fields

• 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: Keyword arguments passed to ArnoldiMethod.partialschur.

• 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 ShiftInvert{T, Tls<:BifurcationKit.AbstractLinearSolver, Teig<:AbstractEigenSolver} <: AbstractEigenSolver

Create an eigensolver based on Shift-Invert strategy. Basically, one computes the eigen-elements of (J - σ⋅I)⁻¹.

Internal fields

• sigma::Any: Shift.

• ls::BifurcationKit.AbstractLinearSolver: Linear solver to compute (J - σ⋅I)⁻¹.

• eig::AbstractEigenSolver: Eigen-solver to compute the eigenvalues of (J - σ⋅I)⁻¹.

struct EigenMassMatrix{Tb, Teig<:AbstractEigenSolver} <: AbstractEigenSolver

Create an eigensolver for DAE, Basically a GEV with mass matrix.

Internal fields

• B::Any: Mass matrix

• eig::AbstractEigenSolver: Eigen-solver

gev(l, A, B, nev; kwargs...)

Solve the generalized eigenvalue problem $A x = \lambda B x$ using l.

Arguments

• l::AbstractEigenSolver: the eigensolver configuration (controls sorting via l.which).
• A, B: the matrix (or operators) pair defining $A x = \lambda B x$.
• nev: number of requested eigenvalues.

Generalized eigen solver based on LinearAlgebra.eigen.

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

Internal fields

• eigensolver::Any: Arpack eigensolver.

• 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

GEigArpack(; kw...) where the keyword arguments are forwarded to the inner EigArpack.

struct GEigKrylovKit{T, Tb} <: BifurcationKit.AbstractMFGEigenSolver

Create an generalised eigen solver based on KrylovKit.jl.

Internal fields

• eigensolver::Any: KryloKit eigensolver.

• B::Any: Mass matrix / operator.

struct GEigArnoldiMethod{T, Tb} <: BifurcationKit.AbstractMFGEigenSolver

Create an generalised eigen solver based on ArnoldiMethod.jl.

Internal fields

• eigensolver::Any: Eigensolver.

• B::Any: Mass matrix / operator.

struct FloquetQaD{E<:AbstractEigenSolver} <: BifurcationKit.AbstractFloquetSolver

Computes Floquet multipliers (eigenvalues of the monodromy matrix) for periodic orbit problems using the Shooting method or Finite Differences (Trapeze method).

Method Description

The "Quick and Dirty" (QaD) method computes Floquet multipliers through sequential matrix products along the periodic orbit. This approach has numerical limitation:

• Precision issues: Accuracy degrades for large or small Floquet exponents when using many time sections, due to accumulated errors from repeated matrix multiplications.

Despite precision limitations, the method is sufficient for bifurcation detection in most cases.

Internal fields

• eigsolver::AbstractEigenSolver: eigensolver used to compute the eigenvalues of the monodromy matrix.
• matrix_free::Bool: whether to use a matrix-free linear operator (automatic when eigsolver is not a direct solver).

Implementation Details

• If eigsolver == DefaultEig(), the full monodromy matrix is explicitly formed and its eigenvalues are computed
• Otherwise, a matrix-free formulation is used, applying the monodromy operator via sequential time evolution

struct FloquetGEV{E<:AbstractEigenSolver, Tb} <: BifurcationKit.AbstractFloquetSolver

Computes Floquet exponents for periodic orbits using a Generalized Eigenvalue Problem (GEV) formulation.

Method Description

This method reformulates the Floquet multiplier computation as a large-dimensional generalized eigenvalue problem. The approach is based on a simplified version of the algorithm described in [1].

Performance characteristics:

• Accuracy: More numerically precise than FloquetQaD, especially for extreme Floquet exponents.
• Speed: Slower due to the large dimension of the eigenvalue problem.
• Recommendation: Use FloquetColl() when possible for better performance with similar accuracy.

Fields

• eigsolver::AbstractEigenSolver: Eigensolver used to solve the generalized eigenvalue problem.
• B::Tb: Mass matrix for the generalized eigenvalue formulation.

Constructor

FloquetGEV(eigls::AbstractEigenSolver, ntot::Int, n::Int; array_zeros = zeros)

Arguments:

• eigls: Eigensolver to use.
• ntot: Total dimension of the generalized eigenvalue problem.
• n: State space dimension.
• array_zeros: Function to allocate zero arrays: defaults to zeros but spzeros can be passed for sparse matrices.

Example

# For a 2D system with 30 time sections and 4 collocation points:eigfloquet = FloquetGEV(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.

struct FloquetColl{E<:AbstractEigenSolver, C} <: BifurcationKit.AbstractFloquetSolver

Computes Floquet exponents for periodic orbits discretized with the orthogonal collocation method.

Method Description

This method uses the condensation of parameters technique, originally described in [1] and implemented in AUTO07p, with improvements from [2,3]. The approach efficiently computes Floquet multipliers by exploiting the structure of the collocation discretization.

Performance characteristics:

• Accuracy: High precision, comparable to FloquetGEV.
• Speed: Faster than FloquetGEV by exploiting collocation structure.
• Recommendation: Preferred method for collocation-based periodic orbit problems.

Fields

• eigsolver::AbstractEigenSolver: Which eigen solver. Defaults to DefaultEig.

• cache::Any: Cache, defaults to nothing. It should be set to COPCACHE. When used with COPBLS, it is automatically set up.

• small_n::Bool: Whether to use optimized COP to compute the Floquet exponents. Note that this does not work well for sparse matrices. Defaults to true.

Constructor

FloquetColl(; eigls::AbstractEigenSolver = DefaultEig(), cache = nothing, small_n = true)

Keyword Arguments:

• eigls: Eigensolver to use (defaults to DefaultEig())
• cache: Optional cache for optimization (set to COPCACHE; automatically configured when using COPBLS)
• small_n: Whether to use optimized algorithm for small state dimensions (defaults to true)

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.

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

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.

Internal fields

• solver::Union{Nothing, BifurcationKit.AbstractLinearSolver}: Linear solver used to invert the full matrix.

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

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

Internal fields

• solver::Union{Nothing, BifurcationKit.AbstractLinearSolver}: Linear solver for the Bordering method. Default: nothing

• tol::Real: Tolerance for checking precision. Default: 1.0e-12

• check_precision::Bool: Check precision of the linear solve? Default: true

• k::Int64: Number of recursions to achieve tolerance. Default: 1

• dot::Any: Inner product used by the solver. Default: VI.inner

• norm::Any: Norm used by the solver. Default: VI.norm

Constructors

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

Reference(s)

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.

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

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

Constructor

MatrixFreeBLS(solver, ::Bool)

Internal fields

• solver::Union{Nothing, BifurcationKit.AbstractLinearSolver}: Linear solver for solving the extended linear system.

• use_bordered_array::Bool: Structure used to hold (x, p). If true, this is achieved using BorderedArray. If false, a Vector is used which is analogous to vcat(x, p).

struct MatrixFreeBLSmap{Tj, Ta, Tb, Tc, Td}

Composite type to save the bordered linear system with expression

┌ ┐ │ J a │ │ b' c │ └ ┘

It then solved using Matrix Free algorithm applied to the full operator and not just J as for MatrixFreeBLS

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

It is interesting for solving the linear system associated with Collocation / Trapeze functionals, for example using BorderingBLS(solver = BK.LSFromBLS(), tol = 1e-9, k = 2, check_precision = true)

Internal fields

• solver::Any: Linear solver used to solve the smaller linear systems.

solve(prob::AbstractBifurcationProblem, ::Newton, options::NewtonPar; normN = norm, callback = (;x, fx, J, residual, step, itlinear, options, x0, residuals; kwargs...) -> true, kwargs...)

This is the Newton-Krylov Solver for F(x, p0) = 0 with Jacobian w.r.t. x written J(x, p0) and initial guess x0. It is important to set the linear solver options.linsolver properly depending on your problem. This linear solver is used to solve $J(x, p_0)u = F(x, p_0)$ in the Newton step (the update then reads $x \leftarrow x - u$). 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 information.

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

Optional Arguments

• normN = norm specifies a norm for the convergence criteria
• callback function passed by the user which is called at the end of each iteration. The default one is the following cb_default((x, fx, J, residual, step, itlinear, options, x0, residuals); k...) = true. Can be used to update a preconditioner 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
  • fx current residual
  • J current jacobian
  • residual current norm of the residual
  • step current newton step
  • itlinear number of iterations to solve the linear system
  • options a copy of the argument options passed to newton
  • residuals the history of residuals
  • kwargs kwargs arguments, contain your initial guess x0
  They can be used as callback = (state;k...) -> state.residual<1 for example.
• kwargs arguments passed to the callback. Useful when newton is called from continuation

Output:

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

solve(
    prob::BifurcationKit.AbstractBifurcationProblem,
    defOp::DeflationOperator{Tp, Tdot, T, vectype},
    options::NewtonPar{T, L, E};
    ...
) -> NonLinearSolution{_A, Tprob, Tres} where {_A, Tprob<:(DeflatedProblem{Tprob, _A, _B, _C, _D, Val{:Custom}} where {Tprob<:BifurcationKit.AbstractBifurcationProblem, _A, _B, _C, _D}), Tres<:(Vector)}
solve(
    prob::BifurcationKit.AbstractBifurcationProblem,
    defOp::DeflationOperator{Tp, Tdot, T, vectype},
    options::NewtonPar{T, L, E},
    _linsolver::BifurcationKit.DeflatedProblemCustomLS;
    kwargs...
) -> NonLinearSolution{_A, Tprob, Tres} where {_A, Tprob<:(DeflatedProblem{Tprob, _A, _B, _C, _D, Val{:Custom}} where {Tprob<:BifurcationKit.AbstractBifurcationProblem, _A, _B, _C, _D}), Tres<:(Vector)}

This is the deflated version of the Krylov-Newton Solver for F(x, p0) = 0.

We refer to the regular solve for more information. It penalises the roots saved in defOp.roots. The other arguments are as for solve. See DeflationOperator for more information on defOp.

Arguments

Compared to solve, the only different arguments are

• defOp::DeflationOperator deflation operator
• linsolver linear solver used to invert the Jacobian of the deflated functional.
  • custom solver DeflatedProblemCustomLS() which 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 Matrix of the deflated problem
  • if passed Val(:fullIterative), then a full matrix free method is used for the deflated problem.

This specific Newton-Krylov method first tries to converge to a solution sol0 close the guess x0. It then attempts to converge from the guess x1 while avoiding the previous converged solution close to sol0. This is very handy for branch switching. The method is based on a deflated Newton-Krylov solver.

Arguments

Compared to newton, the only different arguments are

• defOp::DeflationOperator deflation operator
• linsolver linear solver used to invert the Jacobian of the deflated functional.
  • custom solver DeflatedProblemCustomLS() which 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 Matrix of the deflated problem
  • if passed Val(:fullIterative), then a full matrix free method is used for the deflated problem.

newton(
    br::BifurcationKit.AbstractBranchResult,
    ind_bif::Int64;
    normN,
    options,
    start_with_eigen,
    lens2,
    kwargs...
) -> Any

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.newton_options
• 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
• start_with_eigen = false whether to start the Minimally Augmented problem with information from eigen elements.
• kwargs keywords arguments to be passed to the regular Newton-Krylov solver

newton(
    disc::BifurcationKit.AbstractShootingDiscretization,
    orbitguess,
    options::NewtonPar;
    lens,
    δ,
    kwargs...
) -> NonLinearSolution{_A, Tprob, Tres} where {_A, Tprob<:(BifurcationKit.PeriodicOrbitFunctionalSh{Tdisc, _A, _B, Nothing, Nothing} where {Tdisc<:BifurcationKit.AbstractShootingDiscretization, _A, _B}), Tres<:(Vector)}

This is the Newton-Krylov Solver for computing a periodic orbit using the (Standard / Poincaré) Shooting method. Note that the linear solver has to be appropriately set up in options.

Arguments

Similar to newton except that prob is either a Shooting or a PoincareShooting. These two problems have specific options to be tuned, we refer to their link for more information and to the tutorials.

• disc a problem of type <: AbstractShootingDiscretization encoding the shooting functional G.
• orbitguess a guess for the periodic orbit. See Shooting and See PoincareShooting for information regarding the shape of orbitguess.
• options same as for the regular newton method.

Optional argument

jacobian Specify the choice of the linear algorithm, which must belong to (BifurcationKit.AutoDiffMF(), BifurcationKit.MatrixFree(), BifurcationKit.AutoDiffDense(), BifurcationKit.AutoDiffDenseAnalytical(), BifurcationKit.FiniteDifferences(), BifurcationKit.FiniteDifferencesMF()). This is used to select a way of inverting the jacobian dG

1. For `MatrixFree()`, matrix free jacobian, the jacobian is specified by the user in `prob`. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system
2. For `AutoDiffMF()`, we use Automatic Differentiation (AD) to compute the (matrix-free) derivative of `x -> prob(x, p)` using a directional derivative, also called JVP product. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system
3. 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.
4. For `FiniteDifferences()`, 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.
5. For `AutoDiffDenseAnalytical()`. Same as for `AutoDiffDense` but the jacobian is formed using a mix of AD and analytical formula.
6. For `FiniteDifferencesMF()`, use Finite Differences to compute the matrix-free jacobian of `x -> prob(x, p)` using the `δ = 1e-8` which can be passed as an argument.

newton(
    disc::BifurcationKit.AbstractShootingDiscretization,
    orbitguess,
    defOp::DeflationOperator{Tp, Tdot, T, vectype},
    options::NewtonPar{T, S, E};
    lens,
    kwargs...
) -> NonLinearSolution{_A, Tprob, Tres} where {_A, Tprob<:(DeflatedProblem{Tprob, _A, _B, _C, _D, Val{:Custom}} where {Tprob<:(BifurcationKit.PeriodicOrbitFunctionalSh{Tdisc, _A, _B, Nothing, Nothing} where {Tdisc<:BifurcationKit.AbstractShootingDiscretization, _A, _B}), _A, _B, _C, _D}), Tres<:(Vector)}

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 Shooting or a PoincareShooting.

Optional argument

jacobian Specify the choice of the linear algorithm, which must belong to (BifurcationKit.AutoDiffMF(), BifurcationKit.MatrixFree(), BifurcationKit.AutoDiffDense(), BifurcationKit.AutoDiffDenseAnalytical(), BifurcationKit.FiniteDifferences(), BifurcationKit.FiniteDifferencesMF()). This is used to select a way of inverting the jacobian dG

1. For `MatrixFree()`, matrix free jacobian, the jacobian is specified by the user in `prob`. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system
2. For `AutoDiffMF()`, we use Automatic Differentiation (AD) to compute the (matrix-free) derivative of `x -> prob(x, p)` using a directional derivative, also called JVP product. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system
3. 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.
4. For `FiniteDifferences()`, 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.
5. For `AutoDiffDenseAnalytical()`. Same as for `AutoDiffDense` but the jacobian is formed using a mix of AD and analytical formula.
6. For `FiniteDifferencesMF()`, use Finite Differences to compute the matrix-free jacobian of `x -> prob(x, p)` using the `δ = 1e-8` which can be passed as an argument.

Output:

• solution::NonLinearSolution, see NonLinearSolution

newton(
    trap::Trapeze,
    orbitguess,
    options::NewtonPar;
    kwargs...
) -> NonLinearSolution{_A, Tprob, Tres} where {_A, Tprob<:(BifurcationKit.PeriodicOrbitFunctionalTrap{Tdisc, _A, _B, Nothing, Nothing} where {Tdisc<:Trapeze, _A, _B}), Tres<:(Vector)}

Locate a periodic orbit with a Newton solver applied to the finite-difference functional G of the Trapeze discretization, i.e. solve $G(u) = 0$. The returned solution has u[end] = T equal to the period of the orbit.

Arguments:

• trap a problem of type Trapeze encoding the functional G; its jacobian field selects the linear algebra used, see below.
• orbitguess a guess for the periodic orbit. See Trapeze for more details.
• options same as for the regular newton method.

These methods only differ in the linear algebra used to invert the jacobian dG of the functional G (see Trapeze); the discretization is otherwise the same. The value of jacobian must belong to (BifurcationKit.Dense(), BifurcationKit.AutoDiffDense(), BifurcationKit.FullLU(), BifurcationKit.FullMatrixFree(), BifurcationKit.BorderedLU(), BifurcationKit.BorderedMatrixFree(), BifurcationKit.FullSparseInplace(), BifurcationKit.BorderedSparseInplace(), BifurcationKit.AutoDiffMF()).

• 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 right choice when the sparsity pattern can change.
• For jacobian = FullSparseInplace(), this is the same as for FullLU() but the sparse matrix dG is updated inplace. This method allocates much less and, in some cases, is significantly faster than 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, and it is also updated inplace. This option is useful to study ODEs of small dimension.
• For jacobian = AutoDiffDense(), the jacobian is evaluated using automatic differentiation (ForwardDiff).
• For jacobian = BorderedLU(), we take advantage of the bordered shape of dG and invert it with a bordered linear solver based on a LU decomposition of the cyclic matrix.
• For jacobian = BorderedSparseInplace(), this is the same as for BorderedLU() but the cyclic matrix Jc is updated inplace. This method allocates much less and, in some cases, is significantly faster than 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 (given by options.linsolver) is used to invert dG: note that a preconditioner is very likely required here because of the cyclic shape of dG which negatively affects the convergence properties of GMRES.
• For jacobian = BorderedMatrixFree(), a matrix-free linear solver is used as well but only for Jc (see the docs): options.linsolver is then used to invert Jc. These two matrix-free options thus expose different parts of the jacobian dG in order to apply 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 = AutoDiffMF(), the evaluation map of the differential is derived using automatic differentiation. Thus, unlike the previous two cases, the user does not need to pass a matrix-free differential.

newton(
    trap::Trapeze,
    orbitguess,
    defOp::DeflationOperator{Tp, Tdot, T, vectype},
    options::NewtonPar;
    kwargs...
) -> NonLinearSolution{_A, Tprob, Tres} where {_A, Tprob<:(DeflatedProblem{Tprob, Tp, _A, T, _B, Val{:Custom}} where {Tprob<:(BifurcationKit.PeriodicOrbitFunctionalTrap{Tdisc, _A, _B, Nothing, Nothing} where {Tdisc<:Trapeze, _A, _B}), Tp<:Real, _A, T<:Real, _B}), Tres<:(Vector)}

This function is similar to newton 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 algorithm from converging to the equilibrium point.

newton(
    coll::Collocation,
    orbitguess,
    options::NewtonPar;
    kwargs...
) -> NonLinearSolution{_A, Tprob, Tres} where {_A, Tprob<:(BifurcationKit.PeriodicOrbitFunctionalColl{Tdisc, _A, _B, Nothing, Nothing} where {Tdisc<:Collocation, _A, _B}), Tres<:(Vector)}

This is the Newton solver for computing a periodic orbit using orthogonal collocation method. Note that the linear solver has to be appropriately set up in options.

Arguments

Similar to newton except that coll is a Collocation.

• coll a discretization of type <: Collocation encoding the collocation 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 (BifurcationKit.AutoDiffDense(), BifurcationKit.DenseAnalytical(), BifurcationKit.FullSparse(), BifurcationKit.DenseAnalyticalInplace(), BifurcationKit.FullSparseInplace(), BifurcationKit.AutoDiffMF()). 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.
  • For DenseAnalytical() Same as for AutoDiffDense but the jacobian is formed using a mix of AD and analytical formula.

newton(
    coll::Collocation,
    orbitguess,
    defOp::DeflationOperator,
    options::NewtonPar;
    kwargs...
) -> NonLinearSolution{_A, Tprob, Tres} where {_A, Tprob<:(DeflatedProblem{Tprob, Tp, _A, T, _B, Val{:Custom}} where {Tprob<:(BifurcationKit.PeriodicOrbitFunctionalColl{Tdisc, _A, _B, Nothing, Nothing} where {Tdisc<:Collocation, _A, _B}), Tp<:Real, _A, T<:Real, _B}), Tres<:(Vector)}

This function is similar to newton(::Collocation, orbitguess, options; 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.

struct DotTheta{Tdot, Ta}

Internal fields

• 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 weighted dot product $\langle (u_1,p_1), (u_2,p_2)\rangle_\theta = \theta \Re \langle u_1,u_2\rangle +(1-\theta)p_1p_2$ where $u_i\in\mathbb R^N$. The $\Re$ factor is put to ensure a real valued result despite possible complex valued arguments.

continuation(
    prob::BifurcationKit.AbstractBifurcationProblem,
    alg::BifurcationKit.AbstractContinuationAlgorithm,
    contparams::ContinuationPar;
    linear_algo,
    bothside,
    kwargs...
) -> BifurcationKit.DCResult{Tprob, Tbr, Tit, Tsol, Talg} where {Tprob<:BifurcationKit.BVP.BVPBifProblem, Tbr<:(Vector), Tit<:(ContIterable{BifurcationKit.BoundaryValueProblemCont, Tprob, Talg, T, _A, _B, typeof(LinearAlgebra.norm), typeof(BifurcationKit.finalise_default), typeof(BifurcationKit.cb_default), Nothing} where {Tprob<:BifurcationKit.BVP.BVPBifProblem, Talg<:BifurcationKit.AbstractContinuationAlgorithm, T<:Real, _A, _B}), Tsol<:(Vector), Talg<:DefCont}

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::AbstractBifurcationProblem 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::ContinuationPar parameters for continuation. See ContinuationPar

Optional Arguments:

• plot = false whether to plot the solution/branch/spectrum while computing the branch
• bothside = false compute the branches on the two sides of the initial parameter value p0, merge them and return it.
• normC = norm norm used in the nonlinear solves
• filename to save the computed branch during continuation. The identifier .jld2 will be appended to this filename. This requires using JLD2.
• callback_newton callback for newton iterations. See docs of solve. For example, it can be used to change the preconditioners. It can also be used to limit residuals and parameter steps, see cbMaxNorm and cbMaxNormAndΔp
• finalise_solution = (z, tau, step, contResult; kwargs...) -> true Function called at the end of each continuation step. Can be used to alter the continuation procedure (stop it by returning false), save personal data, plot... The notations are z = BorderedArray(x, p) where x (resp. p) is the current solution (resp. parameter value), tau::BorderedArray is the tangent at z, step::Int is the index of the current continuation step and contResult is the current branch but before saving the current state to it! For advanced use:
  • the state state::ContState of the continuation iterator is passed in kwargs. This can be used for testing whether this is called from bisection for locating bifurcation points / events: in_bisection(state) for example. This allows to escape some personal code in this case.
  Note that you can have a better control over the continuation procedure by using an iterator, see Iterator Interface.
  • the iterator iter::ContIterable of the continuation is passed in kwargs.
• verbosity::Int = 0 controls the amount of information printed during the continuation process. Must belong to {0,1,2,3}. In case contparams.newton_options.verbose = false, the following is valid (otherwise the newton iterations are shown). Each case prints more information than the previous one:
  • case 0: print nothing
  • case 1: print basic information about the continuation: used predictor, step size and parameter values
  • case 2: print newton iterations number, stability of solution, detected bifurcations / events
  • case 3: print information during bisection to locate bifurcations / events
• linear_algo set the linear solver for the continuation algorithm alg. For example, PALC needs a linear solver for an enlarged problem (size n+1 instead of n) and one thus needs to tune the one passed in contparams.newton_options.linsolver. This is a convenient argument to thus change the alg linear solver and is used mostly internally. The proper way is to pass directly to alg the correct linear solver.
• kind::AbstractContinuationKind [Internal] flag to describe continuation kind (equilibrium, codim 2, ...). Default = EquilibriumCont()

Output:

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

continuation(
    prob::BifurcationKit.AbstractBifurcationProblem,
    algdc::DefCont,
    contParams::ContinuationPar;
    verbosity,
    plot,
    linear_algo,
    dot_palc,
    callback_newton,
    filename,
    normC,
    kwcont...
) -> BifurcationKit.DCResult{Tprob, Tbr, Tit, Tsol, Talg} where {Tprob<:BifurcationKit.BVP.BVPBifProblem, Tbr<:(Vector), Tit<:(ContIterable{BifurcationKit.BoundaryValueProblemCont, Tprob, Talg, T, _A, _B, typeof(LinearAlgebra.norm), typeof(BifurcationKit.finalise_default), typeof(BifurcationKit.cb_default), Nothing} where {Tprob<:BifurcationKit.BVP.BVPBifProblem, Talg<:BifurcationKit.AbstractContinuationAlgorithm, T<:Real, _A, _B}), Tsol<:(Vector), Talg<:DefCont}

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,
• callback_newton 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},
• normC = norm norm used in the Newton solves,
• dot_palc = (x, y) -> dot(x, y) / length(x), dot product used to define the weighted dot product (resp. norm) $\|(x, p)\|^2_\theta$ in the constraint $N(x, p)$ (see online docs on PALC). This argument can be used to remove the 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::BifurcationKit.AbstractBranchResult,
    ind_bif,
    lens2::Union{typeof(identity), IndexLens, PropertyLens, ComposedFunction};
    ...
) -> Any
continuation(
    br::BifurcationKit.AbstractBranchResult,
    ind_bif,
    lens2::Union{typeof(identity), IndexLens, PropertyLens, ComposedFunction},
    options_cont::ContinuationPar;
    prob,
    start_with_eigen,
    detect_codim2_bifurcation,
    update_minaug_every_step,
    kwargs...
) -> Any

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:

• linsolve_adjoint solver for (J+iω)˟ ⋅sol = rhs or Jᵗ ⋅sol = rhs

• bdlinsolver bordered linear solver for the constraint equation

• bdlinsolver_adjoint bordered linear solver for the constraint equation with top-left block (J-iω)˟ or Jᵗ. Required in the linear solver for the Minimally Augmented Fold/Hopf functional. This option can be used to pass a dedicated linear solver for example with specific preconditioner.

• update_minaug_every_step update vectors a, b in Minimally Formulation every update_minaug_every_step steps

• detect_codim2_bifurcation ∈ {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. Default value = 2.

• start_with_eigen = false whether to start the Minimally Augmented problem with information from eigen elements. If start_with_eigen = false, then:

  • a::Nothing estimate of null vector of J (resp. J-iω) for Fold (resp. Hopf). If nothing is passed, a random vector is used. In case you do not rely on AbstractArray, you should probably pass this.
  • b::Nothing estimate of null vector of Jᵗ (resp. (J-iω)˟) for Fold (resp. Hopf). If nothing is passed, a random vector is used. In case you do not rely on AbstractArray, you should probably pass this.

• 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(
    prob::BifurcationKit.AbstractBifurcationProblem,
    x0,
    par0,
    x1,
    p1::Real,
    alg,
    lens::Union{typeof(identity), IndexLens, PropertyLens, ComposedFunction},
    contParams::ContinuationPar;
    bothside,
    kwargs...
) -> ContResult

[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 pseudo arc length s=0) with parameter par0 and x1 is the second point with parameter set(par0, lens, p1).

continuation(
    br::BifurcationKit.AbstractResult{Tkind<:Union{BifurcationKit.BoundaryValueProblemCont, BifurcationKit.EquilibriumCont}, Tprob},
    ind_bif::Int64;
    ...
) -> Any
continuation(
    br::BifurcationKit.AbstractResult{Tkind<:Union{BifurcationKit.BoundaryValueProblemCont, BifurcationKit.EquilibriumCont}, Tprob},
    ind_bif::Int64,
    options_cont::ContinuationPar;
    ...
) -> Branch
continuation(
    br::BifurcationKit.AbstractResult{Tkind<:Union{BifurcationKit.BoundaryValueProblemCont, BifurcationKit.EquilibriumCont}, Tprob},
    ind_bif::Int64,
    options_cont::ContinuationPar,
    Teigvec::Type{𝒯eigvec};
    alg,
    δp,
    ampfactor,
    use_normal_form,
    bls,
    bls_block,
    nev,
    scaleζ,
    autodiff,
    start_with_eigen,
    usedeflation,
    verbosedeflation,
    max_iter_deflation,
    perturb,
    plot_solution,
    tol_fold,
    kwargs_deflated_newton,
    kwargs...
) -> Any

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 containing the bifurcation point
• ind_bif: index of the bifurcation point in br from which to branch
• options_cont: continuation parameters for the new branch

Optional arguments

• alg = getalg(br) continuation algorithm to be used, default value: getalg(br)
• δp used to specify a specific value for the parameter on the bifurcated branch which is otherwise determined by options_cont.ds. This allows to use a step larger than options_cont.dsmax.
• ampfactor = 1 factor to alter the amplitude of the bifurcated solution. Useful to magnify the bifurcated solution when the bifurcated branch is very steep. Can also be used to select the upper/lower branch in Pitchfork bifurcations. See also use_normal_form below.
• use_normal_form = true. If use_normal_form = true, the normal form is computed as well as its predictor and a guess is automatically formed. If use_normal_form = false, the parameter value p = p₀ + δp and the guess x = x₀ + ampfactor .* e (where e is a vector of the kernel) are used as initial guess. This is useful in case automatic branch switching does not work.
• nev number of eigenvalues to be computed to get the right eigenvector
• scaleζ = norm pass a norm to normalize vectors during normal form computation
• plot_solution change plot solution method in the problem getprob(br)
• usedeflation = false whether to use nonlinear deflation (see Deflated problems) to help finding the guess on the bifurcated
• verbosedeflation print deflated newton iterations
• max_iter_deflation number of newton steps in deflated newton
• perturb = identity which perturbation function to use during deflated newton
• Teigvec = _getvectortype(br) type of the eigenvector. Useful when br was loaded from a file and this information was lost
• kwargs optional arguments to be passed to continuation, the regular continuation one and to get_normal_form.

continuation(
    discPO::BifurcationKit.AbstractShootingDiscretization,
    orbitguess,
    alg::BifurcationKit.AbstractContinuationAlgorithm,
    contParams::ContinuationPar,
    linear_algo::BifurcationKit.AbstractBorderedLinearSolver;
    δ,
    eigsolver,
    record_from_solution,
    plot_solution,
    kwargs...
) -> Any

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 Shooting or a PoincareShooting. 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 (BifurcationKit.AutoDiffMF(), BifurcationKit.MatrixFree(), BifurcationKit.AutoDiffDense(), BifurcationKit.AutoDiffDenseAnalytical(), BifurcationKit.FiniteDifferences(), BifurcationKit.FiniteDifferencesMF()). This is used to select a way of inverting the jacobian dG

1. For `MatrixFree()`, matrix free jacobian, the jacobian is specified by the user in `prob`. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system
2. For `AutoDiffMF()`, we use Automatic Differentiation (AD) to compute the (matrix-free) derivative of `x -> prob(x, p)` using a directional derivative, also called JVP product. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system
3. 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.
4. For `FiniteDifferences()`, 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.
5. For `AutoDiffDenseAnalytical()`. Same as for `AutoDiffDense` but the jacobian is formed using a mix of AD and analytical formula.
6. For `FiniteDifferencesMF()`, use Finite Differences to compute the matrix-free jacobian of `x -> prob(x, p)` using the `δ = 1e-8` which can be passed as an argument.

continuation(
    disc::BifurcationKit.AbstractBoundaryValueDiscretization,
    orbitguess,
    alg::BifurcationKit.AbstractContinuationAlgorithm,
    _contParams::ContinuationPar;
    linear_algo,
    kwargs...
) -> Any

This is the continuation routine for computing a periodic orbit.

Arguments

Similar to continuation except that prob::AbstractBoundaryValueDiscretization.

Optional argument

• linear_algo::AbstractBorderedLinearSolver

jacobian Specify the choice of the linear algorithm, which must belong to (BifurcationKit.AutoDiffMF(), BifurcationKit.MatrixFree(), BifurcationKit.AutoDiffDense(), BifurcationKit.AutoDiffDenseAnalytical(), BifurcationKit.FiniteDifferences(), BifurcationKit.FiniteDifferencesMF()). This is used to select a way of inverting the jacobian dG

1. For `MatrixFree()`, matrix free jacobian, the jacobian is specified by the user in `prob`. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system
2. For `AutoDiffMF()`, we use Automatic Differentiation (AD) to compute the (matrix-free) derivative of `x -> prob(x, p)` using a directional derivative, also called JVP product. This is to be used with an iterative solver (e.g. GMRES) to solve the linear system
3. 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.
4. For `FiniteDifferences()`, 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.
5. For `AutoDiffDenseAnalytical()`. Same as for `AutoDiffDense` but the jacobian is formed using a mix of AD and analytical formula.
6. For `FiniteDifferencesMF()`, use Finite Differences to compute the matrix-free jacobian of `x -> prob(x, p)` using the `δ = 1e-8` which can be passed as an argument.

continuation(
    br::BifurcationKit.AbstractBranchResult,
    ind_bif::Int64,
    _contParams::ContinuationPar,
    disc::BifurcationKit.AbstractBoundaryValueDiscretization;
    bif_prob,
    detailed,
    start_with_eigen,
    use_normal_form,
    nev,
    kwargs...
) -> Branch

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
• disc discretization used to specify the way toc compute the periodic orbit. It can be Trapeze, Collocation, Shooting or PoincareShooting .

Optional arguments

• alg = getalg(br) continuation algorithm used for the bifurcated branch. It is inherited from the branch br.
• δp used to specify the guess for the parameter on the bifurcated branch which otherwise defaults to contParams.ds. This allows to use an initial step larger than contParams.dsmax.
• ampfactor = 1 multiplicative factor to alter the amplitude of the bifurcated solution. Useful to magnify the bifurcated solution when the bifurcated branch is very steep.
• use_normal_form = true whether to use the normal form in order to compute the predictor. When false, ampfactor and δp are used to make a predictor based on the bifurcating eigenvector. Setting use_normal_form = false can be useful when computing the normal form is not possible for example when higher order derivatives are not available.
• usedeflation = false 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 plot_solution and finalise_solution.

continuation(
    br::BifurcationKit.AbstractResult{BifurcationKit.PeriodicOrbitCont, Tprob<:BifurcationKit.AbstractWrapperPeriodicOrbitProblem},
    ind_bif::Int64,
    _contParams::ContinuationPar;
    alg,
    δp,
    ampfactor,
    usedeflation,
    linear_algo,
    detailed,
    prm,
    use_normal_form,
    autodiff_nf,
    kwargs...
) -> Branch

Branch switching at a bifurcation point of periodic orbits (PO) specified by a br::AbstractBranchResult. The functional for computing the PO is getprob(br). A deflated Newton-Krylov solver can be used to improve the branch switching capabilities.

Arguments

• br branch of periodic orbits
• ind_bif index of the branch point
• _contParams continuation parameters, see continuation

Optional arguments

• δp = _contParams.ds 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 alters the amplitude of the bifurcated solution. Useful to magnify the bifurcated solution when the bifurcated branch is very steep.
• usedeflation = false whether to use nonlinear deflation (see Deflated problems) to help finding the guess on the bifurcated branch
• use_normal_form = true if false, the predictor is based on the couple δp, ampfactor.

For normal form

• detailed = true whether to fully compute the normal form or a very simplified version.
• autodiff_nf = true whether to use autodiff in get_normal_form. This can be used in case automatic differentiation is not working as intented.

For continuation

• linear_algo = nothing, same as for continuation; when nothing, a BorderingBLS is built from _contParams.newton_options.linsolver.
• kwargs keywords arguments used for a call to the regular continuation and the ones specific to periodic orbits (POs).

continuation(
    trap::Trapeze,
    orbitguess,
    alg::BifurcationKit.AbstractContinuationAlgorithm,
    _contParams::ContinuationPar;
    1002,
    record_from_solution,
    linear_algo,
    kwargs...
) -> Any

Convenience wrapper around continuation_po to continue a branch of periodic orbits computed with the Trapeze finite-difference functional.

Arguments

• trap a problem of type Trapeze encoding the functional G.
• orbitguess a guess for a first periodic orbit. See Trapeze for more details.
• alg continuation algorithm.
• contParams same as for the regular continuation method.

Keyword arguments

• linear_algo same as in continuation; it defaults to a bordered linear solver based on contParams.newton_options.linsolver.
• record_from_solution function used to record the solution on the branch; it defaults to (u, p; k...) -> (period = u[end],) so that the period is printed along the branch.

These methods only differ in the linear algebra used to invert the jacobian dG of the functional G (see Trapeze); the discretization is otherwise the same. The value of jacobian must belong to (BifurcationKit.Dense(), BifurcationKit.AutoDiffDense(), BifurcationKit.FullLU(), BifurcationKit.FullMatrixFree(), BifurcationKit.BorderedLU(), BifurcationKit.BorderedMatrixFree(), BifurcationKit.FullSparseInplace(), BifurcationKit.BorderedSparseInplace(), BifurcationKit.AutoDiffMF()).

• 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 right choice when the sparsity pattern can change.
• For jacobian = FullSparseInplace(), this is the same as for FullLU() but the sparse matrix dG is updated inplace. This method allocates much less and, in some cases, is significantly faster than 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, and it is also updated inplace. This option is useful to study ODEs of small dimension.
• For jacobian = AutoDiffDense(), the jacobian is evaluated using automatic differentiation (ForwardDiff).
• For jacobian = BorderedLU(), we take advantage of the bordered shape of dG and invert it with a bordered linear solver based on a LU decomposition of the cyclic matrix.
• For jacobian = BorderedSparseInplace(), this is the same as for BorderedLU() but the cyclic matrix Jc is updated inplace. This method allocates much less and, in some cases, is significantly faster than 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 (given by options.linsolver) is used to invert dG: note that a preconditioner is very likely required here because of the cyclic shape of dG which negatively affects the convergence properties of GMRES.
• For jacobian = BorderedMatrixFree(), a matrix-free linear solver is used as well but only for Jc (see the docs): options.linsolver is then used to invert Jc. These two matrix-free options thus expose different parts of the jacobian dG in order to apply 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 = AutoDiffMF(), the evaluation map of the differential is derived using automatic differentiation. Thus, unlike the previous two cases, the user does not need to pass a matrix-free differential.

Note that by default, the method prints the period of the periodic orbit as function of the parameter. This can be changed by providing your record_from_solution argument.

continuation(
    coll::Collocation,
    orbitguess,
    alg::BifurcationKit.AbstractContinuationAlgorithm,
    _contParams::ContinuationPar,
    linear_algo::BifurcationKit.AbstractBorderedLinearSolver;
    δ,
    eigsolver,
    record_from_solution,
    plot_solution,
    kwargs...
) -> Any

This is the continuation method for computing a periodic orbit using an orthogonal collocation method.

Arguments

Similar to continuation except that prob is a Collocation. 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 FloquetColl

Specific continuation routine for wave problems.

Arguments

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

Iterator over the steps of a continuation procedure. Calling continuation is equivalent to iterating over a ContIterable while saving the states in a ContResult. Constructing the iterator manually allows a finer control of the continuation procedure:

iter = ContIterable(prob, alg, contparams; kwargs...)
for state in iter
    println("Continuation step = ", state.step)
end

The first iteration computes the initial point on the branch, its tangent and its stability. Each subsequent iteration performs one continuation step: Newton correction of the predicted point, step size adaptation and computation of the next predictor. Iteration stops when the number of steps reaches contparams.max_steps, when the parameter exits [p_min, p_max], when the corrector fails or when the user requests it (for example via finalise_solution, see continuation).

Fields

• kind::BifurcationKit.AbstractContinuationKind: type of solution computed during continuation, e.g. EquilibriumCont or PeriodicOrbitCont

• prob::Any: bifurcation problem to be continued, see BifurcationProblem

• alg::Any: continuation algorithm, e.g. PALC or MoorePenrose

• contparams::ContinuationPar{T, S, E} where {T, S, E}: continuation parameters, see ContinuationPar

• plot::Bool: whether to plot the branch during continuation

• event::Any: structure for event detection

• normC::Any: norm used in the nonlinear solver to measure the size of the corrections

• finalise_solution::Any: function called at the end of each continuation step, see continuation

• callback_newton::Any: callback called after each newton iteration, see continuation

• verbosity::UInt8: level of verbosity, must belong to {0,1,2,3}

• filename::String: name of the file where the branch is saved if save_to_file = true in ContinuationPar

Useful functions

• getprob(iter) returns the bifurcation problem
• getalg(iter) returns the continuation algorithm
• getcontparams(iter) returns the continuation parameters
• setparam(iter, p::Real) set parameter with lens iter.prob.lens to p
• is_event_active(iter) whether the event detection is active
• compute_eigenelements(iter) whether to compute eigen elements
• save_eigenvectors(iter) whether to save eigen vectors
• isindomain(iter, p) whether p is in the domain [pmin, pmax]
• is_on_boundary(iter, p) whether p is on the boundary {pmin, pmax}
• length(iter) the maximum number of continuation steps

More information is available on the website

mutable struct ContState{Tv, T, Teigvals, Teigvec, Tcb} <: BifurcationKit.AbstractContinuationState{Tv}

Mutable structure holding the current state of the continuation procedure: the current solution z, its tangent τ, the predictor z_pred and the current step size ds. It is created internally by the iterator ContIterable and updated at each continuation step (Newton correction, eigen-solver, step size control, predictor). It is also passed as a keyword argument to user callbacks like finalise_solution (see continuation).

Internal fields

• z_pred::Any: predictor, i.e. the initial guess for the next Newton correction.

• τ::Any: tangent to the curve at the current solution.

• z::Any: current solution.

• z_old::Any: previous solution on the branch.

• converged::Bool: whether the newton correction converged at the current continuation step.

• itnewton::Int64: Number of newton iteration (in corrector).

• itlinear::Int64: number of linear iteration (in newton corrector).

• step::Int64: current continuation step.

• ds::Any: current step size.

• stopcontinuation::Bool: boolean to stop continuation.

• stepsizecontrol::Bool: perform step size adaptation.

• n_unstable::Tuple{Int64, Int64}: number of unstable eigenvalues (current, previous).

• n_imag::Tuple{Int64, Int64}: number of unstable complex eigenvalues (current, previous).

• convergedEig::Bool: boolean for eigen solver computation.

• eigvals::Any: current eigenvalues.

• eigvecs::Any: current eigenvectors.

• eventValue::Any: current and previous values of the event function.

• in_bisection::Bool: whether the state is in bisection for locating special points.

Useful functions

• copy(state) returns a copy of state.
• copyto!(dest, state) copy state into dest.
• getsolution(state) returns the current solution (x, p).
• gettangent(state) return the tangent at the current solution.
• getpredictor(state) return the predictor at the current solution.
• getx(state) returns the x component of the current solution.
• getp(state) returns the p component of the current solution.
• get_previous_solution(state) returns the previous solution (x, p).
• getpreviousx(state) returns the x component of the previous solution.
• getpreviousp(state) returns the p component of the previous solution.
• is_stable(state) whether the current state is stable.
• in_bisection(state) whether the state is in bisection for locating special points.
• getparams(iter, state) return the current parameter set.

Natural continuation algorithm. The predictor is the constant predictor and the parameter is incremented by ContinuationPar().ds at each continuation step. The corrector is the regular newton algorithm.

Constructor(s)

Natural()

Secant Tangent predictor

Bordered Tangent predictor

Polynomial Tangent predictor

Internal fields

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

Can be used like

PALC(tangent = Polynomial(Bordered(), 2, 6, rand(1)))

Multiple Tangent continuation algorithm.

The predictor is designed [Uecker2014] to avoid spurious branch switching and pass singular points especially in PDE where branch point density can be quite high. It is called pmcont in pde2path.

Internal fields

• 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, α, nb)

Multiple(x0, α, nb)

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

Pseudo-arclength continuation algorithm.

Additional information is available on the website.

Internal fields

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

• θ::Any: θ is a parameter in the arclength constraint. It is very important to tune it. 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. Default: 0.5

• _bothside::Bool: [internal], not yet used. Default: false

• bls::BifurcationKit.AbstractLinearSolver: Bordered linear solver used to invert the jacobian of the bordered problem during newton iterations. It is also used to compute the tangent for the predictor Bordered(). Default: MatrixBLS()

• dotθ::Any: dotθ = 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, ...) or when a specific (FEM) dot product is provided. Default: DotTheta()

struct AutoSwitch{Talg, T} <: BifurcationKit.AbstractContinuationAlgorithm

Continuation algorithm which switches automatically between Natural continuation and PALC (or other if specified) depending on the stiffness of the branch being continued. The formula for switching from PALC to NATURAL is:

(1 - θ) * abs(τ.p) > tol_param

Internal fields

• alg::Any: Continuation algorithm to switch to when Natural is discarded. Typically PALC()

• tol_param::Any: tolerance for switching to PALC(), default value = 1//2

Constructor(s)

AutoSwitch(;alg = PALC(tangent = Bordered()), tol_param = 1//2)

struct MoorePenrose{T, Tls<:BifurcationKit.AbstractLinearSolver} <: BifurcationKit.AbstractContinuationAlgorithm

Moore-Penrose continuation algorithm. The tangent used by the predictor is computed by predictor (default PALC()). The corrector solves the extended system (F(x,p)=0, τ⋅(x-x₀,p-p₀)=0) using a Moore-Penrose inverse.

Available linear solvers for method:

• BifurcationKit.direct (default): solves [J dFdp] ⋅ Δ = rhs directly via \.
• BifurcationKit.pInv: solves using LinearAlgebra.pinv.
• BifurcationKit.iterative: uses a bordered linear solver ls (e.g., MatrixBLS()) for the extended system.

Constructors

MoorePenrose(;predictor = PALC(), method = direct, ls = nothing)

Internal fields

• predictor::Any: Tangent predictor, for example PALC(), Natural(), etc.

• method::MoorePenroseLS: Moore Penrose linear solver. Can be BifurcationKit.direct, BifurcationKit.pInv or BifurcationKit.iterative.

• ls::BifurcationKit.AbstractLinearSolver: (Bordered) linear solver used to invert the jacobian of the (bordered) problem during moore-penrose iterations. It is also used to compute the tangent for the predictor Bordered().

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

Structure which holds the parameters specific to Deflated continuation.

Internal fields

• deflation_operator::Any: Deflation operator, ::DeflationOperator Default: nothing

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

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

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

• max_iter_defop::Int64: maximum number of deflated Newton iterations Default: 5

• perturb_solution::Any: perturb function Default: _perturbSolution

• accept_solution::Any: accept (solution) function Default: _acceptSolution

• update_deflation_op::Any: function to update the deflation operator, ie pushing new solutions Default: _updateDeflationOp

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

struct DiscreteEvent{Tcb, Tl, Tf, Td} <: 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.

Internal fields

• 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 occurring in the first component. You can use labels = ("hopf",) or labels = ("hopf", "fold"). You must have labels::Union{Nothing, NTuple{N, String}}.

• finaliser::Any: Finaliser function

• data::Any: Place to store some personal data

struct ContinuousEvent{Tcb, Tl, T, Tf, Td} <: 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.

Internal fields

• nb::Int64: Number of events, i.e. 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 at 1.0 and at -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.

• finaliser::Any: Finaliser function.

• data::Any: Place to store some personal data.

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

Multiple events can be chained together to form a SetOfEvents. A SetOfEvents is constructed by passing to the constructor ContinuousEvent, DiscreteEvent or other SetOfEvents instances:

SetOfEvents(cb1, cb2, cb3)

Example

 BifurcationKit.SetOfEvents(BK.FoldDetectEvent, BK.BifDetectEvent)

You can pass as many events as you like.

Internal fields

• eventC::Tuple: Continuous event

• eventD::Tuple: Discrete event

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. A PairOfEvents is constructed by passing to the constructor a ContinuousEvent and a DiscreteEvent:

PairOfEvents(contEvent, discreteEvent)

Internal fields

• eventC::BifurcationKit.AbstractContinuousEvent: Continuous event

• eventD::BifurcationKit.AbstractDiscreteEvent: Discrete event

SaveAtEvent(
    positions::Tuple{Vararg{T, N}} where {N, T};
    use_newton
) -> ContinuousEvent{Tcb, Tl, Int64} where {Tcb, Tl}

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

The keyword use_newton controls how the saved state is finalised:

• use_newton = false (default): the event point is simply the last continuation state recorded when the parameter crossed the requested value p; it is saved without further refinement, so its accuracy is that of the continuation step / event location;
• use_newton = true: once the event has been located around a target value p, the event point is refined by running a Newton solve of F(x, p) = 0 at the exact parameter p, starting from the event state event_point.x and using the newton_options of the continuation. If it converges, the state event_point.x and parameter event_point.param are updated to the refined solution and event_point.precision is set to the Newton tolerance. Use this when you need an accurate equilibrium located exactly at the requested parameter values.

`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 PALC(tangent = Bordered()) as continuation algorithm. To use it, pass event = FoldDetectEvent to continuation.

`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 on the website at Detection of bifurcation points of Equilibria.

continuation(
    br::BifurcationKit.AbstractResult{Tkind<:Union{BifurcationKit.BoundaryValueProblemCont, BifurcationKit.EquilibriumCont}, Tprob},
    ind_bif::Int64;
    ...
) -> Any
continuation(
    br::BifurcationKit.AbstractResult{Tkind<:Union{BifurcationKit.BoundaryValueProblemCont, BifurcationKit.EquilibriumCont}, Tprob},
    ind_bif::Int64,
    options_cont::ContinuationPar;
    ...
) -> Branch
continuation(
    br::BifurcationKit.AbstractResult{Tkind<:Union{BifurcationKit.BoundaryValueProblemCont, BifurcationKit.EquilibriumCont}, Tprob},
    ind_bif::Int64,
    options_cont::ContinuationPar,
    Teigvec::Type{𝒯eigvec};
    alg,
    δp,
    ampfactor,
    use_normal_form,
    bls,
    bls_block,
    nev,
    scaleζ,
    autodiff,
    start_with_eigen,
    usedeflation,
    verbosedeflation,
    max_iter_deflation,
    perturb,
    plot_solution,
    tol_fold,
    kwargs_deflated_newton,
    kwargs...
) -> Any

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 containing the bifurcation point
• ind_bif: index of the bifurcation point in br from which to branch
• options_cont: continuation parameters for the new branch

Optional arguments

• alg = getalg(br) continuation algorithm to be used, default value: getalg(br)
• δp used to specify a specific value for the parameter on the bifurcated branch which is otherwise determined by options_cont.ds. This allows to use a step larger than options_cont.dsmax.
• ampfactor = 1 factor to alter the amplitude of the bifurcated solution. Useful to magnify the bifurcated solution when the bifurcated branch is very steep. Can also be used to select the upper/lower branch in Pitchfork bifurcations. See also use_normal_form below.
• use_normal_form = true. If use_normal_form = true, the normal form is computed as well as its predictor and a guess is automatically formed. If use_normal_form = false, the parameter value p = p₀ + δp and the guess x = x₀ + ampfactor .* e (where e is a vector of the kernel) are used as initial guess. This is useful in case automatic branch switching does not work.
• nev number of eigenvalues to be computed to get the right eigenvector
• scaleζ = norm pass a norm to normalize vectors during normal form computation
• plot_solution change plot solution method in the problem getprob(br)
• usedeflation = false whether to use nonlinear deflation (see Deflated problems) to help finding the guess on the bifurcated
• verbosedeflation print deflated newton iterations
• max_iter_deflation number of newton steps in deflated newton
• perturb = identity which perturbation function to use during deflated newton
• Teigvec = _getvectortype(br) type of the eigenvector. Useful when br was loaded from a file and this information was lost
• kwargs optional arguments to be passed to continuation, the regular continuation one and to get_normal_form.

multicontinuation(
    br::BifurcationKit.AbstractBranchResult,
    ind_bif::Int64;
    ...
) -> Union{Nothing, Vector}
multicontinuation(
    br::BifurcationKit.AbstractBranchResult,
    ind_bif::Int64,
    options_cont::ContinuationPar;
    ...
) -> Union{Nothing, Vector}
multicontinuation(
    br::BifurcationKit.AbstractBranchResult,
    ind_bif::Int64,
    options_cont::ContinuationPar,
    Teigvec::Type{𝒯eigvec};
    δp,
    ampfactor,
    nev,
    ζs,
    scaleζ,
    autodiff,
    bls_block,
    start_with_eigen,
    verbosedeflation,
    plot_solution,
    kwargs...
) -> Union{Nothing, Vector}

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
• options_cont options for the call to continuation

Optional arguments

• alg = getalg(br) continuation algorithm to be used, default value: getalg(br)
• δp used to specify a particular guess for the parameter on the bifurcated branch which is otherwise determined by options_cont.ds. This allows to use a step larger than options_cont.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 = false 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
• perturb_guess = 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::BifurcationKit.AbstractBranchResult,
    bpnf::BifurcationKit.NdBranchPoint,
    defOpm::DeflationOperator,
    defOpp::DeflationOperator;
    ...
) -> Vector
multicontinuation(
    br::BifurcationKit.AbstractBranchResult,
    bpnf::BifurcationKit.NdBranchPoint,
    defOpm::DeflationOperator,
    defOpp::DeflationOperator,
    options_cont::ContinuationPar;
    alg,
    δp,
    verbosedeflation,
    max_iter_deflation,
    plot_solution,
    Teigvec,
    kwargs...
) -> Vector

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 non-trivial branches before/after the bifurcation points. The points are located in defOpm.roots and defOpp.roots. Note that we only continue from the second points in the roots vectors, the first one is meant to be the trivial branch.

The rest is as the regular multicontinuation function.

bifurcationdiagram(
    prob::BifurcationKit.AbstractBifurcationProblem,
    alg::BifurcationKit.AbstractContinuationAlgorithm,
    level::Int64,
    options;
    linear_algo,
    kwargs...
) -> BifDiagNode{Tγ, Vector{BifDiagNode}} where Tγ<:BifurcationKit.AbstractBranchResult

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.

bifurcationdiagram!(
    prob::BifurcationKit.AbstractBifurcationProblem,
    node::BifDiagNode,
    maxlevel::Int64,
    options;
    code,
    halfbranch,
    verbosediagram,
    kwargs...
) -> BifDiagNode

Similar to bifurcationdiagram but you pass a previously computed node from which you want to further compute the bifurcated branches. It is usually used with node = get_branch(diagram, code) from a previously computed bifurcation diagram.

Arguments

• node::BifDiagNode a node in the bifurcation diagram
• maxlevel required maximal level of recursion.
• options = (x, p, level; k...) -> 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.
• verbosediagram verbose specific to bifurcation diagram. Print information about the branches as they are being computed.
• kwargs optional arguments as for continuation but also for the different versions listed in Continuation.

get_branch(diagram::BifDiagNode, code) -> BifDiagNode

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

get_branches_from_BP(
    diagram::BifDiagNode,
    indbif::Int64
) -> Vector{BifDiagNode}

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

getperiod(
    ::BifurcationKit.AbstractBoundaryValueDiscretization,
    x
) -> Any
getperiod(
    ::BifurcationKit.AbstractBoundaryValueDiscretization,
    x,
    par
) -> Any

Compute the period of the periodic orbit associated to x.

getperiod(prob::Trapeze, x, p) -> Any

Compute the period of the periodic orbit associated to x.

getperiod(psh::PoincareShooting, x_bar, par) -> Any

Compute the period of the periodic orbit associated to x_bar.

struct SectionSS{Tn} <: 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 Shooting. The hyperplane is defined by a point center and a normal and is defined by

\[0 = < normal, x - center >\]

Internal fields

• normal::Any: Normal to define hyperplane.

• center::Any: Representative point on hyperplane.

Constructor(s)

SectionSS(normals, centers) which only applies to the first normal (resp. center) in the list of normals (resp. centers).

struct SectionPS{Tn, Tc, Tnb, Tcb, Tr} <: BifurcationKit.AbstractSection

This composite type (named for SectionPoincaréShooting) encodes a type of Poincaré sections implemented by hyperplanes. It can be used in conjunction with PoincareShooting. Each hyperplane is defined par a point (one example in centers) and a normal (one example in normals). See [1] for more details.

Internal fields

• M::Int64: number of hyperplanes

• normals::Any: normals to define hyperplanes

• centers::Any: representative point on each hyperplane

• indices::Vector{Int64}: indices to be removed in the operator Ek

• normals_bar::Any: Projected normals

• centers_bar::Any: Projected centers

• radius::Any

Constructor(s)

SectionPS(normals, centers)

Ref(s)

[1] J., M. Net, B. Garcı́a-Archilla, and C. Simó. “Newton–Krylov Continuation of Periodic Orbits for Navier–Stokes Flows.” Journal of Computational Physics 201, no. 1 (November 20, 2004): 13–33. https://doi.org/10.1016/j.jcp.2004.04.018.

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

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.

struct Flow{TF, Tf, Tts, Tff, Td, Tad, Tse, TR01, TR11, TR20, TR30, Tfs, Tcb, Tδ} <: BifurcationKit.AbstractFlow

Structure to encode the flow associated to a Cauchy problem dx/dt = F(x, p).

Internal fields

• 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

• jvp::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

• vjp::Any: The adjoint differential vjpflow of the flow w.r.t. x, (x, p, dx, t) -> vjpflow(x, p, dx, t). One important thing is that we require vjpflow(x, p, dx, t) to return a Named Tuple: (t = t, u = flow(x, p, t), du = vjpflow(x, p, dx, t)), the last component being the value of the derivative of the flow. Default: nothing

• jvpSerial::Any: [Optional] Serial version of dflow. Used internally by parallel multiple shooting. Please use nothing as default. Default: nothing

• R01::Any: [Optional] Derivatives of the flow with respect to the parameter lens. R01(x, pars, t, lens, p) returns ∂ₚφ(x, p, t), the derivative of the flow map with respect to the parameter lens evaluated at p, as a vector of the size of x. It is used by the Poincaré return map and the normal forms. Default: nothing

• R11::Any: [Optional] Derivatives of the flow with respect to the parameter lens. R11(x, pars, dx, t, lens, p) returns ∂ₚ[dφ(x, p, t)⋅dx], the mixed derivative of the JVP with respect to the parameter, as a vector of the size of x. It is used by the Poincaré return map and the normal forms. Default: nothing

• R20::Any: [Optional] Higher-order differentials of the flow with respect to x. R20(x, pars, h1, h2, t) returns d²φ(x, p, t)(h1, h2), the second differential of the flow map applied to h1, h2. Used by the normal forms. Default: nothing

• R30::Any: [Optional] Higher-order differentials of the flow with respect to x. R30(x, pars, h1, h2, h3, t) returns d³φ(x, p, t)(h1, h2, h3), the third differential of the flow map applied to h1, h2, h3. Used by the normal forms. Default: nothing

• flowSerial::Any: [Internal] Serial version of the flow Default: nothing

• callback::Any: [Internal] Store possible callback Default: nothing

• delta::Any: [Internal] Default: 1.0e-8

Simplified constructor(s)

We provide a simple constructor where you only pass the vector field F, the flow ϕ and its differential dϕ:

fl = Flow(F, ϕ, dϕ)

Simplified constructors for DifferentialEquations.jl

These are some simple constructors for which you only have to pass a prob::ODEProblem or prob::EnsembleProblem (for parallel computation) from DifferentialEquations.jl and an ODE time stepper like Tsit5(). Hence, you can do for example

fl = Flow(prob, Tsit5(); kwargs...)

where kwargs is passed to SciMLBase::solve. If your vector field depends on parameters p, you can define a Flow using

fl = Flow(prob, Tsit5(); kwargs...)

Finally, you can pass two ODEProblem where the second one is used to compute the variational equation:

fl = Flow(prob1::ODEProblem, alg1, prob2::ODEProblem, alg2; kwargs...)

get_normal_form(
    prob::BifurcationKit.AbstractBifurcationProblem,
    br::BifurcationKit.AbstractBranchResult,
    id_bif::Int64;
    ...
) -> Any
get_normal_form(
    prob::BifurcationKit.AbstractBifurcationProblem,
    br::BifurcationKit.AbstractBranchResult,
    id_bif::Int64,
    Teigvec::Type{𝒯eigvec};
    nev,
    verbose,
    lens,
    detailed,
    autodiff,
    scaleζ,
    ζs,
    ζs_ad,
    bls,
    bls_adjoint,
    bls_block,
    bls_block_adjoint,
    start_with_eigen
) -> Any

Compute the reduced equation / 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 the jacobian at the bifurcation point. Useful for enforcing the kernel basis used for the normal form.
• lens::Lens specify which parameter to take the partial derivative ∂pF
• scaleζ function to normalize the kernel basis. Indeed, the kernel vectors are normalized using norm, the normal form coefficients can be super small and can imped its analysis. Using scaleζ = norminf can help sometimes.
• autodiff = true whether to use ForwardDiff for the differentiations. Only used for the codim 2 normal forms (e.g. Bogdanov-Takens (BT), Zero-Hopf (ZH), Hopf-Hopf (HH) points).
• detailed = Val(true) whether to compute only a simplified normal form when only basic information is required. This can be useful is cases the computation is "long", for example for a Bogdanov-Takens point.
• bls = MatrixBLS() specify bordered linear solver. Needed to compute the reduced equation Taylor expansion of Branch/BT points. Indeed, it is required to solve L⋅u = rhs where L is the jacobian at the bifurcation point, L is thus singular and we rely on a bordered linear solver to solve this system.
• bls_block = bls specify bordered linear solver when the border has dimension > 1 (1 for bls). (see bls option above).
• start_with_eigen = Val(true) whether to compute the basis of the kernel (the eigenvectors) using the eigensolver (Val(true)) or using a bordered linear system (Val(false)). The latter can be more robust for large scale problems where the eigensolver fails. It requires bls, bls_adjoint and possibly bls_block to be provided.

Available method(s)

You can directly call

get_normal_form(br, ind_bif ; kwargs...)

which is a shortcut for get_normal_form(getprob(br), br, ind_bif ; kwargs...).

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

References

[1] Golubitsky, Martin, and David G Schaeffer. Singularities and Groups in Bifurcation Theory. Springer-Verlag, 1985. http://books.google.com/books?id=rrg-AQAAIAAJ.

[2] Kielhöfer, Hansjörg. Bifurcation Theory: An Introduction with Applications to PDEs. Applied Mathematical Sciences 156. Springer, 2003. https://doi.org/10.1007/978-1-4614-0502-3.

get_normal_form(
    wrap::BifurcationKit.AbstractWrapperPeriodicOrbitProblem,
    br::BifurcationKit.AbstractResult{<:BifurcationKit.PeriodicOrbitCont},
    id_bif::Int64;
    ...
) -> Any
get_normal_form(
    wrap::BifurcationKit.AbstractWrapperPeriodicOrbitProblem,
    br::BifurcationKit.AbstractResult{<:BifurcationKit.PeriodicOrbitCont},
    id_bif::Int64,
    Teigvec::Type{𝒯eigvec};
    nev,
    verbose,
    ζs,
    lens,
    scaleζ,
    autodiff,
    δ,
    k...
) -> Any

Compute the normal form (NF) of bifurcations of periodic orbits. We detail the additional keyword arguments specific to periodic orbits.

Optional arguments

• prm = false compute the normal form using Poincaré return map (PRM). If false, use the Iooss normal form.
• nev = length(eigenvalsfrombif(br, id_bif)),
• verbose = false,
• lens = getlens(br),
• Teigvec = _getvectortype(br) type of the eigenvectors (can be useful for GPU)
• scaleζ = norm, scale the eigenvector
• autodiff = false use autodiff or finite differences in some part of the normal form computation
• detailed = true whether to compute only a simplified normal form when only basic information is required. This can be useful is cases the computation is long.
• δ = getdelta(wrap) delta used for derivatives based on finite differences.

Notes

For collocation, the default method to compute the NF of Period-doubling and Neimark-Sacker bifurcations is Iooss' one [1].

References

[1] Iooss, "Global Characterization of the Normal Form for a Vector Field near a Closed Orbit.", 1988