In previous lessons, we introduced nonlinear sources, phase portraits, and local linearization. Here we go one step further: some nonlinear systems exhibit long-term behavior that is neither convergence to an equilibrium nor divergence to infinity. Instead, trajectories can converge to a closed orbit. Such an isolated closed orbit is called a limit cycle.

We also examine how changing a scalar parameter \( \mu \) or \( r \) can change the number or stability of equilibria. This qualitative change is a bifurcation. We focus on the simplest (normal-form level) bifurcations that are essential in control engineering models with saturation, friction, positive feedback, or self-excitation.

flowchart TD
  A["Nonlinear model xdot = f(x, p)"] --> B["Find equilibria f(x*, p)=0"]
  B --> C["Local test: Jacobian eigenvalues"]
  C --> D["Check global clues in phase plane"]
  D --> E["Closed orbit seen? → \nlimit cycle candidate"]
  D --> F["Equilibria count changes with p? → \nbifurcation candidate"]
  E --> G["Numerical verification (time + phase trajectories)"]
  F --> G
  G --> H["Interpret stability and engineering meaning"]
        

Consider an autonomous nonlinear system in the plane:

\[ \dot{\mathbf{x} } = \mathbf{f}(\mathbf{x}), \quad \mathbf{x} = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}. \]

A trajectory \( \mathbf{x}(t) \) is periodic with period \( T \gt 0 \) if

\[ \mathbf{x}(t+T) = \mathbf{x}(t) \quad \text{for all } t. \]

A periodic orbit is a limit cycle if it is isolated: there is a neighborhood of the closed orbit that contains no other periodic orbits. Stability of a limit cycle is defined using distance to the orbit (not distance to a point).

For a stable limit cycle \( \Gamma \), trajectories starting near \( \Gamma \) satisfy

\[ \operatorname{dist}(\mathbf{x}(t), \Gamma) \to 0 \quad \text{as } t \to \infty. \]

In control engineering, self-oscillation due to nonlinear actuation, relay action, saturation, or negative damping often appears as a limit cycle. Unlike linear undamped oscillators (which have a continuum of closed orbits), a nonlinear limit cycle is isolated and therefore robust to initial-condition perturbations.

A scalar nonlinear system

\[ \dot{x} = f(x; r) \]

can have several equilibria for the same parameter \( r \). Each equilibrium \( x^\star \) satisfies \( f(x^\star; r)=0 \). The local stability test follows from linearization:

\[ \dot{\eta} = f_x(x^\star; r)\,\eta + \mathcal{O}(\eta^2), \quad \eta = x - x^\star. \]

Therefore, if \( f_x(x^\star; r) \lt 0 \), the equilibrium is locally asymptotically stable; if \( f_x(x^\star; r) \gt 0 \), it is unstable. The case \( f_x(x^\star; r)=0 \) is nonhyperbolic and requires higher-order analysis (precisely the setting where bifurcations occur).

Proposition 1 (1D local stability criterion).

Let \( f \in C^1 \) and let \( x^\star \) be an equilibrium. If \( a = f_x(x^\star; r) \neq 0 \), then the equilibrium is locally exponentially stable for \( a \lt 0 \) and unstable for \( a \gt 0 \).

Proof. Write \( x = x^\star + \eta \). By Taylor expansion,

\[ \dot{\eta} = a\eta + \phi(\eta), \quad \lim_{\eta \to 0}\frac{\phi(\eta)}{\eta}=0. \]

For sufficiently small \( |\eta| \), we can bound \( |\phi(\eta)| \le \varepsilon |\eta| \) with \( 0 \lt \varepsilon \lt |a|/2 \). Then

\[ (a-\varepsilon)|\eta| \le \frac{d}{dt}|\eta| \le (a+\varepsilon)|\eta| \]

in the comparison sense. If \( a \lt 0 \), the upper bound decays exponentially; if \( a \gt 0 \), the lower bound grows, proving instability. \(\square\)

The simplest bifurcation analysis is done on normal forms. These are reduced equations that capture the local qualitative behavior near a nonhyperbolic equilibrium.

4.1 Saddle-Node (Fold) Normal Form

\[ \dot{x} = r - x^2. \]

Equilibria satisfy \( x^\star = \pm \sqrt{r} \) for \( r \gt 0 \), one repeated nonhyperbolic equilibrium \( x^\star=0 \) at \( r=0 \), and no real equilibrium for \( r \lt 0 \).

Since \( f_x = -2x \), the branch \( x^\star = +\sqrt{r} \) is stable and \( x^\star = -\sqrt{r} \) is unstable.

4.2 Supercritical Pitchfork Normal Form

\[ \dot{x} = rx - x^3 = x(r - x^2). \]

Equilibria are \( x^\star=0 \) for all \( r \), and additional branches \( x^\star = \pm\sqrt{r} \) for \( r \gt 0 \). The derivative \( f_x = r - 3x^2 \) gives:

\[ f_x(0;r)=r,\qquad f_x(\pm \sqrt{r};r)= -2r \quad (r \gt 0). \]

Hence the origin is stable for \( r \lt 0 \) and unstable for \( r \gt 0 \), while the two nonzero branches are stable for \( r \gt 0 \). This is the canonical symmetry-breaking scenario.

Proposition 2 (classification by derivative signs).

For the two normal forms above, the stability claims follow directly from Proposition 1 because every equilibrium branch away from \( r=0 \) is hyperbolic.

Proof. Substitute each equilibrium into \( f_x \) and inspect the sign. At \( r=0 \), the derivative vanishes and hyperbolic linearization is inconclusive, which is exactly why \( r=0 \) is the bifurcation point. \(\square\)

