Chapter 19 — Differential Equations

19.1 Nature's Native Language

For eighteen chapters you have learned to take a function apart with the derivative and put it back together with the integral. This chapter asks the question those two skills were built to answer: if you know how fast something is changing, can you recover the thing itself?

That question is the whole of this chapter, and the answer reorganizes how you see the world. Almost no law of nature tells you directly what a quantity is. Newton did not write down a formula for where a planet sits; he wrote down how its position changes — its acceleration — and let the mathematics deliver the orbit. Epidemiologists do not write a formula for how many people are sick on day 40; they write down how the number of sick people changes each day, and solve for the curve. The laws come to us as statements about rates. Converting a rate law into the quantity it governs is what solving a differential equation means.

A differential equation (DE) is an equation involving an unknown function and one or more of its derivatives. We study the most common case, the ordinary differential equation (ODE): an unknown function of one variable, related to its own derivatives. A few that will recur:

• $\dfrac{dy}{dx} = ky$ — exponential growth or decay (the simplest, and the most pervasive).
• $\dfrac{dP}{dt} = rP\left(1 - \dfrac{P}{K}\right)$ — logistic growth toward a carrying capacity $K$.
• $\dfrac{d^2 x}{dt^2} = -\omega^2 x$ — the simple harmonic oscillator (a mass on a spring).
• $\dfrac{dS}{dt} = -\beta S I,\quad \dfrac{dI}{dt} = \beta S I - \gamma I,\quad \dfrac{dR}{dt} = \gamma I$ — the SIR epidemic model (a system of ODEs).

The Key Insight. Solving a differential equation is antidifferentiation with feedback. In an ordinary integral $\int f(x)\,dx$ the right-hand side depends only on the independent variable. In a differential equation the rate depends on the unknown function itself — the quantity steers its own change. That feedback is what makes ODEs powerful enough to describe growth, decay, oscillation, and contagion, and what makes most of them impossible to solve with a single antiderivative.

This is Chapter 19, the capstone of Part III. Everything you built — the definite integral (Chapter 13), the Fundamental Theorem of Calculus (Chapter 14), and the integration techniques of Chapters 15 and 16 — converges here, because solving an ODE almost always comes down to computing an integral. We will develop just enough theory to model population growth, disease spread, cooling, and mixing, and we will introduce in full the SIR model of epidemics — an anchor example (see the continuity tracker) that climaxes in the Chapter 39 capstone. The complete theory of differential equations is a separate course, often two; this is your invitation into it.

19.2 Anatomy of a Differential Equation

A first-order ODE has the general form

$$\frac{dy}{dx} = F(x, y),$$

where $F$ is some expression built from the independent variable $x$ and the unknown function $y(x)$. To solve it is to find a function $y(x)$ that makes the equation an identity. Usually there is a whole family of solutions; pinning down a single one requires an initial condition $y(x_0) = y_0$. The equation plus its initial condition is an initial value problem (IVP).

Two vocabulary words organize the entire subject.

Order. The order of an ODE is the highest derivative that appears. $y' = ky$ is first order; $x'' = -\omega^2 x$ is second order. This chapter is almost entirely about first-order equations, with one famous second-order guest (the oscillator).

Linear vs. nonlinear. An ODE is linear if the unknown $y$ and its derivatives appear only to the first power and are not multiplied by each other — the coefficients may depend on the independent variable but not on $y$. So $$y' + p(x)\,y = q(x)$$ is linear, while $y' = y^2$, $y' = \sin y$, and $S' = -\beta S I$ are nonlinear because $y$ enters through a square, a sine, or a product of unknowns. Linearity is the great dividing line: linear equations submit to systematic formulas; most nonlinear equations do not, and we reach for numerical methods.

Check Your Understanding. Classify $\dfrac{dy}{dx} = x^2 y + \cos x$ by order and linearity, and verify whether $y = x^3$ solves $\dfrac{dy}{dx} = \dfrac{3y}{x}$.

Answer

The first equation is first-order linear: $y$ appears to the first power, multiplied only by the function $x^2$, and $\cos x$ is a free term — it fits $y' + p(x)y = q(x)$ with $p(x) = -x^2$, $q(x) = \cos x$. For the second: if $y = x^3$ then $y' = 3x^2$, while $\dfrac{3y}{x} = \dfrac{3x^3}{x} = 3x^2$. The two sides agree, so $y = x^3$ is indeed a solution. (It is one member of the family $y = Cx^3$.)

Before any technique, internalize the one habit that prevents most errors: a candidate solution can always be checked by substitution. Differentiate it, plug into the equation, and confirm both sides match. Unlike most of mathematics, ODEs let you verify the answer even when you could never have guessed it — a fact we will lean on repeatedly.

19.3 Separable Equations

The first solvable class is the most intuitive. A first-order ODE is separable if its right-hand side factors into a function of $x$ times a function of $y$:

$$\frac{dy}{dx} = g(x)\,h(y).$$

The name describes the method. Treat $dy$ and $dx$ as the differentials they secretly are (this is the Leibniz notation paying off), gather everything with $y$ on one side and everything with $x$ on the other, and integrate:

