BrevFeed
Sign in

9. Systems of linear equations and the phase plane

📖 Reading · 11 min
💡 Every code box below is live — edit it and hit Run.

Every equation so far has involved a single unknown function. Many real systems — predator and prey populations, coupled springs, circuits with two components — need two unknown functions evolving together, each one's rate of change depending on both. This closing modeling lesson packages such systems into a single vector equation, solves it, and visualizes the solutions with a picture that needs no formula at all: the phase plane.

Prerequisite note. This lesson leans on a small amount of linear algebra that the calculus series does not develop: 2\times2 matrices as objects that act on vectors, the determinant of a 2\times2 matrix, and above all eigenvalues and eigenvectors — the scalars \lambda and nonzero vectors \vec v satisfying A\vec v=\lambda\vec v. §9 gives you vectors, but not matrices; everything matrix-shaped is introduced inline below and the lesson is self-contained, but if these ideas are new, an hour with any introduction to linear algebra first will make this lesson substantially easier — the companion linear-algebra course's §16 (Foundations) and §19.0 (characteristic polynomial and eigenvalues) cover exactly this. The same applies to §14.2 and §14.5.

Systems as a single vector equation

A linear system of two first-order equations,

x'=ax+by,\qquad y'=cx+dy

packages into a single vector equation using a matrix — a 2\times2 block of coefficients whose action on a vector reproduces both right-hand sides at once:

\vec X'=A\vec X,\qquad\vec X=\begin{pmatrix}x\\y\end{pmatrix},\quad A=\begin{pmatrix}a&b\\c&d\end{pmatrix}

This is y'=ky (§13.0) one dimension up — a scalar equation replaced by a vector one, and (as the solution method below shows) the scalar solution y=Ce^{kt} generalizes almost unchanged.

Solving via eigenvalues and eigenvectors

Guess a solution of the form \vec X=e^{\lambda t}\vec v for a constant vector \vec v and scalar \lambda — the direct vector analogue of §13.5's y=e^{rx} guess. Substituting:

\lambda e^{\lambda t}\vec v=Ae^{\lambda t}\vec v\ \Longrightarrow\ A\vec v=\lambda\vec v

This is exactly the eigenvalue equation: \lambda must be an eigenvalue of A, and \vec v the corresponding eigenvector. Finding \lambda solves \det(A-\lambda I)=0 — a polynomial in \lambda called the characteristic equation, the direct matrix generalization of §13.5's characteristic equation ar^2+br+c=0 (in fact, for a 2\times2 system, expanding \det(A-\lambda I)=0 produces exactly a quadratic in \lambda, with the same three real-distinct / repeated / complex root cases as §13.5).

For two distinct real eigenvalues \lambda_1,\lambda_2 with eigenvectors \vec v_1,\vec v_2:

\vec X(t)=C_1e^{\lambda_1t}\vec v_1+C_2e^{\lambda_2t}\vec v_2

Structurally identical to §13.5's Case 1 — two exponentials, now scaled by vectors instead of by 1, with two constants pinned down by an initial vector condition \vec X(0)=\vec X_0.

The phase plane

Rather than plotting x(t) and y(t) separately against t, the phase plane plots the trajectory \big(x(t),y(t)\big) directly in the xy-plane, with t as an implicit parameter — exactly §6.0's parametric curve, traced by the solution as time advances. The system's right side, \langle ax+by,\,cx+dy\rangle, is itself a vector field (§12.0): at every point (x,y), it gives the instantaneous velocity \langle x',y'\rangle a trajectory passing through that point must have — so a phase-plane sketch is really a slope field (§13.0) one dimension up, now showing a full velocity vector rather than a single slope.

The origin is always an equilibrium (since A\vec0=\vec0, a trajectory starting there never moves), and the sign pattern of the eigenvalues classifies its behavior:

  • Both real, same sign: a node — trajectories flow straight into (both negative, stable) or straight out of (both positive, unstable) the origin along the eigenvector directions.
  • Both real, opposite signs: a saddle — trajectories approach along one eigenvector direction and flee along the other, the dynamical-systems sibling of §10.6's saddle point.
  • Complex, nonzero real part: a spiral — trajectories wind around the origin while growing or shrinking, from Euler's formula (§8.3) turning the complex exponential into a rotating, scaling motion.
  • Purely imaginary: a center — trajectories orbit the origin in closed loops forever, neither approaching nor receding.

Doing it in Python

Solving the system x'=x+y, y'=4x+y, x(0)=1, y(0)=0 via eigenvalues and eigenvectors:

import sympy as sp

A = sp.Matrix([[1, 1], [4, 1]])
eigen_data = A.eigenvects()

for eigenvalue, multiplicity, eigenvectors in eigen_data:
    print(f"eigenvalue = {eigenvalue}, eigenvector = {eigenvectors[0].T}")

Confirming the classification: opposite-sign eigenvalues mean this system has a saddle point at the origin, and solving the full initial value problem directly:

import sympy as sp

t = sp.Symbol('t')
x, y = sp.Function('x'), sp.Function('y')

system = [
    sp.Eq(x(t).diff(t), x(t) + y(t)),
    sp.Eq(y(t).diff(t), 4*x(t) + y(t))
]
solution = sp.dsolve(system, ics={x(0): 1, y(0): 0})
for eq in solution:
    print(eq)

Tracing the phase-plane trajectory numerically — confirming it moves toward the unstable direction (eigenvalue +3) as t grows:

import math

def xt(t):
    return math.exp(3*t)/2 + math.exp(-t)/2

def yt(t):
    return math.exp(3*t) - math.exp(-t)

print(f"{'t':>5} {'x(t)':>10} {'y(t)':>10} {'y/x':>8}")
for t in (0, 0.5, 1, 1.5, 2):
    x_val, y_val = xt(t), yt(t)
    print(f"{t:>5.1f} {x_val:>10.4f} {y_val:>10.4f} {y_val/x_val:>8.4f}")
print("\ny/x approaches 2 -- the trajectory aligns with the (1,2) eigenvector")

Worked example

Classify the equilibrium at the origin for the system x'=x+y, y'=4x+y, and solve the initial value problem x(0)=1, y(0)=0.

A=\begin{pmatrix}1&1\\4&1\end{pmatrix}

Characteristic equation:

\det(A-\lambda I)=(1-\lambda)^2-4=0\ \Longrightarrow\ 1-\lambda=\pm2\ \Longrightarrow\ \lambda=-1,3

Opposite signs (-1 and 3) \Longrightarrow $\boxed{\text{saddle point}}$.

Eigenvectors: for \lambda=3: $(1-3)v_1+v_2=0\Rightarrow v_2=2v_1$, giving \vec v_1=\langle1,2\rangle. For \lambda=-1: (1-(-1))v_1+v_2=0\Rightarrow v_2=-2v_1, giving \vec v_2=\langle1,-2\rangle.

\vec X(t)=C_1e^{3t}\langle1,2\rangle+C_2e^{-t}\langle1,-2\rangle

Apply \vec X(0)=\langle1,0\rangle: C_1+C_2=1 and 2C_1-2C_2=0\Rightarrow C_1=C_2. Combined with C_1+C_2=1: C_1=C_2=\frac12.

\boxed{x(t)=\frac12e^{3t}+\frac12e^{-t},\qquad y(t)=e^{3t}-e^{-t}}

Sanity check. x(0)=\frac12+\frac12=1 ✓, y(0)=1-1=0 ✓. As t\to\infty, the e^{3t} terms dominate both components, and \frac{y}x\to\frac{e^{3t}}{\frac12e^{3t}}=2 — the trajectory bends to align with the unstable eigenvector direction \langle1,2\rangle, exactly the saddle-point behavior described in the concept section: initial data eventually gets swept along the growing eigendirection, regardless of where it started (as long as it has any component along that direction at all). ✓

Your turn

1. Find the eigenvalues of A=\begin{pmatrix}-2&0\\0&-3\end{pmatrix} and classify the equilibrium at the origin (this matrix is already diagonal, so the eigenvalues can be read off directly).

2. Find the eigenvalues of A=\begin{pmatrix}0&1\\-1&0\end{pmatrix} and classify the equilibrium.

3. True or false: a saddle point's trajectories, given enough time, always end up approaching the origin.

Solutions

1. For a diagonal matrix, the eigenvalues are exactly the diagonal entries: \lambda=-2,-3.

Both real and negative (same sign) \Longrightarrow \boxed{\text{a stable node}} — every trajectory flows straight into the origin as t\to\infty.

2. \det(A-\lambda I)=\lambda^2+1=0\Rightarrow\lambda=\pm i — purely imaginary.

\boxed{\text{a center}}

— trajectories orbit the origin in closed loops (in fact, this particular system is exactly x'=y,\,y'=-x, whose solutions are circles, directly related to §9.4's circular parametrization).

3. False. A saddle point has one eigenvalue positive (repelling) and one negative (attracting) — only trajectories starting exactly on the attracting eigenvector's line approach the origin; every other starting point eventually gets swept away along the repelling eigenvector direction, exactly the behavior the worked example's y/x\to2 limit demonstrated. A saddle is called a saddle precisely because it attracts along one direction and repels along another, never both.

Check yourself in code

Find the eigenvalues and eigenvectors of A=\begin{pmatrix}1&1\\4&1\end{pmatrix}.

Print exactly this:

eigenvalue = -1, eigenvector = [-1/2, 1]
eigenvalue = 3, eigenvector = [1/2, 1]

import sympy as sp

A = sp.Matrix([[1, 1], [4, 1]])
eigen_data = A.eigenvects()

for eigenvalue, multiplicity, eigenvectors in eigen_data:
    print(f"eigenvalue = {eigenvalue}, eigenvector = ...")

import sympy as sp

A = sp.Matrix([[1, 1], [4, 1]])
eigen_data = A.eigenvects()

for eigenvalue, multiplicity, eigenvectors in eigen_data:
    print(f"eigenvalue = {eigenvalue}, eigenvector = {list(eigenvectors[0])}")

A linear system \vec X'=A\vec X solves exactly like §13.5's constant-coefficient equation, with the characteristic equation \det(A-\lambda I)=0 producing eigenvalues in place of roots, and eigenvectors supplying the fixed directions each exponential term grows or shrinks along. The phase plane visualizes every solution at once as trajectories through a vector field (§12.0), and the eigenvalues' signs classify the origin's behavior — node, saddle, spiral, or center — without ever needing an explicit formula for x(t) and y(t).

Next: an advanced closing technique — the Laplace transform, which converts differential equations into algebra entirely, handling discontinuous forcing terms that undetermined coefficients could never reach.

← Previous Next →