In two-dimensional autonomous systems, limit cycles can be studied using geometric tools. A complete proof of Poincaré–Bendixson is beyond this lesson, but its statement is fundamental:

Poincaré–Bendixson (informal statement). If a trajectory remains in a closed, bounded region of the plane and the region contains no equilibrium point in its \(\omega\)-limit set, then the \(\omega\)-limit set is a periodic orbit.

This theorem is a key reason why phase-plane analysis is powerful in second-order engineering models. It does not extend directly to dimensions 3 and higher.

A complementary tool is the Bendixson–Dulac criterion, which can rule out periodic orbits. If there exists a continuously differentiable scalar function \( g(x_1,x_2) \) such that

\[ \frac{\partial (g f_1)}{\partial x_1} + \frac{\partial (g f_2)}{\partial x_2} \]

has one strict sign (never zero) in a simply connected region, then no periodic orbit lies entirely in that region.

For engineering interpretation: one can often prove that oscillations are impossible in a region, or narrow the search for self-sustained oscillations to a specific annulus in the phase plane.

The Van der Pol oscillator is a standard self-excited nonlinear system:

\[ \ddot{x} - \mu(1-x^2)\dot{x} + x = 0, \qquad \mu \gt 0. \]

In first-order form with \( x_1 = x \), \( x_2 = \dot{x} \):

\[ \dot{x}_1 = x_2, \qquad \dot{x}_2 = \mu(1-x_1^2)x_2 - x_1. \]

The origin is the only equilibrium. Linearization at the origin gives

\[ \mathbf{A} = \begin{bmatrix} 0 & 1 \\ -1 & \mu \end{bmatrix}, \qquad \lambda^2 - \mu\lambda + 1 = 0. \]

For \( \mu \gt 0 \), the trace is positive, so the origin is unstable. However, trajectories do not diverge to infinity. The nonlinear damping term \( -\mu(1-x^2)\dot{x} \) acts like negative damping for small \( |x| \) and positive damping for large \( |x| \), creating an attracting limit cycle.

A useful energy-like quantity is \( V(x,\dot{x})=\frac{1}{2}(x^2+\dot{x}^2) \). Along trajectories:

\[ \dot{V} = x\dot{x} + \dot{x}\ddot{x} = x\dot{x} + \dot{x}\big(\mu(1-x^2)\dot{x} - x\big) = \mu(1-x^2)\dot{x}^2. \]

Thus \( \dot{V} \gt 0 \) when \( |x| \lt 1 \) (energy injection) and \( \dot{V} \lt 0 \) when \( |x| \gt 1 \) (energy dissipation), which strongly suggests a trapping region and a stable oscillation amplitude.

A Hopf bifurcation is the planar mechanism by which an equilibrium changes stability and a small-amplitude periodic orbit appears (or disappears). At the normal-form level, a supercritical Hopf bifurcation can be written in polar coordinates as

\[ \dot{\rho} = \mu\rho - \rho^3, \qquad \dot{\theta} = \omega, \qquad \omega \neq 0. \]

The radial dynamics is scalar, so the equilibrium radii are \( \rho^\star = 0 \) and \( \rho^\star = \sqrt{\mu} \) for \( \mu \gt 0 \).

Proposition 3 (supercritical Hopf normal-form stability).

For the radial equation above, the origin is stable for \( \mu \lt 0 \) and unstable for \( \mu \gt 0 \). When \( \mu \gt 0 \), the circle \( \rho = \sqrt{\mu} \) is asymptotically stable, corresponding to a stable limit cycle in Cartesian coordinates.

Proof. The radial equation factors as

\[ \dot{\rho} = \rho(\mu - \rho^2). \]

For \( \mu \lt 0 \), all sufficiently small positive radii satisfy \( \dot{\rho} \lt 0 \), so \( \rho \to 0 \). For \( \mu \gt 0 \), we have \( \dot{\rho} \gt 0 \) when \( 0 \lt \rho \lt \sqrt{\mu} \) and \( \dot{\rho} \lt 0 \) when \( \rho \gt \sqrt{\mu} \), hence \( \rho(t) \to \sqrt{\mu} \). Since \( \dot{\theta}=\omega \neq 0 \), the angular motion is periodic, yielding a stable closed orbit. \(\square\)

flowchart TD
  A["Parameter mu < 0"] --> B["Equilibrium stable"]
  C["Parameter mu = 0"] --> D["Pair of imag-axis eigenvalues (critical)"]
  E["Parameter mu > 0"] --> F["Equilibrium unstable + small stable cycle"]
  B --> D
  D --> F
        

For an engineering nonlinear model \( \dot{\mathbf{x} } = \mathbf{f}(\mathbf{x};p) \) with scalar parameter \( p \), a practical analysis sequence is:

1. Compute equilibria by solving \( \mathbf{f}(\mathbf{x}^\star;p)=0 \).
2. Compute the Jacobian \( \mathbf{J}(\mathbf{x}^\star;p) \) and inspect eigenvalues.
3. Track sign changes in determinant/trace (2D) or eigenvalue crossings as \( p \) varies.
4. Use phase portraits and time simulations to detect periodic behavior.
5. Use Poincaré section crossings to estimate period numerically.

In 2D, the local linearized characteristic polynomial is

\[ \lambda^2 - \operatorname{tr}(\mathbf{J})\lambda + \det(\mathbf{J}) = 0. \]

A Hopf-like transition requires (locally) a pair of complex conjugates crossing the imaginary axis, which is associated with

