Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

By the end of this section you will be able to:

Interpret any 2D linear system $\dot{\mathbf{x}} = A\,\mathbf{x}$ as a velocity field on the plane and a flow of curves.
Recognise the five fundamental portraits — saddle, stable node, unstable node, stable spiral, unstable spiral — and the two boundary cases (center, degenerate) by looking at the picture.
Predict the portrait type from just two scalars: $T = \operatorname{tr}A$ and $D = \det A$ , via the discriminant $\Delta = T^2 - 4D$ .
Locate any 2×2 matrix as a single dot in the trace–determinant plane and read off its qualitative behaviour.
Implement the velocity field, integrate trajectories by hand-written RK4, and verify against a closed-form solution.
Batch-classify thousands of systems in one PyTorch call using a vectorised companion-matrix construction.

Why We Need Phase Portraits

"Forget the formula. Look at the picture." Henri Poincaré's decisive move in the 1880s — the move that founded modern dynamical systems.

For a single scalar ODE $\dot x = f(x)$ we already have a graphical habit: draw $f$ against $x$ , mark the zeros (the equilibria), and read off the direction of motion from the sign of $f$ . That picture often tells us everything we want to know without solving the ODE.

For a system of two ODEs $\dot x = f(x, y),\; \dot y = g(x, y)$ , the state lives on the $xy$ -plane. The single- variable trick no longer works because there are two variables and only one time. We need a richer picture — one where every point of the plane has an arrow telling us where the state will move next. That picture is the phase portrait.

The promise of this section

For 2×2 linear systems we will prove that the entire qualitative behaviour is determined by $T = \operatorname{tr}A$ and $D = \det A$ . A 4-number matrix collapses to two scalars, and those two scalars choose one of seven labelled regions on a single 2D map. Once you can read that map, you can sketch the portrait of any linear system in seconds.

State Space and the Velocity Field

A state of our system is a point $\mathbf{x} = (x, y) \in \mathbb{R}^2$ . The plane of all such points is the state space (also called phase space). For a linear system $\dot{\mathbf{x}} = A\mathbf{x}$ , the ODE does not give us where the state IS; it tells us how fast and in which direction it MOVES from wherever it happens to be:

\displaystyle \dot{\mathbf{x}} = A\,\mathbf{x} \quad\Longleftrightarrow\quad \mathbf{v}(\mathbf{x}) = A\,\mathbf{x}.

Think of $\mathbf{v}(\mathbf{x}) = A\mathbf{x}$ as wind blowing across the plane. At every point $\mathbf{x}$ , the wind direction and speed are given by the linear function $A\mathbf{x}$ . Drop a leaf at any starting position; the leaf's path is the trajectory through that seed, and the collection of all such trajectories — the whole "weather pattern" — is the phase portrait.

Three things to notice about linear vector fields

The origin is always a fixed point. Because $A \cdot \mathbf{0} = \mathbf{0}$ , the wind speed at the origin is zero. A leaf placed at the origin never moves.
The field is invariant under scaling. If $\mathbf{x}$ doubles, so does the arrow. The picture is self-similar: zooming in does not change the local geometry.
The eigenvectors are special directions. If $A \mathbf{v} = \lambda \mathbf{v}$ , then along the line through $\mathbf{v}$ the arrow is parallel to the line itself. Solutions that start on an eigenvector stay on it forever — sliding to or from the origin at exponential rate $\lambda$ .

The Five Fundamental Portraits

Even before we look at the algebra, the eye can already classify every linear 2D system into a handful of shapes. Each shape is ruled by the eigenvalues of $A$ :

Portrait	Eigenvalues	What the flow does
Saddle	λ₁ > 0, λ₂ < 0	Two real eigenvalues of opposite sign. Along v₁ trajectories fly outward; along v₂ they fall inward. Everywhere else they bend like hyperbolae. The origin is unstable.
Stable node	λ₁ < λ₂ < 0	Both eigenvalues negative. Every trajectory crashes into the origin, tangent to the SLOWER eigenvector (the one with the less-negative λ).
Unstable node	0 < λ₁ < λ₂	Both eigenvalues positive. Every trajectory shoots away from the origin, tangent to the SLOWER eigenvector.
Stable spiral	α ± iβ, α < 0	Complex-conjugate pair with negative real part. Trajectories rotate at angular frequency β while the radius shrinks like e^(αt). The origin is a stable focus.
Unstable spiral	α ± iβ, α > 0	Same shape rotation, but radius grows. The origin is an unstable focus.

Two further boundary cases sit between these five:

