Ordinary Differential Equations (10): Bifurcation Theory

A lake stays clear for decades, then turns murky in a single season. A power grid hums along stably, then trips into a cascading blackout in seconds. A column under slowly increasing load is straight, straight, straight — and then suddenly buckles.

These are not prediction failures. The universe is doing exactly what dynamical systems theory says it must: cross a bifurcation. When a parameter drifts past a critical value, the topology of phase space rearranges, and what was once impossible becomes inevitable. This chapter classifies these rearrangements. There are only a few, and once you see the catalog, you’ll spot them everywhere.

What You Will Learn#

• What bifurcation means precisely, and why it has to be defined in terms of topology, not quantity
• The four codimension-1 normal forms: saddle-node, transcritical, supercritical pitchfork, subcritical pitchfork
• The Hopf bifurcation: how a stable spiral spawns a limit cycle when a pair of complex eigenvalues crosses the imaginary axis
• Why subcritical bifurcations are catastrophic (they jump and they hyster) while supercritical ones are gentle
• A first taste of global bifurcations — homoclinic orbits, SNIC, and how they open the door to chaos
• The ideas of codimension and universality that explain why nature reuses the same handful of normal forms

Prerequisites: stability and phase-plane analysis from Chapter 8 , and the chaos vocabulary from Chapter 9 .

What is a bifurcation, really?#

The word bifurcation (Latin furca = fork) was coined by Henri Poincaré around 1885. The intuition is geometric. Picture the long-term behaviour of a one-parameter system$\dot{x} = f(x,\mu)$ as a portrait that depends continuously on$\mu$ . For most values of$\mu$ , small perturbations move the portrait around but do not change its essential shape: the same number of fixed points, the same stability assignments, the same cycles. We call such$\mu$ structurally stable.

A bifurcation is the opposite: a value$\mu_c$ at which the portrait changes its shape discontinuously under arbitrarily small perturbation. New equilibria appear out of nowhere, two of them collide and annihilate, a stable focus becomes a limit cycle, or a periodic orbit doubles its period.

The shift is qualitative, not quantitative. A water pipe gradually narrowing is not a bifurcation. A water pipe suddenly bursting is.

A useful mental model#

Think of$f(x,\mu)$ as a landscape that depends on$\mu$ . The fixed points are where the gradient vanishes; their stability is the curvature there. As$\mu$ slides, the landscape morphs continuously — but whenever a hill flattens to a saddle and then flips into a valley, that moment is a bifurcation. Between bifurcations the landscape is just being pushed around.

The four codimension-1 normal forms#

A miracle of bifurcation theory: near any one-parameter bifurcation, the dynamics are locally equivalent (after a smooth change of coordinates) to one of just four canonical equations. These are the normal forms.

Why “codimension-1”? Because each requires tuning exactly one parameter to occur. To meet two of them at once, you need to tune two parameters, etc. Codimension-1 events are the ones you bump into generically when sliding a single dial.

Saddle-node (fold) bifurcation#

Normal form:$\dot{x} = \mu - x^2.$ Setting$\dot{x}=0$ gives$x^* = \pm\sqrt{\mu}$ . So:

The linearisation$f_x = -2x$ tells us the stability immediately: at$+\sqrt{\mu}$ we have$f_x = -2\sqrt{\mu} < 0$ (stable), and at$-\sqrt{\mu}$ we have$f_x = +2\sqrt{\mu} > 0$ (unstable). The two equilibria are born together as$\mu$ increases through 0.

Where it shows up

• Lasers: below threshold pump current, the only stationary state is “no light”. Above threshold, a coherent emitting state appears. The “no-light” state continues to exist but coexists with it.
• Neurons (Class I): below the rheobase current the resting state is the only equilibrium. Above it, the resting state collides with a saddle and disappears — the neuron starts firing.
• Lakes flipping from clear to murky: the clear-water equilibrium disappears in a saddle-node fold as nutrient loading crosses a threshold.