\[ \operatorname{tr}(\mathbf{J}) = 0, \qquad \det(\mathbf{J}) \gt 0 \]

at the critical parameter, plus nondegeneracy conditions (studied formally in advanced nonlinear dynamics courses).

The Python script simulates the Van der Pol oscillator using RK4, estimates the asymptotic period with a Poincaré-section crossing method, and plots equilibrium branches for saddle-node and pitchfork normal forms.

Code: Chapter14_Lesson5.py

# Chapter14_Lesson5.py
# Limit Cycles, Multiple Equilibria, and Basic Bifurcation Notions
# Python implementation: Van der Pol oscillator + 1D bifurcation normal forms
import numpy as np
import matplotlib.pyplot as plt

def rk4_step(f, t, x, h, params):
    k1 = f(t, x, params)
    k2 = f(t + 0.5*h, x + 0.5*h*k1, params)
    k3 = f(t + 0.5*h, x + 0.5*h*k2, params)
    k4 = f(t + h, x + h*k3, params)
    return x + (h/6.0)*(k1 + 2*k2 + 2*k3 + k4)

def vdp_rhs(t, state, params):
    mu = params["mu"]
    x, y = state
    dx = y
    dy = mu*(1.0 - x*x)*y - x
    return np.array([dx, dy], dtype=float)

def simulate_vdp(mu=1.0, x0=(2.0, 0.0), T=80.0, h=0.01):
    n = int(T/h)
    t = np.linspace(0.0, T, n+1)
    x = np.zeros((n+1, 2))
    x[0] = np.array(x0, dtype=float)
    params = {"mu": mu}
    for k in range(n):
        x[k+1] = rk4_step(vdp_rhs, t[k], x[k], h, params)
    return t, x

def estimate_period_from_crossings(t, traj):
    # Upward crossings of y=0 (x-axis): use y sign change from negative to positive
    y = traj[:, 1]
    x = traj[:, 0]
    crossings = []
    for k in range(len(t)-1):
        if y[k] < 0.0 and y[k+1] >= 0.0:
            # Linear interpolation
            alpha = -y[k] / (y[k+1] - y[k] + 1e-14)
            tc = t[k] + alpha*(t[k+1] - t[k])
            xc = x[k] + alpha*(x[k+1] - x[k])
            if xc > 0:  # consistent Poincare section branch
                crossings.append(tc)
    if len(crossings) < 3:
        return np.nan
    # Discard transient crossing intervals
    periods = np.diff(crossings)
    if len(periods) >= 4:
        periods = periods[-4:]
    return float(np.mean(periods))

def classify_pitchfork(r, x):
    # f(x)=r x - x^3 ; stability from f'(x)=r - 3 x^2
    lam = r - 3.0*x*x
    if lam < 0:
        return "stable"
    elif lam > 0:
        return "unstable"
    return "nonhyperbolic"

def classify_saddlenode(r, x):
    # f(x)=r - x^2 ; stability from f'(x)=-2x
    lam = -2.0*x
    if lam < 0:
        return "stable"
    elif lam > 0:
        return "unstable"
    return "nonhyperbolic"

def main():
    # 1) Van der Pol simulation for several mu values
    mu_list = [0.2, 1.0, 3.0]
    results = []
    fig1 = plt.figure(figsize=(7, 5))
    for mu in mu_list:
        t, z = simulate_vdp(mu=mu, x0=(2.0, 0.1), T=80.0, h=0.01)
        # discard transient for plotting
        z_ss = z[int(0.5*len(z)):, :]
        plt.plot(z_ss[:,0], z_ss[:,1], label=f"mu={mu}")
        period = estimate_period_from_crossings(t[int(0.4*len(t)):], z[int(0.4*len(z)):, :])
        amp = float(np.max(np.abs(z_ss[:,0])))
        results.append((mu, amp, period))
    plt.xlabel("x")
    plt.ylabel("x_dot")
    plt.title("Van der Pol Limit Cycles (phase plane)")
    plt.legend()
    plt.grid(True)
    plt.tight_layout()

    # 2) 1D bifurcation normal forms (equilibrium branches)
    r_vals = np.linspace(-2.0, 2.0, 401)
    pf_stable_r, pf_stable_x = [], []
    pf_unstable_r, pf_unstable_x = [], []
    sn_stable_r, sn_stable_x = [], []
    sn_unstable_r, sn_unstable_x = [], []

    # Pitchfork: x*=0 and ±sqrt(r) for r>=0
    for r in r_vals:
        for xeq in [0.0]:
            c = classify_pitchfork(r, xeq)
            (pf_stable_r if c == "stable" else pf_unstable_r).append(r)
            (pf_stable_x if c == "stable" else pf_unstable_x).append(xeq)
        if r >= 0:
            for xeq in [np.sqrt(r), -np.sqrt(r)]:
                c = classify_pitchfork(r, xeq)
                (pf_stable_r if c == "stable" else pf_unstable_r).append(r)
                (pf_stable_x if c == "stable" else pf_unstable_x).append(xeq)

    # Saddle-node: x*=±sqrt(r) for r>=0
    for r in r_vals:
        if r >= 0:
            for xeq in [np.sqrt(r), -np.sqrt(r)]:
                c = classify_saddlenode(r, xeq)
                (sn_stable_r if c == "stable" else sn_unstable_r).append(r)
                (sn_stable_x if c == "stable" else sn_unstable_x).append(xeq)

    fig2 = plt.figure(figsize=(7, 5))
    plt.plot(pf_stable_r, pf_stable_x, linewidth=2, label="Pitchfork stable")
    plt.plot(pf_unstable_r, pf_unstable_x, "--", linewidth=2, label="Pitchfork unstable")
    plt.plot(sn_stable_r, sn_stable_x, linewidth=2, label="Saddle-node stable")
    plt.plot(sn_unstable_r, sn_unstable_x, "--", linewidth=2, label="Saddle-node unstable")
    plt.xlabel("parameter r")
    plt.ylabel("equilibrium x*")
    plt.title("Basic Bifurcation Branches")
    plt.legend()
    plt.grid(True)
    plt.tight_layout()

    # 3) Print numerical summary
    print("Van der Pol asymptotic limit-cycle estimates")
    print("mu\tamplitude(|x|)\tperiod")
    for mu, amp, period in results:
        print(f"{mu:.2f}\t{amp:.4f}\t\t{period:.4f}")

    plt.show()

