Library

options = NewtonPar(tol = 1e-4,...)

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

Arguments (with default values):

• tol = 1e-10: absolute tolerance for F(x)
• maxIter = 50: number of Newton iterations
• verbose = false: display Newton iterations?
• linsolver = DefaultLS(): linear solver, must be <: AbstractLinearSolver
• eigsolver = DefaultEig(): eigen solver, must be <: AbstractEigenSolver

Arguments only used in newtonPALC

• linesearch = false: use line search algorithm
• alpha = 1.0: alpha (damping) parameter for line search algorithm
• almin = 0.001: minimal vslue of the damping alpha

options = ContinuationPar(dsmin = 1e-4,...)

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

Arguments

• dsmin, dsmax are the minimum, maximum arclength allowed value. It controls the density of points in the computed branch of solutions.
• ds is the initial arclength.
• theta is a parameter in the arclength constraint. It is very important to tune it. See the docs of continuation.
• pMin, pMax allowed parameter range for p
• maxSteps maximum number of continuation steps
• newtonOptions::NewtonPar: options for the Newton algorithm
• saveToFile = false: save to file. A name is automatically generated.
• saveSolEveryStep::Int64 = 0 at which continuation steps do we save the current solution`
• plotEveryStep = 3

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

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

Handling bifurcation detection

• precisionStability = 1e-10 lower bound on the real part of the eigenvalues to test for stability of equilibria and periodic orbits
• detectFold = true detect Fold bifurcations? It is a useful option although the detection of Fold is cheap. Indeed, it may happen that there is a lot of Fold points and this can saturate the memory in memory limited devices (e.g. on GPU)
• detectBifurcation::Int ∈ {0, 1, 2, 3} If set to 0, nothing is done. If set to 0, the eigen-elements are computed. If set to 2, bifurcation are detected along the continuation run, but not located precisely. If set to 3, a bisection algorithm is used to locate the bifurcations (slower). The possibility to switch off detection is a useful option. Indeed, it may happen that there are a lot of bifurcation points and this can saturate the memory of memory limited devices (e.g. on GPU)
• dsminBisection dsmin for the bisection algorithm for locating bifurcation points
• nInversion number of sign inversions in bisection algorithm
• maxBisectionSteps maximum number of bisection steps
• tolBisectionEigenvalue 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.
• thetaMin = 1.0e-3 minimum value of theta
• doArcLengthScaling trigger further adaptation of theta

Misc

• finDiffEps::T = 1e-9 ε used in finite differences computations

struct ContResult{Ta, Teigvals, Teigvec, Biftype, Foldtype, Ts, Tfunc, Tpar, Tl<:Setfield.Lens} <: BifurcationKit.BranchResult

Structure which holds the results after a call to continuation.

You can see the propertynames of a result by using propertynames(::ContResult) or by typing br. + TAB where br::ContResult.

Fields

• branch::StructArrays.StructArray{Ta,N,C,I} where I where C<:Union{Tuple, NamedTuple} where N where Ta

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

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

• eig::Array{NamedTuple{(:eigenvals, :eigenvec, :step),Tuple{Teigvals,Teigvec,Int64}},1} where Teigvec where Teigvals

  A vector with eigen-elements at each continuation step.

• foldpoint::Array{Foldtype,1} where Foldtype

  A vector holding the set of detected fold points. See GenericBifPoint for a description of the fields.

• sol::Any

  Vector of solutions sampled along the branch. This is set by the argument saveSolEveryNsteps::Int64 (default 0) in ContinuationPar.

• contparams::ContinuationPar

  The parameters used for the call to continuation which produced this branch.

• type::Symbol

  Type of solutions computed in this branch. Default: :Equilibrium

• functional::Any

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

• params::Any

  Parameters passed to continuation and used in the equation F(x, par) = 0. Default: nothing

• param_lens::Setfield.Lens

  Parameter axis used for computing the branch

• bifpoint::Array{Biftype,1} where Biftype

  A vector holding the set of detectedxbifurcation points (other than fold points). See GenericBifPoint for a description of the fields.

Associated methods

• length(br) number of the continuation steps
• eigenvals(br, ind) returns the eigenvalues for the ind-th continuation step
• eigenvec(br, ind, indev) returns the indev-th eigenvector for the ind-th continuation step
• br[k+1] gives information about the k-th step

DeflationOperator. It is used 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 ; the algorithm works very well despite its simplicity. You can use DeflationOperator to define a function M(u) used to find, with Newton iterations, the zeros of the following function $F(u) \cdot Π_i(dot(u - roots_i, u - roots_i)^{-p} + shift) := F(u) \cdot M(u)$. The fields of the struct DeflationOperator are as follows:

• power p
• dot function, this function has to be bilinear and symmetric for the linear solver to work well
• shift
• roots

The deflation operator is is $M(u) = \prod_{i=1}^{n_{roots}}(shift + norm(u-roots_i)^{-p})$ where $norm(u) = dot(u,u)$.

Given defOp::DeflationOperator, one can access its roots as defOp[n] as a shortcut for defOp.roots[n]. 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)

pb = DeflatedProblem(F, J, M::DeflationOperator)

This creates a deflated problem $M(u) \cdot F(u) = 0$ where M is a DeflationOperator which encodes the penalization term. J is the jacobian of F. Can be used to call newton and continuation.

pb = PeriodicOrbitTrapProblem(F, J, ϕ, xπ, M::Int)

This composite type implements Finite Differences based on a Trapezoidal rule to locate periodic orbits. More details (maths, notations, linear systems) can be found here. The arguments are as follows

• F(x,p) vector field
• J is the 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ᵗ = nothing jacobian tranpose of F (optional), useful for continuation of Fold of periodic orbits. it should not be passed in case the jacobian is a (sparse) matrix as it is computed internally, and it would be computed twice in that case.
• d2F = nothing second derivative of F (optional), useful for continuation of Fold of periodic orbits. It has the definition d2F(x,p,dx1,dx2).`
• ϕ used to set a section for the phase constraint equation
• xπ used in the section for the phase constraint equation
• M::Int number of time slices
• linsolver: = DefaultLS() linear solver for each time slice, i.e. to solve J⋅sol = rhs. This is only needed for the computation of the Floquet multipliers.
• isinplace::Bool whether F and J are inplace functions (Experimental). In this case, the functions F and J must have the following definitions (o, x, p) -> F(o, x, p) and (o, x, p, dx) -> J(o, x, p, dx).
• ongpu::Bool whether the computation takes place on the gpu (Experimental)