The signature of a fold is bistability before annihilation. Just below$\mu_c$ two equilibria coexist; just above$\mu_c$ neither does. The system must go somewhere else, often violently.

Transcritical bifurcation#

Normal form:$\dot{x} = \mu x - x^2.$ Two equilibria always exist:$x^*=0$ and$x^*=\mu$ . They never disappear — they merely swap stability as they cross at$\mu=0$ .

This is the bifurcation you get whenever the system has a “trivial” state ($x=0$ ) that must always remain an equilibrium for symmetry or definitional reasons.

Where it shows up

• Epidemiology: the disease-free equilibrium always exists. When the basic reproduction number$R_0$ crosses 1 (the bifurcation parameter), it loses stability to the endemic equilibrium. The transition is exactly transcritical.
• Population dynamics: extinction is always an equilibrium. As an environmental quality parameter crosses a threshold, the extinction state hands off stability to a positive coexistence state.
• Lasers (alternative model): the off-state is always a fixed point; it loses stability to the lasing state at threshold.

Supercritical pitchfork#

Normal form:$\dot{x} = \mu x - x^3.$ Equilibria: always$x^*=0$ , plus$x^*=\pm\sqrt{\mu}$ when$\mu>0$ . The trivial branch loses stability and two new stable branches are born symmetrically.

This is the universal symmetry-breaking bifurcation. The equation is invariant under$x \to -x$ , so any new equilibrium must come with a partner. Below$\mu_c$ the system sits on the symmetric solution; above$\mu_c$ it must commit to one of two equally valid asymmetric solutions.

Subcritical pitchfork (the dangerous one)#

Normal form:$\dot{x} = \mu x + x^3.$ Now the trivial branch loses stability without any nearby stable branch waiting to catch the system. Below$\mu=0$ we have$x^*=0$ stable plus two unstable branches at$\pm\sqrt{-\mu}$ ; above$\mu=0$ , only the unstable trivial branch remains. Trajectories shoot off to infinity.

In real systems higher-order terms eventually re-stabilise things: adding a$-x^5$ term gives the canonical hysteretic model$\dot{x} = \mu x + x^3 - x^5,$ which has a high-amplitude stable branch coexisting with the trivial state in a window$-\tfrac14 \le \mu \le 0$ . Slowly ramping$\mu$ produces the famous hysteresis loop: the system jumps to large amplitude when$\mu$ crosses 0 from below, and only jumps back down when$\mu$ is pulled past$-\tfrac14$ on the way back.

Why it matters

• Buckled columns (Euler buckling): a slender vertical rod under compression undergoes a supercritical pitchfork at the critical load — it bends a small amount one way or the other, reversibly.
• Snapping shells (von Karman buckling of cylinders): the pitchfork is subcritical. The shell sits straight, sits straight, then collapses with a bang into a heavily-deformed configuration. This is why aerospace engineers calculate buckling loads with safety factors of 3-10.
• Ferromagnets near the Curie point: supercritical pitchfork. Magnetisation grows continuously from zero as temperature drops.
• Climate tipping: glacial$\leftrightarrow$ interglacial transitions are often modelled as subcritical-pitchfork (or fold) bifurcations of a temperature-albedo system. The hysteresis window means a tipped state cannot be “untipped” simply by reverting CO$_2$ to its earlier value.

Summary table#

The Hopf bifurcation: a focus gives birth to a cycle#

The bifurcations above are scalar. The first genuinely two-dimensional bifurcation is the Hopf (Andronov-Hopf, really). It is what allows oscillations to appear.

Normal form (polar coordinates):$\dot{r} = \mu r - r^3, \qquad \dot{\theta} = \omega.$ The radial equation is exactly the supercritical pitchfork. So:

• For$\mu \le 0$ , the only attractor is$r=0$ – a stable spiral.
• For$\mu > 0$ , the origin is an unstable spiral and a stable limit cycle of radius$r=\sqrt{\mu}$ encircles it.

