 codex/past-exam-of-the-mathematics-course-of-the-university-of-cambridge/2011/ia/paper-2-planar-phase.py

"""Python 3.14; numpy 2.3.5 and matplotlib 3.10.7. Emit an opaque PNG to cwd.
The caller supplies MPLCONFIGDIR; no cache location is changed here.
"""
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

def rhs(z):
    x, y = z
    return np.array([x+y, -x+y-2*x*x*y])

def orbit(seed, direction):
    h = direction * .001
    z = np.array(seed, dtype=float)
    pts = [z.copy()]
    for _ in range(16000):
        k1 = rhs(z)
        k2 = rhs(z+h*k1/2)
        k3 = rhs(z+h*k2/2)
        k4 = rhs(z+h*k3)
        z += h*(k1+2*k2+2*k3+k4)/6
        if np.max(np.abs(z)) > 2.3:
            break
        pts.append(z.copy())
    return np.asarray(pts)

fig, ax = plt.subplots(figsize=(8, 6.2), dpi=100, facecolor="white")
grid = np.linspace(-2.12, 2.12, 160)
X, Y = np.meshgrid(grid, grid)
U, V = X+Y, -X+Y-2*X*X*Y
ax.streamplot(grid, grid, U, V, color="#a7adb4", density=1.15,
              linewidth=.8, arrowsize=.9)
ax.plot(grid, -grid, "k--", alpha=.35, lw=1, label="nullclines")
for lo, hi in [(-2.12, -.73), (-.685, .685), (.73, 2.12)]:
    x = np.linspace(lo, hi, 400)
    ax.plot(x, x/(1-2*x*x), "k--", alpha=.35, lw=1)
for p in [np.array([1., -1.]), np.array([-1., 1.])]:
    for vector, direction, color, label in [
            (np.array([1., -3.])/np.sqrt(10), -1, "#2163b5", "stable"),
            (np.array([1., 1.])/np.sqrt(2), 1, "#bb3e36", "unstable")]:
        for sign in [-1, 1]:
            xy = orbit(p + sign*1e-5*vector, direction)
            ax.plot(xy[:, 0], xy[:, 1], color=color, lw=1.9,
                    label=label if np.all(p == [1, -1]) and sign == 1 else None)
            visible = np.where((np.linalg.norm(xy-p, axis=1) > .12) &
                               (np.linalg.norm(xy-p, axis=1) < .28))[0]
            if len(visible) > 30:
                j = visible[len(visible)//2]
                a, b = xy[j-15], xy[j+15]
                if direction < 0:
                    a, b = b, a
                ax.annotate("", xy=b, xytext=a,
                            arrowprops=dict(arrowstyle="->", color=color, lw=1.8))
    ax.scatter(*p, marker="x", color="black", s=65, zorder=5)
    ax.annotate("saddle", p, xytext=(9, 9), textcoords="offset points")
ax.scatter(0, 0, color="black", s=28, zorder=5)
ax.annotate("unstable spiral", (0, 0), xytext=(12, 12),
            textcoords="offset points", bbox=dict(facecolor="white", alpha=.8, edgecolor="none"))
ax.set(xlim=(-2.12, 2.12), ylim=(-2.12, 2.12), xlabel="x", ylabel="y",
       title="Planar flow and saddle separatrices")
ax.set_aspect("equal")
ax.legend(loc="lower left", fontsize=9, framealpha=.95)
fig.tight_layout()
fig.savefig(Path.cwd() / "paper-2-planar-phase.png", dpi=100,
            facecolor="white", transparent=False)
plt.close(fig)