Systems of ODEs

Real systems rarely have just one moving part. Two tanks exchange fluid; the floors of a building sway against each other; a vehicle body bounces on its suspension while its wheels bounce on the tires; a chemical reaction consumes one species while producing another; predators and prey rise and fall together. Each needs several unknowns that depend on each other, governed by a system of differential equations.

This final chapter shows how to write such systems in matrix form, solve linear constant-coefficient systems with eigenvalues and eigenvectors (from MATH-125), interpret the results physically as decay modes and vibration modes, picture behavior with phase portraits, and solve systems numerically, including nonlinear ones, with a few lines of Python.

Systems in matrix form

A linear system of first-order equations with constant coefficients is written

where A is a constant square matrix. Every higher-order equation can be rewritten this way, which is how numerical solvers handle them: introduce the derivatives as new unknowns.

Solving x′ = Ax with eigenvalues

Try \(\mathbf{x} = \mathbf{v}e^{\lambda t}\), a fixed vector times an exponential. Substituting gives \(\lambda\mathbf{v}e^{\lambda t} = A\mathbf{v}e^{\lambda t}\), so

The exponents are the eigenvalues of A and the vectors are its eigenvectors. For an n × n matrix with n independent eigenvectors, the general solution is a combination of these modes:

Each mode is a fixed pattern (the eigenvector) that grows or decays at its own rate (the eigenvalue). The initial condition decides how much of each mode is present.

Vibration modes

For undamped mechanical systems, Newton's law gives a second-order matrix equation \(M\ddot{\mathbf{x}} + K\mathbf{x} = \mathbf{0}\), with a mass matrix M and a stiffness matrix K. Trying \(\mathbf{x} = \boldsymbol{\phi}\cos\omega t\) gives the eigenvalue problem

Each eigenvalue \(\omega^2\) gives a natural frequency, and each eigenvector \(\boldsymbol{\phi}\) a mode shape: the pattern in which all the masses move together at that frequency. A structure with n degrees of freedom has n natural frequencies and n modes. Earthquake engineers design buildings by examining these modes; the first (lowest) mode usually dominates the response.

Phase portraits and stability

For a two-dimensional system, plot \(x_2\) against \(x_1\) instead of each against time. Each solution traces a curve, a trajectory, in this phase plane, and a set of trajectories is a phase portrait. The eigenvalues of A determine its shape around the equilibrium at the origin:

• Both eigenvalues real and negative: a stable node; all trajectories flow into the origin.
• Real with opposite signs: a saddle; almost all trajectories eventually run away. Unstable.
• Complex with negative real part: a stable spiral; trajectories spiral in (a damped oscillation).
• Purely imaginary: a center; closed orbits (an undamped oscillation).
• Any eigenvalue with positive real part: unstable.

Numerical solution

Most real systems, especially nonlinear ones, have no closed-form solution. Numerical solvers step them forward in time with much more accurate versions of Euler's method (ENCM-100). In Python, scipy.integrate.solve_ivp does this for any system written as a function returning \(\mathbf{x}'\).

import numpy as np
from scipy.integrate import solve_ivp

A = np.array([[-0.2, 0.1],
              [ 0.2, -0.2]])

def rates(t, x):
    return A @ x

sol = solve_ivp(rates, [0, 30], [10, 0], dense_output=True, rtol=1e-8)
x1, x2 = sol.sol(10)
print(f"After 10 min: tank 1 = {x1:.3f} kg, tank 2 = {x2:.3f} kg")

Output: After 10 min: tank 1 = 2.948 kg, tank 2 = 3.704 kg, matching the eigenvalue solution. Changing the model to a nonlinear one, or to 50 tanks, changes only the rates function.

A nonlinear system: predators and prey

The Lotka-Volterra model describes two interacting populations, prey x and predators y:

Prey grow on their own and are eaten; predators die off on their own and grow by eating. The same structure appears in ecology, in some chemical reaction networks and in economics. Linearizing near an equilibrium (using the matrix of partial derivatives, the Jacobian, from MATH-209) shows its local behavior.

This is where the course ends and much of engineering analysis begins. The finite element programs that analyze structures, the circuit simulators that design electronics, the process simulators that run chemical plants and the climate and reservoir models that inform policy all solve very large systems of differential equations. Every one of them relies on the ideas here: linearity and superposition, modes and eigenvalues, stability, and careful numerical solution checked against cases you can solve by hand.