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.
• saveSolEveryNsteps::Int64 = 0 at which continuation steps do we save the current solution`
• plotEveryNsteps = 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.
• saveEigEveryNsteps = 1 record eigen vectors every specified steps. Important for memory limited ressource, e.g. GPU.
• saveEigenvectors = true Important 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{T, Teigvals, Teigvec, Biftype, Ts, Tfunc, Tpar, Tl<:Setfield.Lens} <: BifurcationKit.BranchResult

Structure which holds the results after a call to continuation.

• branch::RecursiveArrayTools.VectorOfArray{T,2,Array{Array{T,1},1}} where T

  holds the low-dimensional information about the branch. More precisely, branch[:,i] contains the following information (param, printSolution(u, param), Newton iterations, ds, i) for each continuation step 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{Biftype,1} where Biftype

  A vector holding the set of fold points detected during the computation of the branch.

• stability::Array{Bool,1}

  A Vector{Bool} holding the stability of the computed solution for each continuation step. Hence, the stability stability[k] should match eig[k] which corresponds to branch[k] for a given k

• n_imag::Array{Int64,1}

  A Vector{Int64} holding the number of eigenvalues with positive real part and non zero imaginary part for each continuation step (to detect Hopf bifurcation)

• n_unstable::Array{Int64,1}

  A Vector{Int64} holding the number of eigenvalues with positive real part for each continuation step (to detect stationary bifurcation)

• 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 bifurcation points detected during the computation of the branch. Each entry of the vector contains a tuple with fields:

• • type bifurcation type, :hopf, :bp...,
  • idx is the index in eig (see above) for which the bifurcation occurs.
  • param parameter value at the bifurcation point
  • norm norm of the equilibrium at the bifurcation point
  • printsol = printSolution(x, param)
  • x equilibrium at the bifurcation point
  • tau tangent along the branch at the bifurcation point
  • ind_ev is the eigenvalue index responsible for the bifurcation (if applicable)
  • step is the continuation step at which the bifurcation occurs,
  • status ∈ {:converged, :guess} indicates if the bisection algorithm was successful in detecting the bifurcation point
  • δ = (δ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.
  • precision precision in location of the bifurcation point

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. The arguments are as follows

• F(x,p) vector field
• J jacobian of F. Can be matrix based J(x,p) or Matrix-Free
• Jt = 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)

You can then call pb(orbitguess, p) to compute the functional on a 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 of the limit cycle.

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

$\langle x[1] - x_\pi, \phi\rangle=0.$

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; isparallel = false)

This composite type implements the Standard Simple / Parallel Multiple Standard Shooting method to locate periodic orbits. 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) must return a scalar number where x is a guess for the periodic orbit. Note that the period T of the guess x is always included either as the last component of T = x[end] or as T = x.p. The type of x depends on what is passed to the newton solver.
• isparallel whether the shooting are computed in parallel (threading). Only available through the use of Flows defined by EnsembleProblem.

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.

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

Simplified constructors

• A simpler way to build a 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. isparallel = 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 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; isparallel = false, kwargs...)

or

pb = ShootingProblem(F, p, prob::Union{ODEProblem, EnsembleProblem}, alg, ds, section; isparallel = 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; isparallel = false, kwargs...)

or

pb = ShootingProblem(F, p, prob1::Union{ODEProblem, EnsembleProblem}, alg1, prob2::Union{ODEProblem, EnsembleProblem}, alg2, ds, section; isparallel = 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, isparallel = false)

This composite type implements the Poincaré Shooting method to locate periodic orbits by relying on Poincaré return maps. 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 strict 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.
• δ = 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)
• isparallel = 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

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}

• F::Any

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

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

• flowTimeSol::Any

  Flow which returns the tuple (t, u(t))

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

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

• dfserial::Any

  Serial version of dflow

Simplified constructors

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 = FloquetQaDTrap(eigsolver::AbstractEigenSolver)

This composite type implements the computation of the eigenvalues of the monodromy matrix in the case of periodic orbits problems based on 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 isa DefaultEig, then the monodromy matrix is formed and its eigenvalues are computed. Otherwise, a Matrix-Free version of the monodromy is used.

floquet = FloquetQaDShooting(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, 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, Jt, verbose, ζs, lens, issymmetric, Teigvec)

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.
• Jt = (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, Jt 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 Jt = (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.

	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, p0)u = -F(x, p0)$ 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 for more informations.

Arguments:

• (x, p) -> F(x, p) functional whose zeros are looked for. In particular, it is not inplace,
• dF(x, p) = (x, p) -> J(x, p) compute the jacobian of F at x. It is then passed to options.linsolver. The Jacobian J(x, p) can be a matrix or an out-of-place function.
• 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(Fhandle, 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 mnethod is based on a deflated Newton-Krylov solver.

newton(F, J, br, ind_bif, par, lens; Jt, 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:

• Jt = (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 sparse based jacobian, Jt 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 Jt = (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; 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.

Arguments

Similar as 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.

Output:

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

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

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

Arguments

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

Output:

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

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) -> 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 = (x, p) -> F(x, p) where p is the set of parameters passed to F.
• J = (x, p) -> d_xF(x, p) its associated jacobian. It can be a matrix, a function or a callable struct.
• 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
• 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 -> 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.
• finaliseSolution = (z, tau, step, contResult) -> 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]$. 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)$
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$.

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; Jt, δ, δp, ampfactor, nev, issymmetric, usedeflation, Teigvec, 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

• Jt 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, Jt 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 Jt = (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.
• 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(prob, orbitguess, par, lens, contParams; linearAlgo, printPeriod, kwargs...)

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

Arguments

Similar as 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; Jt, δ, δp, ampfactor, usedeflation, nev, 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)
• Jt 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!!!!!! COMME NEWTON FOLD
• δ = 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

continuation(probPO, orbitguess, p0::Real, _contParams::ContinuationPar; linearPO = :BorderedLU, printSolution = (u,p) -> u[end], linearAlgo = BorderingBLS(), 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 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 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 :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.

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=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. x,p is 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.

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