Systems of Differential Equations

Many engineering models do not move one quantity at a time. Two tanks trade brine, two floors of a building sway against each other, a predator and its prey rise and fall together. Each of these is a set of first-order differential equations that are coupled, meaning the rate of change of one unknown depends on the others. Stack the unknowns into a single vector and the whole tangle collapses into one clean line, x' = A x, and from there the eigenvalue machinery of the previous lessons solves it and even predicts the long-run behavior at a glance. This lesson builds that pipeline from the matrix form up, reads off whether the equilibrium is a node, a saddle, or a spiral, works a coupled-tank problem, and shows how any higher-order ODE folds into the same first-order shape. Every worked example is carried all the way to a finished number.

First-order linear systems and the matrix form x' = A x

A first-order linear system with constant coefficients is a collection of coupled equations such as

x1' = a11 x1 + a12 x2 x2' = a21 x1 + a22 x2,

where each unknown is a function of the same independent variable t and the coefficients are plain numbers. Gather the unknowns into a state vector x = (x1, x2) and the coefficients into a matrix A, and the two lines become a single vector equation:

x' = A x, with A = [[a11, a12], [a21, a22]].

Here x' is the vector of derivatives (x1', x2'). The system is homogeneous because there is no separate forcing term; adding one would give x' = A x + g(t). A homogeneous constant-coefficient system is autonomous, so its only equilibrium (a state where nothing changes, x' = 0) sits at the origin whenever A is nonsingular, because A x = 0 then forces x = 0. Everything below is about how solutions move around that equilibrium.

An n by n system needs n linearly independent solutions, and the general solution is their linear combination with one arbitrary constant apiece. For the 2 by 2 case that means two building-block solutions and two constants, fixed later by two initial conditions.

Solving by eigenvalues and eigenvectors

The single idea that cracks the whole system is to look for a solution that keeps a fixed direction and only grows or shrinks in time. Guess

x = e^(lambda t) v,

where v is a constant nonzero vector and lambda is a number. Differentiating gives x' = lambda e^(lambda t) v, while the right side is A x = e^(lambda t) A v. Setting them equal and cancelling the never-zero factor e^(lambda t) leaves

A v = lambda v, that is (A - lambda I) v = 0.

This is exactly the eigenvalue problem. So lambda must be an eigenvalue of A and v its eigenvector, both found from det(A - lambda I) = 0 as in the eigenvalue lesson. When a 2 by 2 matrix has two distinct real eigenvalues lambda1 and lambda2 with eigenvectors v1 and v2, the general solution is

x(t) = C1 e^(lambda1 t) v1 + C2 e^(lambda2 t) v2.

Each term is a mode: it points along its eigenvector and its size rides the exponential e^(lambda t). A positive eigenvalue grows the mode, a negative one shrinks it, and that single sign is what decides the geometry of the whole picture.

Worked example: Solve x' = A x for A = [[1, 2], [2, 1]] with x(0) = (3, 1), where x1 and x2 are displacements in metres and t is in seconds, then classify the equilibrium. First the eigenvalues: trace = 1 + 1 = 2 and det = (1)(1) - (2)(2) = 1 - 4 = -3, so the characteristic equation is lambda^2 - 2 lambda - 3 = 0, which factors as (lambda - 3)(lambda + 1) = 0, giving lambda = 3 and lambda = -1. For lambda = 3, A - 3I = [[-2, 2], [2, -2]], and the top row -2 v1 + 2 v2 = 0 gives v2 = v1, so v1 = (1, 1). For lambda = -1, A + I = [[2, 2], [2, 2]], and 2 v1 + 2 v2 = 0 gives v2 = -v1, so v2 = (1, -1). The general solution is

x(t) = C1 e^(3t) (1, 1) + C2 e^(-t) (1, -1).

Apply x(0) = (3, 1): C1 + C2 = 3 and C1 - C2 = 1, so C1 = 2 and C2 = 1. Then x1(t) = 2 e^(3t) + e^(-t) and x2(t) = 2 e^(3t) - e^(-t). At t = 0.5 s, e^(1.5) = 4.4817 and e^(-0.5) = 0.6065, so x1(0.5) = 2(4.4817) + 0.6065 = 9.57 m and x2(0.5) = 2(4.4817) - 0.6065 = 8.36 m. Because the eigenvalues 3 and -1 have opposite signs, the origin is a saddle point, and the growing e^(3t) mode confirms the motion runs away.

Classifying the equilibrium: nodes, saddles, and spirals

You often do not need the full solution, only the character of the motion near the origin, and the eigenvalues alone settle it. The complete catalogue for a 2 by 2 system is short enough to memorize.

Eigenvalues of A	Geometry of the origin	Stability
real, distinct, both negative	stable node (sink)	asymptotically stable
real, distinct, both positive	unstable node (source)	unstable
real, opposite signs	saddle point	unstable
complex, negative real part	stable spiral (focus)	asymptotically stable
complex, positive real part	unstable spiral	unstable
purely imaginary, zero real part	center	stable, not asymptotic
repeated real, negative	degenerate or star node	asymptotically stable

Read the pattern this way. Real eigenvalues give straight-line eigendirections and no rotation, so the origin is a node when they share a sign and a saddle when they differ. Complex eigenvalues alpha +/- beta i mean the modes spiral, because the imaginary part beta rotates the vector while the real part alpha grows or decays its length, so the origin is a spiral when alpha is nonzero and a center when alpha = 0. In every case the sign of the real part is the referee of stability: all real parts negative means every trajectory decays to the origin (asymptotically stable), and any positive real part means at least one mode runs away (unstable).

For a 2 by 2 matrix you can skip finding the eigenvalues and read the same verdict from the trace and determinant, since the eigenvalues satisfy lambda^2 - (trace) lambda + (det) = 0 with discriminant (trace)^2 - 4(det).

Trace-determinant test (2 by 2)	Result
det < 0	saddle point, always
det > 0 and (trace)^2 - 4 det > 0	node, stable if trace < 0
det > 0 and (trace)^2 - 4 det < 0	spiral, stable if trace < 0
det > 0 and trace = 0	center

Worked example: Classify the equilibrium of x' = A x for A = [[-1, -4], [1, -1]], with t in seconds. Trace = -1 + (-1) = -2 and det = (-1)(-1) - (-4)(1) = 1 + 4 = 5, so det > 0 and the discriminant (trace)^2 - 4 det = (-2)^2 - 4(5) = 4 - 20 = -16 < 0, which means complex eigenvalues and a spiral. Solving, lambda = (trace +/- sqrt(discriminant)) / 2 = (-2 +/- sqrt(-16)) / 2 = (-2 +/- 4i) / 2 = -1 +/- 2i. The real part alpha = -1 is negative, so the origin is a stable spiral that winds inward. The imaginary part beta = 2 rad/s is the rotation rate, so one loop of the spiral takes a period T = 2 pi / beta = 2 pi / 2 = 3.14 s.

Coupled tanks and mixing problems

Mixing problems are the classic place a system shows up on the board. The rule for each tank is the same conservation statement used for a single tank: the rate of change of dissolved salt equals the rate salt flows in minus the rate it flows out, where each stream carries salt at (flow rate) times (concentration) and concentration is (salt in that tank) divided by (its volume). Two tanks that trade fluid make the two equations depend on each other, and that coupling is the whole point.

Worked example: Two tanks each hold 100 L of brine, connected in a loop that carries 10 L/min from tank 1 to tank 2 and 10 L/min back, so both volumes stay at 100 L. Let x1 and x2 be the kilograms of salt in each tank. Salt enters tank 1 from tank 2 at (10 L/min)(x2 / 100 kg/L) = 0.1 x2 kg/min and leaves to tank 2 at (10 L/min)(x1 / 100 kg/L) = 0.1 x1 kg/min, and the mirror image holds for tank 2, so

x1' = 0.1 x2 - 0.1 x1 x2' = 0.1 x1 - 0.1 x2, that is A = [[-0.1, 0.1], [0.1, -0.1]].

Trace = -0.2 and det = (-0.1)(-0.1) - (0.1)(0.1) = 0.01 - 0.01 = 0, so the characteristic equation lambda^2 + 0.2 lambda = 0 gives lambda = 0 and lambda = -0.2. The eigenvalue 0 has eigenvector (1, 1), the direction of the total salt x1 + x2, which is conserved because the loop is closed; the eigenvalue -0.2 has eigenvector (1, -1), the imbalance that decays away. The general solution is x(t) = C1 (1, 1) + C2 e^(-0.2 t) (1, -1). Start tank 1 with 20 kg and tank 2 with 0 kg: C1 + C2 = 20 and C1 - C2 = 0, so C1 = 10 and C2 = 10, giving x1(t) = 10 + 10 e^(-0.2 t) kg and x2(t) = 10 - 10 e^(-0.2 t) kg. At t = 5 min the exponent is -0.2(5) = -1, so e^(-1) = 0.368 and x1(5) = 10 + 10(0.368) = 13.68 kg (with x2(5) = 6.32 kg, and the two still sum to the conserved 20 kg). As t goes to infinity the exponential dies and each tank settles at 10 kg, the even split of the total.

Converting a higher-order ODE to a first-order system

Any higher-order linear ODE can be rewritten as a first-order system, which is why numerical solvers and stability tools only ever need the first-order form. The recipe is to name each derivative as a new unknown. For a second-order equation y'' + p y' + q y = 0, set

x1 = y and x2 = y'.

Then x1' = y' = x2, and the ODE itself, solved for the top derivative, gives x2' = y'' = -q y - p y' = -q x1 - p x2. In matrix form x' = A x with the companion matrix

A = [[0, 1], [-q, -p]].

The pattern generalizes: an n-th order equation becomes an n by n system whose companion matrix carries the coefficients along its bottom row. The payoff is that the eigenvalues of the companion matrix are exactly the roots of the original characteristic equation, so the two viewpoints agree and you can classify the motion with the same node, saddle, or spiral labels.

Worked example: Convert y'' + 3 y' + 2 y = 0 to a system and solve with y(0) = 1 m and y'(0) = 0, where t is in seconds. Set x1 = y, x2 = y', so x1' = x2 and x2' = -2 x1 - 3 x2, giving A = [[0, 1], [-2, -3]]. Its trace = 0 + (-3) = -3 and det = (0)(-3) - (1)(-2) = 2, so lambda^2 + 3 lambda + 2 = 0 factors as (lambda + 1)(lambda + 2) = 0 and lambda = -1 and lambda = -2, exactly the roots of the original characteristic equation r^2 + 3 r + 2 = 0. Both are negative and real, so the origin is a stable node. In terms of y the general solution is y = C1 e^(-t) + C2 e^(-2t). From y(0) = 1, C1 + C2 = 1; from y'(0) = -C1 - 2 C2 = 0, C1 = -2 C2, so -2 C2 + C2 = 1 gives C2 = -1 and C1 = 2. Then y(t) = 2 e^(-t) - e^(-2t) m, and at t = 1 s, y(1) = 2(0.3679) - 0.1353 = 0.7358 - 0.1353 = 0.60 m.

Exam-day strategy

• Write the system as x' = A x first: stack the unknowns into a vector and pull the constant coefficients into A, then the whole problem becomes an eigenvalue problem for A.
• Find the eigenvalues from lambda^2 - (trace) lambda + (det) = 0, then attach each eigenvector to build x(t) = C1 e^(lambda1 t) v1 + C2 e^(lambda2 t) v2; keep each eigenvalue paired with its own eigenvector.
• Classify without solving when only the behavior is asked: det < 0 is always a saddle, det > 0 with a positive discriminant is a node, and det > 0 with a negative discriminant is a spiral, with the sign of the trace fixing stable versus unstable.
• Remember the rule of thumb on stability: every eigenvalue with negative real part means the origin is asymptotically stable, any positive real part means it is unstable, and purely imaginary eigenvalues give a center that neither grows nor decays.
• For a mixing problem, write each rate as (flow in)(concentration in) minus (flow out)(concentration out), and divide salt by the tank volume for concentration; a closed loop gives a zero eigenvalue that flags a conserved total.
• To handle a higher-order ODE as a system, set x1 = y, x2 = y', and so on, and check your work by confirming the companion matrix eigenvalues match the characteristic roots of the original equation.
• Sanity-check any eigenvalue pair against the trace and determinant: they must sum to the trace and multiply to the determinant, which catches a sign slip before it costs you the item.