$$\frac{dy}{h(y)} = g(x)\,dx \quad\Longrightarrow\quad \int \frac{dy}{h(y)} = \int g(x)\,dx.$$

Each side is now an ordinary integral of the kind you mastered in Chapters 14–16. The Fundamental Theorem does the rest. Solve for $y$ explicitly if you can; otherwise leave the relationship implicit.

Geometric Intuition. Why is "separating $dx$ and $dy$" legitimate and not algebraic vandalism? Picture the solution curve in the plane. At each point it has slope $g(x)h(y)$. Separating variables rewrites this as $\frac{1}{h(y)}\,dy = g(x)\,dx$: it says that as you step along the curve, the accumulated change in the quantity $\int \frac{dy}{h(y)}$ keeps perfect pace with the accumulated change in $\int g(x)\,dx$. Integrating both sides is just totaling those two ledgers and demanding they stay equal. The differentials are not being abused; they are bookkeeping for matched infinitesimal steps, justified rigorously by the chain rule.

Worked Example 19.3.1 — Exponential decay (the master equation)

Solve $\dfrac{dy}{dt} = -ky$ with $k > 0$ and $y(0) = y_0$.

Separate: $\dfrac{dy}{y} = -k\,dt$. Integrate both sides:

$$\ln|y| = -kt + C.$$

Exponentiate: $|y| = e^{-kt + C} = e^{C}e^{-kt}$. Absorbing the sign and the constant $e^C$ into a single constant $A = \pm e^C$ gives $y = A e^{-kt}$. The initial condition fixes $A$: setting $t = 0$ gives $y_0 = A$, so

$$\boxed{\,y(t) = y_0\,e^{-kt}.\,}$$

This single solution is the backbone of half the applications in this chapter:

• Radioactive decay. $k$ is the decay constant. The half-life — the time for $y$ to halve — comes from $\tfrac12 = e^{-k t_{1/2}}$, giving $t_{1/2} = \dfrac{\ln 2}{k}$.
• Drug clearance. $k$ is the elimination rate; $1/k$ is the mean residence time of the drug (built on in the Chapter 15 pharmacokinetics case study).
• Atmospheric pressure with altitude, capacitor discharge, cooling of a small excess temperature — all the same equation wearing different clothes.

Historical Note. The equation $y' = ky$ is sometimes called the Malthusian law after Thomas Malthus, whose 1798 Essay on the Principle of Population argued that populations grow geometrically while food grows arithmetically. The same equation governs compound interest (Jacob Bernoulli, 1683, in the calculation that first uncovered the number $e$) and radioactivity (Rutherford and Soddy, 1902). One differential equation, three centuries, three sciences — a recurring lesson that calculus is field-agnostic.

Worked Example 19.3.2 — A nonlinear separable equation with a surprise

Solve $\dfrac{dy}{dx} = x\,y^2$ with $y(0) = 1$.

Separate and integrate: $\dfrac{dy}{y^2} = x\,dx$ gives $-\dfrac{1}{y} = \dfrac{x^2}{2} + C$. Apply $y(0) = 1$: $-1 = 0 + C$, so $C = -1$. Then

$$-\frac{1}{y} = \frac{x^2}{2} - 1 \quad\Longrightarrow\quad y = \frac{1}{1 - x^2/2} = \frac{2}{2 - x^2}.$$

Check by substitution: $y = 2(2-x^2)^{-1}$, so $y' = 2\cdot 2x\,(2-x^2)^{-2} = 4x/(2-x^2)^2$, while $xy^2 = x\cdot 4/(2-x^2)^2$ — the two agree. The answer is correct, and it is alarming: as $x \to \sqrt{2}$ the denominator vanishes and $y \to \infty$. The solution blows up in finite time even though the equation looks tame. Nonlinearity buys this kind of drama; a linear equation can never reach infinity in finite time.

Warning. A solution to an IVP is only guaranteed on the largest interval containing the initial point on which it stays finite and the equation stays well behaved. Here the solution through $(0,1)$ lives on $(-\sqrt2, \sqrt2)$ and no further — the formula $2/(2-x^2)$ on the far side of the asymptote is a different solution, not a continuation. The existence-and-uniqueness theorem (the Math Major Sidebar below) promises a solution locally, never globally.

Math Major Sidebar — Existence and uniqueness (Picard–Lindelöf). When can we be sure an IVP $y' = F(x,y),\ y(x_0)=y_0$ has exactly one solution? The Picard–Lindelöf theorem says: if $F$ is continuous near $(x_0,y_0)$ and Lipschitz continuous in $y$ (there is a constant $L$ with $|F(x,y_1)-F(x,y_2)| \le L|y_1-y_2|$, which holds whenever $\partial F/\partial y$ is bounded), then a unique solution exists on some interval around $x_0$. Continuity alone (Peano's theorem) guarantees existence but not uniqueness: $y' = \sqrt{|y|},\ y(0)=0$ has both $y\equiv 0$ and $y = x^2/4$ as solutions because $\sqrt{|y|}$ fails the Lipschitz condition at $0$. Uniqueness is what lets a deterministic law have a deterministic future.

