Fold / Hopf Continuation

In this page, we explain how to perform continuation of Fold / Hopf points and detect the associated bifurcations.

For this to work best, it is necessary to have an analytical expression for the jacobian. See the tutorial Temperature model (Simplest example) for more details.

A quite complete example for detection of codim 2 bifurcations of equilibria is CO-oxydation (codim 2) although it is for ODEs.

List of detected codim 2 bifurcation points

Bifurcationsymbol used
Bogdanov-Takensbt
Bautingh
Cuspcusp
Zero-Hopfzh
Hopf-Hopfhh

In a nutshell, all you have to do (see below) is to call continuation(br, ind_bif) to continue the bifurcation point stored in br.specialpoint[ind_bif] and set proper options.

Fold continuation

The continuation of Fold bifurcation points is based on a Minimally Augmented[Govaerts] formulation which is an efficient way to detect singularities. The continuation of Fold points is based on the formulation $G(u,p) = (F(u,p), g(u,p))\in\mathbb R^{n+1}$ where the test function $g$ is solution of

\[\left[\begin{array}{cc} dF(u,p) & w \\ v^{\top} & 0 \end{array}\right]\left[\begin{array}{c} r \\ g(u,p) \end{array}\right]=\left[\begin{array}{c}0_{n} \\1\end{array}\right]\quad\quad (M_f)\]

