A planar autonomous system is written as \( \dot{x} = f(x,y) \), \( \dot{y} = g(x,y) \), or compactly as \( \dot{\mathbf{z} } = \mathbf{F}(\mathbf{z}) \) with \( \mathbf{z}=[x\;y]^\top \in \mathbb{R}^2 \) and \( \mathbf{F}=[f\;g]^\top \).

Throughout, assume \( f,g \) are at least locally Lipschitz so that the initial value problem has a unique maximal solution. Formally, for any initial condition \( \mathbf{z}(t_0)=\mathbf{z}_0 \), the IVP

\[ \dot{\mathbf{z} } = \mathbf{F}(\mathbf{z}), \qquad \mathbf{z}(t_0)=\mathbf{z}_0 \]

has a unique trajectory \( \mathbf{z}(t) \) on some interval \( (t_0-\epsilon, t_0+\epsilon) \).

flowchart TD
  A["Given: z_dot = F(z) in R2"] --> B["Step 1: Find equilibria: solve F(z)=0"]
  B --> C["Step 2: Compute Jacobian J = dF/dz at each equilibrium"]
  C --> D["Step 3: Local type from eigenvalues / trace-det test"]
  D --> E["Step 4: Compute nullclines: f(x,y)=0 and g(x,y)=0"]
  E --> F["Step 5: Determine vector directions in each region"]
  F --> G["Step 6: Integrate trajectories from representative initial points"]
  G --> H["Step 7: Assemble phase portrait and interpret"]
        

An equilibrium point (fixed point) is any \( \mathbf{z}^* \in \mathbb{R}^2 \) satisfying \( \mathbf{F}(\mathbf{z}^*)=\mathbf{0} \), i.e. \( f(x^*,y^*)=0 \) and \( g(x^*,y^*)=0 \).

Proposition (equilibria are invariant). If \( \mathbf{z}(t_0)=\mathbf{z}^* \) and \( \mathbf{F}(\mathbf{z}^*)=\mathbf{0} \), then \( \mathbf{z}(t)=\mathbf{z}^* \) for all times where the solution exists.

Proof. Consider the constant function \( \tilde{\mathbf{z} }(t)=\mathbf{z}^* \). Then \( \dot{\tilde{\mathbf{z} }}(t)=\mathbf{0}=\mathbf{F}(\mathbf{z}^*)=\mathbf{F}(\tilde{\mathbf{z} }(t)) \), so \( \tilde{\mathbf{z} } \) satisfies the ODE and the initial condition. By uniqueness of solutions, \( \mathbf{z}(t)\equiv\tilde{\mathbf{z} }(t) \). □

In applications, solving \( f=0 \), \( g=0 \) can be done analytically (algebraic elimination) or numerically (Newton, continuation, etc.). For 2D systems, geometry via nullclines is especially informative.

The x-nullcline is the set where \( \dot{x}=0 \), i.e. \( f(x,y)=0 \). The y-nullcline is the set where \( \dot{y}=0 \), i.e. \( g(x,y)=0 \). Intersections of nullclines are precisely equilibria.

On the x-nullcline, the vector field has form \( (0,\dot{y}) \) (purely vertical); on the y-nullcline, it has form \( (\dot{x},0) \) (purely horizontal). Between nullclines, the signs of \( f \) and \( g \) determine the quadrant of the local direction.

flowchart LR
  N1["x-nullcline: f(x,y)=0 (x_dot=0)"] --> V["Vectors vertical (up/down)"]
  N2["y-nullcline: g(x,y)=0 (y_dot=0)"] --> H["Vectors horizontal (left/right)"]
  I["Intersection points"] --> E["Equilibria: F(z*)=0"]
  R["Regions between nullclines"] --> S["Sign test: sign(f), sign(g)"]
  S --> Q1["f>0, g>0: move right & up"]
  S --> Q2["f<0, g>0: move left & up"]
  S --> Q3["f<0, g<0: move left & down"]
  S --> Q4["f>0, g<0: move right & down"]
        

A trajectory is the curve traced by \( \mathbf{z}(t) \) in the phase plane. For autonomous systems, the phase portrait is time-translation invariant: if \( \mathbf{z}(t) \) is a solution, then so is \( \mathbf{z}(t+c) \) for any constant \( c \), provided both are defined.

Theorem (no crossing of trajectories). Suppose \( \mathbf{F} \) is locally Lipschitz. If two solutions \( \mathbf{z}_1(t) \) and \( \mathbf{z}_2(t) \) satisfy \( \mathbf{z}_1(t_0)=\mathbf{z}_2(t_0) \) for some \( t_0 \), then \( \mathbf{z}_1(t)=\mathbf{z}_2(t) \) for all times in the common interval of existence. Consequently, distinct trajectories in the phase plane cannot intersect.

Proof. Both trajectories solve the same IVP \( \dot{\mathbf{z} }=\mathbf{F}(\mathbf{z}) \) with the same initial value \( \mathbf{z}(t_0)=\mathbf{z}_0 \). By uniqueness under local Lipschitz continuity, the solution is unique, hence \( \mathbf{z}_1 \equiv \mathbf{z}_2 \). □