19.4 First-Order Linear Equations and the Integrating Factor

Not every useful equation separates. The mixing problems, cooling problems, and circuit equations of this chapter all share the linear form

$$\frac{dy}{dx} + P(x)\,y = Q(x),$$

which separates only when $Q \equiv 0$. We need a different trick — and it is one of the most elegant in elementary calculus.

The obstacle is that $\frac{dy}{dx} + P y$ is not the derivative of anything obvious. The idea is to multiply the whole equation by a cleverly chosen function $\mu(x)$, the integrating factor, so that the left side collapses into a single derivative. We want

$$\mu y' + \mu P y = \frac{d}{dx}\big(\mu y\big) = \mu y' + \mu' y.$$

Matching the two expressions forces $\mu' = \mu P$, a separable equation for $\mu$ itself, whose solution is

$$\mu(x) = e^{\int P(x)\,dx}.$$

With that $\mu$, the equation becomes $\dfrac{d}{dx}\big(\mu y\big) = \mu Q$, and now both sides integrate directly:

$$\mu(x)\,y = \int \mu(x)\,Q(x)\,dx + C \quad\Longrightarrow\quad y = \frac{1}{\mu(x)}\left(\int \mu Q\,dx + C\right).$$

The Key Insight. The integrating factor turns the product rule backward. You recognized in Chapter 15 that $u$-substitution reverses the chain rule and integration by parts reverses the product rule; the integrating factor is a third instance of the same philosophy — reverse-engineer a known differentiation rule to make an unrunnable left side into a perfect derivative. Every solution technique in calculus is some differentiation rule read in reverse.

Worked Example 19.4.1 — A clean integrating factor

Solve $\dfrac{dy}{dx} + 2y = e^{-x}$.

Here $P(x) = 2$ and $Q(x) = e^{-x}$, so $\mu = e^{\int 2\,dx} = e^{2x}$. Multiply through:

$$e^{2x}y' + 2e^{2x}y = e^{2x}e^{-x} = e^{x}, \qquad\text{i.e.}\qquad \frac{d}{dx}\big(e^{2x}y\big) = e^{x}.$$

Integrate: $e^{2x}y = e^{x} + C$, so $y = e^{-x} + C e^{-2x}$. Substituting back confirms it. Notice the structure of the answer: $e^{-x}$ is a particular solution driven by the forcing term $e^{-x}$, and $Ce^{-2x}$ is the homogeneous part that decays away. Every linear ODE solution splits this way — a transient that dies plus a steady response to the input.

Worked Example 19.4.2 — A mixing problem (chemical and environmental engineering)

A tank holds $100$ L of pure water. Brine containing $5$ g of salt per liter flows in at $2$ L/min, and the perfectly stirred mixture drains out at the same $2$ L/min, so the volume stays at $100$ L. How does the salt content evolve, and what is the long-run amount?

Let $y(t)$ be the grams of salt in the tank at time $t$, with $y(0) = 0$. The governing principle is conservation of mass: rate of change of salt = rate in − rate out.

• Rate in: $2\ \tfrac{\text{L}}{\text{min}} \times 5\ \tfrac{\text{g}}{\text{L}} = 10$ g/min.
• Rate out: the concentration in the tank is $\dfrac{y}{100}$ g/L, draining at $2$ L/min, so $2\cdot\dfrac{y}{100} = \dfrac{y}{50}$ g/min.

The model is therefore

$$\frac{dy}{dt} = 10 - \frac{y}{50}, \qquad \text{or} \qquad \frac{dy}{dt} + \frac{1}{50}\,y = 10.$$

This is linear with $P = \tfrac{1}{50}$, $Q = 10$, so $\mu = e^{t/50}$. Then $\dfrac{d}{dt}\big(e^{t/50}y\big) = 10\,e^{t/50}$, and integrating,

$$e^{t/50}y = 500\,e^{t/50} + C \quad\Longrightarrow\quad y(t) = 500 + C e^{-t/50}.$$

The initial condition $y(0) = 0$ gives $C = -500$, so

$$y(t) = 500\big(1 - e^{-t/50}\big).$$

As $t \to \infty$, $y \to 500$ g — exactly the inflow concentration ($5$ g/L) times the tank volume ($100$ L), the equilibrium where what comes in equals what goes out. Mixing problems like this model pollutant flushing in lakes, drug concentration in the bloodstream, and ion exchange in chemical reactors; the equation is always the same conservation statement.

Common Pitfall. Students often compute the "rate out" as $Q\cdot(\text{volume})$ or forget to convert concentration to grams. The clean rule is: rate out $=$ (volumetric outflow rate) $\times$ (concentration in the tank) $= (\text{flow}) \times \dfrac{y}{V}$. And if inflow and outflow rates differ, the volume $V$ is itself a function of time, $V(t) = V_0 + (\text{rate}_{\text{in}} - \text{rate}_{\text{out}})t$, which must go inside the concentration — a trap on exactly the problems that look hardest.

19.5 Direction Fields: Seeing Solutions Without Solving

