# Bifurcations of a first-order equation, from blue sky to hysteresis Take an autonomous equation with a knob on it, $ \frac{dy}{dt} = f(y;\,p), \qquad y(t_0) = y_0 \in \mathbb R , $ where $p$ is a parameter that does not change in time. For each fixed $p$ the phase line tells the whole story: find the equilibria (the roots of $f$), read the sign of $f$ between them, and every solution moves monotonically toward the next equilibrium in that direction or off to infinity. Now turn the knob. Usually nothing interesting happens; the equilibria slide a little and the phase line keeps the same shape. But at special values of $p$ the shape changes: equilibria appear, disappear, or swap stability. That change is a **bifurcation**, and the parameter value where it happens is a **bifurcation value**. At a bifurcation value the flow fails to be *structurally stable*: an arbitrarily small change in $p$ changes the qualitative picture. Why would a parameter be fixed while the state moves? Think of a glass of water that has sat in a room long enough to be at room temperature. Change the thermostat suddenly and the glass takes a while to catch up. Change it *slowly enough* and the glass never visibly lags: at every moment it sits at the equilibrium for the current room temperature. The room temperature is then a slowly varied parameter, and the glass's state is read off the phase line for that value. That is the point of view on this page. We watch a stable equilibrium, move $p$ slowly, and ask what we would observe. Pushed far enough, the theory explains why some systems remember where they have been. Two facts carry the whole page. First, a bifurcation of a smooth scalar equation can only happen where an equilibrium has $f'(y^*) = 0$, since a root with $f' \neq 0$ survives small changes in $p$ (implicit function theorem) and keeps its stability. Second, $f'(y^*)$ is not only the stability test but a *rate*: near a sink a disturbance decays like $e^{f'(y^*)\,t}$, so $ \tau = \frac{1}{|f'(y^*)|} $ is the **relaxation time**. Put the two together and something is visible *before* every bifurcation: as $p$ approaches the bifurcation value, $f'(y^*) \to 0$ and $\tau \to \infty$. The system recovers from disturbances more and more slowly. That is **critical slowing down**, and every applet on this page shows it. > [!info] Reading the applets > Top row, left to right: $f$ plotted *with the state vertical*, the phase line, and solutions from fifteen starting values, all on the same vertical scale so you can read across. Solutions are coloured by the equilibrium they settle on; red ones run off to infinity. Filled circle = sink, open = source, half-filled = semi-stable (filled on the side it attracts from). > > Bottom row: the **bifurcation diagram** (equilibria against the parameter, solid = sink, dashed = source), with a dashed blue marker at the current parameter. Stand on the marker and look along it: what you see is the phase line. Beside it, $\tau = 1/|f'|$ along the stable branches, on a log scale. Tick marks on the parameter axis are the bifurcation values. > > **Kick** displaces each stable state by $\delta$ and times its return to within $\delta/20$; the dashed orange curve is what the linearization predicts. **Ramp** sweeps the parameter slowly up and back with one state carried along (and a little noise, as any real system has); red is the way up, purple the way down. ## The saddle-node, or blue-sky, bifurcation The simplest example is $ \frac{dy}{dt} = f_a(y) = a + y^{2}, \qquad a \in \mathbb R . $ Equilibria solve $a + y^2 = 0$: $ y_{1} = -\sqrt{-a}, \qquad y_{2} = +\sqrt{-a}, \qquad a \le 0 . $ For $a < 0$ the graph of $f_a$ is an upward parabola crossing the axis twice. Since $f_a'(y) = 2y$, the negative root has $f' = -2\sqrt{-a} < 0$ and is a **sink**; the positive root has $f' = +2\sqrt{-a} > 0$ and is a **source**. Raise $a$ and the parabola lifts. At $a = 0$ it touches the axis at a single double root, $y = 0$, with $f'(0) = 0$; the flow is upward on both sides, so the merged point attracts from below and repels above. It is **semi-stable**. For $a > 0$ the parabola clears the axis entirely, there are no equilibria, and every solution increases without bound. In fact it blows up in finite time, since $y' = a + y^2 \ge y^2$. Read in reverse, the picture explains the name: lowering $a$ through zero, a pair of equilibria materialises out of the clear blue sky. **Saddle-node** (from the two-dimensional picture), **fold** and **tangent** bifurcation are other names for the same event. In the $(a, y)$ plane the equilibria trace the parabola $a = -y^2$ lying on its side: the lower half solid (sinks), the upper half dashed (sources), meeting at the nose $(0, 0)$. <div class="applet" data-applet="bifurcation-explorer" data-model="saddle-node" data-models="saddle-node"></div> Things worth trying: - Start at $a = -1$ and drag $a$ toward $0$. The open and filled circles approach each other, and the band of starting values that escape to $+\infty$ (red) grows downward until it is everything above the sink. - Press **a = 0**. The surviving point is half-filled. Solutions from below creep up to it; solutions from above leave. - Press **kick** at $a = -1$, then at $a = -0.01$, then at $a = 0$. The return takes $1.4$, then $10.8$, then $63$ time units. At $a = 0$ the linearization predicts no return at all ($f' = 0$), and the actual return is algebraic: $y' = y^2$ gives $y = y_0/(1 - y_0 t)$, which for $y_0 < 0$ decays like $1/t$. - Press **ramp**. The state follows the sink up toward $0$, runs out of equilibrium, and leaves. On the way back nothing returns it: by then it is above the source. ## The transcritical bifurcation The second normal form does not change how *many* equilibria there are, only which is stable: $ \frac{dy}{dt} = f_a(y) = y\,(a - y), \qquad a \in \mathbb R . $ The equilibria are $y = 0$ and $y = a$, for every $a$. With $f_a'(y) = a - 2y$, $ f_a'(0) = a, \qquad f_a'(a) = -a . $ The two derivatives are negatives of each other, so the two equilibria always have opposite stability. For $a < 0$ the origin is the sink and $y = a$ the source; for $a > 0$ they have swapped. At $a = 0$ they coincide, $f_0(y) = -y^2$, and the merged point is semi-stable (attracting from above). In the bifurcation diagram two straight lines cross at the origin and **exchange stability** as they pass through each other. Nothing is created and nothing is destroyed. <div class="applet" data-applet="bifurcation-explorer" data-model="transcritical" data-models="transcritical"></div> - Watch the relaxation-time panel: $\tau = 1/|a|$ on both sides, diverging at $a = 0$ like $|a|^{-1}$. That is a *faster* divergence than the saddle-node's $|a|^{-1/2}$. - Ramp $a$ up and back. The state rides $y = 0$, then (a little late, because the origin is only weakly unstable just past $a = 0$) peels off onto $y = a$, and comes back down the same way. The loop closes. ## Harvesting: both bifurcations in one pond A logistic population with a harvest is where these two normal forms stop being abstract. Take $ \textbf{Model A (constant quota):}\quad \frac{dP}{dt} = rP\left(1-\frac{P}{N}\right) - H, \qquad \textbf{Model B (constant effort):}\quad \frac{dP}{dt} = rP\left(1-\frac{P}{N}\right) - HP . $ They *are* the two normal forms, after a change of variables. Complete the square in Model A, $rP(1 - P/N) = \tfrac{rN}{4} - \tfrac{r}{N}\left(P - \tfrac N2\right)^2$, and set $ y = -\frac{r}{N}\left(P - \frac N2\right), \qquad a = \frac{r}{N}\left(H - H^{*}\right), \qquad H^{*} = \frac{rN}{4} \quad\Longrightarrow\quad \frac{dy}{dt} = a + y^{2}. $ Model A is a saddle-node in disguise, with the fold at the quota $H^* = rN/4$ (the minus sign flips the picture, so the healthy stock above $N/2$ is the sink at negative $y$). For Model B set $ y = \frac{r}{N}P, \qquad a = r - H \quad\Longrightarrow\quad \frac{dy}{dt} = y\,(a - y), $ a transcritical bifurcation at $H = r$. Same fish, same pond; the way the catch is written decides which normal form you get, and the two fail in entirely different ways. A quota collapses the stock suddenly, in finite time, with no warning in the stock level and no recovery when the quota is lowered again. Fixed effort produces a visible, gradual, reversible decline. <div class="applet" data-applet="harvest-models"></div> The full treatment, including why $P = 0$ being an equilibrium for every $H$ forces Model B to be transcritical, and a careful look at critical slowing down at the fold, is on [[Harvesting and Critical Slowing Down]]. ## The supercritical pitchfork The two examples so far had at most two equilibria because $f$ was quadratic. Cubic nonlinearities make room for three, and for a bifurcation that only happens in the presence of a symmetry: $ \frac{dy}{dt} = f_a(y) = y\,(a - y^{2}), \qquad a \in \mathbb R . $ The right-hand side is odd in $y$, so the equation does not change under $y \mapsto -y$: whatever happens on one side happens on the other. The origin is always an equilibrium, with $f_a'(0) = a$. For $a > 0$ there are two more, $y = \pm\sqrt a$, with $ f_a'(\pm\sqrt a) = a - 3a = -2a < 0 . $ So for $a < 0$ there is one sink, the origin. At $a = 0$ the origin loses stability and, at the same instant and symmetrically, two new sinks grow out of it along the parabola $a = y^2$. The diagram looks like a pitchfork. Near $a = 0$ the origin changes stability as in a transcritical bifurcation, while the new pair arrives together, as a pair does at a fold; but neither of those alone produces this picture, because the symmetry forces the two new branches to appear *at* the origin and to be mirror images. Exactly at $a = 0$ the equation is $y' = -y^3$. The origin is still a sink, but $f'(0) = 0$, and the approach is slow: $ y(t) = \frac{y_0}{\sqrt{1 + 2y_0^{2}\,t}} \sim \frac{1}{\sqrt{2t}} . $ <div class="applet" data-applet="bifurcation-explorer" data-model="pitchfork" data-models="pitchfork"></div> - Kick at $a = 0.5$: back within $\delta/20$ in about $2.5$. Kick at $a = 0$: about $2200$. To get ten times closer at $a = 0$ takes a hundred times longer. - Ramp $a$ up a few times. The state leaves the origin past $a = 0$ and falls to $+\sqrt a$ or $-\sqrt a$, and which one is decided by the noise. The model is symmetric; the outcome of any single run is not. That is **symmetry breaking**. ### A bead on a rotating hoop The pitchfork is easy to see in a mechanical system. A bead slides on a vertical hoop of radius $r$ that spins about its vertical diameter at angular rate $\omega$; let $\phi$ be the bead's angle from the bottom and $b$ a friction coefficient. Gravity pulls the bead down, the centrifugal effect of the spin pushes it outward, and Newton's law along the hoop gives $ m r\,\ddot\phi = -b\,\dot\phi - m g \sin\phi + m r \omega^{2}\sin\phi\cos\phi . $ Measure time in units of $b/(mg)$ and write $ \gamma = \frac{r\omega^{2}}{g}, \qquad \varepsilon = \frac{m^{2} g r}{b^{2}} \qquad\Longrightarrow\qquad \varepsilon\,\phi'' + \phi' = \sin\phi\,(\gamma\cos\phi - 1). $ When friction is strong ($\varepsilon \ll 1$) the inertia term is negligible after a brief transient and the bead obeys the first-order equation $\phi' = \sin\phi\,(\gamma\cos\phi - 1)$. The bottom $\phi = 0$ and the top $\phi = \pi$ are always equilibria, and $\cos\phi^* = 1/\gamma$ adds a symmetric pair once $\gamma > 1$. Expanding near the bottom, $ \phi' = (\gamma - 1)\,\phi - \frac{4\gamma - 1}{6}\,\phi^{3} + O(\phi^{5}), $ which is the pitchfork normal form with $a = \gamma - 1$. Slow spin: gravity wins and the bead sits at the bottom. Past the critical rate $\omega_c = \sqrt{g/r}$ the bottom goes unstable and the bead rides up to $\phi^* = \pm\arccos(1/\gamma)$, with $f'(\phi^*) = (1 - \gamma^2)/\gamma$. The dependent variable is the bead's angle, the parameter is the spin rate, and the two stable states are the two sides of the hoop; a single experiment shows only one of them, and which one is up to whatever small disturbance got there first. <div class="applet" data-applet="bead-hoop"></div> - Set $\gamma = 0.6$ and nudge the bead: it settles back in a few time units. Set $\gamma = 1$ and nudge it: thirty time units later it is still visibly displaced, because near the bottom $\phi' \approx -\phi^3/2$ and the return is algebraic. - Ramp $\gamma$ several times. The bead leaves the bottom on a different side from run to run. - Tick **inertia**. The bead now overshoots and rings, as in the video below, but it comes to rest in the same places. Inertia changes how the bead gets there, not where "there" is: the equilibria and their stability are those of the first-order equation. A video of the experiment, with rings spinning at increasing rates: [AppDynSys: Spinning Hoops](https://www.youtube.com/watch?v=kS9WoYj2AaY). ## Breaking the symmetry: the imperfect pitchfork Real hoops are never perfectly vertical, and real beams never perfectly straight. Add a small constant to the pitchfork, $ \frac{dy}{dt} = h + a\,y - y^{3}, $ and the symmetry $y \mapsto -y$ is gone: $y = 0$ is no longer an equilibrium when $h \ne 0$. Equilibria satisfy $a = (y^3 - h)/y$, and folds sit where $f$ and $f'$ vanish together: $ h + ay - y^{3} = 0, \quad a - 3y^{2} = 0 \quad\Longrightarrow\quad a_{\text{fold}} = 3\left(\frac{|h|}{2}\right)^{2/3}, \qquad y_{\text{fold}} = -\operatorname{sign}(h)\left(\frac{|h|}{2}\right)^{1/3}. $ For $h > 0$ the pitchfork comes apart into two pieces. The upper branch is connected and smooth: raise $a$ and the state glides from near zero up to $+\sqrt a$ with no bifurcation at all. The other side of the old pitchfork appears out of nowhere, a saddle-node pair born at $a_{\text{fold}}$. The imperfection chooses the side in advance. Plot the folds in the $(a, h)$ plane and they trace the **cusp** $ 27\,h^{2} = 4\,a^{3}, $ inside which there are three equilibria (two sinks and a source) and outside which there is one. This two-parameter picture organises everything in the rest of the page: a horizontal line through it (fixed $h$, sweep $a$) is an imperfect pitchfork; a vertical line through it (fixed $a$, sweep $h$) is the S-curve below. <div class="applet" data-applet="bifurcation-explorer" data-model="imperfect" data-models="pitchfork,imperfect"></div> - Set $h = 0.05$ and ramp $a$. The state always ends on the upper branch, because there is no choice to make. - Set $h = 0$: the explorer is the perfect pitchfork again, and the trivial branch $y = 0$ reappears. - Drag the point in the $(a, h)$ plane across the cusp and watch the number of equilibria change from one to three. - Look at the relaxation time. On the connected branch $\tau$ rises near $a = 0$ but stays finite; it diverges only at the fold, on the branch that is about to be born or to die. With an imperfection, slowing down is a warning about the fold, not about the old symmetric point. ## Hysteresis ### A subcritical pitchfork, held in check Flip the sign of the cubic and the pitchfork turns around: $y' = ay + y^3$ has its nonzero branches for $a < 0$, both unstable, and for $a > 0$ every solution off the origin blows up. That **subcritical** pitchfork is dangerous on its own. Add a stabilising quintic, $ \frac{dx}{dt} = f_r(x) = r\,x + x^{3} - x^{5}, \qquad r \in \mathbb R, $ and the runaway is caught. The equilibria are $x = 0$ (with $f_r'(0) = r$) and $ x^{2} = \frac{1 \pm \sqrt{1 + 4r}}{2}, \qquad r \ge -\tfrac14 . $ The outer pair, $x^2 = \tfrac12\left(1 + \sqrt{1 + 4r}\right)$, are sinks. The inner pair, $x^2 = \tfrac12\left(1 - \sqrt{1 + 4r}\right)$, are sources that exist only for $-\tfrac14 \le r \le 0$; they shrink into the origin at $r = 0$ (the subcritical pitchfork) and meet the outer pair at $ r_s = -\tfrac14, \qquad x = \pm \tfrac{1}{\sqrt 2} $ in a pair of folds. Read the diagram as $r$ increases. For $r < r_s$ only the origin exists and it is stable. Two folds at $r_s$ create the outer sinks and the inner sources, out of the blue sky. At $r = 0$ the inner sources close in on the origin and destabilise it. In the window $ -\tfrac14 < r < 0 $ there are **three stable states at once**, the origin and $\pm x_{\text{outer}}$, and which one a system occupies depends on its history. Run the experiment. Start at $r = -\tfrac12$ with $x = 0$ and raise $r$ slowly. The origin stays stable all the way to $r = 0$, where it fails, and the state *jumps* to the outer branch near $x = \pm 1$. Now lower $r$. The outer branch is still there and still stable, so the state stays on it, below $r = 0$, through the window, all the way down to $r_s = -\tfrac14$, where the branch folds away and the state drops back to $0$. Up at one parameter value, down at another. The path up and the path down do not coincide, and the area between them is a memory of where the system has been. This is **hysteresis**. ### The S-curve The cleanest picture of hysteresis is a vertical slice through the cusp: $ \frac{dy}{dt} = h + y - y^{3} \qquad (a = 1,\ \text{sweep } h). $ The equilibria lie on the S-shaped curve $h = y^3 - y$, whose turning points are folds at $ y = \pm\frac{1}{\sqrt3}, \qquad h = \mp\frac{2}{3\sqrt3} \approx \mp 0.385 . $ Between the folds there are two sinks (an upper and a lower branch) separated by a source. Sweep $h$ up from $-0.8$: the state rides the lower branch until it ends at $h = +0.385$, then jumps to the upper branch at $y = 2/\sqrt3$. Sweep back down: it rides the upper branch until *that* ends at $h = -0.385$, then jumps down. The loop has width $4/(3\sqrt3) \approx 0.77$ in $h$. <div class="applet" data-applet="bifurcation-explorer" data-model="quintic" data-models="quintic,s-curve"></div> - Ramp the quintic. The red path climbs at $r = 0$, the purple path drops at $r = -\tfrac14$, and the readout reports the loop. Kick the outer state at $r = -0.24$, just above the fold: it takes a long time to come back. The relaxation-time panel says why: $\tau$ diverges at both ends of the window, at $r = 0$ on the inner branch and at $r = -\tfrac14$ on the outer ones. - Both jumps happen a little *past* the bifurcation value, not at it. At a fold the state is slow (critical slowing down again) and the parameter keeps moving while it decides to leave; at the subcritical pitchfork the origin is only weakly unstable at first. The slower the ramp, the closer the jump comes to the bifurcation value. - Switch to the S-curve and ramp it with **noise: none**, then **larger**. With more noise the jumps come *early*: near a fold the barrier between the two sinks (the source) is close and shallow, and a big enough disturbance carries the state over it before the branch actually ends. This is noise-induced tipping. - On the S-curve, lower $a$ (the second slider, or drag in the $(a, h)$ plane). The folds move toward each other and the loop narrows; at $a = 0$ they meet at the tip of the cusp and the hysteresis is gone. ### Where hysteresis shows up Agar, the jelly made from red algae, melts at about $85\,^\circ\text{C}$ ($185\,^\circ\text{F}$) but does not set again until it has cooled to roughly $32$–$40\,^\circ\text{C}$ ($90$–$104\,^\circ\text{F}$). Between those temperatures it can be solid or liquid, and which one you find depends on whether you arrived from above or from below: a system with two stable branches and a window between two folds. The same structure appears, with different physics, in: 1. **Supercooled water**, which stays liquid below its freezing point until a disturbance tips it: [Supercooled Water – Explained](https://www.youtube.com/watch?v=ph8xusY3GTM) (Veritasium). 2. **Shape-memory alloys**, whose crystal structure switches between two phases at different temperatures on heating and cooling: [Nitinol: The Shape Memory Effect and Superelasticity](https://www.youtube.com/watch?v=wI-qAxKJoSU), [How NASA Reinvented the Wheel – Shape Memory Alloys](https://www.youtube.com/watch?v=2lv6Vs12jLc). 3. **Biological viscoelastic materials**, such as the cornea, whose response to an air puff differs between loading and unloading; the gap is measured clinically as *corneal hysteresis*: [Corneal Biomechanics](https://www.youtube.com/watch?v=TA4ab5GKkYM), [How the Ocular Response Analyzer works](https://www.youtube.com/watch?v=gfHr_XC0cYI), [Corneal Hysteresis in Glaucoma Progression](https://www.youtube.com/watch?v=UnUmXoS3h54). ## Critical slowing down, at every bifurcation Every bifurcation on this page happens where $f'(y^*) = 0$, so every one of them is preceded by a diverging relaxation time. How fast it diverges, and how a disturbance decays *at* the bifurcation value, depends on the leading nonlinear term: | bifurcation | normal form | $\tau$ on the stable branch near $p_c$ | at $p = p_c$ | decay of a disturbance at $p_c$ | |---|---|---|---|---| | saddle-node (fold) | $y' = a + y^2$ | $\dfrac{1}{2\sqrt{-a}} \propto \lvert a\rvert^{-1/2}$ | $y' = y^2$ | $\sim 1/t$, from one side only | | transcritical | $y' = y(a - y)$ | $\dfrac{1}{\lvert a\rvert}$ | $y' = -y^2$ | $\sim 1/t$, from one side only | | supercritical pitchfork | $y' = y(a - y^2)$ | $\dfrac{1}{\lvert a\rvert}$ for $a<0$, $\dfrac{1}{2a}$ for $a>0$ | $y' = -y^3$ | $\sim 1/\sqrt{2t}$, from both sides | | fold of the S-curve or quintic | locally $y' = c\,(p - p_c) - d\,(y - y_c)^2$ | $\propto \lvert p - p_c\rvert^{-1/2}$ | quadratic | $\sim 1/t$, from one side only | Away from $p_c$ the linearization is right, and the return time to within a fraction of the kick is $\tau \ln(\cdot)$. At $p_c$ it is silent, and the leading nonlinear term takes over with algebraic rather than exponential decay. The explorer below holds all six models; kick each one at a few parameter values approaching its bifurcation and compare. <div class="applet" data-applet="bifurcation-explorer" data-model="saddle-node"></div> > [!note] Why this matters > The level of a system near a fold usually says little about how close the fold is: the stock in Model A sits at half its carrying capacity or more right up to the collapse. The *recovery* says a lot. A real system is knocked around constantly, and as $\tau$ grows the knocks fade more slowly, so its fluctuations become larger and more correlated from one observation to the next. Rising variance and rising autocorrelation are the early-warning signals of the tipping-point literature (M. Scheffer et al., "Early-warning signals for critical transitions," *Nature* 461, 53–59, 2009). The imperfect pitchfork adds a caution: the warning belongs to the branch that is about to disappear, and a system that is sitting on the smooth, connected branch shows no such warning, because it has nothing to tip over. --- ## Summary | | saddle-node | transcritical | pitchfork (super) | imperfect pitchfork | subcritical + quintic | S-curve | |---|---|---|---|---|---|---| | normal form | $a + y^2$ | $y(a - y)$ | $y(a - y^2)$ | $h + ay - y^3$ | $rx + x^3 - x^5$ | $h + y - y^3$ | | what happens | a sink and a source collide and vanish | two equilibria cross and trade stability | one sink becomes three equilibria | the pitchfork splits into a smooth branch and a fold | jump at $r = 0$, return at $r = -\tfrac14$ | jumps at $h = \pm 2/(3\sqrt3)$ | | bifurcation value | $a = 0$ | $a = 0$ | $a = 0$ | $a = 3(\lvert h\rvert/2)^{2/3}$ | $r = 0$, $r_s = -\tfrac14$ | $h = \pm 0.385$ | | needs | nothing | an equilibrium that persists for all $p$ | a symmetry $y \mapsto -y$ | a broken symmetry | a subcritical pitchfork and a fold | two folds | | memory? | escape is irreversible | no | no | no | yes (hysteresis) | yes (hysteresis) | ## The applets in code Each applet reduced to the few lines that compute what it shows. Fully working versions, with the checks against the closed forms on this page and the figures, are in the `fourier-codes` repository as `bifurcation_explorer` and `bead_hoop` (`python/`, `matlab/`, `r/`, `mathematica/`). > [!example]- `bifurcation-explorer` — the S-curve: its bifurcation diagram drawn as $h = y^3 - ay$, and a slow ramp of $h$ up and back that traces the hysteresis loop > > > [!info]- Python > > ```python > > import numpy as np, matplotlib.pyplot as plt > > a = 1.0 > > f = lambda y, h: h + a * y - y**3 > > fy = lambda y: a - 3 * y**2 # stability: fy < 0 is a sink > > y = np.linspace(-1.6, 1.6, 1201); h = y**3 - a * y # every equilibrium, exactly > > plt.plot(np.where(fy(y) < 0, h, np.nan), y, "k") # sinks solid > > plt.plot(np.where(fy(y) > 0, h, np.nan), y, "k--") # source dashed > > rng = np.random.default_rng(1); dt, LEG, sig = 0.01, 400.0, 0.01 > > H = np.r_[np.linspace(-0.8, 0.8, int(LEG / dt)), np.linspace(0.8, -0.8, int(LEG / dt))] > > Y = np.empty_like(H); Y[0] = -1.2756 # start on the lower branch > > for k in range(1, len(H)): # Euler–Maruyama, h moving slowly > > Y[k] = Y[k-1] + dt * f(Y[k-1], H[k-1]) + sig * np.sqrt(dt) * rng.standard_normal() > > n = len(H) // 2 > > plt.plot(H[:n], Y[:n], "r", H[n:], Y[n:], "purple"); plt.xlabel("h"); plt.ylabel("y"); plt.show() > > print("jump up near h =", H[:n][np.argmax(Y[:n] > 0)], " down near h =", H[n:][np.argmax(Y[n:] < 0)], > > " folds at ±", 2 / (3 * np.sqrt(3))) > > ``` > > > [!info]- MATLAB / Octave > > ```matlab > > a = 1; f = @(y, h) h + a*y - y.^3; fy = @(y) a - 3*y.^2; > > y = linspace(-1.6, 1.6, 1201); h = y.^3 - a*y; % every equilibrium, exactly > > hs = h; hs(fy(y) > 0) = NaN; hu = h; hu(fy(y) < 0) = NaN; > > plot(hs, y, 'k', hu, y, 'k--'); hold on > > dt = 0.01; LEG = 400; sig = 0.01; m = round(LEG/dt); > > H = [linspace(-0.8, 0.8, m), linspace(0.8, -0.8, m)]; > > Y = zeros(size(H)); Y(1) = -1.2756; % start on the lower branch > > for k = 2:numel(H) % Euler–Maruyama, h moving slowly > > Y(k) = Y(k-1) + dt*f(Y(k-1), H(k-1)) + sig*sqrt(dt)*randn; > > end > > plot(H(1:m), Y(1:m), 'r', H(m+1:end), Y(m+1:end), 'm'); xlabel('h'); ylabel('y') > > ``` > > > [!info]- R > > ```r > > a <- 1; f <- function(y, h) h + a * y - y^3; fy <- function(y) a - 3 * y^2 > > y <- seq(-1.6, 1.6, length.out = 1201); h <- y^3 - a * y # every equilibrium, exactly > > plot(ifelse(fy(y) < 0, h, NA), y, type = "l", xlim = c(-0.8, 0.8), xlab = "h", ylab = "y") > > lines(ifelse(fy(y) > 0, h, NA), y, lty = 2) > > set.seed(1); dt <- 0.01; LEG <- 400; sig <- 0.01; m <- LEG / dt > > H <- c(seq(-0.8, 0.8, length.out = m), seq(0.8, -0.8, length.out = m)) > > Y <- numeric(length(H)); Y[1] <- -1.2756 # start on the lower branch > > for (k in 2:length(H)) # Euler–Maruyama, h moving slowly > > Y[k] <- Y[k - 1] + dt * f(Y[k - 1], H[k - 1]) + sig * sqrt(dt) * rnorm(1) > > lines(H[1:m], Y[1:m], col = "red"); lines(H[-(1:m)], Y[-(1:m)], col = "purple") > > ``` > > > [!info]- Mathematica > > ```mathematica > > a = 1; f[y_, h_] := h + a y - y^3; > > diagram = ParametricPlot[{y^3 - a y, y}, {y, -1.6, 1.6}, > > ColorFunction -> Function[{h, y, u}, If[a - 3 y^2 < 0, Black, LightGray]], ColorFunctionScaling -> False]; > > dt = 0.01; m = 40000; sig = 0.01; SeedRandom[1]; > > H = Join[Subdivide[-0.8, 0.8, m - 1], Subdivide[0.8, -0.8, m - 1]]; > > Y = FoldList[#1 + dt f[#1, #2] + sig Sqrt[dt] RandomVariate[NormalDistribution[]] &, -1.2756, Most[H]]; > > Show[diagram, ListLinePlot[{Transpose[{H[[;; m]], Y[[;; m]]}], Transpose[{H[[m + 1 ;;]], Y[[m + 1 ;;]]}]}, > > PlotStyle -> {Red, Purple}], PlotRange -> {{-0.8, 0.8}, {-1.6, 1.6}}, AxesLabel -> {"h", "y"}] > > ``` > > [!example]- `bead-hoop` — the bead on a rotating hoop: equilibria against $\gamma$, and the bead released near the bottom at a few spin rates > > > [!info]- Python > > ```python > > import numpy as np, matplotlib.pyplot as plt > > F = lambda phi, g: np.sin(phi) * (g * np.cos(phi) - 1) # overdamped: phi' = F(phi) > > gam = np.linspace(1, 3, 200); plt.plot(gam, np.degrees(np.arccos(1 / gam)), "k", gam, -np.degrees(np.arccos(1 / gam)), "k") > > plt.plot([0, 1], [0, 0], "k", [1, 3], [0, 0], "k--") > > for g in [0.6, 1.0, 1.5, 2.5]: # release at 5 degrees, run to t = 60 > > phi, dt = np.radians(5), 0.01 > > for _ in range(6000): phi += dt * F(phi, g) > > plt.plot(g, np.degrees(phi), "ro") > > print(g, np.degrees(phi), np.degrees(np.arccos(1 / g)) if g > 1 else 0.0) > > plt.xlabel("gamma = r w^2 / g"); plt.ylabel("phi (degrees)"); plt.show() > > ``` > > > [!info]- MATLAB / Octave > > ```matlab > > F = @(phi, g) sin(phi) .* (g*cos(phi) - 1); % overdamped: phi' = F(phi) > > gam = linspace(1, 3, 200); plot(gam, acosd(1./gam), 'k', gam, -acosd(1./gam), 'k', [0 1], [0 0], 'k', [1 3], [0 0], 'k--'); hold on > > for g = [0.6 1.0 1.5 2.5] % release at 5 degrees, run to t = 60 > > [~, P] = ode45(@(t, p) F(p, g), [0 60], deg2rad(5)); > > plot(g, rad2deg(P(end)), 'ro') > > end > > xlabel('\gamma = r\omega^2/g'); ylabel('\phi (degrees)') > > ``` > > > [!info]- R > > ```r > > F <- function(phi, g) sin(phi) * (g * cos(phi) - 1) # overdamped: phi' = F(phi) > > gam <- seq(1, 3, length.out = 200) > > plot(gam, acos(1 / gam) * 180 / pi, type = "l", xlim = c(0, 3), ylim = c(-80, 80), xlab = "gamma", ylab = "phi (degrees)") > > lines(gam, -acos(1 / gam) * 180 / pi); lines(c(0, 1), c(0, 0)); lines(c(1, 3), c(0, 0), lty = 2) > > for (g in c(0.6, 1.0, 1.5, 2.5)) { # release at 5 degrees, run to t = 60 > > phi <- 5 * pi / 180; for (k in 1:6000) phi <- phi + 0.01 * F(phi, g) > > points(g, phi * 180 / pi, col = "red", pch = 19) > > } > > ``` > > > [!info]- Mathematica > > ```mathematica > > F[phi_, g_] := Sin[phi] (g Cos[phi] - 1); (* overdamped: phi' = F[phi] *) > > final[g_] := phi[60] /. First @ NDSolve[{phi'[t] == F[phi[t], g], phi[0] == 5 Degree}, phi, {t, 0, 60}]; > > Show[Plot[{ArcCos[1/g]/Degree, -ArcCos[1/g]/Degree}, {g, 1, 3}, PlotStyle -> Black], > > ListPlot[Table[{g, final[g]/Degree}, {g, {0.6, 1.0, 1.5, 2.5}}], PlotStyle -> Red], > > PlotRange -> {{0, 3}, {-80, 80}}, AxesLabel -> {"\[Gamma]", "\[Phi] (deg)"}] > > ``` > > [!quote] Sources > The sequence on this page (blue sky, transcritical, supercritical pitchfork, then hysteresis through $rx + x^3 - x^5$) and the bead on the hoop follow S. H. Strogatz, *Nonlinear Dynamics and Chaos*, Chapter 3 (§3.1–3.6), where the imperfect pitchfork and the cusp are treated as "imperfect bifurcations and catastrophes." The spinning-hoop video is from the AppDynSys channel.