Boundary	Eigenvalues	Behaviour
Center	± iβ (pure imaginary)	Trajectories are perfect closed ellipses around the origin. Energy-like quantities are conserved exactly. Marginally stable.
Degenerate	one λ = 0	An entire line of fixed points (along the kernel of A). Off that line, motion is along the other eigendirection only.

A reading habit

Always look first at where the arrows go near the origin. That tiny region tells you the topological type. The far-field arrows are just a stretched version of the same picture — linear fields look the same at every scale.

From Matrix to Portrait: Trace and Determinant

The eigenvalues of a 2×2 matrix $A = \begin{pmatrix} a & b \\ c & d \end{pmatrix}$ satisfy the characteristic polynomial

\lambda^2 - T\,\lambda + D = 0,

where $T = \operatorname{tr}(A) = a + d$ and $D = \det(A) = ad - bc$ .

Vieta's formulas give us the eigenvalues without ever solving the quadratic:

\lambda_1 + \lambda_2 = T, \qquad \lambda_1 \lambda_2 = D.

From these two identities and the discriminant $\Delta = T^2 - 4D$ , we can decide every question we care about without ever computing the eigenvalues themselves:

Question about the portrait	Answer in (T, D)
Are the eigenvalues real or complex?	Δ > 0 ⇒ real distinct. Δ = 0 ⇒ repeated. Δ < 0 ⇒ complex conjugate.
Do the eigenvalues have the same sign?	Yes iff D > 0 (product is positive).
Is the origin attracting (stable)?	Yes iff Re(λ₁) < 0 and Re(λ₂) < 0 ⇔ T < 0 and D > 0.
Is the origin repelling (unstable)?	Yes iff at least one real part > 0 ⇔ T > 0 or D < 0.
Is it a saddle?	Yes iff D < 0 (eigenvalues opposite sign).
Is it a center?	Yes iff T = 0 and D > 0 (purely imaginary eigenvalues).

The single most useful fact in this chapter

A 2D linear system is stable at the origin iff

T < 0 \quad\text{and}\quad D > 0.

That single sentence is the gateway to control theory, circuit design, and most of mathematical biology. Whenever an engineer asks "is this linearisation stable?", they are checking the sign of two numbers.

The Discriminant Parabola

On the $(T, D)$ plane, the curve $\Delta = T^2 - 4D = 0$ — equivalently $D = T^2/4$ — is a parabola opening upwards. Above the parabola the discriminant is negative and the eigenvalues are complex (spirals and centers). Below the parabola the discriminant is positive and the eigenvalues are real distinct (nodes and saddles). The parabola itself is the boundary — repeated real eigenvalues — where nodes degenerate into a single eigendirection.

The other crucial curves are the two axes:

The positive D-axis (T = 0, D > 0) is the home of centers. Cross it and a stable spiral becomes an unstable spiral.
The T-axis (D = 0) is the boundary between saddles and nodes. Crossing it makes one eigenvalue change sign — the dramatic event of bifurcation in nonlinear theory.

Why this map matters far beyond linear ODEs

In Section 23.7 we will linearise nonlinear systems near each fixed point and read off the local picture from the same map. The (T, D) plane therefore IS the dictionary that translates algebra into geometry for nearly every model you will encounter in physics, biology, and engineering.

Worked Example: A Stable Spiral by Hand

Take the matrix

A = \begin{pmatrix} -\tfrac12 & 1 \\ -1 & -\tfrac12 \end{pmatrix}.

Before any computation, eyeball the structure: the diagonal entries are equal and negative, the off-diagonal entries are equal in magnitude and opposite in sign. That symmetry is the algebraic fingerprint of a rotation plus contraction. Let's verify.

Show the full hand-computation (try it yourself first)

Step 1 — Trace, determinant, discriminant. $T = -\tfrac12 + (-\tfrac12) = -1$ , $D = (-\tfrac12)(-\tfrac12) - (1)(-1) = \tfrac14 + 1 = \tfrac54$ , $\Delta = T^2 - 4D = 1 - 5 = -4 < 0$ . So eigenvalues are complex, T < 0 and D > 0 — by our table, a stable spiral.

Step 2 — Eigenvalues. $\lambda_{1,2} = \tfrac{T \pm \sqrt{\Delta}}{2} = \tfrac{-1 \pm 2i}{2} = -\tfrac12 \pm i$ . Decay rate $\alpha = -\tfrac12$ , angular frequency $\beta = 1$ .

Step 3 — Complex eigenvector for $\lambda = -\tfrac12 + i$ . Solve $(A - \lambda I)\mathbf{v} = 0$ : row 1 reads $-i\,v_1 + v_2 = 0$ , so $v_2 = i\,v_1$ . Take $\mathbf{v} = (1,\, i)^{\top}$ .