In Cartesian form, $\dot{x} = \mu x - \omega y - x(x^2+y^2),$ $\dot{y} = \omega x + \mu y - y(x^2+y^2).$ The Jacobian at the origin is$\bigl(\begin{smallmatrix}\mu & -\omega \\ \omega & \mu\end{smallmatrix}\bigr)$ , with eigenvalues$\lambda = \mu \pm i\omega$ . As$\mu$ crosses zero from below, the complex pair crosses the imaginary axis transversally — this is the Hopf condition.

The Hopf theorem (Hopf 1942). Let$\mathbf{x}_0$ be an equilibrium of$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mu)$ at$\mu=\mu_c$ , and suppose the Jacobian has a complex pair$\lambda(\mu) = \alpha(\mu) \pm i\beta(\mu)$ satisfying

1.$\alpha(\mu_c)=0$ (the pair is on the imaginary axis), 2.$\beta(\mu_c)\ne 0$ (it is genuinely complex, not double-zero), 3.$\alpha'(\mu_c)\ne 0$ (the pair crosses, not merely touches),

and assume the cubic first Lyapunov coefficient$\ell_1$ is non-zero. Then a one-parameter family of periodic orbits emerges from$\mathbf{x}_0$ at$\mu_c$ , with period$\approx 2\pi/\beta(\mu_c)$ . The cycle is stable if$\ell_1 < 0$ (supercritical) and unstable if$\ell_1>0$ (subcritical).

Subcritical Hopf is the cyclic version of the subcritical pitchfork: it is silent until the bifurcation, then suddenly ejects the system to a distant attractor. Aircraft wing flutter, sudden onset of cardiac arrhythmia, and mode switching in lasers are textbook subcritical Hopfs.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

import numpy as np
from scipy.integrate import odeint

def hopf(s, t, mu, omega=1.0):
    x, y = s
    r2 = x*x + y*y
    return [mu*x - omega*y - x*r2,
            omega*x + mu*y - y*r2]

t = np.linspace(0, 60, 6000)
sol = odeint(hopf, [0.05, 0.05], t, args=(0.4,))

Codimension and universality#

A bifurcation has codimension k if it requires tuning k independent parameters to occur generically. Codimension-1 events fill curves in parameter space; codimension-2 events occur at isolated points where two such curves meet.

Common codimension-2 bifurcations:

• Cusp: where two saddle-node curves meet tangentially. The unfolding is the cusp catastrophe$\dot{x} = \mu_1 + \mu_2 x - x^3$ , which contains a hysteresis region bounded by two saddle-node curves meeting at a point.
• Bogdanov-Takens: a double-zero eigenvalue. The unfolding contains saddle-node, Hopf, and homoclinic curves meeting at one point. Whenever you find a Hopf curve and a fold curve approaching each other in a parameter diagram, look for a BT point at the meeting.
• Bautin (generalised Hopf): where the first Lyapunov coefficient passes through zero, marking the boundary between supercritical and subcritical Hopf. Cycles fold over in a saddle-node-of-cycles bifurcation.

The deep reason these classifications exist is the centre-manifold theorem plus the method of normal forms. Near a codimension-k bifurcation, only a handful of variables (the centre directions) carry the slow dynamics; everything else is enslaved to them. After polynomial coordinate changes, the slow dynamics reduces to the universal normal form. This is why bifurcation theory is finite — almost a periodic table.

Global bifurcations: when the topology changes far from any equilibrium#

Local bifurcations rearrange phase space near a single point. Global bifurcations rearrange it on a large scale, typically by reconnecting invariant manifolds.

Homoclinic bifurcation#

A trajectory that leaves a saddle along its unstable manifold and returns along its stable manifold is a homoclinic orbit. It exists only at isolated parameter values (a codimension-1 phenomenon). Near a homoclinic bifurcation, periodic orbits nearby have arbitrarily long periods — the orbit spends ever more time creeping past the saddle. This produces the universal scaling$T \sim -\log|\mu - \mu_c|$ .