This is the rigorous reason “phase curves do not cross” (unlike arbitrary curves). It is also why separatrices of a saddle act as boundaries between regions with qualitatively different behavior.

Let \( \mathbf{z}^* \) be an equilibrium. Define the perturbation \( \mathbf{u}=\mathbf{z}-\mathbf{z}^* \). The first-order Taylor expansion yields

\[ \dot{\mathbf{u} } = \mathbf{J}(\mathbf{z}^*)\,\mathbf{u} + \mathcal{O}(\|\mathbf{u}\|^2), \qquad \mathbf{J}(\mathbf{z}^*) = \begin{bmatrix} \frac{\partial f}{\partial x} & \frac{\partial f}{\partial y} \\ \frac{\partial g}{\partial x} & \frac{\partial g}{\partial y} \end{bmatrix}_{(x^*,y^*)}. \]

The linearized system \( \dot{\mathbf{u} }=\mathbf{J}\mathbf{u} \) has characteristic polynomial

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

Hence eigenvalues are

\[ \lambda_{1,2} = \frac{\operatorname{tr}\mathbf{J} \pm \sqrt{(\operatorname{tr}\mathbf{J})^2 - 4\det\mathbf{J} }}{2}. \]

Define invariants \( \tau=\operatorname{tr}\mathbf{J} \), \( \Delta=\det\mathbf{J} \), and discriminant \( D = \tau^2 - 4\Delta \). The usual planar classification is:

• \( \Delta < 0 \): saddle (one eigenvalue positive, one negative).
• \( \Delta > 0,\; D > 0 \): node (real eigenvalues); stable if \( \tau < 0 \), unstable if \( \tau > 0 \).
• \( \Delta > 0,\; D < 0 \): spiral/focus (complex pair); stable if \( \tau < 0 \), unstable if \( \tau > 0 \).
• \( \Delta > 0,\; \tau = 0,\; D < 0 \): linear center; nonlinear terms decide whether a true center persists.
• \( \Delta = 0 \): degenerate/non-isolated (requires higher-order analysis).

For hyperbolic equilibria (no eigenvalues with zero real part), the linearization predicts the local phase portrait up to a homeomorphism (Hartman–Grobman principle). In this lesson, you use it as a practical guideline: compute \( \mathbf{J} \), evaluate \( (\tau,\Delta,D) \), and classify.

Consider the planar system (parameter \( a > 0 \)):

\[ \dot{x} = x - x^3 - y, \qquad \dot{y} = x + a y. \]

Equilibria satisfy \( y = x - x^3 \) (from \( \dot{x}=0 \)) and \( y = -x/a \) (from \( \dot{y}=0 \)). Eliminating \( y \) gives

\[ x - x^3 = -\frac{x}{a} \quad \Longrightarrow \quad x\left(\frac{a+1}{a} - x^2\right)=0. \]

Hence \( x^*_0=0 \) and \( x^*_{\pm}=\pm\sqrt{(a+1)/a} \). The corresponding \( y^* \) values are \( y^*=-x^*/a \).

The Jacobian is

\[ \mathbf{J}(x,y)= \begin{bmatrix} 1-3x^2 & -1 \\ 1 & a \end{bmatrix}. \]

For the illustrative choice \( a=1 \):

• At \( (0,0) \), \( \tau=2 \), \( \Delta=2 \), and \( D=4-8=-4 < 0 \), so the linearization is an unstable spiral.
• At \( (\pm\sqrt2,\mp\sqrt2) \), \( \Delta < 0 \), so each is a saddle.

The global phase portrait combines these local features with the nullclines’ geometry. In computations below, you will generate a stream plot and integrate representative trajectories to reveal separatrix structure around saddles.

This script uses scipy.optimize.root to locate equilibria from multiple guesses, classifies each equilibrium using the \( (\tau,\Delta,D) \) test, and plots a stream-based phase portrait with trajectories via solve_ivp. Libraries: numpy, scipy, matplotlib.

Chapter14_Lesson2.py

"""
Chapter14_Lesson2.py
Equilibrium points, phase portraits, and trajectories for planar (2D) autonomous systems.

Requires:
  numpy, scipy, matplotlib
Install (example):
  pip install numpy scipy matplotlib
"""

from __future__ import annotations

import numpy as np
from dataclasses import dataclass
from typing import Callable, List, Tuple, Dict

from scipy.integrate import solve_ivp
from scipy.optimize import root
import matplotlib.pyplot as plt


# -----------------------------
# 1) Define a 2D autonomous system
# -----------------------------
@dataclass(frozen=True)
class PlanarSystem:
    """Planar autonomous system: z' = F(z) with z=(x,y)."""
    a: float = 1.0  # parameter for example system

    def f(self, x: float, y: float) -> float:
        # dx/dt
        return x - x**3 - y

    def g(self, x: float, y: float) -> float:
        # dy/dt
        return x + self.a * y

    def F(self, t: float, z: np.ndarray) -> np.ndarray:
        x, y = float(z[0]), float(z[1])
        return np.array([self.f(x, y), self.g(x, y)], dtype=float)