The scheme is as follows. We first consider a partition of $[0,1]$ given by $0<s_0<\cdots<s_m=1$ and one looks for T = x[end] such that

$\left(x_{i} - x_{i-1}\right) - \frac{T\cdot h_i}{2} \left(F(x_{i}) + F(x_{i-1})\right) = 0,\ i=1,\cdots,m-1$

with $u_{0} := u_{m-1}$ and the periodicity condition $u_{m} - u_{1} = 0$ and

where $h_1 = s_i-s_{i-1}$. Finally, the phase of the periodic orbit is constrained by using a section (but you could use your own)

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

Orbit guess

You will see below that you can evaluate the residual of the functional (and other things) by calling pb(orbitguess, p) on an orbit guess orbitguess. Note that orbitguess must be 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]$.

Functional

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

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

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

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

• flow::Flow: implements the flow of the Cauchy problem though the structure Flow.
• ds: vector of time differences for each shooting. Its length is written M. If M==1, then the simple shooting is implemented and the multiple one otherwise.
• section: implements a phase condition. The evaluation section(x, T) must return a scalar number where x is a guess for one point the periodic orbit and T is the period of the guess. The type of x depends on what is passed to the newton solver. See SectionSS for a type of section defined as a hyperplane.
• parallel whether the shooting are computed in parallel (threading). Available through the use of Flows defined by EnsembleProblem.

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

• pb(orbitguess, par) evaluates the functional G on orbitguess
• pb(orbitguess, par, du; δ = 1e-9) evaluates the jacobian dG(orbitguess).du functional at orbitguess on du. The optional argument δ is used to compute a finite difference approximation of the derivative of the section.
• pb(Val(:JacobianMatrixInplace), J, x, par)` compute the jacobian of the functional analytically. This is based on ForwardDiff.jl. Useful mainly for ODEs.
• pb(Val(:JacobianMatrix), x, par) same as above but out-of-place.

You can then call pb(orbitguess, par) to apply the functional to a guess. Note that orbitguess::AbstractVector must be of size M * N + 1 where N is the number of unknowns of the state space and orbitguess[M * N + 1] is an estimate of the period T of the limit cycle. This form of guess is convenient for the use of the linear solvers in IterativeSolvers.jl (for example) which accepts only AbstractVectors. Another accepted guess is of the form BorderedArray(guess, T) where guess[i] is the state of the orbit at the ith time slice. This last form allows for non-vector state space which can be convenient for 2d problems for example, use GMRESKrylovKit for the linear solver in this case.

Simplified constructors

• A simpler way to build the functional is to use

pb = ShootingProblem(F, p, prob::Union{ODEProblem, EnsembleProblem}, alg, centers::AbstractVector; kwargs...)

where F(x,p) is the vector field, p is a parameter (to be passed to the vector field and the flow), prob is an ODEProblem (resp. EnsembleProblem) which is used to create a flow using the ODE solver alg (for example Tsit5()). centers is list of M points close to the periodic orbit, they will be used to build a constraint for the phase. parallel = false is an option to use Parallel simulations (Threading) to simulate the multiple trajectories in the case of multiple shooting. This is efficient when the trajectories are relatively long to compute. Finally, the arguments kwargs are passed to the ODE solver defining the flow. Look at DifferentialEquations.jl for more information. Note that, in this case, the derivative of the flow is computed internally using Finite Differences.

• Another way to create a Shooting problem with more options is the following where in particular, one can provide its own scalar constraint section(x)::Number for the phase

pb = ShootingProblem(F, p, prob::Union{ODEProblem, EnsembleProblem}, alg, M::Int, section; parallel = false, kwargs...)

or

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

• The next way is an elaboration of the previous one

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

or

pb = ShootingProblem(F, p, 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)$.

pb = PoincareShootingProblem(flow::Flow, M, sections; δ = 1e-8, interp_points = 50, parallel = false)

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

• flow::Flow: implements the flow of the Cauchy problem though the structure Flow.
• M: the number of Poincaré sections. If M==1, then the simple shooting is implemented and the multiple one otherwise.
• sections: function or callable struct which implements a Poincaré section condition. The evaluation sections(x) must return a scalar number when M==1. Otherwise, one must implement a function section(out, x) which populates out with the M sections. See SectionPS for type of section defined as a hyperplane.
• δ = 1e-8 used to compute the jacobian of the fonctional by finite differences. If set to 0, an analytical expression of the jacobian is used instead.
• interp_points = 50 number of interpolation point used to define the callback (to compute the hitting of the hyperplan section)
• parallel = false whether the shooting are computed in parallel (threading). Only available through the use of Flows defined by EnsembleProblem.

Simplified constructors

• A simpler way is to create a functional is

pb = PoincareShootingProblem(F, p, prob::ODEProblem, alg, section; kwargs...)

for simple shooting or

pb = PoincareShootingProblem(F, p, prob::Union{ODEProblem, EnsembleProblem}, alg, M::Int, section; kwargs...)

for multiple shooting . Here F(x,p) is the vector field, p is a parameter (to be passed to the vector and the flow), 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 refere to DifferentialEquations.jl for more information.

• Another convenient call is

pb = PoincareShootingProblem(F, p, prob::Union{ODEProblem, EnsembleProblem}, alg, normals::AbstractVector, centers::AbstractVector; δ = 1e-8, kwargs...)

where normals (resp. centers) is a list of normals (resp. centers) which defines a list of hyperplanes $\Sigma_i$. These hyperplanes are used to define partial Poincaré return maps.

Computing the functionals

A functional, hereby called G encodes this shooting problem. You can then call pb(orbitguess, par) to apply the functional to a guess. Note that orbitguess::AbstractVector must be of size M * N where N is the number of unknowns in the state space and M is the number of Poincaré maps. Another accepted guess is such that guess[i] is the state of the orbit on the ith section. This last form allows for non-vector state space which can be convenient for 2d problems for example.

• pb(orbitguess, par) evaluates the functional G on orbitguess
• pb(orbitguess, par, du) evaluates the jacobian dG(orbitguess).du functional at orbitguess on du
• pb(Val(:JacobianMatrixInplace), J, x, par)` compute the jacobian of the functional analytically. This is based on ForwardDiff.jl. Useful mainly for ODEs.
• pb(Val(:JacobianMatrix), x, par) same as above but out-of-place.

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, Tse, Tprob, TprobMono, Tfs}