The Shilnikov theorem says: a homoclinic loop to a saddle-focus with appropriate eigenvalue ratio implies the existence of countably many periodic orbits of all periods in a neighborhood — in other words, chaos.

Heteroclinic bifurcation#

Same idea, but the orbit connects different saddles. Heteroclinic cycles can give rise to slow oscillations with extremely long, almost-pause-like phases near each saddle. They are the standard model for “winnerless competition” in neural circuits.

SNIC (saddle-node on invariant circle)#

A saddle-node bifurcation that happens on a closed invariant curve. Below the bifurcation, the curve is broken at the saddle-node pair; above it, the pair has annihilated and the curve is restored as a limit cycle. The limit cycle is born with infinite period (the system creeps through the place where the equilibria used to be) and the period scales as$T \sim 1/\sqrt{\mu_c - \mu}$ . SNIC is the second standard route from quiescence to firing in neurons (Class I excitability), and a key mechanism in El Niño oscillator models.

The route to chaos: period doubling#

Limit cycles can themselves bifurcate. The most famous route is the period-doubling cascade: a stable cycle of period$T$ loses stability and gives birth to a stable cycle of period$2T$ , which in turn doubles to$4T$ , then$8T$ , then$16T$ , accumulating at a finite parameter value beyond which the dynamics is chaotic.

Mitchell Feigenbaum discovered (1978) that the parameter intervals$\Delta_n = \mu_n - \mu_{n-1}$ between successive doublings shrink geometrically with a universal ratio$\delta = \lim_{n \to \infty} \frac{\Delta_n}{\Delta_{n+1}} = 4.6692016\ldots$ independent of the specific system, as long as it has a smooth quadratic maximum. The same constant governs period doubling in dripping faucets, Rayleigh-Bénard convection, and electronic circuits.

The minimal toy model is the logistic map$x_{n+1} = r x_n(1 - x_n)$ :

Sharkovsky’s theorem (1964) and the famous corollary “period 3 implies chaos” (Li-Yorke 1975) round out the story: any continuous interval map with a period-3 orbit has orbits of every other period, and uncountably many aperiodic orbits.

For a deeper dive into the chaos that lives beyond the cascade, see Chapter 9 .

Numerical detection and continuation#

In practice we rarely have closed-form normal forms. We have a vector field$\mathbf{f}(\mathbf{x},\mu)$ and want to map out its bifurcations as$\mu$ varies. The standard tools are continuation methods:

1. Track an equilibrium branch. Start at a known equilibrium, then use Newton’s method on$\mathbf{f}(\mathbf{x},\mu)=0$ as$\mu$ is incremented. Use pseudo-arclength continuation to handle folds: parametrise the branch by arclength rather than by$\mu$ , so the algorithm can turn corners.
2. Monitor the Jacobian eigenvalues along the branch. A real eigenvalue crossing 0 flags a saddle-node, transcritical, or pitchfork (you tell them apart by symmetry and by the second-order coefficient). A complex pair crossing the imaginary axis flags a Hopf.
3. Switch branches at detected bifurcations using normal-form coefficients to compute the tangent direction of the new branch.

Production tools: AUTO-07p (the gold standard, Fortran/C), MATCONT (MATLAB), PyDSTool and BifurcationKit.jl (Python/Julia). For research-grade bifurcation work these are essentially mandatory; rolling your own continuation algorithm is a lot of work to get right.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

import numpy as np
from scipy.linalg import eigvals