# -----------------------------
# 2) Equilibria and Jacobian
# -----------------------------
def jacobian_fd(sys: PlanarSystem, z: np.ndarray, h: float = 1e-6) -> np.ndarray:
    """Finite-difference Jacobian J = dF/dz at z (2x2)."""
    z = np.asarray(z, dtype=float)
    J = np.zeros((2, 2), dtype=float)
    f0 = sys.F(0.0, z)
    for i in range(2):
        zp = z.copy()
        zp[i] += h
        fp = sys.F(0.0, zp)
        J[:, i] = (fp - f0) / h
    return J


def find_equilibria(
    sys: PlanarSystem,
    guesses: List[Tuple[float, float]],
    tol: float = 1e-8
) -> List[np.ndarray]:
    """
    Find equilibria by running a nonlinear root solver from multiple initial guesses,
    then deduplicating by tolerance.
    """
    sols: List[np.ndarray] = []

    def fun(z):
        return sys.F(0.0, np.asarray(z, dtype=float))

    for gx, gy in guesses:
        res = root(fun, x0=np.array([gx, gy], dtype=float), method="hybr")
        if not res.success:
            continue
        z = res.x.astype(float)
        # Deduplicate
        if all(np.linalg.norm(z - s) > tol for s in sols):
            sols.append(z)
    return sols


# -----------------------------
# 3) Local classification via trace/determinant and eigenvalues
# -----------------------------
@dataclass(frozen=True)
class LinClass:
    kind: str
    trace: float
    det: float
    disc: float
    eig: Tuple[complex, complex]


def classify_2x2(J: np.ndarray, eps: float = 1e-10) -> LinClass:
    """
    Classify planar linearization z' = J z at an equilibrium using invariants:
      tr = trace(J), det = det(J), disc = tr^2 - 4 det.
    """
    tr = float(np.trace(J))
    det = float(np.linalg.det(J))
    disc = float(tr * tr - 4.0 * det)
    eigvals = np.linalg.eigvals(J)
    lam1, lam2 = complex(eigvals[0]), complex(eigvals[1])

    # Robust comparisons with eps
    if det < -eps:
        kind = "saddle (hyperbolic)"
    elif abs(det) <= eps:
        kind = "non-isolated or degenerate (det ~ 0)"
    else:
        # det > 0
        if disc > eps:
            # real eigenvalues
            if tr < -eps:
                kind = "stable node"
            elif tr > eps:
                kind = "unstable node"
            else:
                kind = "improper / star node (tr ~ 0, disc > 0)"
        elif disc < -eps:
            # complex conjugate
            if tr < -eps:
                kind = "stable spiral (focus)"
            elif tr > eps:
                kind = "unstable spiral (focus)"
            else:
                kind = "center (linear); nonlinear terms decide"
        else:
            # disc ~ 0: repeated eigenvalue
            if tr < -eps:
                kind = "stable degenerate node"
            elif tr > eps:
                kind = "unstable degenerate node"
            else:
                kind = "center/degenerate (tr ~ 0, disc ~ 0)"
    return LinClass(kind=kind, trace=tr, det=det, disc=disc, eig=(lam1, lam2))


# -----------------------------
# 4) Phase portrait plotting + trajectories
# -----------------------------
def plot_phase_portrait(
    sys: PlanarSystem,
    equilibria: List[np.ndarray],
    xlim=(-2.5, 2.5),
    ylim=(-2.5, 2.5),
    ngrid: int = 25,
    tmax: float = 12.0,
    initials: List[Tuple[float, float]] | None = None,
    title: str = "Phase portrait"
) -> None:
    if initials is None:
        initials = [(-2, -2), (-2, 0), (-2, 2), (0.5, -2), (0.5, 2), (2, -2), (2, 0), (2, 2)]

    # Vector field grid
    xs = np.linspace(xlim[0], xlim[1], ngrid)
    ys = np.linspace(ylim[0], ylim[1], ngrid)
    X, Y = np.meshgrid(xs, ys)
    U = np.zeros_like(X)
    V = np.zeros_like(Y)
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            U[i, j] = sys.f(X[i, j], Y[i, j])
            V[i, j] = sys.g(X[i, j], Y[i, j])

    plt.figure()
    plt.streamplot(X, Y, U, V, density=1.0)

    # Plot equilibria
    for z in equilibria:
        plt.plot([z[0]], [z[1]], marker="o")
        J = jacobian_fd(sys, z)
        cls = classify_2x2(J)
        plt.text(z[0] + 0.05, z[1] + 0.05, cls.kind, fontsize=8)

    # Plot trajectories
    def rhs(t, z):
        return sys.F(t, z)

    for x0, y0 in initials:
        sol = solve_ivp(rhs, (0.0, tmax), y0=np.array([x0, y0], dtype=float),
                        rtol=1e-7, atol=1e-9, max_step=0.05)
        plt.plot(sol.y[0], sol.y[1], linewidth=1.0)

    plt.xlim(*xlim)
    plt.ylim(*ylim)
    plt.xlabel("x")
    plt.ylabel("y")
    plt.title(title)
    plt.grid(True)
    plt.tight_layout()
    plt.show()