if __name__ == "__main__":
    main()

These implementations reproduce the RK4 simulation and branch classification in strongly typed compiled languages. They are useful for embedded simulation prototypes, high-performance batches, or integration into control software stacks.

Code: Chapter14_Lesson5.cpp

// Chapter14_Lesson5.cpp
// Limit Cycles, Multiple Equilibria, and Basic Bifurcation Notions
// C++17: RK4 simulation of Van der Pol oscillator and equilibrium classification
#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>
#include <fstream>
#include <string>

struct State {
    double x;
    double y;
};

State vdp_rhs(const State& s, double mu) {
    State ds;
    ds.x = s.y;
    ds.y = mu * (1.0 - s.x * s.x) * s.y - s.x;
    return ds;
}

State add(const State& a, const State& b, double scale=1.0) {
    return {a.x + scale * b.x, a.y + scale * b.y};
}

State rk4_step(const State& s, double h, double mu) {
    State k1 = vdp_rhs(s, mu);
    State k2 = vdp_rhs(add(s, k1, 0.5*h), mu);
    State k3 = vdp_rhs(add(s, k2, 0.5*h), mu);
    State k4 = vdp_rhs(add(s, k3, h), mu);
    State out;
    out.x = s.x + (h/6.0) * (k1.x + 2.0*k2.x + 2.0*k3.x + k4.x);
    out.y = s.y + (h/6.0) * (k1.y + 2.0*k2.y + 2.0*k3.y + k4.y);
    return out;
}

double estimate_period(const std::vector<double>& t, const std::vector<State>& z) {
    std::vector<double> crossings;
    for (size_t k = 0; k + 1 < z.size(); ++k) {
        if (z[k].y < 0.0 && z[k+1].y >= 0.0) {
            double denom = (z[k+1].y - z[k].y);
            if (std::abs(denom) < 1e-14) continue;
            double alpha = -z[k].y / denom;
            double tc = t[k] + alpha * (t[k+1] - t[k]);
            double xc = z[k].x + alpha * (z[k+1].x - z[k].x);
            if (xc > 0.0) crossings.push_back(tc);
        }
    }
    if (crossings.size() < 3) return std::nan("");
    size_t start = (crossings.size() > 5) ? crossings.size() - 5 : 1;
    double sum = 0.0;
    int cnt = 0;
    for (size_t i = start; i < crossings.size(); ++i) {
        sum += (crossings[i] - crossings[i-1]);
        cnt++;
    }
    return (cnt > 0) ? sum / cnt : std::nan("");
}

std::string stability_pitchfork(double r, double xeq) {
    double lambda = r - 3.0 * xeq * xeq;
    if (lambda < 0.0) return "stable";
    if (lambda > 0.0) return "unstable";
    return "nonhyperbolic";
}

std::string stability_saddlenode(double xeq) {
    double lambda = -2.0 * xeq;
    if (lambda < 0.0) return "stable";
    if (lambda > 0.0) return "unstable";
    return "nonhyperbolic";
}

int main() {
    const double mu = 1.0;
    const double h = 0.01;
    const double T = 80.0;
    const int N = static_cast<int>(T / h);

    std::vector<double> t(N + 1);
    std::vector<State> z(N + 1);
    z[0] = {2.0, 0.1};
    t[0] = 0.0;

    for (int k = 0; k < N; ++k) {
        t[k+1] = t[k] + h;
        z[k+1] = rk4_step(z[k], h, mu);
    }

    // Write phase-portrait data to CSV for plotting externally
    std::ofstream csv("Chapter14_Lesson5_vdp_phase.csv");
    csv << "t,x,y\n";
    for (int k = 0; k <= N; ++k) {
        csv << t[k] << "," << z[k].x << "," << z[k].y << "\n";
    }
    csv.close();

    // Estimate asymptotic period from latter part
    std::vector<double> t_tail(t.begin() + N/2, t.end());
    std::vector<State> z_tail(z.begin() + N/2, z.end());
    double period = estimate_period(t_tail, z_tail);

    double amp = 0.0;
    for (size_t k = N/2; k < z.size(); ++k) amp = std::max(amp, std::abs(z[k].x));

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "Van der Pol (mu=" << mu << ") limit-cycle estimate\n";
    std::cout << "Amplitude |x|_max ~= " << amp << "\n";
    std::cout << "Period ~= " << period << "\n\n";

    // Print equilibrium branches for two normal forms (sampled values)
    std::cout << "Pitchfork normal form: xdot = r*x - x^3\n";
    for (double r : {-1.0, 0.0, 1.0}) {
        std::cout << "r = " << r << ": ";
        std::cout << "x*=0 (" << stability_pitchfork(r, 0.0) << ")";
        if (r >= 0.0) {
            double x1 = std::sqrt(r), x2 = -std::sqrt(r);
            std::cout << ", x*=+" << x1 << " (" << stability_pitchfork(r, x1) << ")";
            std::cout << ", x*=" << x2 << " (" << stability_pitchfork(r, x2) << ")";
        }
        std::cout << "\n";
    }

    std::cout << "\nSaddle-node normal form: xdot = r - x^2\n";
    for (double r : {-1.0, 0.0, 1.0}) {
        std::cout << "r = " << r << ": ";
        if (r < 0.0) {
            std::cout << "no real equilibria\n";
        } else if (r == 0.0) {
            std::cout << "x*=0 (nonhyperbolic)\n";
        } else {
            double x1 = std::sqrt(r), x2 = -std::sqrt(r);
            std::cout << "x*=+" << x1 << " (" << stability_saddlenode(x1) << "), ";
            std::cout << "x*=" << x2 << " (" << stability_saddlenode(x2) << ")\n";
        }
    }

    std::cout << "\nCSV written: Chapter14_Lesson5_vdp_phase.csv\n";
    return 0;
}