Most differential equations cannot be solved in closed form. But every first-order equation $y' = F(x,y)$ hands you something for free: at every point $(x,y)$ of the plane it tells you the slope a solution curve must have there. Draw a short segment of that slope at a grid of points and you get a direction field (or slope field) — a picture of all solutions at once, with not a single integral computed.

A solution curve is then any curve that flows tangent to the segments everywhere, like a streamline following arrows in a current. Pick a starting point (an initial condition) and trace the arrows: that traces the solution through that point.

# Direction field for y' = x - y, with one solution curve overlaid.
import numpy as np
import matplotlib.pyplot as plt

X, Y = np.meshgrid(np.linspace(-2, 2, 20), np.linspace(-2, 2, 20))
slope = X - Y                       # y' = F(x, y) = x - y
U = np.ones_like(X)                 # run = 1 in x
V = slope                           # rise = slope
N = np.hypot(U, V)                  # normalize arrows to equal length
plt.quiver(X, Y, U / N, V / N, color="gray", scale=30)
plt.xlabel("x"); plt.ylabel("y")
plt.title("Direction field for y' = x - y")
# The exact solution through (−2, 2): y = x − 1 + 3 e^{−(x+2)}
xs = np.linspace(-2, 2, 200)
plt.plot(xs, xs - 1 + 3 * np.exp(-(xs + 2)), "C0", lw=2, label="y = x − 1 + 3e^{−(x+2)}")
plt.legend(); plt.grid(alpha=0.3); plt.show()
# The curve hugs the arrows everywhere and is pulled toward the line y = x − 1.

The field for $y' = x - y$ reveals its long-term behavior at a glance: away from the diagonal line $y = x - 1$ the arrows steer solutions toward that line (above it slopes are negative, below it positive), so every solution is funneled onto the asymptote $y \approx x - 1$. We learned that without solving anything — the qualitative story lives in the field. For genuinely unsolvable equations this is often all we can know analytically, and it is frequently enough.

Real-World Application — Phase analysis of a thermostat (engineering). A room with a thermostat obeys $T' = -k(T - T_{\text{set}}) + (\text{disturbances})$. Engineers rarely solve this exactly; they read its direction field (and its multidimensional cousin, the phase portrait) to confirm the system is stable — that every nearby trajectory is pulled back to the set point rather than oscillating away. Qualitative ODE analysis, not closed-form solutions, is how control engineers certify that autopilots, HVAC systems, and cruise controls converge.

19.6 Euler's Method: Walking Along the Arrows

A direction field invites an obvious idea: if I know the slope at my current point, let me take a small step in that direction, recompute the slope, and step again. Iterating that idea is a numerical solution. It is Euler's method, the simplest and most important entry point to numerical analysis.

Deriving it. Start from the IVP $y' = F(x,y),\ y(x_0) = y_0$. The derivative is a limit of difference quotients, so for a small step $h$,

$$\frac{y(x+h) - y(x)}{h} \approx y'(x) = F(x, y(x)) \quad\Longrightarrow\quad y(x+h) \approx y(x) + h\,F(x, y(x)).$$

That approximation is exactly the linearization (tangent-line approximation) of Chapter 11: follow the tangent line for one step instead of the true curve. Turn it into an iteration with $x_n = x_0 + nh$ and $y_n \approx y(x_n)$:

$$\boxed{\,y_{n+1} = y_n + h\,F(x_n, y_n), \qquad x_{n+1} = x_n + h.\,}$$

Each step advances along the local slope. The price of simplicity is accuracy: the local error per step is $O(h^2)$ (the discarded second-order term of the Taylor expansion), and these accumulate over the $\sim 1/h$ steps needed to cross a fixed interval into a global error of $O(h)$. Halving the step size roughly halves the error — first-order accuracy, which is slow.

Worked Example 19.6.1 — Euler vs. the exact solution

Solve $y' = -2y,\ y(0) = 1$ on $[0, 5]$ with $h = 0.1$. The exact solution is $y(t) = e^{-2t}$, so we can measure the error.

# Euler's method on y' = −2y, compared against the exact solution e^{−2t}.
import numpy as np
import matplotlib.pyplot as plt
from typing import Callable

def euler(F: Callable[[float, float], float],
          x0: float, y0: float, h: float, n: int):
    xs = np.empty(n + 1); ys = np.empty(n + 1)
    xs[0], ys[0] = x0, y0
    for i in range(n):
        ys[i + 1] = ys[i] + h * F(xs[i], ys[i])   # step along the local slope
        xs[i + 1] = xs[i] + h
    return xs, ys

xs, ys = euler(lambda x, y: -2 * y, 0.0, 1.0, 0.1, 50)
exact = np.exp(-2 * xs)
plt.plot(xs, ys, "o-", ms=3, label="Euler (h = 0.1)")
plt.plot(xs, exact, "r-", label="Exact e^{−2t}")
plt.xlabel("t"); plt.ylabel("y"); plt.legend(); plt.grid(True); plt.show()
print(f"At t = 5:  Euler = {ys[-1]:.6f},  Exact = {exact[-1]:.6f}")
# At t = 5:  Euler = 0.000037,  Exact = 0.000045

