7. Numerical Continuation

The previous chapters used SimulationExperiment.run("auto-07p" | "bifurcationkit.jl" | "pyrates-bifurcation") as a black box. This chapter explains what those backends actually do and which numerical machinery TVBO routes your YAML system through.

Equilibria as roots

Continuation tracks equilibria of \(\dot x = f(x, a)\), i.e. solutions of

\[ f(x, a) = 0 \]

as the parameter \(a\) varies. For each \(a\) this is a nonlinear root-finding problem in \(x\). The basic tool is Newton’s method:

\[ x^{(k+1)} = x^{(k)} - J(x^{(k)}, a)^{-1} f(x^{(k)}, a), \qquad J = \partial f / \partial x. \]

Newton converges quadratically when started close enough to a root and when \(J\) is non-singular. Both conditions fail at fold points, which is exactly where bifurcations occur.

Pseudo-arclength continuation

To get past folds (where \(J\) becomes singular and a naive parameter sweep fails), continuation methods reformulate the problem so that the arclength \(s\) along the solution branch, not the parameter \(a\), is the independent variable. The augmented system is:

\[ \begin{cases} f(x, a) = 0 \\[4pt] \dot x(s)^\top \,(x - x_{\text{prev}}) + \dot a(s)\,(a - a_{\text{prev}}) - \Delta s = 0 \end{cases} \]

The second equation is the pseudo-arclength constraint: it pins the projection of the new point onto the previous tangent vector to a fixed step \(\Delta s\). The augmented Jacobian remains non-singular through the fold, allowing the solver to turn around smoothly.

Predictor-corrector

Each continuation step has two stages:

  1. Predictor: extrapolate along the tangent of the previous solution. \(\;(x_p, a_p) = (x_{\text{prev}}, a_{\text{prev}}) + \Delta s\,(\dot x, \dot a)\).
  2. Corrector: Newton-iterate on the augmented system until convergence.

This is the pseudo-arclength predictor-corrector algorithm at the heart of AUTO-07p, MATCONT, BifurcationKit.jl, and (indirectly via AUTO) PyCoBi.

Detecting bifurcations

Each backend monitors test functions along the branch:

Bifurcation Test function Meaning
Fold (saddle-node, LP) \(\det J = 0\) one real eigenvalue crosses 0
Hopf (HB) \(\mathrm{Re}\,\lambda_{1,2}(J) = 0\), \(\mathrm{Im} \neq 0\) complex pair crosses imaginary axis
Branch point (BP) rank drop in \(\big[\,J\;\partial f/\partial a\,\big]\) branches intersect
Period-doubling (PD) Floquet \(\mu = -1\) period of a limit cycle doubles
Torus / Neimark-Sacker (TR) Floquet \(|\mu| = 1\), \(\mu \neq \pm 1\) quasiperiodic torus emerges

These appear as labelled markers on every BifurcationResult.plot() output in this tutorial. The BIF_STYLES registry in tvbo/analysis/bifurcation.py defines the marker shape and colour for each type, consistently across all backends.

Branch switching

When a BP is detected, continuation can switch onto the new branch by restarting the Newton corrector with the eigenvector of the singular Jacobian as the search direction. In TVBO this is requested declaratively via branches: in the Continuation YAML:

branches:
  - name: po_from_hopf
    source_point: "hopf:all"
    bothside: true

source_point accepts:

  • "hopf:all": switch onto the periodic orbit from every detected Hopf.
  • "hopf:<index>": only the \(i\)-th Hopf.
  • "bp:<index>": switch at a specific branch point (e.g. for the pitchfork).
  • "lp:<index>": restart from a fold.

Backend landscape

Backend Language Strengths Caveats
AUTO-07p Fortran mature, robust, periodic-orbit continuation, two-parameter requires Fortran toolchain; clunky scripting
PyCoBi / PyRates Python wrapper around AUTO-07p declarative model spec requires NDIM >= 2
BifurcationKit.jl Julia modern, GPU-friendly, large-scale, BifurcationKit’s Newton-Krylov solvers needs Julia + juliacall
MatCont MATLAB normal-form coefficients, codim-2 detection not free, not auto-callable from TVBO

TVBO’s SimulationExperiment.run(backend) selects between them at runtime from the same YAML. The continuation methodology is identical across backends; only the numerics, the bifurcation-detection thresholds, and the ergonomics differ.

Two-parameter continuation

Once a codim-1 bifurcation point is located, it can be itself continued in a second parameter, tracing out a curve of fold or Hopf points in \((a_1, a_2)\) space. This is the gateway to codim-2 bifurcations (cusps CP, Bautin GH, Bogdanov-Takens BT, zero-Hopf ZH, double-Hopf HH), each of which appears as a labelled point in BIF_STYLES. Two-parameter continuation is supported by all three TVBO backends but its YAML syntax is left as a follow-up topic.

Reading the bifurcation diagram

Every plot in this tutorial follows the same conventions:

  • Solid line: stable branch (real part of all eigenvalues \(< 0\)).
  • Dashed line: unstable branch.
  • Marker: labelled bifurcation point (see BIF_STYLES).
  • Tube / wireframe: limit cycle in 3-D, with min-max envelope.

These are produced uniformly across all backends because TVBO normalises the result into a single BifurcationResult data class before plotting.