Step 4 — Real basis of solutions. Compute $e^{\lambda t}\mathbf{v}$ :

e^{(-1/2 + i)t}\,\begin{pmatrix}1\\ i\end{pmatrix} = e^{-t/2}(\cos t + i \sin t)\,\begin{pmatrix}1\\ i\end{pmatrix} = e^{-t/2}\!\begin{pmatrix}\cos t + i\sin t \\ -\sin t + i\cos t\end{pmatrix}.

Real and imaginary parts give two independent real solutions:

\mathbf{x}_R(t) = e^{-t/2}\!\begin{pmatrix}\cos t \\ -\sin t\end{pmatrix},\quad \mathbf{x}_I(t) = e^{-t/2}\!\begin{pmatrix}\sin t \\ \phantom{-}\cos t\end{pmatrix}.

Step 5 — Match the seed $\mathbf{x}(0) = (2, 0)^\top$ . Write $\mathbf{x}(t) = C_1\,\mathbf{x}_R(t) + C_2\,\mathbf{x}_I(t)$ . At $t = 0$ : $\mathbf{x}_R(0) = (1, 0)$ , $\mathbf{x}_I(0) = (0, 1)$ , so $C_1 = 2$ , $C_2 = 0$ . Therefore

\boxed{\;\mathbf{x}(t) = e^{-t/2}\begin{pmatrix} 2\cos t \\ -2\sin t \end{pmatrix}.\;}

Step 6 — Numerical samples.

t	x(t)	y(t)	Radius e^(-t/2) · 2
0	2.000	0.000	2.000
π/2	0.000	-0.910	0.910
π	-0.416	0.000	0.416
3π/2	0.000	0.189	0.189
2π	0.086	0.000	0.086

After one full revolution $(t = 2\pi)$ the radius has shrunk by a factor of $e^{-\pi} \approx 0.0432$ — a ~23-times decay per loop. That number, the logarithmic decrement, is exactly how engineers measure damping from oscilloscope traces.

Interactive Phase Portrait Explorer

Below, every slider edits one entry of the 2×2 matrix $A$ . The vector field updates instantly, eigenvalues are recomputed, the classification label changes, and seed trajectories are re-integrated by RK4. Click anywhere in the canvas to drop a new seed — useful for studying how nearby orbits differ.

Loading interactive phase portrait…

Things to try

Move slider b from -2 to +2 with all others fixed. Watch the spiral first tighten, then become a node, then become a saddle as the determinant crosses zero.
Set a = d = 0 and b = -c. You get exact centers — closed elliptical orbits at every radius.
Set a = d = -1 and b = c = 0. Repeated real eigenvalue, but every direction is an eigenvector — a star node.

The Trace–Determinant Plane

The previous explorer lets you change the matrix and watch the portrait. The next visualiser inverts that: it lets you navigate the qualitative-behaviour map directly. Drag the red dot over the shaded regions and the canvas on the right shows the corresponding portrait. The yellow parabola $D = T^2/4$ is the discriminant boundary.

Loading trace–determinant classifier…

Crossing boundaries = bifurcation

Walking the dot from the stable-spiral region across the positive D-axis (T = 0) takes a damped oscillator and turns it into a pure oscillator and then into an unstable spiral. That three-step crossing — stable → marginal → unstable — is the Hopf bifurcation, one of the universal routes a dynamical system can take into self-sustained oscillation (think heartbeat onset, laser threshold, predator-prey cycles).

Python: Drawing the Velocity Field

Three lines of NumPy — meshgrid, the matrix–vector product, and $\texttt{quiver}$ — produce every phase- portrait poster you have ever seen. The function below is the skeleton you can reuse for any 2×2 linear system in this course.

From matrix A to a quiver plot in eight expressive lines

🐍velocity_field.py

Explanation(11)

Code(56)

1Import NumPy and matplotlib

NumPy gives us vectorized arithmetic on the entire grid in one shot — no Python for-loops over points. matplotlib's quiver routine draws one arrow per grid cell, which is exactly what a velocity field is.

4Function signature

The function takes the 2×2 system matrix A and the plotting window. We will reuse it for every linear system in this section. Everything you see on a phase-portrait poster — saddles, spirals, nodes — is generated by this exact six-line core.

20Build a regular grid of base points

np.linspace(-limit, limit, density) is a 1-D array of evenly spaced sample positions. np.meshgrid expands that 1-D array into a 2-D coordinate grid: X0[i,j] holds the x-component of cell (i,j), and X1[i,j] holds the y-component. After this line, every cell of the picture has a known (x, y) location.