• F::Any

  The vector field (x, p) -> F(x, p) associated to a Cauchy problem, Default: nothing

• flow::Any

  The flow (or semigroup) associated to the Cauchy problem (x, p, t) -> flow(x, p, t). Only the last time point must be returned. Default: nothing

• flowTimeSol::Any

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

• flowFull::Any

  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, mainly used for plotting on the user side. Please use nothing as default. Default: nothing

• dflow::Any

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

• dfSerial::Any

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

• prob::Any

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

• probMono::Any

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

• flowSerial::Any

  [Internal] Serial version of the flow Default: nothing

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

There 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(F, prob, Tsit5(); kwargs...)

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

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

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

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

floquet = FloquetQaD(eigsolver::AbstractEigenSolver)

This composite type implements the computation of the eigenvalues of the monodromy matrix in the case of periodic orbits problems (based on the Shooting method or Finite Differences (Trapeze method)), also called the Floquet multipliers. The method, dubbed Quick and Dirty (QaD), is not numerically very precise for large / small Floquet exponents. It allows, nevertheless, to detect bifurcations. The arguments are as follows:

• eigsolver::AbstractEigenSolver solver used to compute the eigenvalues.

If eigsolver == DefaultEig(), then the monodromy matrix is formed and its eigenvalues are computed. Otherwise, a Matrix-Free version of the monodromy is used.

guessFromHopf(br, ind_hopf, eigsolver::AbstractEigenSolver, M, amplitude; phase = 0)

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

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

The arguments are

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

computeNormalForm(F, dF, d2F, d3F, br, id_bif; δ, nev, Jᵗ, verbose, ζs, lens, issymmetric, Teigvec, scaleζ)

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

Arguments

• F, dF, d2F, d3F vector field (x, p) -> F(x, p) and its derivatives w.r.t. x.
• br result from a call to continuation
• ind_bif index of the bifurcation point in br.bifpoint

Optional arguments

• δ used to compute ∂pF with finite differences
• nev number of eigenvalues used to compute the spectral projection. This number has to be adjusted when used with iterative methods.
• Jᵗ = (x,p) -> ... 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 avoid recomputing the jacobian as it would be if you pass Jᵗ = (x, p) -> transpose(dF(x, p)).
• verbose whether to display information
• ζs list of vectors spanning the kernel of dF at the bifurcation point. Useful to enforce the basis for the normal form.
• lens::Lens specify which parameter to take the partial derivative ∂pF
• issymmetric whether the Jacobian is Symmetric, avoid computing the left eigenvectors.
• scaleζ function to normalise the kernel basis. Indeed, when used with large vectors and norm, it results in ζs and the normal form coeffocient being super small.

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

	newton(F, J, x0, p0, options::NewtonPar; normN = norm, callback = (x, f, J, res, iteration, itlinear, optionsN; kwargs...) -> true, kwargs...)

This is the Newton-Krylov Solver for F(x, p0) = 0 with Jacobian w.r.t. x written J(x, p0) and initial guess x0. The function normN allows to specify a norm for the convergence criteria. It is important to set the linear solver options.linsolver properly depending on your problem. This linear solver is used to solve $J(x, p_0)u = -F(x, p_0)$ in the Newton step. You can for example use linsolver = DefaultLS() which is the operator backslash: it works well for Sparse / Dense matrices. See Linear solvers (LS) for more informations.

Arguments:

• F is a function with input arguments (x, p) returning a vector r that represents the functional and for type stability, the types of x and r should match. In particular, it is not inplace.
• J is the jacobian of F at (x, p). It can assume two forms. Either J is a function and J(x, p) returns a ::AbstractMatrix. In this case, the default arguments of NewtonPar will make newton work. Or J is a function and J(x, p) returns a function taking one argument dx and returns dr of the same type of dx. In our notation, dr = J * dx. In this case, the default parameters of NewtonPar will not work and you have to use a Matrix Free linear solver, for example GMRESIterativeSolvers.
• x0 initial guess
• p0 set of parameters to be passed to F and J
• options variable holding the internal parameters used by the newton method
• callback function passed by the user which is called at the end of each iteration. Can be used to update a preconditionner for example. The arguments passed to the callback are as follows
  • x current solution
  • f current residual
  • J current jacobian
  • res current norm of the residual
  • iteration current newton iteration
  • itlinear number of iterations to solve the linear system
  • optionsN a copy of the argument options passed to newton
  • kwargs kwargs arguments, contain your initial guess x0