def main() -> None:
    sys = PlanarSystem(a=1.0)

    # Try a grid of initial guesses to find equilibria
    guesses = []
    for gx in np.linspace(-2.0, 2.0, 9):
        for gy in np.linspace(-2.0, 2.0, 9):
            guesses.append((float(gx), float(gy)))

    eqs = find_equilibria(sys, guesses)
    print("Equilibria (approx):")
    for z in eqs:
        J = jacobian_fd(sys, z)
        cls = classify_2x2(J)
        print(f"  z* = [{z[0]: .6f}, {z[1]: .6f}]  trace={cls.trace: .6f} det={cls.det: .6f} -> {cls.kind}")

    plot_phase_portrait(sys, eqs, title="Example planar nonlinear system: x' = x - x^3 - y,  y' = x + a y")


if __name__ == "__main__":
    main()
      

This C++ program implements the same example system, classifies equilibria from the analytic Jacobian, and generates trajectory samples using a fixed-step RK4 integrator. It writes trajectories.csv for plotting elsewhere (MATLAB/Python/Excel). This is deliberately “from scratch” and avoids external libraries.

Chapter14_Lesson2.cpp

/*
Chapter14_Lesson2.cpp
Equilibrium points, local classification, and trajectory simulation (RK4) for a 2D autonomous system.

Build (example):
  g++ -O2 -std=c++17 Chapter14_Lesson2.cpp -o Chapter14_Lesson2

Run:
  ./Chapter14_Lesson2
Outputs:
  trajectories.csv   (columns: traj_id,t,x,y)
*/

#include 
#include 
#include 
#include 
#include 
#include 

struct Vec2 {
  double x;
  double y;
};

struct LinClass {
  std::string kind;
  double tr;
  double det;
  double disc;
};

static Vec2 F(double /*t*/, const Vec2& z, double a) {
  // Example system:
  // x' = x - x^3 - y
  // y' = x + a y
  return Vec2{z.x - z.x*z.x*z.x - z.y, z.x + a*z.y};
}

static void jacobian(const Vec2& z, double a, double J[2][2]) {
  // Analytical Jacobian for the example system
  // f_x = 1 - 3 x^2, f_y = -1
  // g_x = 1,        g_y = a
  J[0][0] = 1.0 - 3.0*z.x*z.x;
  J[0][1] = -1.0;
  J[1][0] = 1.0;
  J[1][1] = a;
}

static LinClass classify2x2(const double J[2][2], double eps=1e-12) {
  const double tr  = J[0][0] + J[1][1];
  const double det = J[0][0]*J[1][1] - J[0][1]*J[1][0];
  const double disc = tr*tr - 4.0*det;

  std::string kind;
  if (det < -eps) {
    kind = "saddle (hyperbolic)";
  } else if (std::abs(det) <= eps) {
    kind = "degenerate (det ~ 0)";
  } else {
    if (disc > eps) {
      if (tr < -eps) kind = "stable node";
      else if (tr > eps) kind = "unstable node";
      else kind = "improper/star node";
    } else if (disc < -eps) {
      if (tr < -eps) kind = "stable spiral (focus)";
      else if (tr > eps) kind = "unstable spiral (focus)";
      else kind = "center (linear); nonlinear decides";
    } else {
      if (tr < -eps) kind = "stable degenerate node";
      else if (tr > eps) kind = "unstable degenerate node";
      else kind = "degenerate/center";
    }
  }
  return LinClass{kind, tr, det, disc};
}

static Vec2 rk4_step(double t, const Vec2& z, double h, double a) {
  auto add = [](const Vec2& u, const Vec2& v) { return Vec2{u.x+v.x, u.y+v.y}; };
  auto mul = [](double c, const Vec2& v) { return Vec2{c*v.x, c*v.y}; };

  Vec2 k1 = F(t, z, a);
  Vec2 k2 = F(t + 0.5*h, add(z, mul(0.5*h, k1)), a);
  Vec2 k3 = F(t + 0.5*h, add(z, mul(0.5*h, k2)), a);
  Vec2 k4 = F(t + h, add(z, mul(h, k3)), a);

  Vec2 incr = add(add(mul(1.0, k1), mul(2.0, k2)), add(mul(2.0, k3), mul(1.0, k4)));
  return add(z, mul(h/6.0, incr));
}

int main() {
  const double a = 1.0;

  // Analytical equilibria for a>0:
  // y = x - x^3 and y = -x/a -> x[(a+1)/a - x^2] = 0
  const double xp = std::sqrt((a+1.0)/a);

  std::vector eqs = {
      {0.0, 0.0},
      {+xp, -xp},
      {-xp, +xp}
  };

  std::cout << "Equilibria and local classification (linearization):\n";
  for (const auto& z : eqs) {
    double J[2][2];
    jacobian(z, a, J);
    LinClass cls = classify2x2(J);
    std::cout << "  z* = (" << z.x << ", " << z.y << ")"
              << "  trace=" << cls.tr << " det=" << cls.det
              << " -> " << cls.kind << "\n";
  }

  // Trajectory simulation (RK4) and CSV output
  std::vector initials = {
      {-2.0, -2.0}, {-2.0, 0.0}, {-2.0, 2.0},
      { 0.5, -2.0}, { 0.5, 2.0},
      { 2.0, -2.0}, { 2.0, 0.0}, { 2.0, 2.0}
  };

  const double h = 0.01;
  const double tmax = 12.0;
  const int steps = static_cast(tmax / h);

  std::ofstream out("trajectories.csv");
  out << "traj_id,t,x,y\n";

  for (size_t k = 0; k < initials.size(); ++k) {
    Vec2 z = initials[k];
    double t = 0.0;
    for (int i = 0; i <= steps; ++i) {
      out << k << "," << t << "," << z.x << "," << z.y << "\n";
      z = rk4_step(t, z, h, a);
      t += h;
    }
  }

  std::cout << "\nWrote trajectories.csv (plot it in Python/MATLAB/Excel).\n";
  return 0;
}
      

