CDIMA: Stability and the Hopf Point

Module 1: Low-Dimensional Dynamics

This page continues CDIMA: Nullclines and Fixed Points. You decide the stability of the fixed point from the Jacobian and check the Hopf condition \(b_c = 3a/5 - 25/a\) from Hopf Bifurcation (Strogatz 2024, sec. 8.3).

Prediction

Exercise

For \(a = 10\), which of \(b = 3\) and \(b = 4\) gives a stable fixed point? What should the eigenvalues look like at \(b = 3.5\)?

Solution

\(b_c = 3 \cdot 10/5 - 25/10 = 3.5\). The fixed point is stable for \(b > b_c\), so \(b = 4\) gives a stable spiral and \(b = 3\) an unstable one. At \(b = b_c\) the eigenvalues are a purely imaginary pair: the real part crosses zero while the imaginary part stays nonzero.

Implementation: the Jacobian

The Jacobian collects the partial derivatives of the right-hand side:

\[ J = \begin{bmatrix} \partial f/\partial x & \partial f/\partial y \\ \partial g/\partial x & \partial g/\partial y \end{bmatrix}. \]

As in one dimension, a stable equilibrium attracts nearby trajectories and an unstable one repels them. In two dimensions the eigenvalues of \(J\) at the fixed point decide which. You can compute \(J\) analytically or numerically. The numerical route uses central differences and works for any model.

Exercise

Complete the central difference Jacobian. Column \(j\) is the derivative of the right-hand side with respect to variable \(j\).

def numerical_jacobian(func, point, args, h=1e-6):
    jac = np.zeros((2, 2))
    for j in range(2):
        e = np.zeros(2)
        e[j] = h
        jac[:, j] = # TODO: (f(point + e) - f(point - e)) / (2h)
    return jac

Solution

def numerical_jacobian(func, point, args, h=1e-6):
    jac = np.zeros((2, 2))
    for j in range(2):
        e = np.zeros(2)
        e[j] = h
        jac[:, j] = (np.asarray(func(0, point + e, *args))
                     - np.asarray(func(0, point - e, *args))) / (2 * h)
    return jac

amlab.odes_2d.cdima_analysis.numerical_jacobian is the reference version for CDIMA.

Then compute the eigenvalues with np.linalg.eigvals:

• both real parts negative: stable (node or spiral);
• at least one real part positive: unstable.

Validation

Compare the numerical Jacobian with the analytical one and check where the real part of the eigenvalues changes sign.

import numpy as np

from amlab.odes_2d.cdima_analysis import (
    cdima_jacobian,
    hopf_b,
    locate_hopf_numerically,
    numerical_jacobian,
)

a = 10.0
for b in [3.0, 3.5, 4.0]:
    ev = np.linalg.eigvals(numerical_jacobian(a, b))
    diff = np.max(np.abs(numerical_jacobian(a, b) - cdima_jacobian(a, b)))
    print(f"b = {b}: eigenvalues {np.round(ev, 4)}, |J_num - J_exact| = {diff:.1e}")

for a in [8.0, 10.0, 15.0]:
    print(f"a = {a}: numerical Hopf b = {locate_hopf_numerically(a):.4f}, formula 3a/5 - 25/a = {hopf_b(a):.4f}")

b = 3.0: eigenvalues [0.1+2.4474j 0.1-2.4474j], |J_num - J_exact| = 4.9e-10
b = 3.5: eigenvalues [-0.+2.6458j -0.-2.6458j], |J_num - J_exact| = 4.9e-10
b = 4.0: eigenvalues [-0.1+2.8267j -0.1-2.8267j], |J_num - J_exact| = 4.9e-10
a = 8.0: numerical Hopf b = 1.6750, formula 3a/5 - 25/a = 1.6750
a = 10.0: numerical Hopf b = 3.5000, formula 3a/5 - 25/a = 3.5000
a = 15.0: numerical Hopf b = 7.3333, formula 3a/5 - 25/a = 7.3333

At \(b = 3.5\) the eigenvalues are \(\pm i\sqrt 7 \approx \pm 2.6458i\), as predicted, and the numerical Hopf point agrees with the formula to the accuracy of the root finder.

Experiment

Exercise

Run the phase plane from CDIMA: Phase Plane with \(a = 10\) and \(b = 4, 3.6, 3.4, 3\). Where does the limit cycle appear? How does its size change as \(b\) decreases below \(b_c\)?

Solution

For \(b = 4\) and \(b = 3.6\) the oscillations die out (slowly at \(b = 3.6\), where the real part of the eigenvalues is only \(-0.02\)). The limit cycle appears once \(b < b_c = 3.5\). Its amplitude grows as \(b\) decreases: about 0.43 at \(b = 3.4\) and 0.86 at \(b = 3\).

from amlab.odes_2d.cdima_analysis import limit_cycle_amplitude

for b in [4.0, 3.6, 3.4, 3.0]:
    print(f"b = {b}: amplitude of x after transients = {limit_cycle_amplitude(10.0, b, t_end=400.0):.3f}")

b = 4.0: amplitude of x after transients = 0.000
b = 3.6: amplitude of x after transients = 0.000
b = 3.4: amplitude of x after transients = 0.430
b = 3.0: amplitude of x after transients = 0.862

Interpretation

For \(b > b_c\) the fixed point is a stable spiral and the reaction settles to constant concentrations. For \(b < b_c\) it is an unstable spiral surrounded by a stable limit cycle: the concentrations oscillate. The formula gives the whole stability boundary in the \((a, b)\) plane. Note that \(b_c > 0\) requires \(a > \sqrt{125/3} \approx 6.45\); for smaller \(a\) the fixed point is stable for every \(b > 0\).

Limitations

The linear analysis says nothing about the size of the limit cycle, nor whether the Hopf bifurcation is supercritical. You answered those questions with simulations in Hopf Bifurcation.

Resources

• Strogatz (2024), Section 8.3, for the theory of the CDIMA model and its phase plane.
• SciPy documentation for solve_ivp.
• Reference implementation: amlab/odes_2d.

Next

You now have the full phase-plane workflow: trajectories, nullclines, fixed points and their stability. Next you animate the trajectories, then reproduce the same plots for the Van der Pol oscillator and the FitzHugh-Nagumo model.

References

Strogatz, Steven H. 2024. Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. Chapman; Hall/CRC.