---
title: "5. Hysteresis (Double Saddle-Node)"
sidebar: usage
jupyter: python3
---

A small constant tilt added to the pitchfork unfolds it into a hysteresis
curve with **two** saddle-node bifurcations:

$$
\dot x = a + x - x^3.
$$

The bifurcation diagram in $(a, x)$ is the classic S-curve with two folds
(vertical tangents) at $a_{\pm} = \pm\,2/(3\sqrt 3) \approx \pm 0.385$,
framing a *bistable* region in which the system has two coexisting stable
states separated by an unstable middle branch. Reproduces `Hysteresis.png`.

```{python}
from tvbo import Dynamics, Continuation, SimulationExperiment
import matplotlib.pyplot as plt
import numpy as np

HY_YAML = """
name: Hysteresis
description: Double saddle-node / hysteresis $\\dot x = a + x - x^3$.
parameters:
  a:
    name: a
    value: 0.0
state_variables:
  x:
    name: x
    domain: { lo: -2.0, hi: 2.0 }
    equation:
      lhs: Derivative(x, t)
      rhs: a + x - x**3
    initial_value: 1.0
"""

CONT_YAML = """
name: hy_in_a
dynamics: Hysteresis
free_parameters:
  - name: a
    domain: { lo: -2.0, hi: 2.0 }
max_steps: 600
ds: 0.005
bothside: true
"""

dyn = Dynamics.from_string(HY_YAML)
cont = Continuation.from_string(CONT_YAML)
exp = SimulationExperiment(dynamics=dyn, continuations=[cont])
```

## BifurcationKit.jl

```{python}
result_jl = exp.run("bifurcationkit.jl")
result_jl.continuations["hy_in_a"].plot(VOI="x")
plt.gca().set_xlim(-1, 1); plt.gca().set_ylim(-2, 2); plt.show()
```

## PyRates / PyCoBi

PyCoBi requires `NDIM >= 2`; the adapter raises a clear error for this 1-D
system. Use `auto-07p` or `bifurcationkit.jl` instead.

```{python}
try:
    result_py = exp.run("pyrates-bifurcation")
    result_py.continuations["hy_in_a"].plot(VOI="x")
except NotImplementedError as e:
    print(e)
```

## AUTO-07p (NumCont)

Direct analogue of the standalone `DoubleSaddleNode.py` example. AUTO-07p
time-integrates from the model defaults to land on the upper stable branch,
then sweeps $a$ in both directions so both saddle-nodes are picked up in a
single merged diagram.

```{python}
result_auto = exp.run("auto-07p")
result_auto.continuations["hy_in_a"].plot(VOI="x")
plt.gca().set_xlim(-1, 1); plt.gca().set_ylim(-2, 2); plt.show()
```

## Forward and backward sweep — visualising the hysteresis loop

Numerically integrate the system while *slowly* sweeping $a$ first up
($-1 \to 1$) then back down ($1 \to -1$). Trajectories *track* the upper
branch on the way up and the lower branch on the way down — the system
*remembers* its history.

```{python}
def sweep(dyn, a_start, a_end, n_steps=4000, x0=1.0):
    rhs = lambda x, a: a + x - x**3
    a_traj = np.linspace(a_start, a_end, n_steps)
    x = np.empty(n_steps); x[0] = x0
    dt = 0.05
    for k in range(1, n_steps):
        x[k] = x[k-1] + dt * rhs(x[k-1], a_traj[k-1])
    return a_traj, x

a_up,   x_up   = sweep(dyn, -1.0, +1.0, x0=-1.0)
a_down, x_down = sweep(dyn, +1.0, -1.0, x0=+1.0)

fig, ax = plt.subplots(figsize=(6, 4))
ax.plot(a_up, x_up, color="C0", lw=2, label=r"sweep $a:-1\to 1$ (start low)")
ax.plot(a_down, x_down, color="C3", lw=2, label=r"sweep $a:+1\to -1$ (start high)")
af = 2.0/(3.0*np.sqrt(3.0))
ax.axvline(+af, color="0.5", ls="--", lw=0.8)
ax.axvline(-af, color="0.5", ls="--", lw=0.8)
ax.set_xlabel("$a$"); ax.set_ylabel("$x$")
ax.set_xlim(-1, 1); ax.set_ylim(-2, 2)
ax.legend(loc="upper left", fontsize=9, handlelength=1.4)
plt.tight_layout(); plt.show()
```

The system *jumps* discontinuously to the opposite branch precisely when
each saddle-node is crossed. The enclosed area of the loop is the
*hysteresis area* — a quantitative measure of memory.