• kwargs arguments passed to the callback. Useful when newton is called from continuation

Output:

• solution:
• history of residuals
• flag of convergence
• number of iterations

Simplified calls

When J is not passed, the jacobian matrix is then computed with finite differences (beware of large systems of equations!). The call is as follows:

newton(F, x0, p0, options::NewtonPar; kwargs...)

You can also pass functions which do not have parameters x -> F(x), x -> J(x) as follows

newton(F, J, x0, options::NewtonPar;  kwargs...)

or

newton(F, x0, options::NewtonPar;  kwargs...)

Example

julia> F(x, p) = x.^3 .- 1
julia> Jac(x, p) = spdiagm(0 => 3 .* x.^2) # sparse jacobian
julia> x0 = rand(1_000)
julia> opts = NewtonPar()
julia> sol, hist, flag, _ = newton(F, Jac, x0, nothing, opts, normN = x->norm(x, Inf))

function newton(F, J, x0::vectype, p0, options:: NewtonPar{T}, defOp::DeflationOperator{T, Tf, vectype}, linsolver = DeflatedLinearSolver(); kwargs...) where {T, Tf, vectype}

This is the deflated version of the Krylov-Newton Solver for F(x, p0) = 0 with Jacobian J(x, p0). We refer to newton for more information. It penalises the roots saved in defOp.roots. The other arguments are as for newton. See DeflationOperator for more information.

Arguments

Compared to newton, the only different arguments are

• defOp deflation operator
• linsolver linear solver used to invert the Jacobian of the deflated functional. We have a custom solver DeflatedLinearSolver() with requires solving two linear systems J⋅x = rhs. For other linear solvers, a matrix free method is used for the deflated functional.

Output:

• solution:
• history of residuals
• flag of convergence
• number of iterations

Simplified call

When J is not passed. It then computed with finite differences. The call is as follows:

newton(F, x0, p0, options, defOp; kwargs...)

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

newton(F, J, br, ind_bif, par, lens; Jᵗ, d2F, normN, options, kwargs...)

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

• F = (x, p) -> F(x, p) where p is a set of parameters.
• J = (x, p) -> d_xF(x, p) associated jacobian
• br results returned after a call to continuation
• ind_bif bifurcation index in br
• par parameters used for the vector field
• lens parameter axis used to locate the Fold/Hopf point.
• options::NewtonPar

Optional arguments:

• Jᵗ = (x, p) -> transpose(d_xF(x, p)) 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 Matrix / Sparse based jacobian, Jᵗ should not be passed as it is computed internally more efficiently, i.e. it avoid recomputing the jacobian as it would be if you pass Jᵗ = (x, p) -> transpose(dF(x, p))
• d2F = (x, p, v1, v2) -> d2F(x, p, v1, v2) a bilinear operator representing the hessian of F. It has to provide an expression for d2F(x,p)[v1,v2].
• normN = norm
• options You can pass newton parameters different from the ones stored in br by using this argument options.
• kwargs keywords arguments to be passed to the regular Newton-Krylov solver

newton(prob, orbitguess, par, options; linearPO, δ, kwargs...)

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

Arguments

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

• prob a problem of type <: AbstractShootingProblem encoding the shooting functional G.
• orbitguess a guess for the periodic orbit. See ShootingProblem and See PoincareShootingProblem for information regarding the shape of orbitguess.
• par parameters to be passed to the functional
• options same as for the regular newton method.

Optional argument

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

Output:

• solution
• history of residuals
• flag of convergence
• number of iterations

newton(prob, orbitguess, par, options, defOp; linearPO, kwargs...)

This is the deflated Newton-Krylov Solver for computing a periodic orbit using a (Standard / Poincaré) Shooting method.

Arguments

Similar to newton except that prob is either a ShootingProblem or a PoincareShootingProblem.

Optional argument

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

Output:

• solution
• history of residuals
• flag of convergence
• number of iterations

newton(probPO, orbitguess, par, options; linearPO, kwargs...)

This is the Newton-Krylov Solver for computing a periodic orbit using a functional G based on Finite Differences and a Trapezoidal rule.

Arguments:

• prob a problem of type PeriodicOrbitTrapProblem encoding the functional G
• orbitguess a guess for the periodic orbit where orbitguess[end] is an estimate of the period of the orbit. It should be a vector of size N * M + 1 where M is the number of time slices, N is the dimension of the phase space. This must be compatible with the numbers N,M in prob.
• par parameters to be passed to the functional
• options same as for the regular newton method
• linearPO = :BorderedLU. Specify the choice of the linear algorithm, which must belong to [:FullLU, :FullSparseInplace, :BorderedLU, :FullMatrixFree, :BorderedMatrixFree, :FullSparseInplace]. This is used to select a way of inverting the jacobian dG of the functional G.
  • For :FullLU, we use the default linear solver based on a sparse matrix representation of dG. This matrix is assembled at each newton iteration. This is the default algorithm.
  • For :FullSparseInplace, this is the same as for :FullLU but the sparse matrix dG is updated inplace. This method allocates much less. In some cases, this is significantly faster than using :FullLU. Note that this method can only be used if the sparsity pattern of the jacobian is always the same.
  • For :Dense, same as above but the matrix dG is dense. It is also updated inplace. This option is useful to study ODE of small dimension.
  • For :BorderedLU, we take advantage of the bordered shape of the linear solver and use a LU decomposition to invert dG using a bordered linear solver. This is the default algorithm.
  • For :BorderedSparseInplace, this is the same as for :BorderedLU but the cyclic matrix dG is updated inplace. This method allocates much less. In some cases, this is significantly faster than using :FullLU. Note that this method can only be used if the sparsity pattern of the jacobian is always the same.
  • For :FullMatrixFree, a matrix free linear solver is used for dG: note that a preconditioner is very likely required here because of the cyclic shape of dG which affects negatively the convergence properties of GMRES.
  • For :BorderedMatrixFree, a matrix free linear solver is used but for Jc only (see docs): it means that options.linsolver is used to invert Jc. These two Matrix-Free options thus expose different part of the jacobian dG in order to use specific preconditioners. For example, an ILU preconditioner on Jc could remove the constraints in dG and lead to poor convergence. Of course, for these last two methods, a preconditioner is likely to be required.