24Apply A to every grid point at once

Because A is 2×2, the matrix product A·x simplifies to two scalar formulas: v₀ = A[0,0]·x + A[0,1]·y and v₁ = A[1,0]·x + A[1,1]·y. We write those formulas using NumPy's elementwise multiplication, so a single line evaluates the velocity at all density² grid points simultaneously. This is the entire computational content of a phase portrait — one matrix–vector product, broadcast.

28Open a square axes if none was given

Phase portraits MUST be drawn with equal aspect ratio. Otherwise a circle would look like an ellipse and the eye would mis-read curvature. We default to a 5-inch square.

32Quiver: one arrow per grid point

ax.quiver(X, Y, U, V, C) draws an arrow at each (X[i,j], Y[i,j]) pointing in the direction (U[i,j], V[i,j]) and coloured by C[i,j]. Here we pass np.hypot(V0, V1) as the color channel, so faster regions appear yellow and slow regions appear deep purple — the velocity magnitude becomes a heat-map underneath the arrows.

35pivot='middle', scale=30

pivot='middle' centers each arrow on its base point so the cell looks balanced. scale controls how many world units one inch of arrow represents — bigger scale = shorter arrows. We pick 30 so the longest arrows fit inside one grid cell.

38Equal aspect, axes through origin

set_aspect('equal') is non-negotiable for phase pictures. axhline/axvline at 0 mark the coordinate axes — useful because eigenvectors of diagonal matrices lie exactly on them.

47Choose a stable-spiral matrix

A = [[-0.5, 1], [-1, -0.5]]. We will hand-classify it in a moment. Trace = -0.5 + (-0.5) = -1 (negative → contraction overall). Determinant = (-0.5)·(-0.5) - (1)·(-1) = 0.25 + 1 = 1.25 (positive → both eigenvalues on the same side). Discriminant T² - 4D = 1 - 5 = -4 (negative → complex eigenvalues). Negative-trace + complex-eigenvalues = stable spiral.

51Verify with NumPy

np.trace(A) returns -1. np.linalg.det(A) returns 1.25. np.linalg.eigvals(A) returns the complex pair -0.5 + 1.0j and -0.5 - 1.0j. The real part -0.5 is the radial decay rate, the imaginary part 1.0 is the angular frequency of the spiral. So solutions spiral inward and complete one full revolution every 2π ≈ 6.28 time units, while the radius shrinks by a factor of e^(-0.5·2π) ≈ 0.043 per revolution.

55Render the field

Call our function and add a title. The output is a square of inward-curling arrows that hint at the spiral structure even before any trajectory is drawn. The arrows alone tell the eye the answer; the trajectories will only confirm it.

45 lines without explanation

1import numpy as np
2import matplotlib.pyplot as plt
3
4def draw_velocity_field(A, limit=4.0, density=21, ax=None):
5    """
6    Draw the velocity field  v(x) = A x  on a square window.
7
8    Parameters
9    ----------
10    A : np.ndarray, shape (2, 2)
11        The coefficient matrix of the linear system  x' = A x.
12    limit : float
13        Half-width of the window. Plot covers [-limit, limit]^2.
14    density : int
15        Number of arrows per side. Total arrows = density^2.
16    ax : matplotlib axes (optional). A new one is made if None.
17
18    The arrow length is proportional to the local speed |A x|.
19    """
20    # 1. Build a regular grid of base points (x0, x1).
21    xs = np.linspace(-limit, limit, density)
22    X0, X1 = np.meshgrid(xs, xs)
23
24    # 2. Apply A elementwise: v = A x  computed at every grid point.
25    V0 = A[0, 0] * X0 + A[0, 1] * X1
26    V1 = A[1, 0] * X0 + A[1, 1] * X1
27
28    if ax is None:
29        _, ax = plt.subplots(figsize=(5, 5))
30
31    # 3. quiver draws an arrow at each (X0, X1) pointing along (V0, V1).
32    ax.quiver(X0, X1, V0, V1, np.hypot(V0, V1),
33              cmap="viridis", pivot="middle",
34              scale=30, width=0.003)
35
36    ax.set_aspect("equal")
37    ax.set_xlim(-limit, limit)
38    ax.set_ylim(-limit, limit)
39    ax.set_xlabel("x")
40    ax.set_ylabel("y")
41    ax.axhline(0, color="gray", lw=0.5)
42    ax.axvline(0, color="gray", lw=0.5)
43    return ax
44
45
46# Stable spiral: trace = -1 (stable), det = 1.25 > T^2/4 (complex)
47A = np.array([[-0.5,  1.0],
48              [-1.0, -0.5]])
49
50print("trace =", np.trace(A))                 # -1.0
51print("det   =", np.linalg.det(A))            # 1.25
52print("eigs  =", np.linalg.eigvals(A))        # [-0.5+1.j  -0.5-1.j]
53
54ax = draw_velocity_field(A)
55ax.set_title("Stable spiral velocity field")
56plt.show()