Code: Chapter14_Lesson5.java

// Chapter14_Lesson5.java
// Limit Cycles, Multiple Equilibria, and Basic Bifurcation Notions
// Java implementation: RK4 simulation of Van der Pol + basic bifurcation classifications
import java.io.FileWriter;
import java.io.IOException;
import java.util.ArrayList;
import java.util.List;

public class Chapter14_Lesson5 {
    static class State {
        double x, y;
        State(double x, double y) { this.x = x; this.y = y; }
    }

    static State rhsVanDerPol(State s, double mu) {
        double dx = s.y;
        double dy = mu * (1.0 - s.x * s.x) * s.y - s.x;
        return new State(dx, dy);
    }

    static State add(State a, State b, double scale) {
        return new State(a.x + scale * b.x, a.y + scale * b.y);
    }

    static State rk4Step(State s, double h, double mu) {
        State k1 = rhsVanDerPol(s, mu);
        State k2 = rhsVanDerPol(add(s, k1, 0.5 * h), mu);
        State k3 = rhsVanDerPol(add(s, k2, 0.5 * h), mu);
        State k4 = rhsVanDerPol(add(s, k3, h), mu);
        return new State(
            s.x + (h / 6.0) * (k1.x + 2*k2.x + 2*k3.x + k4.x),
            s.y + (h / 6.0) * (k1.y + 2*k2.y + 2*k3.y + k4.y)
        );
    }

    static String pitchforkStability(double r, double xeq) {
        double lambda = r - 3.0 * xeq * xeq;
        if (lambda < 0) return "stable";
        if (lambda > 0) return "unstable";
        return "nonhyperbolic";
    }

    static String saddleNodeStability(double xeq) {
        double lambda = -2.0 * xeq;
        if (lambda < 0) return "stable";
        if (lambda > 0) return "unstable";
        return "nonhyperbolic";
    }

    static double estimatePeriod(List<Double> t, List<State> z) {
        List<Double> crossings = new ArrayList<>();
        for (int k = 0; k < z.size() - 1; k++) {
            if (z.get(k).y < 0.0 && z.get(k+1).y >= 0.0) {
                double denom = z.get(k+1).y - z.get(k).y;
                if (Math.abs(denom) < 1e-14) continue;
                double alpha = -z.get(k).y / denom;
                double tc = t.get(k) + alpha * (t.get(k+1) - t.get(k));
                double xc = z.get(k).x + alpha * (z.get(k+1).x - z.get(k).x);
                if (xc > 0) crossings.add(tc);
            }
        }
        if (crossings.size() < 3) return Double.NaN;
        int start = Math.max(1, crossings.size() - 5);
        double sum = 0.0;
        int count = 0;
        for (int i = start; i < crossings.size(); i++) {
            sum += crossings.get(i) - crossings.get(i - 1);
            count++;
        }
        return (count > 0) ? sum / count : Double.NaN;
    }

    public static void main(String[] args) throws IOException {
        double mu = 1.0;
        double h = 0.01;
        double T = 80.0;
        int N = (int)(T / h);

        List<Double> t = new ArrayList<>(N + 1);
        List<State> z = new ArrayList<>(N + 1);
        t.add(0.0);
        z.add(new State(2.0, 0.1));

        for (int k = 0; k < N; k++) {
            t.add(t.get(k) + h);
            z.add(rk4Step(z.get(k), h, mu));
        }

        double amp = 0.0;
        for (int k = N/2; k < z.size(); k++) {
            amp = Math.max(amp, Math.abs(z.get(k).x));
        }

        List<Double> tTail = t.subList(N/2, t.size());
        List<State> zTail = z.subList(N/2, z.size());
        double period = estimatePeriod(tTail, zTail);

        try (FileWriter fw = new FileWriter("Chapter14_Lesson5_vdp_phase.csv")) {
            fw.write("t,x,y\n");
            for (int k = 0; k < z.size(); k++) {
                fw.write(t.get(k) + "," + z.get(k).x + "," + z.get(k).y + "\n");
            }
        }

        System.out.println("Van der Pol limit-cycle summary");
        System.out.printf("mu = %.3f%n", mu);
        System.out.printf("Amplitude |x|_max ~= %.6f%n", amp);
        System.out.printf("Period ~= %.6f%n%n", period);

        double[] rVals = {-1.0, 0.0, 1.0};
        System.out.println("Pitchfork normal form: xdot = r*x - x^3");
        for (double r : rVals) {
            System.out.print("r = " + r + ": x*=0 (" + pitchforkStability(r, 0.0) + ")");
            if (r >= 0.0) {
                double x1 = Math.sqrt(r), x2 = -Math.sqrt(r);
                System.out.print(", x*=+" + x1 + " (" + pitchforkStability(r, x1) + ")");
                System.out.print(", x*=" + x2 + " (" + pitchforkStability(r, x2) + ")");
            }
            System.out.println();
        }

        System.out.println("\nSaddle-node normal form: xdot = r - x^2");
        for (double r : rVals) {
            System.out.print("r = " + r + ": ");
            if (r < 0.0) {
                System.out.println("no real equilibria");
            } else if (r == 0.0) {
                System.out.println("x*=0 (nonhyperbolic)");
            } else {
                double x1 = Math.sqrt(r), x2 = -Math.sqrt(r);
                System.out.println("x*=+" + x1 + " (" + saddleNodeStability(x1) + "), x*=" + x2 + " (" + saddleNodeStability(x2) + ")");
            }
        }

        System.out.println("\nCSV written: Chapter14_Lesson5_vdp_phase.csv");
    }
}