Output:

• solution
• history of residuals
• flag of convergence
• number of iterations

newton(probPO::PeriodicOrbitTrapProblem, orbitguess, options::NewtonPar, defOp::DeflationOperator{T, Tf, vectype}, linearPO = :BorderedLU; kwargs...) where {T, Tf, vectype}

This function is similar to newton(probPO, orbitguess, options, linearPO; kwargs...) except that it uses deflation in order to find periodic orbits different from the ones stored in defOp. We refer to the mentioned method for a full description of the arguments. The current method can be used in the vicinity of a Hopf bifurcation to prevent the Newton-Krylov algorithm from converging to the equilibrium point.

continuation(F, J, x0, par, lens::Lens, contParams::ContinuationPar; plot = false, normC = norm, dotPALC = (x,y) -> dot(x,y) / length(x), printSolution = norm, plotSolution = (x, p; kwargs...)->nothing, finaliseSolution = (z, tau, step, contResult; kwargs...) -> true, callbackN = (x, f, J, res, iteration, itlinear, options; kwargs...) -> true, linearAlgo = BorderingBLS(), tangentAlgo = SecantPred(), verbosity = 0)

Compute the continuation curve associated to the functional F and its jacobian J.

Arguments:

• F is a function with input arguments (x, p), where p is the set of parameters passed to F, and returning a vector r that represents the functional. For type stability, the types of x and r should match. In particular, it is not inplace,
• J is the 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.
• x0 initial guess,
• par initial set of parameters,
• lens::Lens specifies which parameter axis among par is used for continuation. For example, if par = (α = 1.0, β = 1), we can perform continuation w.r.t. α by using lens = (@lens _.α). If you have an array par = [ 1.0, 2.0] and want to perform continuation w.r.t. the first variable, you can use lens = (@lens _[1]). For more information, we refer to SetField.jl.
• contParams parameters for continuation. See ContinuationPar for more information about the options

Optional Arguments:

• plot = false whether to plot the solution while computing
• printSolution = (x, p) -> norm(x) function used to plot in the continuation curve. It is also used in the way results are saved. It could be norm or (x, p) -> x[1]. This is also useful when saving several huge vectors is not possible for memory reasons (for example on GPU...). This function can return pretty much everything but you should keep it small. For example, you can do (x, p) -> (x1 = x[1], x2 = x[2], nrm = norm(x)) or simply (x, p) -> (sum(x), 1). This will be stored in contres.branch (see below).
• plotSolution = (x, p; kwargs...) -> nothing function implementing the plot of the solution.
• finaliseSolution = (z, tau, step, contResult; kwargs...) -> true Function called at the end of each continuation step. Can be used to alter the continuation procedure (stop it by returning false), saving personal data, plotting... The notations are $z=(x,p)$, tau is the tangent at z (see below), step is the index of the current continuation step and ContResult is the current branch. Note that you can have a better control over the continuation procedure by using an iterator, see Iterator Interface.
• callbackN callback for newton iterations. see docs for newton. Can be used to change preconditioners
• tangentAlgo = SecantPred() controls the algorithm used to predict the tangents along the curve of solutions or the corrector. Can be NaturalPred, SecantPred or BorderedPred. See below for more information.
• linearAlgo = BorderingBLS(). Used to control the way the extended linear system associated to the continuation problem is solved. Can be MatrixBLS, BorderingBLS or MatrixFreeBLS.
• verbosity::Int controls the amount of information printed during the continuation process. Must belong to {0,1,2,3}
• normC = norm norm used in the different Newton solves
• dotPALC = (x, y) -> dot(x, y) / length(x), dot product used to define the weighted dot product (resp. norm) $\|(x, p)\|^2_\theta$ in the constraint $N(x, p)$ (see below). This arguement can be used to remove the factor 1/length(x) for example in problems where the dimension of the state space changes (mesh adaptation, ...)
• filename name of a file to save the computed branch during continuation. The identifier .jld2 will be appended to this filename

Outputs:

• contres::ContResult composite type which contains the computed branch. See ContResult for more information.
• u::BorderedArray the last solution computed on the branch

Simplified call:

You can also use the following call for which the jacobian matrix (beware of large systems of equations!) is computed internally using Finite Differences

continuation(Fhandle, x0, par, lens, contParams::ContinuationPar; kwargs...)

Method

Bordered system of equations

In what follows, we abuse of notations, p refers to the scalar value of the parameter we perform continuation with. Hence, it should be p = get(par, lens).

The pseudo-arclength continuation method solves the equation $F(x, p) = 0$ (of dimension N) together with the pseudo-arclength constraint $N(x, p) = \frac{\theta}{length(x)} \langle x - x_0, dx_0\rangle + (1 - \theta)\cdot(p - p_0)\cdot dp_0 - ds = 0$ and $\theta\in[0,1]$ (see Keller, Herbert B. Lectures on Numerical Methods in Bifurcation Problems. Springer, 1988). In practice, a curve $\gamma$ of solutions is sought and is parametrised by $s$: $\gamma(s) = (x(s), p(s))$ is a curve of solutions to $F(x, p)$. This formulation allows to pass turning points (where the implicit theorem fails). In the previous formula, $(x_0, p_0)$ is a solution for a given $s_0$, $\tau_0\equiv(dx_0, dp_0)$ is the tangent to the curve $\gamma$ at $s_0$. Hence, to compute the curve of solutions, we need to solve an equation of dimension N+1 which is called a Bordered system.