Python: Integrating Trajectories with RK4

Arrows alone are not enough — students want to see a leaf ride the wind. We add a Runge–Kutta-4 integrator that walks a seed point forward and backward in time, then verify the result against the closed-form solution we derived by hand above.

Hand-written RK4 produces every trajectory in the explorer

🐍rk4_trajectories.py

Explanation(11)

Code(54)

1Just NumPy

We do not need SciPy here. Writing RK4 by hand makes every assumption explicit, and it is the single most important ODE integrator to internalize: solve_ivp default in SciPy and ode45 in MATLAB are both adaptive-step versions of this same scheme.

3RK4 single step

Classical Runge–Kutta evaluates the right-hand side at four carefully chosen points (start, two midpoints, end) and combines them with weights (1, 2, 2, 1)/6. Local truncation error is O(h⁵); global error O(h⁴). Three orders of magnitude more accurate than Euler at the same step size — well worth the extra three function calls.

11Wrap the system as f(t, x) = A @ x

Our ODE is autonomous (no explicit time dependence) and linear (no x² or sin x). f(t, x) just returns the matrix product A · x using NumPy's @ operator. The integrator treats f as a black box, so the same RK4 routine would work for nonlinear systems too — we just gave it a linear one.

13Forward time-stepping loop

Allocate the output array, set the first row to x0, then march. At each step we replace x with the RK4 result and write it to row n+1. The accumulator fwd holds the entire forward trajectory, one row per time-stamp.

21Backward by flipping the sign of dt

A linear ODE is reversible: if x(t) is a solution, then x(-t) is the solution of the same system with starting point x(0) and run backwards. We integrate by feeding -dt into the same RK4 routine — no separate code path. Backward orbits show where points CAME FROM, which is essential for visualizing saddles whose unstable manifold is invisible if you only run forward.

28Concatenate and reverse

bwd[::-1] reverses the backward array so its earliest time is at index 0. We then drop the duplicate seed at the boundary and stack. The output is the COMPLETE orbit through x0 — past and future — in chronological order.

32Build the stable-spiral system

Same matrix as the velocity-field figure. Reusing it means the trajectory plot can be overlaid on top of the same arrows for the final figure students see in textbooks.

36Five carefully chosen seeds

We spread seeds at radius ~2.5 in five different directions. For a spiral this is enough to suggest the whole flow pattern; for a saddle you would also want seeds near the stable and unstable eigendirections. Choosing seeds well is half the art of drawing a phase portrait.

44List comprehension over seeds

orbits is now a list of five (N, 2) arrays. Each one is a complete bidirectional trajectory. We can plot them all with a single matplotlib loop.

49Cross-check against the closed form

For seed (2, 0) the analytic solution is x(t) = e^(-0.5 t) · (2 cos t, -2 sin t). At t = π that gives e^(-π/2)·(-2, 0) ≈ (-0.4158, 0). The RK4 sample at the corresponding index should match this to ~10⁻⁸. If it doesn't, the integrator has drifted and we need a smaller dt. Doing this check on EVERY new integrator is a non-negotiable habit.

53Index into the bidirectional orbit

The bidirectional orbit has length 2·Nf + 1 = 2001 with the seed at the middle (index 1000). Going forward, time = index_offset · dt, so t = π corresponds to offset ≈ π/0.02 ≈ 157. Hence index 1000 + 157 = 1157 is the closest sample to t = π — RK4 returns approximately (-0.4158, 0).

43 lines without explanation

