Part IV · Calculus Chapter 31 of 60
Differential equations
The laws of nature rarely say how much of something there will be; they say how fast it changes. Welcome to a lab where processes are grown from such laws: cooling coffee, a spring, wolves and hares, an epidemic. And we’ll find out why a uniquely determined future isn’t the same as a predictable one.
Builds on: 30 · Series
You will learn
- read a differential equation as a law of change and see its solutions in a direction field
- solve equations numerically by Euler’s method and understand where the error comes from
- solve second-order linear equations through the characteristic equation and explain the epidemic threshold through R₀
The last chapter ended with an observation: the laws of nature say not “how much” but “how it changes”. A cup of coffee cools the faster the hotter it is. A pendulum accelerates the harder the further it’s displaced. An epidemic grows the faster the more people are ill and the more there are who haven’t been ill yet. In each such law the unknown function stands next to its own derivative, and simply integrating the rate, as in chapter 28, doesn’t work: the rate depends on what we’re looking for.
This chapter is a lab. Each bench has its own model and its own knobs: cooling and growth, a spring, wolves and hares, an epidemic. At each we’ll ask “what happens if…”, look at the answer and explain it. The lab rules, without which all the experiments would be meaningless, are written up at the end, together with an experiment that almost breaks them.
A law instead of an answer
Let’s start with coffee. In 1701 an unsigned note appeared in the Philosophical Transactions of the Royal Society, later recognised as Newton’s work: a body cools at a rate proportional to the difference between its temperature and the temperature of the air. If the room is at $20°$ and the coffee’s temperature after $t$ minutes is $T(t)$, the law is written like this:
$$T'(t) = -k\,\bigl(T(t) - 20\bigr), \qquad k > 0.$$The minus is there because hot coffee cools; $k$ depends on the cup: for thin porcelain it’s larger, for a vacuum flask it’s tiny. The unknown here isn’t a number but a whole function $T(t)$, and the equation links it to its own derivative.
A differential equation is an equation in which the unknown is a function and its derivatives appear. The highest order of derivative is called the order of the equation. A solution is a function that turns the equation into a true equality for all $t$ in some interval. An equation together with an initial condition $y(t_0) = y_0$ is called an initial value problem, or a Cauchy problem; for a second-order equation the initial velocity $y'(t_0)$ is added as well.
Here the unknown function depends on one variable, and such equations are called ordinary. Equations where the function depends on several variables, like temperature depending on time and the point in the room, are called partial differential equations, and they’re harder. We’ve already met the simplest ordinary equation in chapter 29: the rate of growth is proportional to the quantity itself.
A function $y(t)$ defined on an interval satisfies the equation $y' = ky$ if and only if $y = Ce^{kt}$ for some constant $C$. The initial value problem $y(t_0) = y_0$ has exactly one solution: $y = y_0e^{k(t - t_0)}$.
The functions $Ce^{kt}$ fit: $(Ce^{kt})' = kCe^{kt}$. Conversely, let $y' = ky$. The trick is to divide $y$ by the expected answer $e^{kt}$ and check that the quotient doesn’t change. By the product rule $\left(y e^{-kt}\right)' = y'e^{-kt} - ky\,e^{-kt} = (y' - ky)\,e^{-kt} = 0$. A function with zero derivative on an interval is constant (a consequence of Lagrange’s theorem from chapter 27), so $y e^{-kt} = C$ and $y = Ce^{kt}$. The condition $y(t_0) = y_0$ gives $C = y_0e^{-kt_0}$ and allows no other $C$.
Bench one: the direction field
Most equations of the form $y' = f(t, y)$ can’t be solved by a formula. But even without a formula, the equation knows everything about its solutions. At each point $(t, y)$ it reports the slope with which a solution passes through that point. Draw short dashes with these slopes at the nodes of a grid, and you get a picture of a flow, with the solutions as streamlines.
The direction field of the equation $y' = f(t, y)$ is a collection of short segments: at each point $(t, y)$, a segment with slope $f(t, y)$. The graph of a solution is called an integral curve; at each of its points it touches a dash of the field.
For the equation $y' = t - y$ all the integral curves hug the line $y = t - 1$. It’s a solution itself: its slope is $1$, and $t - (t - 1) = 1$. The other solutions are $y = t - 1 + Ce^{-t}$ (check by substituting), and the extra $Ce^{-t}$ melts away. For the logistic equation $y' = y(1 - y)$ you can see two horizontal lines, $y = 0$ and $y = 1$: the slope on them is zero, and a solution that starts there goes nowhere. Such constant solutions are called equilibria. Neighbouring curves are drawn in towards the line $y = 1$ and run away from the line $y = 0$: one equilibrium is stable, the other isn’t.
Bench two: step by step
The field suggests how to find a solution without a formula: stand at the starting point, go a little way along the dash, look at the new dash and step again. Leonhard Euler described this recipe in his Institutions of Integral Calculus (1768).
Euler’s method is an approximate solution of the initial value problem $y' = f(t, y)$, $y(t_0) = y_0$ in steps of length $h$: each step is taken along the tangent given by the field at the current point.
Euler’s method makes mistakes: it follows the tangent, while the solution bends. How big the mistakes are can be seen from the same example. At $n = 10$ the error is $0.125$, at $n = 100$ it’s $0.0135$, at $n = 1000$ it’s $0.00136$: a step ten times smaller, an error ten times smaller. The reason lies in the Taylor series from the last chapter. Over one step the tangent departs from the curve by roughly $\frac{y''h^2}{2}$, and there are $\frac1h$ steps to the target, so an error of order $h$ builds up. For $y' = y$ this can be checked exactly: $\ln\left(1 + h\right)^{1/h} = \frac{1}{h}\left(h - \frac{h^2}{2} + \dots\right) = 1 - \frac h2 + \dots$, so $\left(1 + h\right)^{1/h} \approx e\left(1 - \frac h2\right)$, and the error is close to $\frac{e}{2}h \approx 1.36h$.
Methods whose error falls faster take several slopes at each step and average them. The best known of them, the classical Runge–Kutta method (Carl Runge, 1895; Martin Kutta, 1901), gives $e$ in ten steps with an error of $2 \cdot 10^{-6}$, and halving the step cuts the error roughly $16$-fold. Variants of it solve the equations in almost all the widgets in this chapter.
The last example is worth remembering. For the spring $y'' = -y$, Euler’s method slightly increases the energy at every step, and instead of everlasting oscillations you get a spiral winding outwards: the amplitude grows by a factor of $\sqrt{1 + h^2}$ per step. A numerical method can violate conservation laws that the exact solution respects, which is why special methods that don’t pump in energy have been devised for long calculations of planetary orbits.
Bench three: growth with a ceiling
The law $y' = ky$ with $k > 0$ promises endless growth, but any population runs out of food and space. In 1838 the Belgian mathematician Pierre Verhulst proposed braking the growth with a factor that becomes zero when the population reaches its limit.
The equation $y' = ry\left(1 - \frac yK\right)$ with $r, K > 0$ is called the logistic equation, and the graphs of its solutions are called logistic curves.
If $0 < y_0 < K$, the solution of the initial value problem $y' = ry\left(1 - \frac yK\right)$, $y(0) = y_0$ is unique and equals
$$y(t) = \frac{K}{1 + \left(\frac{K}{y_0} - 1\right)e^{-rt}}.$$The trick is to turn the unknown upside down: for $u = \frac1y$ the equation becomes linear. Let $y$ be a solution, and for now let $y > 0$. Then $u' = -\frac{y'}{y^2} = -\frac{r}{y}\left(1 - \frac{y}{K}\right) = -ru + \frac rK$. For $v = u - \frac1K$ this gives $v' = u' = -r\left(u - \frac1K\right) = -rv$, and by the theorem on exponential growth $v = \left(\frac1{y_0} - \frac1K\right)e^{-rt}$. Going back: $y = \frac1u = \frac{1}{\frac1K + \left(\frac1{y_0} - \frac1K\right)e^{-rt}}$; multiplying numerator and denominator by $K$, we get the formula in the statement. All the steps are reversible: for this function $u = \frac1y$ satisfies the equation $u' = -ru + \frac rK$, and so the function itself satisfies the logistic one. The denominator of the formula is positive and finite for all $t$, so $y$ never becomes zero, and the proviso “for now let $y > 0$” always holds. So any solution with the same initial condition coincides with the one we found.
The curve looks like a stretched letter S. At first the growth is almost exponential, then it slows down and hugs $K$. The population grows fastest at $y = \frac K2$: the right-hand side $ry\left(1 - \frac yK\right)$ is a parabola in $y$ with its vertex in the middle. Logistic curves do a decent job of describing both the spread of new technologies and the growth of bacterial colonies, and a discrete relative of this equation, the logistic map, will open the door to chaos in chapter 59. The field of this equation is on bench one.
Bench four: the spring
A mass $m$ on a spring. The spring pulls it towards equilibrium with a force $-\kappa x$, where $x$ is the displacement and $\kappa$ the stiffness (Hooke’s law). Air resistance brakes it with a force $-cx'$. By Newton’s second law $mx'' = -\kappa x - cx'$. Divide by $m$ and write $\frac{\kappa}{m} = \omega^2$, $\frac{c}{m} = 2\gamma$:
$$x'' + 2\gamma x' + \omega^2 x = 0.$$This is a second-order equation, and an initial value problem needs two numbers: where the mass starts and how fast it’s moving. Without friction, at $\gamma = 0$, the solutions are familiar from mechanics: $x = A\cos\omega t + B\sin\omega t$, oscillations with period $\frac{2\pi}{\omega}$. A $0.5$ kg mass on a spring of stiffness $50$ N/m gives $\omega = 10$ rad/s and a period of about $0.63$ s. And with friction? Euler found a trick that works for all linear equations with constant coefficients: look for a solution of the form $x = e^{\lambda t}$. Substitute: $\lambda^2 e^{\lambda t} + 2\gamma\lambda e^{\lambda t} + \omega^2 e^{\lambda t} = 0$, and after dividing by $e^{\lambda t}$ a quadratic equation is left.
For the equation $x'' + px' + qx = 0$ with constant coefficients, the equation $\lambda^2 + p\lambda + q = 0$ is called its characteristic equation. It’s the characteristic polynomial from chapter 38: if the equation is written as a system for $x$ and $v = x'$, its matrix $\left(\begin{smallmatrix} 0 & 1 \\ -q & -p \end{smallmatrix}\right)$ has exactly this polynomial.
Let $\lambda_1$, $\lambda_2$ be the roots of the characteristic equation for $x'' + px' + qx = 0$ with real $p$ and $q$. All the real solutions of this equation are as follows:
- $x = Ae^{\lambda_1 t} + Be^{\lambda_2 t}$ if the roots are real and distinct;
- $x = (A + Bt)\,e^{\lambda t}$ if the roots coincide: $\lambda_1 = \lambda_2 = \lambda$;
- $x = e^{\alpha t}(A\cos\beta t + B\sin\beta t)$ if the roots are complex: $\lambda_{1,2} = \alpha \pm i\beta$, $\beta \ne 0$.
Here $A$ and $B$ are any real numbers, and for any $x(0)$, $x'(0)$ they’re determined uniquely.
The idea: break the second-order equation into two first-order equations, each of which we already know how to solve.
By Vieta’s formulas (chapter 10) $p = -(\lambda_1 + \lambda_2)$ and $q = \lambda_1\lambda_2$. So $x'' + px' + qx = x'' - (\lambda_1 + \lambda_2)x' + \lambda_1\lambda_2 x$. If we write $w = x' - \lambda_2 x$, then $w' - \lambda_1 w = (x'' - \lambda_2 x') - \lambda_1(x' - \lambda_2 x)$, exactly the same expression. So for a solution $x$ the function $w$ satisfies the equation $w' = \lambda_1 w$, and by the theorem on exponential growth $w = Ce^{\lambda_1 t}$. If the roots are complex, the functions have to be complex-valued too. The proof of that theorem works for them as well: by Euler’s formula from chapter 29 $e^{(a + ib)t} = e^{at}(\cos bt + i\sin bt)$, and differentiating directly gives $\left(e^{\lambda t}\right)' = \lambda e^{\lambda t}$ for complex $\lambda$.
Now solve $x' - \lambda_2 x = Ce^{\lambda_1 t}$. Multiply by $e^{-\lambda_2 t}$: on the left we get the derivative of a product, $\left(xe^{-\lambda_2 t}\right)' = Ce^{(\lambda_1 - \lambda_2)t}$. If $\lambda_1 \ne \lambda_2$, an antiderivative of the right side is $\frac{C}{\lambda_1 - \lambda_2}e^{(\lambda_1 - \lambda_2)t} + D$, and $x = Ae^{\lambda_1 t} + De^{\lambda_2 t}$. If $\lambda_1 = \lambda_2 = \lambda$, the right side is just $C$, and $xe^{-\lambda t} = Ct + D$, that is, $x = (D + Ct)e^{\lambda t}$. Conversely, all such functions are solutions: it’s enough to run these calculations backwards.
With complex roots $\lambda_2 = \bar\lambda_1$, and a real solution $x$ coincides with its real part. The real part of the function $Ae^{(\alpha + i\beta)t}$ is $e^{\alpha t}(a\cos\beta t + b\sin\beta t)$ with some real $a$, $b$, and the second term has the same form. Finally, the constants: with distinct roots $x(0) = A + B$, $x'(0) = \lambda_1 A + \lambda_2 B$, a system with determinant $\lambda_2 - \lambda_1 \ne 0$, which has a unique solution; with a double root $x(0) = A$, $x'(0) = \lambda A + B$.
For the spring these are three regimes. For $\gamma < \omega$ the roots are complex, and the mass oscillates with frequency $\sqrt{\omega^2 - \gamma^2}$ inside the decaying envelope $e^{-\gamma t}$. For $\gamma > \omega$ the roots are real and negative: released from a displaced position, the mass slides back to equilibrium without ever crossing it. The boundary, $\gamma = \omega$, is critical damping, with the double root $-\omega$ and the solution $(A + Bt)e^{-\omega t}$. That’s the very factor $t$ that chapter 38 talked about: with a double root the matrix of the system isn’t diagonalisable. With critical damping the mass returns fastest without swinging, which is why needle gauges and door closers are tuned close to this regime.
Find the solution of the initial value problem $x'' + 2x' + 5x = 0$, $x(0) = 0$, $x'(0) = 2$.
The roots of the characteristic equation $\lambda^2 + 2\lambda + 5 = 0$ are $-1 \pm 2i$, so $x = e^{-t}(A\cos 2t + B\sin 2t)$. From $x(0) = 0$ we get $A = 0$. Then $x' = e^{-t}(-B\sin 2t + 2B\cos 2t)$, and $x'(0) = 2B = 2$, $B = 1$. Answer: $x = e^{-t}\sin 2t$, oscillations with period $\pi$ inside the envelope $\pm e^{-t}$.
You can get your hand in with these equations, from growth and decay to the initial value problem for a spring, in the trainer.
Bench five: wolves and hares
After the First World War the Italian biologist Umberto D’Ancona was going through the statistics of Adriatic fish markets and noticed something odd. In the war years, when little fishing was done, the share of predatory fish in the catches had grown. Why had the slackening of fishing helped the predators more than their prey? D’Ancona asked his father-in-law, the mathematician Vito Volterra, and in 1926 Volterra built a model. A year earlier the American Alfred Lotka had published the same model independently.
Two unknowns are handy to draw not as two graphs but as a single point $(x, y)$ moving across the plane over time.
The phase plane of a system of two equations is the plane whose points are the states $(x, y)$. A solution is shown as a moving point, its path is called a phase trajectory, and the picture of all trajectories a phase portrait. A state in which all the rates are zero is called an equilibrium point.
The predator–prey trajectories close up: the hares breed, the wolves fatten up, the hares get fewer, the wolves go hungry, and it all repeats. That they close isn’t an impression from the picture but a theorem.
For any solution of the Lotka–Volterra system with $x, y > 0$, the quantity $V(x, y) = \delta x - \gamma\ln x + \beta y - \alpha\ln y$ doesn’t change over time.
Differentiate $V$ along a solution by the chain rule and substitute the equations: $V' = \left(\delta - \frac{\gamma}{x}\right)x' + \left(\beta - \frac{\alpha}{y}\right)y' = \left(\delta - \frac\gamma x\right)x(\alpha - \beta y) + \left(\beta - \frac\alpha y\right)y(\delta x - \gamma)$. Multiply out: $(\delta x - \gamma)(\alpha - \beta y) + (\beta y - \alpha)(\delta x - \gamma) = (\delta x - \gamma)\bigl((\alpha - \beta y) + (\beta y - \alpha)\bigr) = 0$. The derivative is zero on an interval, so $V$ is constant.
The function $\delta x - \gamma\ln x$ decreases up to $x = \frac\gamma\delta$ and then increases, and $\beta y - \alpha\ln y$ behaves the same way around $y = \frac\alpha\beta$. So $V$ is shaped like a bowl with its bottom at the equilibrium, and its level curves are closed curves around the bottom. A trajectory has to stay on one such curve, and the point goes round and round it. Near the equilibrium the period is close to $\frac{2\pi}{\sqrt{\alpha\gamma}}$; in our example that’s $\frac{2\pi}{0.6} \approx 10.5$ years. Similar ten-year cycles show up in the famous data on lynx and snowshoe hare pelts bought by the Hudson’s Bay Company in Canada, although real ecology is of course more complicated than two equations.
The average values of $x$ and $y$ over a full cycle equal the coordinates of the equilibrium: $\bar x = \frac\gamma\delta$, $\bar y = \frac\alpha\beta$.
Divide the first equation by $x$: $(\ln x)' = \frac{x'}{x} = \alpha - \beta y$. Integrate over a cycle of length $T$: on the left $\ln x(T) - \ln x(0) = 0$, because over a cycle the point has come back, and on the right $\alpha T - \beta\int_0^T y\,dt$. So the average $\bar y = \frac1T\int_0^T y\,dt = \frac\alpha\beta$. In exactly the same way, $(\ln y)' = \delta x - \gamma$ gives $\bar x = \frac\gamma\delta$.
Now D’Ancona’s riddle is solved in two lines. Fishing takes the same fraction $f$ per unit time from both populations: $\alpha$ turns into $\alpha - f$, and $\gamma$ into $\gamma + f$. The average number of prey, $\frac{\gamma + f}{\delta}$, grows with fishing, and the average number of predators, $\frac{\alpha - f}{\beta}$, falls. The war cut fishing, $f$ decreased, and there were more predators. This conclusion is called Volterra’s principle.
The same bench has a pendulum. Its equation $\theta'' = -\frac gL\sin\theta$ is non-linear, but for small angles $\sin\theta \approx \theta$ (the series from the last chapter), and we get a spring with $\omega = \sqrt{g/L}$ and period $2\pi\sqrt{L/g}$. For a one-metre pendulum that’s $2.006$ s. For large swings the approximation gets worse: at a displacement of $90°$ the period is $18\,\%$ longer. The phase portrait of the pendulum shows both oscillations (closed loops) and rotation over the top (wavy lines), and they’re separated by a special trajectory leading to the upside-down position. The loops close for the same reason as in Volterra’s case: the quantity $\frac{(\theta')^2}{2} - \frac gL\cos\theta$, that is, the energy per unit mass divided by $L^2$, doesn’t change along a solution.
Bench six: an epidemic
In 1927 William Kermack and Anderson McKendrick proposed a model that’s still used today. The population is divided into three groups: the susceptible $S$, who haven’t been ill yet, the infectious $I$, and the removed $R$, who have recovered and can’t be infected again. All three quantities are fractions of the population, $S + I + R = 1$.
The basic reproduction number $R_0 = \frac\beta\gamma$ is how many people one infected person infects on average over the course of the illness if everyone around is susceptible.
The number of infectious people grows at a moment when $R_0 S > 1$ and falls when $R_0 S < 1$. In particular, an outbreak starts only if $R_0 S(0) > 1$.
Take $I$ out of the second equation: $I' = I(\beta S - \gamma) = \gamma I\,(R_0 S - 1)$. The factor $\gamma I$ is positive as long as there are infectious people, and the sign of $I'$ is the sign of $R_0 S - 1$. And since $S' = -\beta SI \le 0$, the share of susceptibles only decreases: if $R_0 S < 1$ at the start, it stays that way, and the number of infectious people falls from day one.
From this comes the main formula of vaccination. If a fraction $p$ of the population is vaccinated in advance, then $S(0) \approx 1 - p$, and there’ll be no outbreak if $R_0(1 - p) < 1$, that is, if $p > 1 - \frac1{R_0}$. For $R_0 = 2.5$ that’s $60\,\%$. For measles $R_0$ is usually estimated at $12$–$18$, and the herd immunity threshold comes out at about $92$–$94\,\%$. But if an epidemic is already under way, it doesn’t stop at the threshold: how many people end up falling ill is answered by the next theorem.
Let $R(0) = 0$, $S(0) = s_0$ and $I(0) = 1 - s_0 > 0$ at the start. Then $S(t)$ tends to a number $s_\infty > 0$ for which $s_\infty = s_0\,e^{-R_0(1 - s_\infty)}$, and $I(t) \to 0$.
The idea: eliminate time and follow the epidemic as a point in the $(S, I)$ plane.
For a disease with $R_0 = 4$, what fraction of the population has to be vaccinated in advance for an outbreak not to start? Give the answer as a fraction.
We need $R_0(1 - p) < 1$, that is, $1 - p < \frac14$ and $p > \frac34$. The threshold is $\frac34$, or $75\,\%$ of the population.
The lab rules
So far every initial value problem has had exactly one solution. That’s the lab’s unspoken rule: run an experiment twice from the same initial conditions and you get the same result. Is it always so? No.
Take the equation $y' = 3\sqrt[3]{y^2}$ with the condition $y(0) = 0$. It’s satisfied by $y = 0$, and by $y = t^3$ too: $(t^3)' = 3t^2 = 3\sqrt[3]{t^6}$. What’s more, a function that sits at zero for as long as it likes and then, at $t = c$, sets off along the curve $(t - c)^3$ fits as well. The solution isn’t unique, and a lab with such a law is unpredictable in the most literal sense. Another trouble: for $y' = y^2$ with $y(0) = 1$ the solution $y = \frac{1}{1 - t}$ runs off to infinity at $t = 1$, although the right-hand side is a harmless parabola. The solution exists only on a finite interval. It all depends on how sharply $f(t, y)$ can change as $y$ changes.
A function $f(t, y)$ satisfies a Lipschitz condition in $y$ with constant $K$ if $|f(t, y_1) - f(t, y_2)| \le K|y_1 - y_2|$ for all the $t$, $y_1$, $y_2$ under consideration. This holds, for example, if $f$ has a derivative in $y$ bounded in absolute value by $K$: that’s Lagrange’s theorem. For $3\sqrt[3]{y^2}$ the slope near zero is infinite, and there’s no Lipschitz condition.
Let the function $f(t, y)$ be continuous on the rectangle $|t - t_0| \le a$, $|y - y_0| \le b$ and satisfy a Lipschitz condition in $y$ there. Then on some interval $|t - t_0| \le h$ the initial value problem $y' = f(t, y)$, $y(t_0) = y_0$ has a solution, and only one.
The idea: replace the equation by an equivalent integral one and get the solution by successive approximations, each substituted into the right-hand side. A continuous function $y$ solves the initial value problem if and only if $y(t) = y_0 + \int_{t_0}^{t} f(s, y(s))\,ds$: one way it’s the Newton–Leibniz formula, the other way it’s differentiating an integral with respect to its upper limit (chapter 28). Let $|f| \le M$ on the rectangle; take $h = \min\left(a, \frac bM\right)$ and consider $t \ge t_0$ (to the left everything is the same). Any solution with $t_0 \le t \le t_0 + h$ doesn’t leave the rectangle: while it’s inside, $|y'| \le M$, and in time $h$ it can’t get further than $Mh \le b$ from $y_0$.
Uniqueness. Let $y$ and $z$ be two solutions, $u(t) = |y(t) - z(t)|$, and $U$ the largest value of $u$ on $[t_0, t_0 + h]$. Subtract the integral equations and apply the Lipschitz condition: $u(t) \le K\int_{t_0}^t u(s)\,ds$. Hence $u(t) \le KU(t - t_0)$; substitute this into the integral again: $u(t) \le K^2U\frac{(t - t_0)^2}{2}$, and by induction $u(t) \le U\frac{(K(t - t_0))^n}{n!}$ for every $n$. The right side is the general term of the convergent series for $Ue^{K(t - t_0)}$, and it tends to zero (chapter 30). So $u = 0$: the solutions coincide.
Existence. Put $y_0(t) = y_0$ and $y_{n+1}(t) = y_0 + \int_{t_0}^t f(s, y_n(s))\,ds$. From $|f| \le M$ and the choice of $h$ it follows that all the approximations stay in the rectangle: $|y_{n+1}(t) - y_0| \le M(t - t_0) \le b$. To begin with, $|y_1(t) - y_0(t)| \le M(t - t_0)$, and after that each step adds a factor $K$ by the Lipschitz condition, plus an integral: $|y_{n+1}(t) - y_n(t)| \le K\int_{t_0}^t |y_n(s) - y_{n-1}(s)|\,ds$. By induction $|y_{n+1}(t) - y_n(t)| \le \frac{M K^n (t - t_0)^{n+1}}{(n + 1)!}$. So the series $y_0 + (y_1 - y_0) + (y_2 - y_1) + \dots$ converges for each $t$ by the comparison test, and the approximations tend to some function $y(t)$, differing from it by less than the tail of the numerical series $\sum\frac{MK^nh^{n+1}}{(n+1)!}$, the same for all $t$. It remains to pass to the limit in the equation $y_{n+1} = y_0 + \int f(s, y_n)$: the integral differs from $\int f(s, y)$ by at most $Kh$ times this tail. For this we also need to know that the limit function is continuous. It is, because the convergence is equally fast for all $t$; this is called uniform convergence, and the corresponding theorem will be proved in chapter 53. Together with it, the proof is complete.
Picard’s iterations can be seen in action. For $y' = y$, $y(0) = 1$ we get $y_1 = 1 + t$, $y_2 = 1 + t + \frac{t^2}{2}$, $y_3 = 1 + t + \frac{t^2}{2} + \frac{t^3}{6}$: the Taylor polynomials of the exponential from the last chapter, which converge to $e^t$.
Can two different integral curves of the equation $y' = t - y$ cross?
The function $t - y$ satisfies a Lipschitz condition with $K = 1$. If two curves had a common point $(t_0, y_0)$, both would be solutions of the same initial value problem, and by Picard’s theorem the solution is unique. So integral curves neither cross nor touch, which is exactly what you can see on bench one.
Picard’s theorem is the mathematical form of determinism. In 1814 Pierre-Simon Laplace wrote of an imaginary intellect that would know the positions and velocities of all the particles in the world and would therefore see the future as clearly as the past. By Picard’s theorem such an intellect really could compute the future for any system with smooth laws. But only if it knew the initial conditions absolutely exactly. The last bench shows what that proviso is worth.
Bench seven: the butterfly
In 1961 the meteorologist Edward Lorenz of the Massachusetts Institute of Technology was running a simplified model of the atmosphere on a computer. To repeat a run from the middle, he typed in an intermediate value from a printout where it had been rounded: $0.506$ instead of $0.506127$. The difference was about one ten-thousandth, less than a typical measurement error. At first the repeat run tracked the first one step for step, and then the model weather went off in a completely different direction. In 1963 Lorenz published an even simpler system of three equations that shows the same thing:
$$x' = \sigma(y - x), \qquad y' = x(\rho - z) - y, \qquad z' = xy - \beta z,$$with the numbers $\sigma = 10$, $\rho = 28$, $\beta = \frac83$. The right-hand sides are polynomials, the Lipschitz condition holds, and Picard’s theorem guarantees that every initial state determines the future uniquely. See what that means in practice.
Each tenfold reduction of the initial error puts off the parting by only about $2.5$ units of time. The error grows exponentially, like a solution of the equation $y' = ky$, only here $k$ is set not by an outside law but by the system itself. Each new correct digit in the initial data adds the same amount of time to the forecast horizon: the precision of measurements has to be increased exponentially, while the forecast gets longer only linearly.
The first part of the answer. The equations that govern the atmosphere are deterministic: by Picard’s theorem the present determines the future uniquely. But uniquely doesn’t mean predictably. In systems like the Lorenz system a small inaccuracy in the initial data grows exponentially, and each new correct digit in the measurements extends the forecast only by a fixed amount. Why this happens, what the Lorenz strange attractor looks like, and where the limit of weather predictability of about two weeks comes from is in chapter 59.
Where next
On the wolves-and-hares bench the solutions ran along the level curves of the function $V(x, y)$, like hikers who go round a mountain along a contour line, neither climbing nor descending. Functions of several variables have turned up at every step of this chapter: $f(t, y)$ in the direction field, $V(x, y)$ in Volterra’s work, SIR with its three shares. But we only looked at them out of the corner of our eye. Picture a mountain landscape: height is a function of two coordinates. How do you measure the steepness of a slope if it depends on the direction? Which way should you go to climb fastest? Where does such a function have peaks, saddles and hollows? The answer comes in the chapter on functions of several variables.