The MATLAB script includes both numerical simulation (ode45) and a programmatic Simulink model construction for the Van der Pol oscillator. The Mathematica notebook uses NDSolveValue and symbolic/numeric plotting for limit-cycle and normal-form visualization.

Code: Chapter14_Lesson5.m

% Chapter14_Lesson5.m
% Limit Cycles, Multiple Equilibria, and Basic Bifurcation Notions
% MATLAB + basic Simulink automation (if Simulink is installed)

clear; clc; close all;

%% Part A: Van der Pol oscillator (ode45) and limit-cycle period estimation
mu = 1.0;
f = @(t,z) [z(2); mu*(1 - z(1)^2)*z(2) - z(1)];
tspan = [0 80];
z0 = [2; 0.1];
opts = odeset('RelTol',1e-8,'AbsTol',1e-10);
[t,z] = ode45(f, tspan, z0, opts);

% phase portrait (transient removed)
idx0 = floor(numel(t)/2);
figure('Name','Van der Pol Phase Portrait');
plot(z(idx0:end,1), z(idx0:end,2), 'LineWidth', 1.5); grid on;
xlabel('x'); ylabel('x_dot');
title('Van der Pol Limit Cycle (mu = 1)');

% Poincare crossing estimate: y = 0 upward, x > 0
crossings = [];
for k = idx0:(numel(t)-1)
    if z(k,2) < 0 && z(k+1,2) >= 0
        alpha = -z(k,2)/(z(k+1,2)-z(k,2) + eps);
        tc = t(k) + alpha*(t(k+1)-t(k));
        xc = z(k,1) + alpha*(z(k+1,1)-z(k,1));
        if xc > 0
            crossings(end+1) = tc; %#ok<SAGROW>
        end
    end
end
if numel(crossings) >= 3
    P = mean(diff(crossings(max(1,end-4):end)));
else
    P = NaN;
end
A = max(abs(z(idx0:end,1)));
fprintf('Van der Pol limit-cycle estimate: amplitude = %.6f, period = %.6f\n', A, P);

%% Part B: Basic 1D bifurcation normal forms
r = linspace(-2, 2, 801);
pf_stable = []; pf_unstable = [];
sn_stable = []; sn_unstable = [];

for rv = r
    % Pitchfork: x*=0 always
    lam0 = rv;
    if lam0 < 0
        pf_stable = [pf_stable; rv 0]; %#ok<AGROW>
    else
        pf_unstable = [pf_unstable; rv 0]; %#ok<AGROW>
    end
    if rv >= 0
        xe = [sqrt(rv), -sqrt(rv)];
        for xeq = xe
            lam = rv - 3*xeq^2;
            if lam < 0
                pf_stable = [pf_stable; rv xeq]; %#ok<AGROW>
            else
                pf_unstable = [pf_unstable; rv xeq]; %#ok<AGROW>
            end
        end
    end

    % Saddle-node: x*= +-sqrt(r) for r>=0
    if rv >= 0
        xplus = sqrt(rv); xminus = -sqrt(rv);
        if -2*xplus < 0, sn_stable = [sn_stable; rv xplus]; else, sn_unstable = [sn_unstable; rv xplus]; end %#ok<AGROW>
        if -2*xminus < 0, sn_stable = [sn_stable; rv xminus]; else, sn_unstable = [sn_unstable; rv xminus]; end %#ok<AGROW>
    end
end

figure('Name','Bifurcation Branches');
hold on; grid on;
if ~isempty(pf_stable),   plot(pf_stable(:,1), pf_stable(:,2),  'LineWidth', 1.6); end
if ~isempty(pf_unstable), plot(pf_unstable(:,1), pf_unstable(:,2),'--','LineWidth', 1.2); end
if ~isempty(sn_stable),   plot(sn_stable(:,1), sn_stable(:,2),   'LineWidth', 1.6); end
if ~isempty(sn_unstable), plot(sn_unstable(:,1), sn_unstable(:,2),'--','LineWidth', 1.2); end
xlabel('r'); ylabel('x^*');
title('Pitchfork and Saddle-Node Equilibrium Branches');
legend('Pitchfork stable','Pitchfork unstable','Saddle-node stable','Saddle-node unstable');