Java implementation mirrors the C++ structure: analytic Jacobian classification plus RK4 trajectory generation into trajectories_java.csv using standard libraries only.

Chapter14_Lesson2.java

/*
Chapter14_Lesson2.java
Equilibria, local classification, and RK4 simulation for a planar autonomous system.

Compile:
  javac Chapter14_Lesson2.java
Run:
  java Chapter14_Lesson2
Outputs:
  trajectories_java.csv
*/

import java.io.FileWriter;
import java.io.IOException;
import java.util.ArrayList;
import java.util.List;

public class Chapter14_Lesson2 {

    static class Vec2 {
        double x, y;
        Vec2(double x, double y) { this.x = x; this.y = y; }
    }

    static class LinClass {
        String kind;
        double tr, det, disc;
        LinClass(String kind, double tr, double det, double disc) {
            this.kind = kind; this.tr = tr; this.det = det; this.disc = disc;
        }
    }

    static Vec2 F(double t, Vec2 z, double a) {
        // Example system:
        // x' = x - x^3 - y
        // y' = x + a y
        double dx = z.x - z.x*z.x*z.x - z.y;
        double dy = z.x + a*z.y;
        return new Vec2(dx, dy);
    }

    static double[][] jacobian(Vec2 z, double a) {
        // Analytical Jacobian
        double[][] J = new double[2][2];
        J[0][0] = 1.0 - 3.0*z.x*z.x;
        J[0][1] = -1.0;
        J[1][0] = 1.0;
        J[1][1] = a;
        return J;
    }

    static LinClass classify2x2(double[][] J, double eps) {
        double tr  = J[0][0] + J[1][1];
        double det = J[0][0]*J[1][1] - J[0][1]*J[1][0];
        double disc = tr*tr - 4.0*det;

        String kind;
        if (det < -eps) {
            kind = "saddle (hyperbolic)";
        } else if (Math.abs(det) <= eps) {
            kind = "degenerate (det ~ 0)";
        } else {
            if (disc > eps) {
                if (tr < -eps) kind = "stable node";
                else if (tr > eps) kind = "unstable node";
                else kind = "improper/star node";
            } else if (disc < -eps) {
                if (tr < -eps) kind = "stable spiral (focus)";
                else if (tr > eps) kind = "unstable spiral (focus)";
                else kind = "center (linear); nonlinear decides";
            } else {
                if (tr < -eps) kind = "stable degenerate node";
                else if (tr > eps) kind = "unstable degenerate node";
                else kind = "degenerate/center";
            }
        }
        return new LinClass(kind, tr, det, disc);
    }

    static Vec2 rk4Step(double t, Vec2 z, double h, double a) {
        Vec2 k1 = F(t, z, a);
        Vec2 k2 = F(t + 0.5*h, new Vec2(z.x + 0.5*h*k1.x, z.y + 0.5*h*k1.y), a);
        Vec2 k3 = F(t + 0.5*h, new Vec2(z.x + 0.5*h*k2.x, z.y + 0.5*h*k2.y), a);
        Vec2 k4 = F(t + h, new Vec2(z.x + h*k3.x, z.y + h*k3.y), a);

        double dx = (k1.x + 2.0*k2.x + 2.0*k3.x + k4.x) * (h/6.0);
        double dy = (k1.y + 2.0*k2.y + 2.0*k3.y + k4.y) * (h/6.0);
        return new Vec2(z.x + dx, z.y + dy);
    }