The parameter ds is adjusted internally depending on the number of Newton iterations and other factors. See the function stepSizeControl for more information. An important parameter to adjust the magnitude of this adaptation is the parameter a in the struct ContinuationPar.

Algorithm

The algorithm works as follows:

1. Start from a known solution $(x_0, p_0)$ with tangent to the curve of solutions: $(dx_0 ,dp_0)$
2. Predictor: set $(x_1, p_1) = (x_0, p_0) + ds\cdot (dx_0, dp_0)$. Note that a different predictor can be used.
3. Corrector: solve $F(x, p)=0,\ N(x, p)=0$ with a (Bordered) Newton Solver with initial guess $(x_1, p_1)$.
   • if Newton in 3. did not converge, update ds/2 ⟶ ds in $N$ and go to 1.
4. New tangent: Compute a new tangent (see below) $(dx_1, dp_1)$ and update $N$ with it. Set $(x_0, p_0, dx_0, dp_0) = (x_1, p_1, dx_1, dp_1)$ and return to step 2

Natural continuation

We speak of natural continuation when we do not consider the constraint $N(x, p)=0$. Knowing $(x_0, p_0)$, we use $x_0$ as a guess for solving $F(x, p_1)=0$ with $p_1$ close to $p_0$. Again, this fails at Turning points but it can be faster to compute than the constrained case. This is set by the option tangentAlgo = NaturalPred() in continuation.

Tangent computation (step 4)

There are various ways to compute $(dx_1, dp_1)$. The first one is called secant and is parametrised by the option tangentAlgo = SecantPred(). It is computed by $(dx_1, dp_1) = (z_1, p_1) - (z_0, p_0)$ and normalised by the norm $\|(x, p)\|^2_\theta = \frac{\theta}{length(x)} \langle x,x\rangle + (1 - \theta)\cdot p^2$. Another method is to compute $(dx_1, dp_1)$ by solving solving the bordered linear system $\begin{bmatrix} F_x & F_p ; \ \frac{\theta}{length(x)}dx_0 & (1-\theta)dp_0\end{bmatrix}\begin{bmatrix}dx_1 ; dp_1\end{bmatrix} =\begin{bmatrix}0 ; 1\end{bmatrix}$ ; it is set by the option tangentAlgo = BorderedPred().

Bordered linear solver

When solving the Bordered system $F(x, p) = 0,\ N(x, p)=0$, one faces the issue of solving the Bordered linear system $\begin{bmatrix} J & a ; b^T & c\end{bmatrix}\begin{bmatrix}X ; y\end{bmatrix} =\begin{bmatrix}R ; n\end{bmatrix}$. This can be solved in many ways via bordering (which requires two Jacobian inverses), by forming the bordered matrix (which works well for sparse matrices) or by using a full Matrix Free formulation. The choice of method is set by the argument linearAlgo. Have a look at the struct linearBorderedSolver for more information.

Linear Algebra

Let us discuss here more about the norm and dot product. First, the option normC gives a norm that is used to evaluate the residual in the following way: $max(normC(F(x,p)), \|N(x,p)\|)<tol$. It is thus used as a stopping criterion for a Newton algorithm. The dot product (resp. norm) used in $N$ and in the (iterative) linear solvers is LinearAlgebra.dot (resp. LinearAlgebra.norm). It can be changed by importing these functions and redefining it. Not that by default, the $L^2$ norm is used. These details are important because of the constraint $N$ which incorporates the factor length. For some custom composite type implementing a Vector space, the dot product could already incorporates the length factor in which case you should either redefine the dot product or change $\theta$.

Step size control

As explained above, each time the corrector phased failed, the step size $ds$ is halved. This has the disavantage of having lost Newton iterations (which costs time) and impose small steps (which can be slow as well). To prevent this, the step size is controlled internally with the idea of having a constant number of Newton iterations per point. This is in part controlled by the aggressiveness factor a in ContinuationPar. Further tuning is performed by using doArcLengthScaling=true in ContinuationPar. This adjusts internally $\theta$ so that the relative contributions of $x$ and $p$ are balanced in the constraint $N$.

continuation(Fhandle, Jhandle, x0, par0, x1, p1, lens, contParams; linearAlgo, kwargs...)

This function is the analog of continuation when the two first points on the branch are passed (instead of a single one). Hence x0 is the first point on the branch (with palc s=0) with parameter par0 and x1 is the second point with parameter set(par0, lens, p1).

continuation(F, dF, d2F, d3F, br, ind_bif, optionsCont; Jᵗ, δ, δp, ampfactor, nev, issymmetric, usedeflation, Teigvec, scaleζ, verbosedeflation, maxIterDeflation, perturb, kwargs...)

Automatic branch switching at branch points based on a computation of the normal form. More information is provided in Branch switching. An example of use is provided in A generalized Bratu–Gelfand problem in two dimensions.

Arguments

• F, dF, d2F, d3F: function (x, p) -> F(x, p) and its differentials (x, p, dx) -> d1F(x, p, dx), (x, p, dx1, dx2) -> d2F(x, p, dx1, dx2)...
• br branch result from a call to continuation
• ind_bif index of the bifurcation point in br from which you want to branch from
• optionsCont options for the call to continuation

Optional arguments