1import numpy as np
2
3def rk4_step(f, y, t, h):
4    """One classical Runge–Kutta-4 step for y' = f(t, y)."""
5    k1 = f(t,         y)
6    k2 = f(t + h/2,   y + h * k1 / 2)
7    k3 = f(t + h/2,   y + h * k2 / 2)
8    k4 = f(t + h,     y + h * k3)
9    return y + h * (k1 + 2*k2 + 2*k3 + k4) / 6
10
11
12def integrate_orbit(A, x0, t_end=20.0, dt=0.02):
13    """Integrate  x'(t) = A x(t)  forward AND backward from x0."""
14    f  = lambda t, x: A @ x
15    Nf = int(t_end / dt)
16    fwd = np.zeros((Nf + 1, 2))
17    fwd[0] = x0
18    x = x0.astype(float)
19    for n in range(Nf):
20        x = rk4_step(f, x, n * dt, dt)
21        fwd[n + 1] = x
22
23    # backward by replacing dt with -dt
24    bwd = np.zeros((Nf + 1, 2))
25    bwd[0] = x0
26    x = x0.astype(float)
27    for n in range(Nf):
28        x = rk4_step(f, x, -n * dt, -dt)
29        bwd[n + 1] = x
30
31    return np.vstack([bwd[::-1], fwd[1:]])   # full bidirectional orbit
32
33
34# The same stable-spiral matrix as before
35A = np.array([[-0.5,  1.0],
36              [-1.0, -0.5]])
37
38# Pick five seeds and orbit them all
39seeds = np.array([
40    [ 2.0,  0.0],
41    [-2.0,  0.0],
42    [ 0.0,  2.5],
43    [-1.5, -2.5],
44    [ 2.5, -1.5],
45])
46orbits = [integrate_orbit(A, s) for s in seeds]
47
48# Check the analytic prediction for seed (2, 0)
49# x(t) = e^(-0.5 t) (2 cos t, -2 sin t)
50t_check = np.pi
51predicted = np.exp(-0.5 * t_check) * np.array([2 * np.cos(t_check),
52                                              -2 * np.sin(t_check)])
53print("at t = pi  predicted =", predicted)            # ~ [-0.4158, 0.0]
54print("at t = pi  RK4       =", orbits[0][1000 + 157])  # index ≈ t_end/dt at t=pi

Pure NumPy is fine here

For a 2D linear system the right-hand side is one matrix–vector multiply. We do NOT need SciPy's adaptive solvers. The point of writing RK4 by hand is to expose the FOUR-stage structure that gives the method its O(h⁴) accuracy. After this section you should be able to integrate any linear system to ~10⁻⁸ precision with twenty lines of NumPy.

PyTorch: Batched Classification of Many Systems

A single 2×2 system is easy. But what if you have ten thousand? Imagine sweeping over candidate controller gains in a robotics loop, or scanning a (T, D) grid to colour the bifurcation diagram of a model. We want to classify them all at once — and PyTorch is the natural language for that. Eigenvalue decomposition broadcasts over leading batch dimensions automatically.

One torch.linalg.eigvals call classifies a 41×41 portrait grid

🐍batched_classify.py

Explanation(16)

Code(71)

1Import PyTorch

We use PyTorch for two reasons. First, torch.linalg.eigvals supports BATCHES of matrices in one call, so we can analyse thousands of systems at once. Second, the same code runs on GPU with a single .cuda(). For a parameter sweep — "which (T, D) points give stable spirals?" — vectorization on GPU turns minutes into milliseconds.

3@torch.no_grad()

We are doing pure forward analysis — no training, no gradients. Disabling autograd skips graph construction and saves both memory and time. If we later want to OPTIMIZE T or D (say, design a circuit whose damping is exactly some target), we would remove this decorator.

21Flatten the input grids

Whether the caller passes 1-D arrays or 2-D meshes, we treat every point uniformly by flattening to (B,). B is the batch dimension — the integer count of (T, D) samples to classify simultaneously.

25Allocate the batch of companion matrices

We construct B copies of the 2×2 companion matrix A = [[0, -D], [1, T]]. Its characteristic polynomial is λ² - T·λ + D, so trace(A) = T and det(A) = D by construction. Every (T, D) pair on the plane corresponds to one canonical matrix here. Other matrices with the same (T, D) are SIMILAR to this one and have the SAME phase-portrait classification.

26Fill A in vectorized form

We assign to slices of the (B, 2, 2) tensor instead of looping. A[:, 0, 1] = -D sets the top-right entry of every matrix to its corresponding -D in one operation. PyTorch dispatches this to a single optimized kernel.

31Eigenvalues of every matrix in one call

torch.linalg.eigvals(A) where A has shape (B, 2, 2) returns a (B, 2) complex tensor — the two eigenvalues of each matrix. Internally LAPACK is called once per matrix, but we never write a Python loop. For B = 1681 (41×41 grid), this is ~1 ms on a modern CPU and submillisecond on GPU.

32Discriminant

T² - 4D is the discriminant of λ² - T·λ + D = 0. Its sign decides whether eigenvalues are real (disc > 0), complex conjugate (disc < 0), or repeated (disc = 0). We will USE this scalar test instead of inspecting eigs.imag because it is numerically more stable near the parabola boundary.