Euler slightly undershoots here because the true solution is convex (curving upward), so each straight tangent step lands a hair below the curve. Shrink $h$ and the gap shrinks proportionally — try $h = 0.05$ and watch the error roughly halve.

Computational Note. Plain Euler is a teaching tool, not a production tool. It is too inaccurate (first order) and can go unstable: for $y' = -2y$ with $h > 1$, the factor $1 + hF = 1 - 2h$ exceeds $1$ in magnitude and the numerical solution oscillates and explodes even though the true solution decays smoothly. Real software uses higher-order, adaptive methods. The progression is: Heun's method (improved Euler) averages the slopes at the start and end of a step for $O(h^2)$ accuracy; Runge–Kutta 4 (RK4) blends four slope evaluations per step for $O(h^4)$ accuracy and is the workhorse default; adaptive methods (RK45) adjust $h$ on the fly to hit a target error. In Python you reach for scipy.integrate.solve_ivp, which uses RK45 by default — we use it for the SIR model below.

19.7 The Logistic Equation: Growth With a Ceiling

The exponential model $P' = rP$ predicts a population doubling forever. No real population does — food runs out, predators arrive, crowding spreads disease. The fix, due to Pierre François Verhulst in 1838, is to make the per-capita growth rate fall as the population approaches a ceiling. That ceiling is the carrying capacity $K$, and the model is the logistic equation:

$$\frac{dP}{dt} = rP\left(1 - \frac{P}{K}\right).$$

Read the right side as a story. When $P$ is tiny relative to $K$, the factor $(1 - P/K) \approx 1$ and growth is nearly exponential, $P' \approx rP$. As $P$ climbs toward $K$, that factor shrinks toward $0$ and growth stalls. If $P$ ever exceeds $K$, the factor goes negative and the population declines back toward $K$. The capacity $K$ is a stable equilibrium; $P = 0$ is an unstable one.

Solving it. The equation is separable, and the left side is a partial-fractions exercise straight out of Chapter 16:

$$\frac{dP}{P\left(1 - P/K\right)} = r\,dt, \qquad \frac{1}{P(1 - P/K)} = \frac{1}{P} + \frac{1/K}{1 - P/K}.$$

Integrating each term gives $\ln P - \ln\!\left(1 - \dfrac{P}{K}\right) = rt + C$, hence $\ln\dfrac{P}{1 - P/K} = rt + C$. Exponentiate and solve for $P$:

$$\boxed{\,P(t) = \frac{K}{1 + A e^{-rt}}, \qquad A = \frac{K - P_0}{P_0}.\,}$$

The constant $A$ is fixed by the initial population $P_0$. The graph is the famous S-curve (sigmoid): near-exponential takeoff, an inflection point at $P = K/2$ where growth is fastest, then a graceful leveling at $P = K$.

Geometric Intuition. The logistic curve is the exponential's reflection coming to meet it. Early on it tracks the runaway exponential almost exactly; then the ceiling bends it over. The single most telling feature is the inflection point at half capacity, $P = K/2$: that is where the population grows fastest in absolute terms (you can confirm $\frac{d}{dP}\big[rP(1-P/K)\big] = 0$ at $P = K/2$). Before half-capacity, growth is accelerating; after, it is decelerating. Knowing where the inflection sits lets epidemiologists and product managers alike recognize, in real data, whether the steepest growth is still ahead or already behind.

# Logistic solution vs. its Euler approximation and the runaway exponential.
import numpy as np
import matplotlib.pyplot as plt

r, K, P0 = 0.8, 1000.0, 10.0
t = np.linspace(0, 20, 400)
A = (K - P0) / P0
P_logistic = K / (1 + A * np.exp(-r * t))         # exact solution
P_exponential = P0 * np.exp(r * t)                # Malthusian, for contrast
plt.plot(t, P_logistic, label="Logistic (carrying capacity K = 1000)")
plt.plot(t, np.minimum(P_exponential, 1500), "--", label="Exponential (no ceiling)")
plt.axhline(K, color="gray", ls=":", label="K"); plt.axhline(K / 2, color="C2", ls=":")
plt.xlabel("t"); plt.ylabel("Population"); plt.ylim(0, 1300); plt.legend(); plt.grid(True); plt.show()
# The two curves agree while P ≪ K, then the logistic bends over at P = K/2 = 500 and saturates.

The logistic equation is one of the most widely reused models in all of applied mathematics: population biology (any species in a bounded habitat), tumor growth (Gompertz and logistic models), technology adoption (the Bass diffusion curve), and the early phase of an epidemic before recovery dominates — which is exactly the bridge to the SIR model.

Check Your Understanding. A bacterial culture follows the logistic equation with $r = 1.2$/hr and carrying capacity $K = 10^6$ cells, starting from $P_0 = 100$ cells. Without solving for $t$, find the population at which the culture is growing fastest, and roughly what fraction of capacity that is.

Answer