    public static void main(String[] args) throws IOException {
        final double a = 1.0;

        // Analytical equilibria for a>0
        double xp = Math.sqrt((a+1.0)/a);

        List eqs = new ArrayList<>();
        eqs.add(new Vec2(0.0, 0.0));
        eqs.add(new Vec2(+xp, -xp));
        eqs.add(new Vec2(-xp, +xp));

        System.out.println("Equilibria and local classification (linearization):");
        for (Vec2 z : eqs) {
            double[][] J = jacobian(z, a);
            LinClass cls = classify2x2(J, 1e-12);
            System.out.printf("  z*=(%.6f, %.6f)  trace=%.6f det=%.6f -> %s%n",
                    z.x, z.y, cls.tr, cls.det, cls.kind);
        }

        // Simulate trajectories and write CSV
        List initials = List.of(
                new Vec2(-2.0, -2.0), new Vec2(-2.0, 0.0), new Vec2(-2.0, 2.0),
                new Vec2( 0.5, -2.0), new Vec2( 0.5, 2.0),
                new Vec2( 2.0, -2.0), new Vec2( 2.0, 0.0), new Vec2( 2.0, 2.0)
        );

        double h = 0.01;
        double tmax = 12.0;
        int steps = (int)Math.round(tmax / h);

        try (FileWriter fw = new FileWriter("trajectories_java.csv")) {
            fw.write("traj_id,t,x,y\n");
            for (int k = 0; k < initials.size(); k++) {
                Vec2 z = new Vec2(initials.get(k).x, initials.get(k).y);
                double t = 0.0;
                for (int i = 0; i <= steps; i++) {
                    fw.write(k + "," + t + "," + z.x + "," + z.y + "\n");
                    z = rk4Step(t, z, h, a);
                    t += h;
                }
            }
        }

        System.out.println("\nWrote trajectories_java.csv");
    }
}
      

The MATLAB script draws a vector field with quiver, integrates trajectories with ode45, and prints local classifications using the trace–determinant test. A helper function builds a simple Simulink model programmatically (two Integrators + a MATLAB Function block). You can enable it by uncommenting the call.

Chapter14_Lesson2.m

% Chapter14_Lesson2.m
% Equilibrium points, phase portrait, and trajectories for a planar nonlinear system.
%
% Requires: MATLAB (Optimization Toolbox for fsolve is helpful but not strictly required).
% If you do not have fsolve, skip the equilibrium solver section and use the analytic equilibria.

function Chapter14_Lesson2()
    a = 1.0;

    % System definition: z = [x; y], z' = [f; g]
    f = @(x,y) x - x.^3 - y;
    g = @(x,y) x + a*y;
    F = @(t,z) [f(z(1), z(2)); g(z(1), z(2))];

    % ---- 1) Equilibria (analytic for this example) ----
    xp = sqrt((a+1)/a);
    eqs = [0, 0; +xp, -xp; -xp, +xp];

    % ---- 2) Local classification via Jacobian invariants ----
    fprintf('Equilibria and local classification (linearization):\n');
    for k = 1:size(eqs,1)
        x = eqs(k,1); y = eqs(k,2);
        J = [1 - 3*x^2, -1; 1, a];
        tr = trace(J);
        detJ = det(J);
        disc = tr^2 - 4*detJ;

        kind = classify2x2(tr, detJ, disc);
        fprintf('  z*=(% .6f, % .6f)  trace=% .6f det=% .6f -> %s\n', x, y, tr, detJ, kind);
    end

    % ---- 3) Phase portrait (vector field) ----
    figure; hold on; grid on;
    xlim([-2.5 2.5]); ylim([-2.5 2.5]);
    xlabel('x'); ylabel('y');
    title('Phase portrait: x''=x-x^3-y, y''=x+a y');

    [X,Y] = meshgrid(linspace(-2.5,2.5,25), linspace(-2.5,2.5,25));
    U = f(X,Y);
    V = g(X,Y);
    quiver(X,Y,U,V,'AutoScale','on');

    % Plot equilibria
    plot(eqs(:,1), eqs(:,2), 'o', 'MarkerSize', 7);

    % ---- 4) Trajectories via ode45 ----
    initials = [-2 -2; -2 0; -2 2; 0.5 -2; 0.5 2; 2 -2; 2 0; 2 2];
    for k = 1:size(initials,1)
        z0 = initials(k,:)';
        [t, z] = ode45(F, [0 12], z0);
        plot(z(:,1), z(:,2), 'LineWidth', 1.0);
    end
    hold off;

    % ---- 5) Simulink (programmatic build) ----
    % This builds a simple Simulink model with two Integrators and a MATLAB Function block.
    % Uncomment to generate and save an .slx model:
    % build_simulink_model(a);
end

function kind = classify2x2(tr, detJ, disc)
    eps = 1e-12;
    if detJ < -eps
        kind = 'saddle (hyperbolic)';
        return;
    end
    if abs(detJ) <= eps
        kind = 'degenerate (det ~ 0)';
        return;
    end

    if disc > eps
        if tr < -eps
            kind = 'stable node';
        elseif tr > eps
            kind = 'unstable node';
        else
            kind = 'improper/star node';
        end
    elseif disc < -eps
        if tr < -eps
            kind = 'stable spiral (focus)';
        elseif tr > eps
            kind = 'unstable spiral (focus)';
        else
            kind = 'center (linear); nonlinear decides';
        end
    else
        if tr < -eps
            kind = 'stable degenerate node';
        elseif tr > eps
            kind = 'unstable degenerate node';
        else
            kind = 'degenerate/center';
        end
    end
end