%% Part C: Simulink (optional) - programmatically build a Van der Pol model
% This section creates a simple Simulink model with two integrators if Simulink exists.
if license('test','Simulink') && exist('new_system','file')
    mdl = 'Chapter14_Lesson5_SimulinkModel';
    if bdIsLoaded(mdl), close_system(mdl,0); end
    if exist([mdl '.slx'],'file'), delete([mdl '.slx']); end
    new_system(mdl); open_system(mdl);

    add_block('simulink/Sources/Constant', [mdl '/mu'], 'Position', [30 40 60 60], 'Value', '1');
    add_block('simulink/Continuous/Integrator', [mdl '/Int_x'], 'Position', [420 90 450 120]);
    add_block('simulink/Continuous/Integrator', [mdl '/Int_y'], 'Position', [420 180 450 210]);

    add_block('simulink/Math Operations/Product', [mdl '/x2'], 'Position', [120 140 150 170], 'Inputs', '**');
    add_block('simulink/Math Operations/Sum', [mdl '/1_minus_x2'], 'Position', [180 140 205 170], 'Inputs', '+-');
    add_block('simulink/Sources/Constant', [mdl '/one'], 'Position', [120 90 150 110], 'Value', '1');

    add_block('simulink/Math Operations/Product', [mdl '/mu_term'], 'Position', [240 120 270 150], 'Inputs', '**');
    add_block('simulink/Math Operations/Product', [mdl '/mu_term_times_y'], 'Position', [300 120 330 150], 'Inputs', '**');
    add_block('simulink/Math Operations/Gain', [mdl '/minus_x'], 'Position', [300 200 340 230], 'Gain', '-1');
    add_block('simulink/Math Operations/Sum', [mdl '/ydot_sum'], 'Position', [360 155 385 185], 'Inputs', '++');

    add_block('simulink/Sinks/Scope', [mdl '/Scope'], 'Position', [520 90 550 120]);

    % Connect x and y states
    add_line(mdl, 'Int_x/1', 'x2/1');
    add_line(mdl, 'Int_x/1', 'x2/2');
    add_line(mdl, 'one/1', '1_minus_x2/1');
    add_line(mdl, 'x2/1', '1_minus_x2/2');
    add_line(mdl, 'mu/1', 'mu_term/1');
    add_line(mdl, '1_minus_x2/1', 'mu_term/2');
    add_line(mdl, 'mu_term/1', 'mu_term_times_y/1');
    add_line(mdl, 'Int_y/1', 'mu_term_times_y/2');
    add_line(mdl, 'mu_term_times_y/1', 'ydot_sum/1');
    add_line(mdl, 'Int_x/1', 'minus_x/1');
    add_line(mdl, 'minus_x/1', 'ydot_sum/2');

    % State equations: xdot = y, ydot = mu(1-x^2)y - x
    add_line(mdl, 'Int_y/1', 'Int_x/1', 'autorouting', 'on');
    add_line(mdl, 'ydot_sum/1', 'Int_y/1', 'autorouting', 'on');
    add_line(mdl, 'Int_x/1', 'Scope/1', 'autorouting', 'on');

    set_param(mdl, 'StopTime', '40');
    save_system(mdl);
    fprintf('Simulink model created: %s.slx\n', mdl);
else
    fprintf('Simulink not available. MATLAB ODE implementation executed.\n');
end

Code: Chapter14_Lesson5.nb

(* Content-type: application/vnd.wolfram.mathematica *)
Notebook[{
 Cell["Chapter14_Lesson5.nb - Limit Cycles, Multiple Equilibria, and Basic Bifurcation Notions", "Title"],
 Cell["Van der Pol oscillator, 1D bifurcation normal forms, and a basic Hopf normal form check.", "Text"],

 Cell["Van der Pol simulation (NDSolveValue)", "Section"],
 Cell[BoxData[
  ToBoxes[
   mu = 1;
   sol = NDSolveValue[
     {
      x'[t] == y[t],
      y'[t] == mu (1 - x[t]^2) y[t] - x[t],
      x[0] == 2, y[0] == 0.1
     },
     {x, y},
     {t, 0, 80},
     MaxStepFraction -> 1/500
   ];
  ]
 ], "Input"],

 Cell[BoxData[
  ToBoxes[
   ParametricPlot[
    Evaluate[{sol[[1]][t], sol[[2]][t]}],
    {t, 40, 80},
    AxesLabel -> {"x", "xdot"},
    PlotLabel -> "Van der Pol Limit Cycle (phase plane)",
    PlotRange -> All
   ]
  ]
 ], "Input"],

 Cell["Estimate period using upward crossings of y[t]==0 with x[t]>0", "Text"],
 Cell[BoxData[
  ToBoxes[
   pts = Table[{tt, sol[[1]][tt], sol[[2]][tt]}, {tt, 40, 80, 0.01}];
   crossTimes = Reap[
      Do[
       If[pts[[k, 3]] < 0 && pts[[k + 1, 3]] >= 0,
        alpha = -pts[[k, 3]]/(pts[[k + 1, 3]] - pts[[k, 3]] + 10^-14);
        tc = pts[[k, 1]] + alpha (pts[[k + 1, 1]] - pts[[k, 1]]);
        xc = pts[[k, 2]] + alpha (pts[[k + 1, 2]] - pts[[k, 2]]);
        If[xc > 0, Sow[tc]];
       ],
       {k, 1, Length[pts] - 1}
      ]
    ][[2, 1]];
   periodEstimate = Mean[Differences[crossTimes[[-5 ;; -1]]]];
   ampEstimate = Max[Abs[pts[[All, 2]]]];
   {ampEstimate, periodEstimate}
  ]
 ], "Input"],

 Cell["1D bifurcation normal forms", "Section"],
 Cell[BoxData[
  ToBoxes[
   pitchforkStable = Join[
     Table[{r, 0}, {r, -2, -10^-3, 0.01}],
     Flatten[Table[{ {r, Sqrt[r]}, {r, -Sqrt[r]} }, {r, 0, 2, 0.01}], 1]
   ];
   pitchforkUnstable = Table[{r, 0}, {r, 0, 2, 0.01}];
   saddleStable = Table[{r, Sqrt[r]}, {r, 0, 2, 0.01}];
   saddleUnstable = Table[{r, -Sqrt[r]}, {r, 0, 2, 0.01}];
   Show[
    ListLinePlot[{pitchforkStable, saddleStable}, PlotRange -> All],
    ListLinePlot[{pitchforkUnstable, saddleUnstable}, PlotStyle -> Dashed],
    AxesLabel -> {"r", "x*"},
    PlotLabel -> "Pitchfork and Saddle-Node Branches"
   ]
  ]
 ], "Input"],

 Cell["Supercritical Hopf normal form in polar coordinates", "Section"],
 Cell[BoxData[
  ToBoxes[
   Clear[rho, th, muH, om];
   muH = 0.5; om = 1.0;
   hopf = NDSolveValue[
    {
     rho'[t] == muH rho[t] - rho[t]^3,
     th'[t] == om,
     rho[0] == 0.2, th[0] == 0
    },
    {rho, th}, {t, 0, 40}
   ];
   ParametricPlot[
    Evaluate[{hopf[[1]][t] Cos[hopf[[2]][t]], hopf[[1]][t] Sin[hopf[[2]][t]]}],
    {t, 0, 40},
    PlotLabel -> "Supercritical Hopf Normal Form Orbit",
    AxesLabel -> {"x", "y"}
   ]
  ]
 ], "Input"]
},
WindowSize -> {1000, 700},
WindowTitle -> "Chapter14_Lesson5"
]