33Boolean mask for complex eigenvalues

is_complex is a Boolean tensor of shape (B,). Each entry is True iff the discriminant is significantly negative. We use a -1e-9 threshold instead of 0 so floating-point noise on the parabola itself does not flip classification spuriously.

35Allocate the integer label array

kind starts as all-6 (degenerate). We then overwrite specific entries with the correct label based on Boolean masks. This is the canonical PyTorch idiom for vectorized branching — no Python if/else loop over batch entries.

38Saddle mask

If det(A) < 0, the two eigenvalues are real and of OPPOSITE SIGN — by Vieta's formulas, det = λ₁·λ₂ < 0 forces one root negative and one positive. That is the saddle. We tag those entries with kind = 0.

41Real-distinct positive-det subset

real_pos_det = (D > 0) AND (not complex). For these, both eigenvalues are real with the same sign. The sign is determined by trace T: positive trace = positive sum = both positive (unstable node); negative trace = both negative (stable node). Two more masked assignments tag the nodes.

46Spirals and centers

When eigenvalues are complex (α ± βi) the REAL PART α is exactly T/2 (because sum of conjugate roots = 2 Re = T). Sign of T determines whether the spiral winds inward or outward; T = 0 means α = 0 — neutral, neither expanding nor shrinking — which is the center, the famous case of perfectly closed elliptical orbits.

52Return eigenvalues and labels

We return BOTH the raw eigenvalue tensor and the integer label tensor. The labels are what a phase-portrait drawing routine consumes; the eigenvalues are kept around for the legend (so we can write λ = -0.5 ± 1.0i next to a stable-spiral point).

56Build a 41×41 sample grid

torch.meshgrid with indexing='xy' returns Ts[i,j] = j-th value of the T-axis and Ds[i,j] = i-th value of the D-axis — the standard matplotlib convention. After flattening this is 1681 (T, D) points covering the [-3, 3]² square.

62Count each kind

A quick sanity check that our classifier covers the whole square. (kind == k).sum() counts entries of a particular type. The total must equal kind.numel() = 1681. If you see a surprise — say, an unexpected number of degenerate points — it is almost always a tolerance-threshold issue near the axes.

71Probe the center of the grid

Index 20·41 + 20 corresponds to (T, D) = (0, 0). That is a doubly-degenerate point (zero matrix); eigenvalues are exactly (0, 0). PyTorch returns 0+0j twice. The function reports kind = 6 (degenerate) — the only point in the plane where every direction is invariant and the "flow" is just standing still.

55 lines without explanation

1import torch
2
3@torch.no_grad()
4def classify_grid(t_grid, d_grid):
5    """
6    For every (T, D) in the meshgrid, build the companion matrix
7        A = [[0, -D], [1, T]]
8    which has trace T and determinant D, take its matrix exponential
9    over one time unit, and return:
10
11      - eigenvalues  (B, 2) complex
12      - kind         (B,)    integer label
13                              0 saddle
14                              1 stable node
15                              2 unstable node
16                              3 stable spiral
17                              4 unstable spiral
18                              5 center
19                              6 degenerate
20    """
21    T = torch.as_tensor(t_grid, dtype=torch.float64).flatten()
22    D = torch.as_tensor(d_grid, dtype=torch.float64).flatten()
23    B = T.shape[0]
24
25    # Stack a batch of companion matrices
26    A = torch.zeros(B, 2, 2, dtype=torch.float64)
27    A[:, 0, 1] = -D
28    A[:, 1, 0] = 1.0
29    A[:, 1, 1] = T
30
31    # Eigenvalues (complex) — used only for labelling
32    eigs = torch.linalg.eigvals(A)              # (B, 2) complex
33    disc = (T * T - 4.0 * D)
34    is_complex = disc < -1e-9
35
36    kind = torch.full((B,), 6, dtype=torch.long)  # default degenerate
37
38    # Saddle: D < 0  (one positive, one negative eigenvalue)
39    kind[D < -1e-6] = 0
40
41    # Stable / unstable node: D > 0, disc >= 0
42    real_pos_det = (D > 1e-6) & (~is_complex)
43    kind[real_pos_det & (T < 0)] = 1
44    kind[real_pos_det & (T > 0)] = 2
45
46    # Stable / unstable spiral: complex eigenvalues, T != 0
47    kind[is_complex & (T < -1e-6)] = 3
48    kind[is_complex & (T >  1e-6)] = 4
49    kind[is_complex & (T.abs() < 1e-6)] = 5      # center
50
51    return eigs, kind
52
53
54# Sample a 41x41 (T, D) grid
55Ts, Ds = torch.meshgrid(
56    torch.linspace(-3.0, 3.0, 41, dtype=torch.float64),
57    torch.linspace(-3.0, 3.0, 41, dtype=torch.float64),
58    indexing="xy",
59)
60eigs, kind = classify_grid(Ts, Ds)
61print("Total points:", kind.numel())
62print("Saddles      :", (kind == 0).sum().item())
63print("Stable nodes :", (kind == 1).sum().item())
64print("Unstable     :", (kind == 2).sum().item())
65print("Stable spir. :", (kind == 3).sum().item())
66print("Unstable spir:", (kind == 4).sum().item())
67print("Centers      :", (kind == 5).sum().item())
68
69# Sample one specific point and verify
70i = 20 * 41 + 20  # row=20, col=20  → (T, D) = (0, 0)
71print("eig at (0,0):", eigs[i])