def detect_bifurcations(f, fx, x_branch, mu_grid):
    """
    Track eigenvalue crossings along a branch of equilibria.

    Parameters
    ----------
    f, fx        : RHS and its Jacobian, both functions of (x, mu).
    x_branch     : array of equilibria, one per mu in mu_grid.
    mu_grid      : parameter values.

    Returns
    -------
    list of dicts, one per detected eigenvalue crossing.
    """
    events = []
    prev = None
    for mu, x in zip(mu_grid, x_branch):
        eigs = eigvals(fx(x, mu))
        if prev is not None:
            # Real crossings -> fold/transcritical/pitchfork
            for a, b in zip(np.sort(prev.real), np.sort(eigs.real)):
                if a * b < 0 and abs(a.imag) + abs(b.imag) < 1e-8:
                    events.append({'kind': 'real_crossing', 'mu': mu})
            # Complex-pair crossings -> Hopf
            prev_pair = prev[np.iscomplex(prev)]
            curr_pair = eigs[np.iscomplex(eigs)]
            if prev_pair.size and curr_pair.size:
                if prev_pair[0].real * curr_pair[0].real < 0:
                    events.append({'kind': 'hopf', 'mu': mu,
                                   'omega': abs(curr_pair[0].imag)})
        prev = eigs
    return events

Why this matters#

The deepest message of bifurcation theory is that smooth causes can produce abrupt effects, but only through a small number of canonical mechanisms. When you suspect a system is approaching a tipping point, you can ask concrete diagnostic questions:

• What kind of bifurcation is approaching? Increasing variance and slow recovery from perturbations (critical slowing down) signal a fold or pitchfork. Growing oscillations signal a Hopf.
• Is it super- or sub-critical? If the system is bistable as you approach, it is sub-critical, and the post-bifurcation jump will be large and possibly irreversible.
• Is there hysteresis? If yes, do not assume that reversing the parameter will restore the original state.
• Are there early warning signals? For folds, the recovery rate from small perturbations decays towards zero as$\mu \to \mu_c$ , with universal scaling$\sim |\mu - \mu_c|^{1/2}$ . This is now used in ecology, climate science, and even epidemiology to forecast tipping events.

The taxonomy is small. The phenomena are everywhere.

Exercises#

1. Imperfect pitchfork. Show that$\dot{x} = h + \mu x - x^3$ has a saddle-node bifurcation curve in the$(\mu, h)$ plane that meets at a cusp at the origin. Sketch the bifurcation diagrams for$h>0$ ,$h=0$ , and$h<0$ .
2. Bifurcations of$\dot{x} = \mu - x - e^{-x}$ . Find all equilibria implicitly via$\mu = x + e^{-x}$ . Show there is a saddle-node at$\mu_c = 1$ ,$x_c = 0$ , and identify the stability of each branch.
3. Logistic period-doubling numerically. Compute the doubling parameters$r_n$ for the logistic map up to$n=6$ . Use$\delta_n := (r_n - r_{n-1})/(r_{n+1} - r_n)$ to estimate the Feigenbaum constant.
4. Hopf in a predator-prey model. For$\dot{x} = x(1-x) - \tfrac{xy}{a+x},\;\dot{y} = -dy + \tfrac{xy}{a+x}$ , find the parameter combinations producing a Hopf bifurcation of the coexistence equilibrium. Decide super- vs sub-critical numerically.
5. Subcritical pitchfork with$x^5$ saturation. For$\dot{x} = \mu x + x^3 - x^5$ , derive the locations of the saddle-node folds at$\mu = -1/4$ and the resulting hysteresis interval. Plot the loop traversed by quasi-statically ramping$\mu$ up and back.

References#

• Strogatz, S. H. (2015). Nonlinear Dynamics and Chaos, 2nd ed. Westview / CRC. The single best entry point.
• Kuznetsov, Y. A. (2004). Elements of Applied Bifurcation Theory, 3rd ed. Springer. The reference for codimension-1 and -2 normal forms.
• Guckenheimer, J. & Holmes, P. (1983). Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. Springer.
• May, R. M. (1976). “Simple mathematical models with very complicated dynamics.” Nature 261, 459-467.
• Scheffer, M. et al. (2009). “Early-warning signals for critical transitions.” Nature 461, 53-59.
• Doedel, E. & Oldeman, B. (2012). AUTO-07p Continuation Software for ODEs.

Previous Chapter: Chapter 9: Chaos Theory and the Lorenz System

Next Chapter: Chapter 11: Numerical Methods for Differential Equations

This is Part 10 of the 18-part series on Ordinary Differential Equations.