where $w,v$ are chosen in order to have a non-singular matrix $(M_f)$. More precisely, $v$ (resp. $w$) should be close to a null vector of dF(u,p) (resp. dF(u,p)'). During continuation, the vectors $w,v$ are updated so that the matrix $(M_f)$ remains non-singular ; this is controlled with the argument updateMinAugEveryStep (see below).

note that there are very simplified calls for this, see Newton refinement below. In particular, you don't need to set up the Fold Minimally Augmented problem yourself. This is done in the background.

Linear Method

You can pass the bordered linear solver to solve $(M_f)$ using the option bdlinsolver (see below). Note that the choice bdlinsolver = BorderingBLS() can lead to singular systems. Indeed, in this case, $(M_f)$ is solved by inverting dF(u,p) which is singular at Fold points.

Detection of codim 2 bifurcation points

You can detect the following codim 2 bifurcation points by using the option detectCodim2Bifurcation in the method continuation (see Codim 2 continuation). Under the hood, the detection of these bifurcations is done by using Event detection as explained in Event Handling.

  • the detection of Cusp (Cusp) is done by the detection of Fold bifurcation points along the curve of Folds by monitoring the parameter component of the tangent.
  • the detection of Bogdanov-Takens (BT) is performed using the test function[Bindel] $\psi_{BT}(p) = \langle w(p),v(p)\rangle$
  • the detection of Zero-Hopf (ZH) is performed by monitoring the number of eigenvalues $\lambda$ such that $\Re\lambda > \min\limits_{\nu\in\Sigma(dF)}|\Re\nu|$ and $\Im\lambda > \epsilon$ where $\epsilon$ is the Newton tolerance.

Hopf continuation

The continuation of Fold bifurcation points is based on a Minimally Augmented (see [Govaerts] p. 87) formulation which is an efficient way to detect singularities. The continuation of Hopf points is based on the formulation $G(u,\omega,p) = (F(u,\omega,p), g(u,\omega,p))\in\mathbb R^{n+2}$ where the test function $g$ is solution of

\[\left[\begin{array}{cc} dF(u,p)-i\omega I_n & w \\ v^{\top} & 0 \end{array}\right]\left[\begin{array}{c} r \\ g(u,\omega,p) \end{array}\right]=\left[\begin{array}{c} 0_{n} \\ 1 \end{array}\right]\quad\quad (M_h)\]

where $w,v$ are chosen in order to have a non-singular matrix $(M_h)$. More precisely, $w$ (resp. $v$) should be a left (resp. right) approximate null vector of $dF(u,p)-i\omega I_n$. During continuation, the vectors $w,v$ are updated so that the matrix $(M_h)$ remains non-singular ; this is controlled with the argument updateMinAugEveryStep (see below).

note that there are very simplified calls to this, see Newton refinement below. In particular, you don't need to set up the Hopf Minimally Augmented problem yourself. This is done in the background.

Linear Method

You can pass the bordered linear solver to solve $(M_h)$ using the option bdlinsolver (see below). Note that the choice bdlinsolver = BorderingBLS() can lead to singular systems. Indeed, in this case, $(M_h)$ is solved by inverting dF(u,p)-iω I_n which is singular at Hopf points.

Detection of codim 2 bifurcation points

You can detect the following codim 2 bifurcation points by using the option detectCodim2Bifurcation in the method continuation (see Codim 2 continuation). Under the hood, the detection of these bifurcations is done by using Event detection as explained in Event Handling.

  • the detection of Bogdanov-Takens (BT) is performed using the test function[Bindel] $\psi_{BT}(p) = \langle w(p),v(p)\rangle$
  • the detection of Bautin (GH) is based on the test function $\psi_{GH}(p) = \Re(l_1(p))$ where $l_1$ is the Lyapunov coefficient defined in Simple Hopf point.
  • the detection of Zero-Hopf (ZH) is performed by monitoring the eigenvalues.
  • the detection of Hopf-Hopf (HH) is performed by monitoring the eigenvalues.

The continuation of Hopf points is stopped at BT and when $\omega<100\epsilon$ where $\epsilon$ is the newton tolerance.

Newton refinement

Once a Fold / Hopf point has been detected after a call to br = continuation(...), it can be refined using newton iterations. Let us say that ind_bif is the index in br.specialpoint of a Fold / Hopf point. This guess can be refined as follows:

outfold = newton(br::AbstractBranchResult, ind_bif::Int;  
	normN = norm, options = br.contparams.newtonOptions, 
	bdlinsolver = BorderingBLS(options.linsolver),
	startWithEigen = false, kwargs...)

For the options parameters, we refer to Newton.

It is important to note that for improved performances, a function implementing the expression of the hessian should be provided. This is by far the fastest. Reader interested in this advanced usage should look at the code example/chan.jl of the tutorial Temperature model (Simplest example).

Codim 2 continuation

To compute the codim 2 curve of Fold / Hopf points, one can call continuation with the following options

BifurcationKit.continuationFunction
continuation(br, ind_bif, lens2)
continuation(br, ind_bif, lens2, options_cont; startWithEigen, detectCodim2Bifurcation, kwargs...)

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

  • br results returned after a call to continuation
  • ind_bif bifurcation index in br
  • lens2 second parameter used for the continuation, the first one is the one used to compute br, e.g. getLens(br)
  • options_cont = br.contparams arguments to be passed to the regular continuation

Optional arguments:

  • bdlinsolver bordered linear solver for the constraint equation
  • updateMinAugEveryStep update vectors a,b in Minimally Formulation every updateMinAugEveryStep steps
  • startWithEigen = false whether to start the Minimally Augmented problem with information from eigen elements
  • detectCodim2Bifurcation ∈ {0,1,2} whether to detect Bogdanov-Takens, Bautin and Cusp. If equals 1 non precise detection is used. If equals 2, a bisection method is used to locate the bifurcations.
  • kwargs keywords arguments to be passed to the regular continuation

where the parameters are as above except that you have to pass the branch br from the result of a call to continuation with detection of bifurcations enabled and index is the index of Hopf point in br you want to refine.

ODE problems

For ODE problems, it is more efficient to pass the Bordered Linear Solver using the option bdlinsolver = MatrixBLS()

where the options are as above except with have an additional parameter axis lens2 which is used to locate the bifurcation points.

See Temperature model (Simplest example) for an example of use.

Advanced use

Here, we expose the solvers that are used to perform newton refinement or codim 2 continuation in case the above methods fails. This is useful in case it is too involved to expose the linear solver options. An example of advanced use is the continuation of Folds of periodic orbits, see Continuation of Fold of periodic orbits.

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

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

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

Optional arguments:

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

Simplified call

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

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

The parameters / options are as usual except that you have to pass the branch br from the result of a call to continuation with detection of bifurcations enabled and index is the index of bifurcation point in br you want to refine. You can pass newton parameters different from the ones stored in br by using the argument options.

Jacobian transpose

The adjoint of the jacobian J is computed internally when Jᵗ = nothing by using transpose(J) which works fine when J is an AbstractArray. In this case, do not pass the jacobian adjoint like Jᵗ = (x, p) -> transpose(d_xF(x, p)) otherwise the jacobian will be computed twice!

ODE problems

For ODE problems, it is more efficient to pass the Bordered Linear Solver using the option bdlinsolver = MatrixBLS()

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

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

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

Optional arguments:

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

Simplified call:

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

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

The parameters / options are as usual except that you have to pass the branch br from the result of a call to continuation with detection of bifurcations enabled and index is the index of bifurcation point in br you want to refine. You can pass newton parameters different from the ones stored in br by using the argument options.

Jacobian transpose

The adjoint of the jacobian J is computed internally when Jᵗ = nothing by using transpose(J) which works fine when J is an AbstractArray. In this case, do not pass the jacobian adjoint like Jᵗ = (x, p) -> transpose(d_xF(x, p)) otherwise the jacobian will be computed twice!

ODE problems

For ODE problems, it is more efficient to pass the Bordered Linear Solver using the option bdlinsolver = MatrixBLS()

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

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

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

Optional arguments:

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

Simplified call

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

where the parameters are as above except that you have to pass the branch br from the result of a call to continuation with detection of bifurcations enabled and index is the index of Fold point in br that you want to continue.

Jacobian transpose

The adjoint of the jacobian J is computed internally when Jᵗ = nothing by using transpose(J) which works fine when J is an AbstractArray. In this case, do not pass the jacobian adjoint like Jᵗ = (x, p) -> transpose(d_xF(x, p)) otherwise the jacobian would be computed twice!

ODE problems

For ODE problems, it is more efficient to pass the Bordered Linear Solver using the option bdlinsolver = MatrixBLS()

Detection of Bogdanov-Takens and Cusp bifurcations

In order to trigger the detection, pass detectEvent = 1,2 in options_cont.

BifurcationKit.continuationHopfFunction

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

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

Optional arguments:

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

Simplified call:

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

where the parameters are as above except that you have to pass the branch br from the result of a call to continuation with detection of bifurcations enabled and index is the index of Hopf point in br that you want to refine.

ODE problems

For ODE problems, it is more efficient to pass the Bordered Linear Solver using the option bdlinsolver = MatrixBLS()

Jacobian transpose

The adjoint of the jacobian J is computed internally when Jᵗ = nothing by using transpose(J) which works fine when J is an AbstractArray. In this case, do not pass the jacobian adjoint like Jᵗ = (x, p) -> transpose(d_xF(x, p)) otherwise the jacobian would be computed twice!

Detection of Bogdanov-Takens and Bautin bifurcations

In order to trigger the detection, pass detectEvent = 1,2 in options_cont. Note that you need to provide d3F in prob.

References

  • Govaerts

    Govaerts, Willy J. F. Numerical Methods for Bifurcations of Dynamical Equilibria. Philadelphia, Pa: Society for Industrial and Applied Mathematics, 2000.

  • Blank

    Blank, H. J. de, Yu. A. Kuznetsov, M. J. Pekkér, and D. W. M. Veldman. “Degenerate Bogdanov–Takens Bifurcations in a One-Dimensional Transport Model of a Fusion Plasma.” Physica D: Nonlinear Phenomena 331 (September 15, 2016): 13–26. https://doi.org/10.1016/j.physd.2016.05.008.

  • Bindel

    Bindel, D., M. Friedman, W. Govaerts, J. Hughes, and Yu.A. Kuznetsov. “Numerical Computation of Bifurcations in Large Equilibrium Systems in Matlab.” Journal of Computational and Applied Mathematics 261 (May 2014): 232–48. https://doi.org/10.1016/j.cam.2013.10.034.