The growth rate $P' = rP(1 - P/K)$ is a downward parabola in $P$, maximized at the vertex $P = K/2 = 5\times 10^5$ cells — exactly half of carrying capacity, the inflection point of the S-curve. This is independent of $r$ and $P_0$: every logistic process grows fastest at $50\%$ saturation. (The maximum growth rate there is $rK/4 = 1.2 \times 10^6/4 = 3\times 10^5$ cells/hr.)

19.8 Newton's Law of Cooling

A hot object placed in cooler surroundings loses heat at a rate proportional to how much hotter it is than its environment. With $T(t)$ the object's temperature and $T_{\text{env}}$ the (constant) ambient temperature,

$$\frac{dT}{dt} = -k\,(T - T_{\text{env}}), \qquad k > 0.$$

This is both separable and linear; either method (or a substitution $u = T - T_{\text{env}}$ reducing it to $u' = -ku$) gives

$$T(t) = T_{\text{env}} + \big(T_0 - T_{\text{env}}\big)\,e^{-kt}.$$

The temperature relaxes exponentially toward the ambient, fast at first (large temperature gap) and ever slower as it closes in — the gap $T - T_{\text{env}}$ is itself just exponential decay. The same equation describes a warming object (then $T_0 < T_{\text{env}}$ and the exponential climbs); only the sign of the initial gap changes.

Worked Example 19.8.1 — Estimating time of death (forensics)

A body is discovered at 6:00 PM with core temperature $30^\circ$C in a room held at $22^\circ$C. Two hours later, at 8:00 PM, it has cooled to $28^\circ$C. Assuming the body was at the normal $37^\circ$C at the moment of death, estimate when death occurred.

Step 1 — find the cooling constant $k$ from the two measurements. Put the time origin at 6:00 PM, so $T(0) = 30$ and $T(2) = 28$. With $T(t) = 22 + (30-22)e^{-kt} = 22 + 8e^{-kt}$,

$$28 = 22 + 8e^{-2k} \;\Longrightarrow\; e^{-2k} = \frac{6}{8} = \frac34 \;\Longrightarrow\; k = \frac{\ln(4/3)}{2} \approx 0.1438\ \text{hr}^{-1}.$$

Step 2 — re-anchor time to the moment of death. Now reset the clock so $t = 0$ is death, where $T_0 = 37$. The cooling law reads $T(t) = 22 + (37-22)e^{-kt} = 22 + 15e^{-kt}$. We know the body was $30^\circ$C at the discovery; let $t_1$ be the elapsed time from death to discovery:

$$30 = 22 + 15\,e^{-k t_1} \;\Longrightarrow\; e^{-k t_1} = \frac{8}{15} \;\Longrightarrow\; t_1 = \frac{\ln(15/8)}{k} = \frac{\ln(1.875)}{0.1438} \approx 4.37\ \text{hr}.$$

So death occurred about $4.4$ hours before the 6:00 PM discovery — roughly 1:38 PM. (Real forensic estimates use multi-compartment cooling models because a body is not a uniform lump, but the calculus is exactly this. Newton's law of cooling also dates the cooling of coffee, the quenching of steel, and the thermal lag of buildings — the same exponential approach to equilibrium each time.)

19.9 The SIR Model of Epidemics

We now develop the anchor example of this chapter and one of the most consequential differential equation models of the modern era. The SIR model, introduced by Kermack and McKendrick in 1927, partitions a population into three compartments and tracks the flow between them. It is the mathematical skeleton beneath every epidemic forecast you saw in 2020, and it returns in full as the Chapter 39 capstone, where you will calibrate it to data and explore interventions.

The compartments. Let $N$ be the total population, split at each instant into:

• $S(t)$ — susceptible: never infected, can still catch the disease;
• $I(t)$ — infectious: currently infected and able to transmit;
• $R(t)$ — recovered/removed: immune or deceased, out of the chain of transmission.

The dynamics. Two processes move people between compartments. Infection converts susceptibles to infectious at a rate proportional to the number of contacts between the two groups, $\beta S I$ (the mass-action assumption: encounters are proportional to the product of the populations, just like reacting chemicals). Recovery converts infectious to recovered at a constant per-capita rate $\gamma$. The system is

$$\frac{dS}{dt} = -\beta S I, \qquad \frac{dI}{dt} = \beta S I - \gamma I, \qquad \frac{dR}{dt} = \gamma I.$$

Add the three equations and the right side telescopes to zero: $\dfrac{d}{dt}(S + I + R) = 0$, so $S + I + R = N$ is conserved — nobody leaves the population, they only change status. The two parameters carry the epidemiology: $\beta$ is the transmission rate per susceptible–infectious encounter, and $\gamma$ is the recovery rate, with $1/\gamma$ the average infectious period.

This is a nonlinear system — the $\beta S I$ terms are products of unknowns — and it has no closed-form solution. We solve it numerically, which is precisely why we built Euler's method and why solve_ivp exists.

The basic reproduction number $R_0$

Before computing anything, ask the question that decides whether an outbreak even happens: how many new infections does one infectious person cause in an otherwise fully susceptible population? At the start of an outbreak $S \approx N$, so the infectious equation reads

$$\frac{dI}{dt} = \beta S I - \gamma I \approx (\beta N - \gamma)\,I.$$

That is exponential growth or decay depending on the sign of $\beta N - \gamma$. The infection takes off when $\beta N > \gamma$ and fizzles when $\beta N < \gamma$. The dimensionless ratio at the crossover is the basic reproduction number:

$$\boxed{\,R_0 = \frac{\beta N}{\gamma}.\,}$$

$R_0 > 1$ means each case more than replaces itself and the epidemic grows; $R_0 < 1$ means it dies out. This one number, distilled from a differential equation, is the most quoted quantity in all of epidemiology:

Real-World Application — Vaccination policy and herd immunity (public health). The same $R_0$ that predicts whether an epidemic grows also dictates how many people must be immune to stop it. If a fraction $p$ of the population is immune, the effective reproduction number is $R_0(1-p)$; transmission collapses once $R_0(1-p) < 1$, i.e. once $p > 1 - \dfrac{1}{R_0}$. This herd immunity threshold is why measles ($R_0 \approx 15$) requires $\approx 1 - 1/15 \approx 93\%$ vaccination coverage to eliminate, while a milder disease with $R_0 = 1.3$ needs only $\approx 23\%$. A single ODE-derived inequality sets vaccination targets for entire nations.

Solving the SIR system numerically

We integrate the system with scipy.integrate.solve_ivp (RK45 under the hood — the production-grade successor to the Euler method we hand-coded in §19.6).

# SIR epidemic model solved with an adaptive Runge–Kutta integrator.
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

def sir(t, y, beta, gamma):
    S, I, R = y
    dS = -beta * S * I
    dI = beta * S * I - gamma * I
    dR = gamma * I
    return [dS, dI, dR]

N = 1000
y0 = [N - 1, 1, 0]                 # one infectious individual seeds the outbreak
gamma = 0.1                        # average infectious period 1/gamma = 10 days
R0_target = 3.0
beta = R0_target * gamma / N       # so that R0 = beta*N/gamma = 3.0  → beta = 3e-4

sol = solve_ivp(sir, [0, 160], y0, args=(beta, gamma),
                dense_output=True, max_step=1.0)
t = np.linspace(0, 160, 400)
S, I, R = sol.sol(t)

plt.figure(figsize=(9, 5))
plt.plot(t, S, label="Susceptible")
plt.plot(t, I, label="Infectious")
plt.plot(t, R, label="Recovered")
plt.xlabel("Days"); plt.ylabel("People"); plt.title(f"SIR dynamics  (R0 = {R0_target})")
plt.legend(); plt.grid(True); plt.show()

peak_day = t[np.argmax(I)]
print(f"Peak infectious: {I.max():.0f} people on day {peak_day:.0f}")
print(f"Final susceptible (never infected): {S[-1]:.0f} of {N}")
# Peak infectious: ~300 people around day 40; final susceptible ~60 — about 94% were infected.

The output shows the signature epidemic curve: $I(t)$ climbs roughly exponentially while susceptibles are plentiful, peaks when $S$ falls to the threshold $S = \gamma/\beta = N/R_0$ (the point where $dI/dt = 0$), then declines as the susceptible pool runs dry. Crucially, the epidemic ends with susceptibles to spare — not everyone gets infected, because transmission stalls once $S$ drops below $N/R_0$. "Flattening the curve," the slogan of 2020, means reducing $\beta$ (through distancing and masks) to lower and delay that peak so it never overwhelms hospital capacity.

Common Pitfall. It is tempting to read $R_0 > 1$ as "everyone eventually gets infected." Not so. The SIR model predicts a final size strictly less than $100\%$: the outbreak burns out when susceptibles drop below $N/R_0$, leaving a residual susceptible fraction. The larger $R_0$, the smaller that residue — but it is never zero. Confusing the threshold for growth ($R_0 > 1$) with the final attack rate is the most common misreading of epidemic models.

Extensions you will meet again. The SIR skeleton supports a whole family of richer models: SEIR inserts an exposed (latent, not-yet-infectious) compartment; age-structured models give each age group its own parameters; time-varying $\beta(t)$ encodes lockdowns and seasonality; spatial models couple many local SIR systems across a transport network. Every one of them is built from the same differential-equation machinery in this chapter, integrated by the same numerical methods. In Chapter 39 you will take this exact model, fit $\beta$ and $\gamma$ to real case data, simulate interventions, and write up the result as the centerpiece of your modeling portfolio.

19.10 A Gallery of ODE Models

The reach of first- and second-order ODEs across the sciences is hard to overstate. Four more classics, each a single equation that captures a whole phenomenon.

Predator and prey — the Lotka–Volterra system. Two coupled nonlinear equations model a predator population $y$ feeding on prey $x$:

$$\frac{dx}{dt} = \alpha x - \beta x y, \qquad \frac{dy}{dt} = \delta x y - \gamma y.$$

Prey grow ($\alpha x$) but are eaten ($-\beta x y$); predators die off ($-\gamma y$) but thrive on prey ($\delta x y$). The solutions are closed loops — perpetual oscillations in which prey booms feed predator booms that crash the prey that starve the predators that let prey rebound. Developed independently by Lotka (1925) and Volterra (1926), it reproduced the puzzling population cycles of Canadian lynx and snowshoe hare recorded in fur-trapping records.

The harmonic oscillator — Newton's second law. A mass $m$ on a spring of stiffness $k$ obeys $F = ma$ with the restoring force $F = -kx$:

$$m\ddot{x} = -kx \quad\Longrightarrow\quad \ddot{x} + \omega^2 x = 0, \qquad \omega = \sqrt{k/m}.$$

This second-order linear equation has solution $x(t) = A\cos(\omega t + \phi)$ — simple harmonic motion, the eternal sinusoid behind pendulums, AC circuits, sound waves, and the orbit of every vibrating molecule. Adding a damping term $c\dot x$ produces decaying oscillations; adding a forcing term produces resonance.

The RC circuit. A capacitor charging through a resistor satisfies $\dot V + \dfrac{V}{RC} = \dfrac{V_{\text{in}}}{RC}$ — first-order linear, identical in form to Newton's cooling law (built on in the Chapter 16 case study). The capacitor voltage relaxes exponentially toward $V_{\text{in}}$ with time constant $RC$.

Continuous compound interest. Money growing at instantaneous rate $r$ obeys $\dfrac{dM}{dt} = rM$, giving $M(t) = M_0 e^{rt}$ — the exponential equation of §19.3 again, now denominated in dollars (the Chapter 7 case study computed $e$ from exactly this limit).

Add to Your Modeling Portfolio. This chapter contributes a dynamical model to your portfolio — a differential equation that evolves your system in time, plus a numerical solver to run it. Biology: implement the SIR model (or logistic growth) for your population, choose biologically plausible $\beta, \gamma$ (or $r, K$), compute $R_0$ (or the carrying capacity), and simulate with solve_ivp. Save this — Chapter 39 builds its capstone directly on it. Economics: model a market with continuous compounding or a logistic product-adoption curve; estimate the growth rate and saturation level, and simulate the path to equilibrium. Physics: model your system's motion with Newton's second law (oscillator, projectile with drag, or cooling), write it as an ODE, and integrate it numerically; compare to the analytic solution where one exists. Data Science: treat gradient descent as the discrete Euler integration of the gradient flow ODE $\dot{\theta} = -\nabla L(\theta)$ — the learning rate is the step size $h$. Implement it and show that too large a step makes Euler (and training) diverge, the §19.6 stability lesson in disguise.

19.11 Summary and the Road Ahead

Differential equations are how the sciences write down the laws of change, and you now have the core toolkit to read and solve them.

• Separable equations $y' = g(x)h(y)$: separate the variables, integrate both sides, apply the initial condition.
• First-order linear equations $y' + P(x)y = Q(x)$: multiply by the integrating factor $\mu = e^{\int P\,dx}$ to collapse the left side into $(\mu y)'$, then integrate.
• Direction fields: visualize all solutions at once without solving; read off stability and long-term behavior.
• Euler's method $y_{n+1} = y_n + hF(x_n, y_n)$: the linear-approximation idea iterated into a numerical solver; first-order accurate, the gateway to RK4 and solve_ivp.
• Exponential growth/decay $y' = ky \Rightarrow y = y_0 e^{kt}$: half-lives, drug clearance, compound interest.
• Logistic growth $P' = rP(1-P/K) \Rightarrow$ the S-curve saturating at $K$, fastest at $K/2$.
• Newton's cooling $T' = -k(T - T_{\text{env}})$: exponential approach to ambient; forensic time-of-death.
• The SIR model $S' = -\beta S I,\ I' = \beta S I - \gamma I,\ R' = \gamma I$ with $R_0 = \beta N/\gamma$: the epidemic anchor, climaxing in Chapter 39.

Three threads of this book run straight through the chapter. Calculus is the mathematics of change in its purest form here — we have been doing nothing but turning rate laws into the quantities they govern. Hand computation built the understanding (you separated, integrated, and verified by substitution) while machine computation delivered the power (Euler and solve_ivp cracked the nonlinear SIR system no formula can touch). And every solution rested on the Fundamental Theorem of Calculus (Chapter 14): solving a differential equation is, at bottom, antidifferentiation with feedback.

This concludes Part III: Integration. Across seven chapters you built the definite integral, proved the Fundamental Theorem, mastered the techniques of integration, tamed improper integrals, applied integration to geometry and physics, and now solved the equations that govern dynamic systems. Part IV turns to sequences and series (Chapters 20–24): the art of adding infinitely many numbers, Taylor series, and the route to Euler's identity $e^{i\pi} + 1 = 0$. From there, Parts V through VIII generalize everything you know to curves, surfaces, and several variables, culminating in Stokes' theorem (Chapter 37) and the universal Fundamental Theorem (Chapter 38).

Differential equations are the corner of calculus that most directly shapes the modern world — from epidemic forecasts to climate models to the training dynamics of neural networks. You have just learned to speak nature's native language. Turn the page and learn to add infinitely many numbers.

Continue to: Exercises · Quiz · Case Study 1 · Case Study 2 · Key Takeaways · Further Reading