function build_simulink_model(a)
    model = 'Chapter14_Lesson2_Simulink';
    if bdIsLoaded(model); close_system(model, 0); end
    new_system(model); open_system(model);

    % Add blocks
    add_block('simulink/Sources/Constant', [model '/x0'], 'Value', '0.5');
    add_block('simulink/Sources/Constant', [model '/y0'], 'Value', '2.0');
    add_block('simulink/Continuous/Integrator', [model '/Int_x'], 'InitialCondition', 'x0');
    add_block('simulink/Continuous/Integrator', [model '/Int_y'], 'InitialCondition', 'y0');

    % MATLAB Function block computing derivatives
    add_block('simulink/User-Defined Functions/MATLAB Function', [model '/Dynamics']);
    set_param([model '/Dynamics'], 'Script', sprintf([ ...
        'function dz = fcn(z)\n' ...
        '%% z = [x; y]\n' ...
        'x = z(1); y = z(2);\n' ...
        'dx = x - x^3 - y;\n' ...
        'dy = x + %.15g*y;\n' ...
        'dz = [dx; dy];\n' ...
        'end\n'], a));

    % Mux/Demux
    add_block('simulink/Signal Routing/Mux', [model '/Mux'], 'Inputs', '2');
    add_block('simulink/Signal Routing/Demux', [model '/Demux'], 'Outputs', '2');

    % Scopes
    add_block('simulink/Sinks/XY Graph', [model '/XY']);
    add_block('simulink/Sinks/Scope', [model '/Scope']);

    % Wiring
    add_line(model, 'x0/1', 'Int_x/1');
    add_line(model, 'y0/1', 'Int_y/1');

    add_line(model, 'Int_x/1', 'Mux/1');
    add_line(model, 'Int_y/1', 'Mux/2');
    add_line(model, 'Mux/1', 'Dynamics/1');

    add_line(model, 'Dynamics/1', 'Demux/1');
    add_line(model, 'Demux/1', 'Int_x/1', 'autorouting', 'on');
    add_line(model, 'Demux/2', 'Int_y/1', 'autorouting', 'on');

    add_line(model, 'Int_x/1', 'XY/1');
    add_line(model, 'Int_y/1', 'XY/2');

    add_line(model, 'Mux/1', 'Scope/1');

    % Layout niceness
    set_param(model, 'StopTime', '12');
    save_system(model);
    fprintf('Saved Simulink model: %s.slx\n', model);
end
      

The notebook below defines the system, computes equilibria, classifies via invariants, and generates a vector/stream plot with sample trajectories using NDSolveValue.

Chapter14_Lesson2.nb

(* Chapter14_Lesson2.nb
   Wolfram Mathematica notebook (plain-text Notebook expression) implementing
   equilibria, classification, and phase portrait for a planar system.
*)

Notebook[{
  Cell["Chapter 14 - Lesson 2: Equilibria, Phase Portraits, Trajectories (2D Systems)", "Title"],

  Cell["System definition", "Section"],
  Cell[BoxData[
    ToBoxes[
      Row[{
        "a = 1.0; ",
        "f[x_, y_] := x - x^3 - y; ",
        "g[x_, y_] := x + a y;"
      }]
    ]
  ], "Input"],

  Cell["Equilibria (analytic for this example)", "Section"],
  Cell[BoxData[
    ToBoxes[
      Row[{
        "xp = Sqrt[(a + 1)/a]; ",
        "eqs = { {0, 0}, {xp, -xp}, {-xp, xp} };"
      }]
    ]
  ], "Input"],

  Cell["Local classification using trace/determinant", "Section"],
  Cell[BoxData[
    ToBoxes[
      Row[{
        "J[x_, y_] := { {1 - 3 x^2, -1}, {1, a} }; ",
        "classify2x2[tr_, det_, disc_] := Which[",
        "det < 0, \"saddle\", ",
        "det == 0, \"degenerate\", ",
        "disc > 0 && tr < 0, \"stable node\", ",
        "disc > 0 && tr > 0, \"unstable node\", ",
        "disc < 0 && tr < 0, \"stable spiral\", ",
        "disc < 0 && tr > 0, \"unstable spiral\", ",
        "disc < 0 && tr == 0, \"center (linear)\", ",
        "True, \"other/degenerate\"",
        "];"
      }]
    ]
  ], "Input"],

  Cell[BoxData[
    ToBoxes[
      "Table[{p, tr = Tr[J @@ p], det = Det[J @@ p], disc = tr^2 - 4 det, classify2x2[tr, det, disc]}, {p, eqs}]"
    ]
  ], "Input"],

  Cell["Phase portrait (vector field + equilibria)", "Section"],
  Cell[BoxData[
    ToBoxes[
      Row[{
        "vf = VectorPlot[{f[x, y], g[x, y]}, {x, -2.5, 2.5}, {y, -2.5, 2.5}, ",
        "VectorPoints -> 25, StreamPoints -> Fine, PlotRange -> All]; ",
        "pts = Graphics[{PointSize[0.02], Point[eqs]}]; ",
        "Show[vf, pts]"
      }]
    ]
  ], "Input"],

  Cell["Sample trajectories via NDSolve", "Section"],
  Cell[BoxData[
    ToBoxes[
      Row[{
        "sol[z0_] := NDSolveValue[{x'[t] == f[x[t], y[t]], y'[t] == g[x[t], y[t]], ",
        "x[0] == z0[[1]], y[0] == z0[[2]]}, {x, y}, {t, 0, 12}];"
      }]
    ]
  ], "Input"],

  Cell[BoxData[
    ToBoxes[
      Row[{
        "inits = { {-2, -2}, {-2, 0}, {-2, 2}, {0.5, -2}, {0.5, 2}, {2, -2}, {2, 0}, {2, 2} }; ",
        "traj = ParametricPlot[Evaluate[Table[{x[t], y[t]} /. sol[z0], {z0, inits}]], {t, 0, 12}, ",
        "PlotRange -> { {-2.5, 2.5}, {-2.5, 2.5} }]; ",
        "Show[vf, pts, traj]"
      }]
    ]
  ], "Input"]
}]
      

