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 failures of prediction. They are the universe doing exactly what dynamical systems theory says it must do: cross a bifurcation. When a parameter drifts past a critical value, the topology of phase space rearranges itself, and what was once impossible becomes inevitable. This chapter is about classifying those rearrangements. There turn out to be only a handful of them, and once you see the catalogue you start spotting 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 .

1. What is a bifurcation, really?

The word bifurcation (Latin furca = fork) was coined by Henri Poincare 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.

2. The four codimension-1 normal forms

A miracle of bifurcation theory: near any one-parameter bifurcation, the dynamics is 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.

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

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

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

2.4 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

3. 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
12

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,))
# sol now spirals out and lands on the limit cycle r = sqrt(0.4) ~ 0.632

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.

5. 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 Nino oscillator models.

6. 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-Benard 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 .

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

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