First-order systems, matrix form, phase planes, eigenvalue analysis, stability classification, and models of interacting populations, coupled mechanics, and economic dynamics
Until now, every differential equation in this book has described a single unknown function changing over time. A population grows. A temperature falls. A charge builds on a capacitor. Each of these situations is governed by one equation in one unknown, and solving it yields one function that describes the system completely.
But most real systems do not work this way. Most real systems involve several quantities changing simultaneously, each one affecting the others. A predator population rises when prey are plentiful, and falls when prey become scarce. The prey population, in turn, declines when predators are numerous and recovers when predators retreat. Neither population can be understood in isolation from the other, and no single differential equation can model both at once.
This chapter introduces systems of differential equations, collections of equations that together describe how several quantities change and interact. The key ideas are: how to write a system from a verbal description or a higher-order equation, how to express a linear system in matrix form, how to find equilibrium points, how to sketch and interpret phase plane portraits, how to solve linear systems using eigenvalues and eigenvectors, and how to classify the stability of equilibrium points.
Throughout, the emphasis is on understanding interaction and feedback. A system of differential equations is not merely several separate equations stacked together. It is a model of mutual dependence, in which the behavior of each quantity is shaped by the current state of all the others.
A rabbit population and a fox population share a meadow. When rabbits are abundant, foxes thrive and their numbers grow. As foxes multiply, rabbits are hunted more heavily and their population falls. When rabbits become scarce, foxes begin to starve and decline. When foxes are few, rabbits recover and their numbers rise again. The cycle repeats.
How should a mathematician write equations for a system like this? How can those equations be solved, or at least analyzed? Will the populations oscillate indefinitely? Will one species drive the other to extinction? Will the system settle into a stable balance? By the end of this chapter, you will have the tools to answer all of these questions.
Systems of differential equations appear throughout every scientific and applied field that uses mathematical modeling. In physics, Newton's laws for coupled mechanical systems produce systems of second-order equations that reduce to first-order systems. In engineering, electrical circuits with multiple loops, feedback control systems, and mechanical systems with multiple degrees of freedom are all modeled by systems of differential equations.
In biology, models of competing species, predator-prey dynamics, epidemic spread through compartments, enzyme kinetics, and neural activity all involve systems. In economics and finance, models of interacting markets, coupled asset prices, and macroeconomic systems use differential equation systems. In chemistry, reaction networks with multiple species interacting follow systems of differential equations. In data science and control systems, state-space representations, which are systems of differential equations in matrix form, are the standard language of dynamic models.
A student who understands systems of differential equations holds a key that opens the mathematics underlying a remarkable range of scientific and applied fields.
In Chapters 3 through 8, every differential equation described a single unknown function of one variable. The equation dy/dt = ky describes how a single quantity y changes over time, and its solution is a single function y(t). This framework is powerful, but it describes only isolated change: a single population in an unlimited environment, a single temperature object cooling toward a fixed ambient temperature, a single electrical component in isolation.
Real systems are rarely isolated. They involve multiple quantities that interact, depend on each other, and change together. To see why a single equation is insufficient, consider three situations.
| Situation | What changes | Why a single equation fails |
|---|---|---|
| Epidemic | Susceptible, infected, and recovered groups each change over time | Each group's rate of change depends on the sizes of the other groups |
| Coupled tanks | Concentration in Tank A and Tank B | Each tank's concentration changes based on the concentration in the other |
| Two-mass spring | Position and velocity of each mass | Each mass is subject to forces from both springs, so its motion affects and is affected by the other mass |
None of these can be described by a single differential equation in a single unknown. Each one requires a system of differential equations, a collection of equations that together track how several quantities change and interact simultaneously.
A system of differential equations is a collection of two or more differential equations involving two or more unknown functions, where the rate of change of each function may depend on any or all of the unknown functions.
The simplest systems involve first-order equations. Higher-order equations can often be reduced to first-order systems, which is one reason first-order systems are so central to applied mathematics.
Before writing equations, it helps to be deliberate about the components of a system model. Every system model involves four ingredients: state variables, rate equations, initial conditions, and parameters.
| Ingredient | Description | Example (predator-prey) |
|---|---|---|
| State variables | The unknown quantities that change over time | Prey population x(t), predator population y(t) |
| Rate equations | One differential equation per state variable, expressing how it changes | x' = ax − bxy and y' = −cy + dxy |
| Initial conditions | The values of all state variables at the starting time | x(0) = 200 rabbits, y(0) = 15 foxes |
| Parameters | Constants that characterize the system | Growth rate a, predation rate b, death rate c, conversion rate d |
A system of n first-order equations requires exactly n initial conditions to determine a unique solution. If your system has two equations in two unknowns, you need both x(0) and y(0). Never apply only one initial condition to a two-dimensional system and expect a unique answer.
One important source of systems is higher-order equations. A second-order equation in one unknown can always be rewritten as a first-order system in two unknowns. The idea is to introduce new variables for each derivative and then write a chain of first-order equations connecting them.
For example, the mass-spring equation from Chapter 7 is mx'' + bx' + kx = 0. Let x₁ = x (position) and x₂ = x' (velocity). Then x₁' = x₂ and x₂' = (1/m)(−bx₂ − kx₁). This is a first-order system in two unknowns. The same idea works for any higher-order equation.
x₁' = y' = x₂
x₂' = y'' = −3y' − 2y = −3x₂ − 2x₁
The system is:
x₁' = x₂
x₂' = −2x₁ − 3x₂
A first-order system in two unknowns x₁(t) and x₂(t) has the general form:
x₁' = f₁(t, x₁, x₂)
x₂' = f₂(t, x₁, x₂)
Each rate of change is a function that may depend on time t and on both current state variables. For n unknowns, the system has n such equations. This chapter focuses on linear systems, where each fᵢ is a linear combination of the state variables.
When writing a system from a verbal description, define each variable clearly, with its name, meaning, and units, before writing any equations. Skipping this step leads to equations that do not match the situation, and errors that are very difficult to trace. A well-labeled variable definition is not bookkeeping, it is part of the mathematics.
A' = A(0.4 − 0.02B)
The net per-capita growth rate of Species B is 0.3 − 0.01A, so:
B' = B(0.3 − 0.01A)
When a system is linear and homogeneous, it can be written in a compact and powerful matrix form. This is one of the places where linear algebra connects most directly to differential equations.
Consider a linear homogeneous system in two unknowns:
x₁' = a₁₁x₁ + a₁₂x₂
x₂' = a₂₁x₁ + a₂₂x₂
Define the state vector x = (x₁, x₂)ᵀ and the coefficient matrix A with entries aᵢⱼ. Then the entire system collapses to:
x' = Ax
This single matrix equation replaces the entire system. The notation is powerful: it shows that a linear system of differential equations has the same structural form as the scalar equation y' = ay, except that y is replaced by a vector x and the constant a is replaced by the matrix A. For n equations in n unknowns, the same format applies with an n × n matrix.
A = | 3 -1 | x = | x₁ |
| 1 2 | | x₂ |
The matrix equation is:
| x₁' | | 3 -1 | | x₁ |
| x₂' | = | 1 2 | | x₂ |
Before finding explicit solutions to a system, it is often valuable to ask whether there are constant solutions, states where the system simply stays put. These are called equilibrium points, and they represent the system at rest.
A point x* = (x₁*, x₂*) is an equilibrium point of the system x' = f(x) if f(x*) = 0. In other words, at an equilibrium point, every rate of change is zero. The system stays there if it starts there.
For a linear homogeneous system x' = Ax, the origin x* = 0 is always an equilibrium point. If det(A) ≠ 0 (that is, if A is invertible), then the origin is the only equilibrium point.
Case 1: x = 0 and y = 0 → Equilibrium: (0, 0)
Case 2: x = 0 and 3 − 2y = 0 → y = 3/2 → Equilibrium: (0, 3/2)
Case 3: 2 − x = 0 and y = 0 → x = 2 → Equilibrium: (2, 0)
Case 4: 2 − x − y = 0 and
3 − x − 2y = 0
Subtract: 1 − y = 0 → y = 1
Then x = 2 − y = 1 → Equilibrium: (1, 1)
When a system has two state variables, its behavior can be visualized in the phase plane, a two-dimensional coordinate plane where the horizontal axis represents x₁ and the vertical axis represents x₂. Rather than plotting x₁ and x₂ separately as functions of t, the phase plane shows the relationship between them directly, tracing a path called a trajectory as t increases.
Every point in the phase plane represents a possible state of the system: a particular combination of values of x₁ and x₂ at some moment in time. As time moves forward, the state moves from point to point, tracing a curve. The collection of all trajectories is called the phase portrait.
Draw a coordinate plane with x₁ on the horizontal axis and x₂ on the vertical axis. At several grid points (−2, −2), (−2, 0), (−2, 2), (0, −2), (0, 2), (2, −2), (2, 0), (2, 2), and so on, compute (x₁', x₂') from the system's equations and draw a short arrow at each grid point pointing in the direction of (x₁', x₂'). The result is a vector field. Mark equilibrium points with filled circles where all arrows vanish to zero. Sketch several trajectories by following the arrows from various starting points. Label the direction of increasing t along each trajectory.
Trajectories in the phase plane do not cross each other, except possibly at an equilibrium point. This is because at any given point in the phase plane, the direction of motion is uniquely determined by the system. Crossing would imply two different directions at the same point, which contradicts the uniqueness of solutions to initial value problems.
| What the phase portrait reveals | Geometric signature |
|---|---|
| Stable equilibria | Nearby trajectories converge toward the equilibrium as t → ∞ |
| Unstable equilibria | Nearby trajectories move away from the equilibrium |
| Saddle points | Some trajectories approach; others move away, depending on direction |
| Centers | Trajectories form closed loops, indicating perpetual oscillation |
| Spirals | Trajectories spiral toward or away from the equilibrium |
A phase plane plot and a time series plot give different information. The phase plane shows how x₁ and x₂ relate to each other, with t implicit. A time series shows each variable as a function of t explicitly. Both are useful, and a complete analysis often uses both. A closed loop in the phase plane corresponds to periodic oscillation in a time series graph.
For the scalar equation x' = ax, the solution is x(t) = Ceᵃᵗ. For a system x' = Ax, the analogous guess is a vector-valued exponential: x(t) = eλtv, where λ is a scalar and v is a constant vector. Substituting into x' = Ax:
λ e^(λt) v = A e^(λt) v
Divide by e^(λt) (never zero): Av = λv
This is precisely the eigenvalue equation. The scalar λ must be an eigenvalue of A, and the vector v must be a corresponding eigenvector. The eigenvalues of a 2×2 matrix A are the roots of the characteristic equation:
det(A − λI) = 0
Which expands to: λ² − tr(A)·λ + det(A) = 0
where tr(A) = a₁₁ + a₂₂ is the trace and det(A) = a₁₁a₂₂ − a₁₂a₂₁ is the determinant.
When the characteristic equation has two distinct real roots λ₁ and λ₂, with corresponding eigenvectors v₁ and v₂, the general solution is:
x(t) = c₁ e^(λ₁t) v₁ + c₂ e^(λ₂t) v₂
The constants c₁ and c₂ are determined by initial conditions.
det(A − λI) = (4−λ)(3−λ) − (1)(2)
= λ² − 7λ + 12 − 2
= λ² − 7λ + 10
= (λ−2)(λ−5) = 0
Eigenvalues: λ₁ = 2, λ₂ = 5
Step 2. Eigenvectors.
For λ₁ = 2: Solve (A − 2I)v = 0:
| 2 1 | |v₁| |0|
| 2 1 | |v₂| = |0|
From row 1: 2v₁ + v₂ = 0, so v₂ = −2v₁.
Choose v₁ = 1: v₁ = (1, −2)ᵀ
For λ₂ = 5: Solve (A − 5I)v = 0:
|−1 1 | |v₁| |0|
| 2 −2 | |v₂| = |0|
From row 1: −v₁ + v₂ = 0, so v₂ = v₁.
Choose v₁ = 1: v₂ = (1, 1)ᵀ
Step 3. General solution.
x(t) = c₁ e^(2t) | 1 | + c₂ e^(5t) | 1 |
|-2 | | 1 |
When the characteristic equation yields a repeated eigenvalue λ₁ = λ₂ = λ, the situation is more delicate. If A has two linearly independent eigenvectors for λ, the general solution uses both directly. But for most 2×2 matrices with a repeated eigenvalue, there is only one independent eigenvector v, and a generalized eigenvector w is needed. The vector w satisfies (A − λI)w = v, and the general solution takes the form:
x(t) = c₁ e^(λt) v + c₂ e^(λt)(tv + w)
The factor t in the second term closely mirrors the repeated-root case in second-order equations, where a factor of t multiplies the exponential when characteristic roots repeat.
det(A − λI) = (3−λ)² = 0 → λ = 3 (repeated)
Eigenvector: Solve (A − 3I)v = 0:
|0 1| |v₁| |0|
|0 0| |v₂| = |0|
v₂ = 0, v₁ free. Choose v = (1, 0)ᵀ.
Generalized eigenvector: Solve (A − 3I)w = v:
|0 1| |w₁| |1|
|0 0| |w₂| = |0|
From row 1: w₂ = 1. Let w₁ = 0. So w = (0, 1)ᵀ.
General solution:
x(t) = c₁ e^(3t) |1| + c₂ e^(3t) (t|1| + |0|)
|0| ( |0| |1|)
= c₁ e^(3t)|1| + c₂ e^(3t)| t |
|0| | 1 |
When the characteristic equation yields complex roots, the system produces oscillatory behavior: spirals or closed loops in the phase plane. For a real matrix A, complex eigenvalues always come in conjugate pairs: λ = α ± βi, where α is the real part and β is the imaginary part.
If v = p + iq is an eigenvector for λ = α + βi, the two real-valued solutions are:
x₁(t) = e^(αt)(cos(βt) p − sin(βt) q)
x₂(t) = e^(αt)(sin(βt) p + cos(βt) q)
The general solution is x(t) = c₁x₁(t) + c₂x₂(t). The exponential factor e^(αt) controls growth or decay. The trigonometric factors produce rotation in the phase plane. The sign of α determines whether the system spirals in, spirals out, or moves along a closed orbit.
det(A − λI) = (−1−λ)² + 4
= λ² + 2λ + 1 + 4
= λ² + 2λ + 5 = 0
λ = (−2 ± √(4−20))/2 = −1 ± 2i
So α = −1 and β = 2.
Eigenvector for λ = −1 + 2i.
Solve (A − (−1+2i)I)v = 0:
|−2i 2 | |v₁| |0|
|−2 −2i | |v₂| = |0|
From row 1: −2iv₁ + 2v₂ = 0, so v₂ = iv₁.
Let v₁ = 1:
v = |1| = |1| + i|0|
|i| |0| |1|
So p = (1, 0)ᵀ and q = (0, 1)ᵀ.
Real-valued solutions.
x₁(t) = e^(−t)(cos(2t)|1| − sin(2t)|0|) = e^(−t)| cos(2t)|
|0| |1| |−sin(2t)|
x₂(t) = e^(−t)(sin(2t)|1| + cos(2t)|0|) = e^(−t)| sin(2t)|
|0| |1| | cos(2t)|
General solution.
x(t) = c₁ e^(−t)| cos(2t)| + c₂ e^(−t)| sin(2t)|
|−sin(2t)| | cos(2t)|
When eigenvalues are complex, the procedure produces two real-valued solution vectors by applying Euler's formula to separate real and imaginary parts. Do not write the complex exponential e^((α+βi)t) directly as part of the final answer. The real-valued general solution is what is needed to satisfy real initial conditions and to make physical sense.
The eigenvalues of A determine not just the form of the solution but also the qualitative behavior of trajectories near the equilibrium at the origin. This classification is one of the most useful tools in applied systems analysis, because it tells you whether a system at rest will stay there, return there after a small disturbance, or drift away indefinitely.
| Eigenvalue type | Equilibrium name | Stability | Phase portrait behavior |
|---|---|---|---|
| Two real, both negative | Stable node | Stable | Trajectories converge to origin along eigenvector directions |
| Two real, both positive | Unstable node | Unstable | Trajectories diverge from origin |
| Two real, opposite signs | Saddle point | Unstable | Approach along one direction, escape along another |
| Repeated, negative | Stable star / degenerate node | Stable | Trajectories converge, possibly along one direction |
| Repeated, positive | Unstable star / degenerate node | Unstable | Trajectories diverge, possibly along one direction |
| Complex, α < 0 | Stable spiral | Stable | Trajectories spiral inward toward origin |
| Complex, α > 0 | Unstable spiral | Unstable | Trajectories spiral outward from origin |
| Complex, α = 0 | Center | Neutral | Closed loops around origin, perpetual oscillation |
Draw six separate coordinate planes labeled (a) through (f). In (a), sketch a stable node: trajectories converge to the origin along two straight-line directions corresponding to the two eigenvectors. In (b), an unstable node: all trajectories diverge. In (c), a saddle point: trajectories approach along one axis and escape along the other, forming a hyperbolic pattern with the two eigenvector lines as asymptotes. In (d), a stable spiral: trajectories rotate and converge, tightening toward the origin. In (e), an unstable spiral: trajectories rotate and expand outward. In (f), a center: closed elliptical trajectories encircle the origin indefinitely. Label each with the corresponding eigenvalue condition.
Stability is determined by the sign of the real part of the eigenvalues. For complex eigenvalues λ = α ± βi, stability depends on α, not on β. A center (α = 0) is neutrally stable: nearby trajectories neither approach nor escape the equilibrium. This is distinct from asymptotic stability, where trajectories actually converge.
det(A − λI) = (1−λ)(−1−λ) + 3
= −1 − λ + λ + λ² + 3
= λ² + 2
Setting equal to zero: λ² = −2, so λ = ±i√2.
The eigenvalues are purely imaginary: α = 0, β = √2.
Systems of differential equations arise naturally in population biology when two or more species compete for the same resources. The most widely studied model is the competitive Lotka-Volterra system, named for Alfred Lotka and Vito Volterra, who developed these equations independently in the 1920s.
Let x(t) and y(t) denote the populations of two competing species. In the absence of competition, each grows logistically toward its own carrying capacity. With competition, each species is additionally inhibited by the presence of the other:
x' = r₁x(1 − x/K₁ − α₁₂y/K₁)
y' = r₂y(1 − y/K₂ − α₂₁x/K₂)
Here r₁ and r₂ are intrinsic growth rates, K₁ and K₂ are carrying capacities, and α₁₂ and α₂₁ are competition coefficients measuring how strongly each species inhibits the other. The competition term α₁₂y/K₁ represents the reduction in Species x's effective resources due to the presence of Species y.
Setting both right-hand sides to zero yields four possible equilibria: (0, 0) where both species are extinct; (K₁, 0) where Species x is at carrying capacity and y is absent; (0, K₂) where Species y is at carrying capacity and x is absent; and a potential coexistence equilibrium inside the positive quadrant, if one exists.
The stability of each equilibrium depends on the competition coefficients. When competition is strong in one direction, one species can drive the other to extinction, a phenomenon called competitive exclusion. When both species inhibit each other more weakly than they inhibit themselves (that is, when α₁₂α₂₁ < 1), stable coexistence is possible.
(0, 0): Both species extinct — trivial equilibrium.
(0, 3/2): Species y at carrying capacity; species x absent.
(2, 0): Species x at carrying capacity; species y absent.
(1, 1): Both species coexist at positive population levels.
Whether the coexistence equilibrium (1, 1) is stable requires eigenvalue analysis of the Jacobian matrix at that point. The outcome, stable coexistence or competitive exclusion, depends on the relative strengths of the competition coefficients.
In competing species models, both species are harmed by the interaction. Predator-prey systems are different: one species benefits from the encounter (the predator), while the other is harmed (the prey). This asymmetry produces a distinctive oscillatory dynamics that has fascinated mathematicians and biologists for a century.
Let x(t) = prey population, y(t) = predator population, with a, b, c, d > 0:
x' = ax − bxy (prey grow exponentially without predators; predation reduces growth at rate bxy)
y' = −cy + dxy (predators decline without prey; consuming prey increases growth at rate dxy)
The equilibria are (0, 0) and (c/d, a/b). At the coexistence equilibrium (c/d, a/b), linearization produces purely imaginary eigenvalues, which means the equilibrium is a center in the linearized model. Trajectories in the full nonlinear model are closed loops, producing perpetual oscillation in both populations.
Draw two side-by-side figures. In the first (phase plane), place prey x on the horizontal axis and predators y on the vertical axis. Mark the coexistence equilibrium (c/d, a/b). Sketch several closed, roughly elliptical trajectories surrounding the equilibrium, with arrows showing counterclockwise motion: starting from a point with high prey and low predators, moving to high predators, then to low prey, then to low predators, then back to high prey. In the second (time series), draw time on the horizontal axis and plot x(t) and y(t) as oscillating curves. Show that the predator peak lags behind the prey peak by approximately a quarter period, reflecting the ecological reality that predator growth responds to prey abundance with a delay.
The populations cycle through a characteristic sequence. When prey are abundant, predators thrive and their numbers grow. As predators multiply, they consume more prey, and the prey population falls. As prey become scarce, predators can no longer sustain themselves and their numbers decline. With fewer predators, prey have a chance to recover. When prey are once again plentiful, the cycle begins again.
This periodic oscillation of predator and prey is one of the most celebrated predictions of mathematical biology. It has been observed in real ecosystems, most famously in data on Canadian lynx and snowshoe hare populations collected over more than a century.
Coexistence equilibrium: (c/d, a/b) = (0.4/0.2, 0.6/0.3) = (2, 2)
Consider two masses m₁ and m₂ connected by springs. Let x₁(t) and x₂(t) denote the displacements of the two masses from equilibrium. Suppose Mass 1 is attached to a wall by a spring with constant k₁, and the two masses are connected by a spring with constant k₂. Applying Newton's second law to each mass:
m₁x₁'' = −k₁x₁ + k₂(x₂ − x₁)
m₂x₂'' = −k₂(x₂ − x₁)
Introducing velocities v₁ = x₁' and v₂ = x₂', this becomes a first-order system in four unknowns (x₁, v₁, x₂, v₂). Its solution involves the eigenvalues of a 4×4 matrix. The system exhibits normal modes, specific oscillation patterns in which all parts of the system move at the same frequency. When driven at one of these natural frequencies, resonance can occur, with profound importance in structural engineering and the design of buildings, bridges, and mechanical components.
In macroeconomics, systems of differential equations describe how interconnected quantities such as national income, investment, and interest rates evolve over time. A simplified model of two interacting economies might specify:
Y₁' = α(C₁ + I₁ − Y₁) + β(Y₂ − Y₁)
Y₂' = α(C₂ + I₂ − Y₂) + β(Y₁ − Y₂)
Here Y₁ and Y₂ are national incomes, C and I represent consumption and investment in each economy, α is an adjustment speed, and β represents trade interaction. This is a linear system in Y₁ and Y₂. The equilibrium is where both economies reach steady state, and the eigenvalues of the coefficient matrix determine whether the economies converge to that steady state, oscillate around it, or diverge. Systems of this type are used in international economics and monetary policy analysis, where coupling between national economies has important implications for coordinated policy decisions.
Correction: The phase plane plots x₂ against x₁ with time as an implicit parameter. A time series plots x₁ or x₂ against t explicitly. A closed loop in the phase plane corresponds to periodic oscillation in a time series. Learn to move between both representations.
Correction: Any nonzero scalar multiple of an eigenvector is also an eigenvector. Choose a convenient normalization (such as setting one component to 1). The choice affects the constants c₁ and c₂ in the general solution but not the solution itself, because those constants are determined by initial conditions.
Correction: When eigenvalues are complex, apply Euler's formula to separate real and imaginary parts and write two real-valued solution vectors. The real-valued general solution is required to satisfy real initial conditions and to have physical meaning.
Correction: The characteristic equation for a 2×2 matrix is λ² − tr(A)λ + det(A) = 0, with a minus sign on the trace. The trace appears with a negative sign because det(A − λI) introduces a (−λ)² = λ² term and two diagonal entries each contributing −λ. Always verify by direct computation of det(A − λI).
Correction: Stability is determined by the sign of the real part of the eigenvalue, not its magnitude. λ = 3 + 0.1i is unstable because its real part is 3 > 0. λ = −1 is stable because its real part is −1 < 0. The imaginary part affects oscillation frequency, not stability.
Correction: A system of n first-order equations requires n initial conditions. For a two-dimensional system, both x₁(0) and x₂(0) must be specified to determine the unique solution. Applying only one initial condition leaves at least one constant undetermined.
Correction: The linearized system (obtained from the Jacobian at an equilibrium) approximates the nonlinear system only near that equilibrium. Conclusions about stability apply locally. A nonlinear system may have several equilibria with different stability types, and global behavior can differ markedly from local behavior.
Correction: An equilibrium point is a state where the system can remain at rest. Stability is a separate question: whether nearby trajectories return to or move away from the equilibrium after a small perturbation. A system can have many equilibria with different stability types, and simply knowing an equilibrium exists says nothing about whether the system will reach it.
Systems of differential equations mark a significant expansion of the subject. Instead of one unknown function, there are many. Instead of one equation, there is a coordinated collection. The mathematical tools, eigenvalues, phase planes, stability analysis, have grown correspondingly more powerful.
But systems also reveal a challenge that will become increasingly important in the chapters ahead: not all systems can be solved exactly. When the system is nonlinear, when the coefficients are not constant, or when the forcing is complicated, explicit solution formulas may be unavailable or impractical. Two major approaches address this challenge.
Laplace transforms, covered in Chapter 10, convert differential equations, including systems, into algebraic equations in a transformed variable. This is especially powerful when initial conditions are given from the start, when the system has discontinuous inputs such as switches or impulses, or when the system is being analyzed as a control system or circuit. For a system x' = Ax + f(t), the Laplace transform converts the equation into sX(s) − x(0) = AX(s) + F(s), which is solved by linear algebra and then inverted to recover x(t).
Numerical methods, covered in Chapter 11, compute approximate solutions by stepping through time from the initial condition. Euler's method, the improved Euler method, and Runge-Kutta methods all extend naturally to systems: at each time step, the update rule is applied to all n variables simultaneously. This makes numerical methods especially important for nonlinear systems, coupled systems with many components, and systems in scientific computing, data science, and engineering simulation.
The skills built in this chapter, writing systems in standard form, recognizing equilibria, interpreting phase portraits, applying eigenvalue analysis, and understanding stability, transfer directly into both Laplace and numerical work.
Chapter 10 introduces the Laplace transform, beginning with its definition and basic properties, then extending to initial value problems, step functions, impulse inputs, and circuit and control system applications. The eigenvalue structure of systems connects naturally to the pole structure of Laplace-domain transfer functions.
Chapter 11 introduces numerical methods, beginning with Euler's method for a single equation and extending immediately to systems of equations. A system of n equations is handled by applying the update formula to all n variables simultaneously, making numerical methods as natural for systems as for scalar equations.
Each problem shows a student's work. Identify the error and write the correct solution.
What is wrong? The student found one eigenvalue and one eigenvector but missed the second. The characteristic equation gives two eigenvalues. Find them, find both eigenvectors, and write the complete general solution.
Why is this form problematic for a real physical system? Describe the correct procedure: extract real and imaginary parts of the eigenvector, apply the real-solution formula, and write the general solution using real-valued functions only.
Is this classification correct? Explain the difference between neutral stability (center, trajectories bounded but not converging) and asymptotic stability (trajectories actually converge to the equilibrium as t → ∞). A center is neutrally stable, not asymptotically stable.
Identify the sign error. The characteristic equation is λ² − tr(A)λ + det(A) = 0. Here tr(A) = 7 and det(A) = 12 − 2 = 10, so the correct equation is λ² − 7λ + 10 = 0.
What error has the student made? Both c₁v₁ and c₂v₂ are present at t = 0. The correct equation is x(0) = c₁v₁ + c₂v₂, giving two scalar equations in two unknowns. Both initial conditions x₁(0) and x₂(0) must be used.
Chapter 9 expanded the subject of differential equations from single equations to systems of interacting quantities. The eigenvalue method provides elegant general solutions when the coefficient matrix is constant and the system is linear. But two important situations remain difficult: when the input to the system switches on or off at a specific time (a step function or an impulsive load), and when it is desirable to incorporate initial conditions directly into the solution process rather than solving for constants after the fact.
Laplace transforms address both challenges. The Laplace transform converts a differential equation (or system of differential equations) into an algebraic equation in a new variable s. Initial conditions appear immediately in the algebra. Discontinuous inputs are handled naturally through step functions and their transforms. The solution is recovered by inverting the transform back into time behavior.
Chapter 10 builds these ideas from the ground up: the transform definition, its basic properties, transform tables, inverse transforms using partial fractions, step functions, impulse inputs, and applications to circuits and control systems. The eigenvalue and phase plane thinking developed in Chapter 9 will connect naturally to the pole structure and stability analysis that arises in Laplace-domain control theory.