Why "companion matrix"?

Every 2×2 matrix is similar (via change of basis) to its companion matrix $A = \begin{pmatrix} 0 & -D \\ 1 & T \end{pmatrix}$ — the unique matrix whose entries are read directly off the characteristic polynomial $\lambda^2 - T\lambda + D$ . Two matrices with the same trace AND the same determinant give phase portraits of the same topological type, so for classification we lose nothing by sampling the companion family.

Where This Shows Up

Domain	What x = (x₁, x₂) is	Why the portrait matters
Mass-spring-damper	(position, velocity)	Stable spiral = underdamped ringing. Stable node = overdamped slump. Repeated root on the parabola = critical damping (the engineer's target).
RLC circuit	(capacitor charge, loop current)	Exactly the same picture as the mechanical oscillator. The trace = -R/L is the damping rate; det = 1/(LC) is the natural frequency squared.
Predator–prey (linearised)	(prey deviation, predator deviation)	Lotka–Volterra near coexistence gives pure imaginary eigenvalues — a center. The neutrally closed orbits are why both species oscillate indefinitely in the textbook model.
Chemical reactor	(temperature deviation, concentration deviation)	Saddles indicate runaway risk. Stable spirals indicate damped chemical oscillations. Hopf crossings predict the onset of sustained chemical clocks.
Neural ODEs	(hidden state coordinates)	A learnable A trained to be stable (T<0, D>0) gives memory cells that decay. Pushing onto the imaginary axis gives oscillator-like memory used in continuous RNNs.

Common Pitfalls

Forgetting unequal aspect ratio. A non-square axes box turns circles into ellipses and lies about spirals versus centers. Always $\texttt{set\_aspect("equal")}$ .
Reading the eigenvector direction off the wrong row. From $(A - \lambda I)\mathbf{v} = 0$ , if row 1 is degenerate, USE row 2. Otherwise you get the zero vector.
Confusing "stable node" tangency rules. Trajectories enter a stable node tangent to the slowest eigenvector (least-negative eigenvalue) — not the fastest. The fastest direction is the one along which trajectories first shed radius.
Calling a center a spiral. Numerical integrators on a center system will drift over many revolutions because floating-point noise injects a small spurious $\alpha \neq 0$ . To check if you have a true center, look at the matrix (T = 0?), not the integrator output (which always lies a little).
Drawing only forward trajectories near a saddle. The unstable manifold is the one going OUT; the stable manifold is the one coming IN. Without backward-time integration you will miss the stable manifold entirely.

Summary

The phase portrait of $\dot{\mathbf{x}} = A\mathbf{x}$ is the velocity field $A\mathbf{x}$ plus the integral curves it produces.
All 2D linear portraits fall into five generic types — saddle, two kinds of node, two kinds of spiral — and two boundary cases (center, degenerate).
The classification is governed by two scalars, $T = \operatorname{tr}A$ and $D = \det A$ . The discriminant parabola $\Delta = T^2 - 4D = 0$ separates the real-eigenvalue regions (below) from the complex ones (above).
$T < 0 \text{ and } D > 0$ is the gateway condition for stability — the half-line every control engineer aims at.
Computationally, the velocity field is one broadcast $A\mathbf{x}$ ; trajectories follow from any first-order ODE integrator (we used RK4); batched classification is one PyTorch eigenvalue call.

Looking forward: In Section 23.4 we will read off the trajectory formula when eigenvalues are real and distinct, and in Section 23.5 we will handle the complex-eigenvalue spirals algebraically. Section 23.7 then explains why these linear portraits remain useful for nonlinear systems near each fixed point.