• Jᵗ associated jacobian transpose, 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 avoid recomputing the jacobian as it would be if you pass Jᵗ = (x, p) -> transpose(dF(x, p)).
• δ used internally to compute derivatives w.r.t the parameter p.
• δp used to specify a particular guess for the parameter on the bifurcated branch which is otherwise determined by optionsCont.ds. This allows to use a step larger than optionsCont.dsmax.
• ampfactor = 1 factor which alter the amplitude of the bifurcated solution. Useful to magnify the bifurcated solution when the bifurcated branch is very steep.
• nev number of eigenvalues to be computed to get the right eigenvector
• issymmetric whether the Jacobian is Symmetric, avoid computing the left eigenvectors in the computation of the reduced equation.
• usedeflation = true whether to use nonlinear deflation (see Deflated problems) to help finding the guess on the bifurcated branch
• kwargs optional arguments to be passed to continuation, the regular continuation one.

continuation(F, J, par, lens, contParams, defOp; verbosity, maxBranches, seekEveryStep, showplot, tangentAlgo, linearAlgo, dotPALC, printSolution, 367, plotSolution, perturbSolution, 370, callbackN, acceptSolution, normN)

The 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 tangentAlgo.

Arguments:

• F is a function with input arguments (x, p), where p is the set of parameters passed to F, and returning a vector r that represents the functional. For type stability, the types of x and r should match. In particular, it is not inplace,
• J is the 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.
• par initial set of parameters,
• lens::Lens specifies which parameter axis among par is used for continuation. For example, if par = (α = 1.0, β = 1), we can perform continuation w.r.t. α by using lens = (@lens _.α). If you have an array par = [ 1.0, 2.0] and want to perform continuation w.r.t. the first variable, you can use lens = (@lens _[1]). For more information, we refer to SetField.jl,
• contParams parameters for continuation. See ContinuationPar for more information about the options,
• defOp::DeflationOperator a Deflation Operator (see DeflationOperator) which contains the set of solution guesses for the parameter par.

Optional Arguments:

• seekEveryStep::Int = 1 we look for additional solution, using deflated newton, every seekEveryStep step,
• maxBranches::Int = 10 maximum number of branches considered,
• showplot = false whether to plot the solution while computing,
• printSolution = (x, p) -> norm(x) function used to plot in the continuation curve. It is also used in the way results are saved. It could be norm or (x, p) -> x[1]. This is also useful when saving several huge vectors is not possible for memory reasons (for example on GPU...),
• plotSolution = (x, p; kwargs...) -> nothing function implementing the plot of the solution,
• callbackN callback for newton iterations. see docs for newton. Can be used to change preconditioners or affect the newton iterations. In the deflation part of the algorithm, when seeking for new branches, the callback is passed the keyword argument fromDeflatedNewton = true to tell the user can it is not in the continuation part (regular newton) of the algorithm,
• tangentAlgo = NaturalPred() controls the algorithm used to predict the tangents along the curve of solutions or the corrector. Can be NaturalPred, SecantPred or BorderedPred,
• verbosity::Int controls the amount of information printed during the continuation process. Must belong to {0,⋯,5},
• normN = norm norm used in the different Newton solves,
• dotPALC = (x, y) -> dot(x, y) / length(x), dot product used to define the weighted dot product (resp. norm) $\|(x, p)\|^2_\theta$ in the constraint $N(x, p)$ (see below). This arguement can be used to remove the factor 1/length(x) for example in problems where the dimension of the state space changes (mesh adaptation, ...),
• perturbSolution = (x, p, id) -> x .+ (1 .+ 0.001 * rand(size(x)...)), perturbation applied to the solution when trying to fimnd new solutions using Deflated Newton.

Outputs:

• contres::Vector{ContResult} composite type which contains the computed branches. See ContResult for more information,
• the solutions at the last parameter value,
• current parameter value.

continuation(prob, orbitguess, par, lens, _contParams; linearAlgo, kwargs...)

This is the continuation routine for computing a periodic orbit using a (Standard / Poincaré) Shooting method.

Arguments

Similar to continuation except that prob is either a ShootingProblem or a PoincareShootingProblem.

• printPeriod boolean to print the period of the solution. This is useful for prob::PoincareShootingProblem as this information is not easily available.

continuation(F, dF, d2F, d3F, br, ind_bif, _contParams, prob; Jᵗ, δ, δp, ampfactor, usedeflation, nev, updateSectionEveryStep, kwargs...)

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

Arguments

• F, dF, d2F, d3F: function (x,p) -> F(x,p) and its differentials (x,p,dx) -> d1F(x,p,dx), (x,p,dx1,dx2) -> d2F(x,p,dx1,dx2)... These are used to compute the Hopf normal form.
• br branch result from a call to continuation
• ind_hopf index of the bifurcation point in br
• contParams parameters for the call to continuation
• prob problem used to specify the way the periodic orbit is computed. It can be PeriodicOrbitTrapProblem, ShootingProblem or PoincareShootingProblem .

Optional arguments

• linearPO linear algorithm used for the computation of periodic orbits when prob is PeriodicOrbitTrapProblem)
• Jᵗ is the jacobian adjoint, used for computation of the eigen-elements of the jacobian adjoint, needed to compute the spectral projector for the Hopf normal form. For matrix-free methods, transpose is not readily available and the user must provide a dedicated method. In the case of Matrix / Sparse based jacobian, Jᵗ should not be passed as it is computed internally more efficiently, i.e. it avoid recomputing the jacobian as it would be if you pass Jᵗ = (x, p) -> transpose(dF(x, p))
• δ = 1e-8 used for finite differences
• δp used to specify a particular guess for the parameter on the bifurcated branch which is otherwise determined by contParams.ds. This allows to use a step larger than contParams.dsmax.
• ampfactor = 1 factor which alter the amplitude of the bifurcated solution. Useful to magnify the bifurcated solution when the bifurcated branch is very steep.
• usedeflation = true whether to use nonlinear deflation (see Deflated problems) to help finding the guess on the bifurcated branch
• updateSectionEveryStep = 0 updates the section every updateSectionEveryStep step during continuation
• linearPO specify the way the jacobian is computed.
• all kwargs from continuation

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