Problem 1 (Saddle-node equilibria and stability). For \( \dot{x}=r-x^2 \), determine the number of equilibria for \( r \lt 0 \), \( r=0 \), and \( r \gt 0 \), and classify their stability.

Solution. Solve \( r-x^2=0 \).

• For \( r \lt 0 \): no real equilibrium.
• For \( r=0 \): one nonhyperbolic equilibrium \( x^\star=0 \).
• For \( r \gt 0 \): two equilibria \( x^\star=\pm\sqrt{r} \).

Since \( f_x=-2x \), the equilibrium \( +\sqrt{r} \) is stable and \( -\sqrt{r} \) is unstable.

Problem 2 (Pitchfork local behavior). For \( \dot{x}=rx-x^3 \), prove that the origin changes stability at \( r=0 \) and identify the stable equilibria for \( r \gt 0 \).

Solution. Equilibria satisfy \( x(r-x^2)=0 \), so \( x^\star=0 \) always, and \( x^\star=\pm\sqrt{r} \) for \( r \gt 0 \). The derivative is \( f_x=r-3x^2 \). Then \( f_x(0)=r \) (stable if \( r \lt 0 \), unstable if \( r \gt 0 \)) and \( f_x(\pm\sqrt{r})=-2r \lt 0 \) for \( r \gt 0 \), so the two nonzero branches are stable.

Problem 3 (Energy sign in Van der Pol). Let \( V=\frac{1}{2}(x^2+\dot{x}^2) \) for \( \ddot{x}-\mu(1-x^2)\dot{x}+x=0 \). Compute \( \dot{V} \) and interpret the sign.

Solution.

\[ \dot{V} = x\dot{x} + \dot{x}\ddot{x} = x\dot{x} + \dot{x}\big(\mu(1-x^2)\dot{x} - x\big) = \mu(1-x^2)\dot{x}^2. \]

Therefore \( \dot{V} \gt 0 \) for \( |x| \lt 1 \) (small oscillations grow) and \( \dot{V} \lt 0 \) for \( |x| \gt 1 \) (large oscillations decay). This sign change is the mechanism behind amplitude regulation and the attracting limit cycle.

Problem 4 (Hopf normal-form limit cycle radius). For \( \dot{\rho}=\mu\rho-\rho^3 \), prove that when \( \mu \gt 0 \), the circle \( \rho=\sqrt{\mu} \) is asymptotically stable.

Solution. The scalar radial dynamics satisfies \( \dot{\rho}=\rho(\mu-\rho^2) \). If \( 0 \lt \rho \lt \sqrt{\mu} \), then \( \dot{\rho}\gt 0 \). If \( \rho \gt \sqrt{\mu} \), then \( \dot{\rho}\lt 0 \). Hence all positive radii near the branch converge to \( \sqrt{\mu} \). With \( \dot{\theta}=\omega\neq 0 \), this is a stable periodic orbit in the plane.

Problem 5 (Numerical period estimation by section crossings). Suppose a simulated 2D trajectory crosses a Poincaré section at times \( t_1,t_2,\dots,t_N \) after transients. Show a consistent estimator of the period and explain why transients must be discarded.

Solution. Use

\[ \hat{T} = \frac{1}{N-1}\sum_{k=1}^{N-1}(t_{k+1}-t_k). \]

During transients, the trajectory has not yet converged to the periodic orbit, so the inter-crossing intervals are not constant. After convergence, section crossings occur once per period (for a well-defined section and orientation), making the estimator consistent.

We introduced limit cycles as isolated periodic orbits, showed how multiple equilibria arise in nonlinear scalar systems, and formalized basic bifurcation notions through saddle-node, pitchfork, and Hopf normal forms. We also connected these concepts to engineering analysis using Jacobians, phase-plane reasoning, and numerical period estimation. The multi-language implementations in this lesson provide a reusable computational foundation for later chapters on numerical simulation and complex dynamics.