The problems below focus on equilibrium computation, nullcline reasoning, and local classification. Assume smoothness/Lipschitz conditions where needed.

Problem 1 (Equilibria + classification): Consider \( \dot{x} = x - y \), \( \dot{y} = x + y - x^3 \). (a) Find all equilibria. (b) Classify each equilibrium using the Jacobian invariants \( \tau,\Delta,D \).

Solution: Equilibria satisfy \( x-y=0 \Rightarrow y=x \) and \( x+y-x^3=0 \Rightarrow x+x-x^3=0 \Rightarrow x(2-x^2)=0 \). Thus \( x^*=0 \) or \( x^*=\pm\sqrt2 \), with \( y^*=x^* \).

Jacobian:

\[ \mathbf{J}(x,y)= \begin{bmatrix} 1 & -1 \\ 1-3x^2 & 1 \end{bmatrix}. \]

At \( (0,0) \): \( \tau=2 \), \( \Delta = 1\cdot 1 - (-1)(1)=2 \), \( D=4-8=-4 < 0 \) so it is an unstable spiral. At \( (\pm\sqrt2,\pm\sqrt2) \): \( 1-3x^2 = 1-6=-5 \), hence \( \Delta = 1\cdot 1 - (-1)(-5)=1-5=-4 < 0 \), so each is a saddle.

Problem 2 (Non-intersection): Prove that trajectories of a planar autonomous system \( \dot{\mathbf{z} }=\mathbf{F}(\mathbf{z}) \) cannot intersect if \( \mathbf{F} \) is locally Lipschitz.

Solution: If two trajectories intersect at time \( t_0 \), then they share the same initial state \( \mathbf{z}(t_0)=\mathbf{z}_0 \). Both satisfy the same IVP; uniqueness implies the two solutions are identical on the common interval. Therefore distinct trajectories cannot cross. □

Problem 3 (Linear phase portrait): Consider the linear system \( \dot{x} = 2x + y \), \( \dot{y} = -3x - 4y \). (a) Classify the origin. (b) Determine whether solutions converge as \( t → \infty \).

Solution: The Jacobian is the constant matrix \( \mathbf{A}=\begin{bmatrix}2&1\\-3&-4\end{bmatrix} \). Compute invariants: \( \tau = 2 + (-4) = -2 \), \( \Delta = 2\cdot(-4) - 1\cdot(-3) = -8 + 3 = -5 < 0 \). Since \( \Delta < 0 \), the origin is a saddle, hence unstable. Because one eigenvalue is positive, some trajectories diverge as \( t → \infty \).

Problem 4 (Center vs. saddle in a nonlinear oscillator form): Consider \( \dot{x} = y \), \( \dot{y} = -x + x^3 \). (a) Find equilibria. (b) Classify each using the Jacobian. (c) Show the origin is a linear center.

Solution: Equilibria satisfy \( y=0 \) and \( -x + x^3=0 \Rightarrow x(x^2-1)=0 \), so \( (0,0) \), \( (1,0) \), \( (-1,0) \). Jacobian:

\[ \mathbf{J}(x,y)= \begin{bmatrix} 0 & 1 \\ -1+3x^2 & 0 \end{bmatrix}. \]

At \( (0,0) \): \( \tau=0 \), \( \Delta = 1 \), \( D = -4 < 0 \), hence purely imaginary eigenvalues (linear center). At \( (\pm 1,0) \): \( -1+3x^2=2 \) so \( \Delta = -2 < 0 \), hence saddles. Part (c) follows from the eigenvalues \( \lambda = \pm i \) at the origin for the linearization.

Problem 5 (Nullclines and direction): For the system in Section 6 with \( a=1 \), determine the sign of \( \dot{x} \) and \( \dot{y} \) in the region \( x > 0 \), \( y > 0 \), and far from nullclines (e.g., take a representative point \( (x,y)=(2,2) \)).

Solution: Evaluate: \( \dot{x}=x-x^3-y = 2-8-2=-8 < 0 \) and \( \dot{y}=x+y = 4 > 0 \). Thus trajectories locally move left and up at that point. Repeating this test in each region separated by nullclines yields a qualitative vector-field sketch without integration.

You can now analyze planar nonlinear dynamics by: (i) computing equilibria (intersections of nullclines), (ii) using existence/uniqueness to reason about trajectories and non-crossing, (iii) classifying equilibria locally via the Jacobian and the \( (\tau,\Delta,D) \) test, and (iv) assembling these ingredients into a phase portrait by combining nullcline geometry with representative integrated trajectories. This qualitative viewpoint is essential for understanding when linear intuition works and when it fails in nonlinear dynamics.