continuation(prob, orbitguess, par, lens, _contParams; linearPO, printSolution, linearAlgo, updateSectionEveryStep, kwargs...)

This is the continuation routine for computing a periodic orbit using a functional G based on Finite Differences and a Trapezoidal rule.

Arguments

• prob::PeriodicOrbitTrapProblem encodes the functional G
• orbitguess a guess for the periodic orbit where orbitguess[end] is an estimate of the period of the orbit. It could be a vector of size N * M + 1 where M is the number of time slices, N is the dimension of the phase space. This must be compatible with the numbers N, M in prob.
• p0 set of parameters passed to the vector field
• contParams same as for the regular continuation method
• linearAlgo same as in continuation
• linearPO = :BorderedLU. Same as newton when applied to a PeriodicOrbitTrapProblem. More precisely:
  • For :FullLU, we use the default linear solver on a sparse matrix representation of dG. This matrix is assembled at each newton iteration.
  • For :FullSparseInplace, this is the same as :FullLU but the sparse matrix dG is updated inplace. This method thus allocates much less. In some cases, this is significantly faster than using :FullLU. Note that this method can only be used if the sparsity pattern of the jacobian is always the same.
  • For :Dense, same as above but the matrix dG is dense. It is also updated inplace. This option is useful to study ODE of small dimension.
  • For :BorderedLU, we take advantage of the bordered shape of the linear solver and use LU decomposition to invert dG using a bordered linear solver. This is the default algorithm.
  • For :FullMatrixFree, a matrix free linear solver is used for dG: note that a preconditioner is very likely required here because of the cyclic shape of dG which affects negatively the convergence properties of GMRES.
  • For :BorderedMatrixFree, a matrix free linear solver is used but for Jc only (see docs): it means that options.linsolver is used to invert Jc. These two Matrix-Free options thus expose different part of the jacobian dG in order to use specific preconditioners. For example, an ILU preconditioner on Jc could remove the constraints in dG and lead to poor convergence. Of course, for these last two methods, a preconditioner is likely to be required.
• updateSectionEveryStep = 1 updates the section every when mod(step, updateSectionEveryStep) == 1 during continuation

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 printSolution argument.

bifurcationdiagram(F, dF, d2F, d3F, x0, par0, lens, level, options; usedeflation, kwargs...)

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

Arguments

• F, dF, d2F, d3F functional and its derivatives
• x0 initial guess
• par0 parameter values at x0
• lens lens to select the parameter axis
• 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 as for continuation but also for the different versions listed in Continuation.

Simplified call:

We also provide the call

bifurcationdiagram(F, dF, d2F, d3F, br::ContResult, level::Int, options; usedeflation = false, kwargs...)

where br is a branch computed after a call to continuation from which we want to compute the bifurcating branches recursively.

bifurcationdiagram!(F, dF, d2F, d3F, node, level, options; code, usedeflation, kwargs...)

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

getBranch(tree, code)

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

getBranchesFromBP(tree, indbif)

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

struct GenericBifPoint{T, Tp, Tv} <: BifurcationKit.BifurcationPoint

Structure to record a generic bifurcation point which was only detected by a change in the number of stable eigenvalues.

• type::Symbol

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

• idx::Int64

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

• param::Any

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

• norm::Any

  Norm of the equilibrium at the bifurcation point Default: 0.0

• printsol::Any

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

• x::Any

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

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

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

• ind_ev::Int64

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

• step::Int64

  Continuation step at which the bifurcation occurs Default: 0

• status::Symbol

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

• precision::Any

  Precision in the location of the bifurcation point Default: -1

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

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

getPeriod(sh, x)
getPeriod(sh, x, par)

Compute the period of the periodic orbit associated to x.

getPeriod(psh, x_bar, par)

Compute the period of the periodic orbit associated to x_bar.

getPeriod(prob, x, p)

Compute the period of the periodic orbit associated to x.

getAmplitude(prob, x, p; ratio)

Compute the amplitude of the periodic orbit associated to x. The keyword argument ratio = 1 is used as follows. If length(x) = ratio * n, the call returns the amplitude over x[1:n].

getAmplitude(prob, x, p; ratio)

Compute the amplitude of the periodic orbit associated to x. The keyword argument ratio = 1 is used as follows. If length(x) = 1 + ratio * n, the call returns the amplitude over x[1:n].

getMaximum(prob, x, p; ratio)

Compute the maximum of the periodic orbit associated to x. The keyword argument ratio = 1 is used as follows. If length(x) = ratio * n, the call returns the amplitude over x[1:n].

getMaximum(prob, x, p; ratio)

Compute the maximum of the periodic orbit associated to x. The keyword argument ratio = 1 is used as follows. If length(x) = 1 + ratio * n, the call returns the amplitude over x[1:n].

struct SectionSS{Tn, Tc} <: BifurcationKit.AbstractSection

This composite type (named for Section Standard Shooting) encodes a type of section implemented by a hyperplane. It can be used in conjunction with ShootingProblem. The hyperplane is defined by a point center and a normal.

• normal::Any

  Normal to define hyperplane

• center::Any

  Representative point on hyperplane

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

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

• M::Int64

• normals::Any

• centers::Any

• indices::Array{Int64,1}

• normals_bar::Any

• centers_bar::Any

Constructor(s)

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