Mathematics RU

Part IX · Horizons Chapter 59 of 60

Chaos and fractals

A one-line formula, Lorenz’s three equations, a game of dice, and behaviour that no computer can predict. A lab notebook of nine experiments: why the weather can’t be forecast a month ahead, and where a coastline gets a fractional dimension.

Years 1–2 60 minutes

Builds on: 58 · Non-Euclidean geometries

You will learn

  • explore iterations of a map with a cobweb diagram and check the stability of a fixed point by the derivative
  • explain why chaos eats the digits of the initial data, and estimate the forecast horizon from the Lyapunov exponent
  • find the fractal dimension by box counting and from self-similarity

6Why can’t the weather be forecast a month ahead?

The last chapter ended with Jacques Hadamard’s billiards. On a surface of negative curvature two paths sent out in almost the same direction diverge exponentially, and a tiny error in the shot makes the ball’s path unpredictable. The rules of motion are simple and exact, yet the distant future is hidden. This chapter is about the fact that not only exotic surfaces behave this way, but also a calculator formula, rabbits on an island, a double pendulum and the weather. There’s no randomness in the rules: they simply blow up any inaccuracy, however small, at tremendous speed.

The course has run into this before. In chapter 31 Edward Lorenz rounded the number $0.506127$ to $0.506$ in a weather computation and got different weather. There we found the first half of the answer to the sixth big question: the equations of the atmosphere determine the future uniquely, but uniquely doesn’t mean predictably. Now we’ll assemble the whole answer: why the error grows exponentially, what shape the motion lives on, and where the forecast limit of about two weeks comes from.

The chapter is set out like a lab notebook. Nine experiments, each with a question, an apparatus (a widget) and an explanation. Let’s start with the simplest apparatus imaginable: a single calculator formula.

Experiment 1. Rabbits year by year

Question: what happens to the size of a population if breeding is limited by food?

In chapter 31 growth with a ceiling was described by Verhulst’s logistic equation, and its solutions were calm S-shaped curves. But many animals breed not continuously but once a year. Then it’s convenient to record the population year by year: $x_n$ is the share of the largest population the island can feed, in year number $n$. The more rabbits, the more offspring, but the less grass for each. The simplest law with both effects is a product.

Next year’s population, as a share of the largest possible: a number from $0$ to $1$. The breeding rate. To keep the values from jumping out of the interval $[0, 1]$, take $0 < r \le 4$: the largest value of $rx(1 - x)$ is attained at $x = \frac12$ and equals $\frac r4$. The more rabbits now, the more offspring. Free room: the closer the population is to the ceiling, the less food for each. Example: $r = 2.5$, $x_0 = 0.2$. Then $x_1 = 2.5 \cdot 0.2 \cdot 0.8 = 0.4$, $x_2 = 2.5 \cdot 0.4 \cdot 0.6 = 0.6$, $x_3 = 2.5 \cdot 0.6 \cdot 0.4 = 0.6$, and from then on always $0.6$.

A rule that gives the next state from the current one is called a map, and applying it repeatedly is called iteration. In 1976 the biologist Robert May published a review in Nature, “Simple mathematical models with very complicated dynamics”. Its main hero was exactly this map, and May urged that students be introduced to it as early as possible, so they’d get used to the idea that a simple equation needn’t behave simply.

Let $f$ be a map of an interval to itself. The orbit of a point $x_0$ is the sequence $x_0,\ x_1 = f(x_0),\ x_2 = f(x_1),\ \dots$

A point that the map leaves in place, $f(x^*) = x^*$, is called a fixed point, like the point on the map from chapter 55 that lies exactly over the place it depicts. An orbit started at a fixed point stays put.

An orbit is conveniently drawn as a cobweb. Put the point $x_0$ on the axis and go up to the graph: the height is $x_1$. To make $x_1$ the new argument, go horizontally to the diagonal $y = x$ and then vertically to the graph again. The fixed points are where the graph crosses the diagonal. The logistic map has two: $x = 0$ and $x^* = 1 - \frac1r$ (solve $x = rx(1 - x)$ and cancel $x$).

Drag $x_0$ along the axis and change $r$. At $r = 2.8$ the cobweb spirals in to the fixed point, at $3.2$ it goes round a rectangle, at $3.5$ round two nested ones, at $4$ it fills almost the whole square.

Observation: at $r = 2.8$ the cobweb spirals inwards, to the point $x^* = \frac{9}{14} \approx 0.643$. The slope of the graph at that point is $-0.8$. At $r = 3.2$ the fixed point $x^* = 0.6875$ hasn’t gone anywhere, but the orbit no longer approaches it: the slope has become $-1.2$. It looks as if the slope decides everything, and so it does.

Let a function $f$ be continuously differentiable near a fixed point $x^*$. If $|f'(x^*)| < 1$, then $x^*$ has a neighbourhood from which every orbit converges to $x^*$, with $|x_n - x^*| \le q^n |x_0 - x^*|$ for some $q < 1$. If $|f'(x^*)| > 1$, then an orbit started arbitrarily close to $x^*$, but not at it, is bound to leave some neighbourhood of $x^*$.

The contraction factor: a number less than one that bounds $|f'|$ near the fixed point. Each step reduces the distance by a factor of at least $\frac1q$. The initial distance to the fixed point. At $r = 2.8$ the slope at the fixed point is $-0.8$. If $|f'| \le 0.9$ near it, then in 20 steps the distance shrinks by a factor of at least $0.9^{-20} \approx 8.2$, in 50 by a factor of $194$.

The trick is to link one step of the orbit with the slope of the graph through the mean value theorem, and then just repeat.

The drawing shows the graph of $f$, the diagonal, and the fixed point $x^*$ where they cross. The condition for being fixed: $f(x^*) = x^*$. The slope of the tangent at $x^*$ is $k = f'(x^*)$, and $|k| < 1$. The derivative is continuous, so near $x^*$ it differs little from $k$: there’s an interval $[x^* - \rho, x^* + \rho]$ (shaded) and a number $q < 1$ such that $|f'(x)| \le q$ on the whole interval. Any $q$ between $|k|$ and $1$ will do, if $\rho$ is taken small enough. Take a start $x_0$ in this interval. Write $d_0 = |x_0 - x^*|$. Take a step: $x_1 = f(x_0)$. By the mean value theorem $f(x_0) - f(x^*) = f'(\xi)(x_0 - x^*)$ for some point $\xi$ between $x_0$ and $x^*$: the chord is parallel to the tangent at $\xi$. The point $\xi$ lies in our interval, so $|f'(\xi)| \le q$, and $f(x^*) = x^*$. So $d_1 = |x_1 - x^*| \le q\,d_0$. Since $d_1 \le q d_0 < d_0 \le \rho$, the point $x_1$ is in the interval again, and the same argument gives $d_2 \le q d_1 \le q^2 d_0$. By induction $d_n \le q^n d_0$ for all $n$. In the drawing the cobweb shrinks towards $x^*$. Since $q < 1$, $q^n \to 0$, and $x_n \to x^*$. If instead $|f'(x^*)| > 1$, then near $x^*$ we have $|f'| \ge Q > 1$, and the same mean value theorem gives $|x_{n+1} - x^*| \ge Q|x_n - x^*|$ while the point stays in the neighbourhood. The distance grows like $Q^n$ and can’t stay below the radius of the neighbourhood for ever: the orbit leaves it.

This is a local relative of Banach’s theorem: on the shaded interval $f$ is a contraction with factor $q$, and the fixed point attracts everything around it. Only here the contraction isn’t given in advance but comes from the derivative.

For the logistic map $f'(x) = r(1 - 2x)$, and at the point $x^* = 1 - \frac1r$ we get $f'(x^*) = r\left(1 - 2 + \frac2r\right) = 2 - r$. The condition $|2 - r| < 1$ means $1 < r < 3$. The sign of the slope explains the shape of the cobweb too: for $1 < r < 2$ the slope is positive, and the values approach $x^*$ in a staircase from one side; for $2 < r < 3$ it’s negative, and the population oscillates around $x^*$, one year more, one year less, dying down.

Find the non-zero fixed point of the map $x \mapsto 2.8\,x(1 - x)$. Give the answer as a fraction.

$x^* = 1 - \frac{1}{2.8} = 1 - \frac{5}{14} = \frac{9}{14} \approx 0.643$. Check: $2.8 \cdot \frac{9}{14} \cdot \frac{5}{14} = \frac{126}{196} = \frac{9}{14}$.

Experiment 2. Doubling

Question: what takes the place of a fixed point when it stops attracting?

At $r = 3.2$ the orbit settles on two values, $0.5130$ and $0.7995$, and jumps between them: a year of plenty, a year of hunger. Such a pair is called a cycle of period 2.

A cycle of period $p$ is an orbit that returns to its starting point after exactly $p$ steps and not sooner: $f^p(x_0) = x_0$, where $f^p$ is the map applied $p$ times. Each point of the cycle is a fixed point of the map $f^p$.

This gives a way to find a cycle: solve $f(f(x)) = x$. This is a quartic equation, but we know two of its roots in advance, the fixed points $0$ and $1 - \frac1r$, since they certainly come back after two steps. After dividing by $x\left(x - 1 + \frac1r\right)$ we’re left with the quadratic equation $r^2x^2 - r(r + 1)x + (r + 1) = 0$, with roots

$$x_{1,2} = \frac{(r + 1) \pm \sqrt{(r - 3)(r + 1)}}{2r}.$$

The number under the root is positive only for $r > 3$: the cycle is born exactly when the fixed point loses stability. At $r = 3.2$ the formula gives $\frac{4.2 \pm \sqrt{0.84}}{6.4}$, that is, $0.5130$ and $0.7995$, as in the widget.

The period-2 cycle of the logistic map is stable (attracts neighbouring orbits) for $3 < r < 1 + \sqrt6 \approx 3.449$.

The points $p$ and $s$ of the cycle are fixed points of the map $g = f \circ f$, so it’s enough to apply the stability test to $g$. By the chain rule $g'(p) = f'(f(p))\,f'(p) = f'(s)\,f'(p)$, and at the point $s$ we get the same product. Let’s compute it: $f'(p)f'(s) = r^2(1 - 2p)(1 - 2s) = r^2\bigl(1 - 2(p + s) + 4ps\bigr)$. By Vieta’s formulas for the quadratic above, $p + s = \frac{r + 1}{r}$ and $ps = \frac{r + 1}{r^2}$. Substituting: $r^2 - 2r(r + 1) + 4(r + 1) = -r^2 + 2r + 4$. At $r = 3$ this is $1$, and as $r$ grows the function $-r^2 + 2r + 4$ decreases. For $r > 3$ the condition $\left|-r^2 + 2r + 4\right| < 1$ reduces to $-r^2 + 2r + 4 > -1$, that is, to $r^2 - 2r - 5 < 0$, and the upper limit is the root $r = 1 + \sqrt6$, where the product equals $-1$.

Beyond $1 + \sqrt6$ the story repeats: the period-2 cycle loses stability, and a cycle of period 4 is born. Then 8, 16, 32… The doublings come faster and faster, and near $r \approx 3.5699$ the period becomes infinite. The whole story fits into one picture: for each $r$, plot vertically the values the orbit visits after it has forgotten its start.

A qualitative change of behaviour as a parameter changes smoothly (for example, the birth of a cycle of twice the period) is called a bifurcation (from the Latin bifurcus, forked).

Tap the diagram or drag round a piece with the mouse: a magnifier zooms in. The ticks $r_1, r_2, \dots$ are the doubling points. Look into the period-3 window near $r \approx 3.83$ and zoom in on its middle branch: there’s a reduced, upside-down copy of the whole diagram.

Let’s write down the doubling points. We derived the first two ourselves; the rest were found numerically.

$n$$r_n$: a cycle of period $2^n$ is born$\dfrac{r_{n} - r_{n-1}}{r_{n+1} - r_n}$
1$3$
2$1 + \sqrt6 \approx 3.449490$$4.751$
3$3.544090$$4.656$
4$3.564407$$4.668$
5$3.568759$$4.669$
$\infty$$3.569946$

The intervals between doublings shrink each time by roughly the same factor, and this factor tends to a limit.

The factor by which each next interval between doublings is shorter than the one before. The length of the interval on which the cycle of period $2^{n-1}$ lives. The length of the next interval, for the cycle of twice the period. From the table: $\frac{3.449490 - 3}{3.544090 - 3.449490} = \frac{0.449490}{0.094600} \approx 4.751$, then $\frac{0.094600}{0.020317} \approx 4.656$, then $4.668$ and $4.669$. The intervals decrease almost like a geometric progression with ratio $\frac1\delta$, so their sum is finite and the doublings accumulate at the point $r_\infty \approx 3.5699$.

Experiment 3. The window with the number three

Question: what does the word “chaos” mean, and how did it get into mathematics?

Beyond the point $r_\infty$ the diagram darkens: orbits no longer settle on cycles. But chaos is interspersed with windows of order. The widest opens at $r = 1 + \sqrt8 \approx 3.8284$: the orbit suddenly settles on a cycle of period 3 (at $r = 3.83$ it’s $0.1561$, $0.5047$, $0.9574$), and inside the window there’s a doubling cascade of its own: 3, 6, 12, 24…

In 1975 Tien-Yien Li and James Yorke published a paper with a title that became a slogan: “Period three implies chaos”. That’s how the word “chaos” entered mathematics. Their main observation is strikingly general: it holds for any continuous function on the line.

If a continuous function $f\colon \mathbb R \to \mathbb R$ has a cycle of period 3, then it has cycles of all periods $n = 1, 2, 3, \dots$

The idea: make two intervals from the period-3 cycle that the function stretches over each other, and build a route of any length from them. We’ll need the intermediate value theorem: a continuous function on an interval takes all values between its values at the ends.

Let $a \to b \to c \to a$ be a cycle of period 3, with $a$ the smallest of the three points. In the drawing $a < b < c$, so $f(a) = b$, $f(b) = c$, $f(c) = a$. If the order is different, $a < c < b$, the argument is the same, only the intervals taken are $[c, b]$ and $[a, c]$. Take two intervals: $I_0 = [a, b]$ and $I_1 = [b, c]$. They have one point in common, $b$. At the ends of $I_0$ the function takes the values $f(a) = b$ and $f(b) = c$, so by the intermediate value theorem $f(I_0) \supset [b, c] = I_1$: the graph over $I_0$ passes through all heights from $b$ to $c$. At the ends of $I_1$ the values are $f(b) = c$ and $f(c) = a$, so $f(I_1) \supset [a, c] = I_0 \cup I_1$. The interval $I_1$ covers both itself and its neighbour. Lemma: if $f(I) \supset J$ for intervals $I$ and $J$, then inside $I$ there’s an interval $K$ with $f(K) = J$. Indeed, take a point of $I$ where $f$ equals one end of $J$, and the nearest point to it where $f$ equals the other end; between them $f$ doesn’t leave $J$, otherwise it would have crossed one of the ends sooner. Now for $n \ge 4$ follow the route: $n - 1$ times from $I_1$ to $I_1$, then to $I_0$, and from $I_0$ back to $I_1$. Every move is allowed: $I_1$ covers $I_1$ and $I_0$, and $I_0$ covers $I_1$. Go along the route from the end and use the lemma to build nested intervals. We get an interval $K_0 \subset I_1$ whose points follow the route: their images from the first to the $(n-2)$th lie in $I_1$, the $(n-1)$th in $I_0$, and $f^n(K_0) = I_1 \supset K_0$. The function $f^n(x) - x$ has values of opposite sign, or zero, at the ends of $K_0$, so by the intermediate value theorem $K_0$ contains a point $x$ with $f^n(x) = x$. It remains to check that the period of $x$ is exactly $n$, not a divisor $m < n$. If $f^m(x) = x$, then $f^{n-1}(x) = f^{m-1}(x)$, which is a point of $I_1$, since $m - 1 \le n - 2$ (for $m = 1$ it’s $x \in I_1$ itself). But by construction $f^{n-1}(x) \in I_0$, so $f^{n-1}(x) = b$, the only common point. Then $x = f^n(x) = f(b) = c$ and $f(x) = a$, and $a$ doesn’t lie in $I_1$, although by the route $f(x) \in I_1$. A contradiction. Periods $n \ge 4$ have been obtained. Period 3 is given. Period 1: $f(I_1) \supset I_1$, so $f(x) - x$ changes sign at the ends of $I_1$, and it contains a fixed point. Period 2 comes from the route $I_0 \to I_1 \to I_0$ by the same device. So $f$ has cycles of all periods.

Li and Yorke didn’t know that 11 years before them the Kyiv mathematician Oleksandr Sharkovsky had proved a much stronger theorem. He arranged the natural numbers in a special order that begins $3, 5, 7, 9, \dots$ and ends $\dots, 8, 4, 2, 1$, and showed that if a continuous function on the line has a cycle of some period, it has cycles of all the periods to the right of it in this order. Three comes first, so everything follows from it. Sharkovsky’s paper came out in 1964 in the Ukrainian Mathematical Journal and long remained unknown in the West.

Infinitely many cycles isn’t yet chaos in the everyday sense: almost all of them are unstable, and they don’t attract an orbit started at random. Real disorder begins where there are no attracting cycles at all. The purest such case is $r = 4$.

Experiment 4. The calculator loses digits

Question: can’t we just compute the orbit more precisely and predict it a hundred steps ahead?

Take $r = 4$ and start from Lorenz’s number $0.506127$. Compute the orbit in four ways: almost exactly (with four hundred binary digits), with ordinary computer numbers of 16 significant digits, with single-precision numbers of seven digits, and, like Lorenz, from a start rounded to $0.506$.

The graph shows the decimal logarithm of each method’s error. All three lines rise with the same slope, about $0.3$ per step, until they hit zero. Move the step slider and watch the correct digits melt away.

Observation: 16 digits last 46 steps, seven digits 16, three just 7. Every three or four steps eat one decimal digit, and twice the precision buys one extra step. This isn’t a failing of the computer; it’s a property of the map itself, and at $r = 4$ it can be seen completely. Back in 1947 Stanisław Ulam and John von Neumann noticed that $4x(1 - x)$ is angle doubling seen from the side.

If $x = \frac{1 - \cos\varphi}{2}$, then $4x(1 - x) = \frac{1 - \cos 2\varphi}{2}$. Hence the orbit of the map $x \mapsto 4x(1 - x)$ starting at $x_0 = \frac{1 - \cos\varphi_0}{2}$ is $x_n = \frac{1 - \cos(2^n\varphi_0)}{2}$.

The idea: a number $x$ in the interval $[0, 1]$ is the shadow of a point of a circle on its diameter.

Take a circle of radius 1 and take its horizontal diameter as the interval $[0, 1]$. The point $P$ stands at the angle $\varphi$, measured from the left end. Its shadow on the diameter is at distance $1 - \cos\varphi$ from the left end, and the whole diameter is 2, so on our scale it’s $x = \frac{1 - \cos\varphi}{2}$. Double the angle: the point $P'$ stands at the angle $2\varphi$. For the same reason its shadow is $x' = \frac{1 - \cos 2\varphi}{2}$. Link $x'$ with $x$. The definition of $x$ gives $\cos\varphi = 1 - 2x$, and by the double-angle formula $\cos 2\varphi = 2\cos^2\varphi - 1$. So $x' = \frac{1 - (2\cos^2\varphi - 1)}{2} = 1 - \cos^2\varphi = 1 - (1 - 2x)^2 = 4x - 4x^2 = 4x(1 - x)$. One step of the map $x \mapsto 4x(1 - x)$ doubles the angle. Doing this $n$ times, we get the angle $2^n\varphi_0$ and $x_n = \frac{1 - \cos(2^n\varphi_0)}{2}$.

Now you can see where the digits go. Write the angle as a fraction of a full turn in binary: $\varphi_0 = 2\pi \cdot 0.b_1b_2b_3\dots_2$. Doubling shifts the point one place to the right, and whole turns don’t matter, so after $n$ steps what’s left is $0.b_{n+1}b_{n+2}\dots_2$. Each step throws away the leading binary digit and brings forward the next, previously invisible one. A double-precision number stores 53 binary digits, and after about that many steps digits that weren’t in it come to the front: they were invented by rounding error. Hence exactly the slope on the graph: $\log_{10} 2 \approx 0.301$ decimal digits per step.

The rate at which errors multiply is measured by a single number. A small error $\varepsilon$ becomes roughly $|f'(x)|\,\varepsilon$ after a step, and $|f'(x_0)|\,|f'(x_1)| \cdots |f'(x_{n-1})|\,\varepsilon$ after $n$ steps. The logarithm of this product is a sum, and it’s natural to take its average.

The average rate of exponential growth of small errors: in $n$ steps an error grows by a factor of about $e^{\lambda n}$. Averaging over the steps of the orbit. By how much, in logarithms, the error is stretched at step $k$. For $r = 4$ angle doubling gives $\lambda = \ln 2 \approx 0.693$: the error doubles each step on average. Numerically along an orbit one gets $0.6931$. At $r = 2.8$ the exponent is negative, $\ln 0.8 \approx -0.223$: errors die out. At $r = 3.7$ it’s about $0.355$, and in the period-3 window at $r = 3.83$ it’s negative again, about $-0.37$.

For multidimensional systems, like Lorenz’s, the definition is the same, only instead of the derivative one takes the stretching of a small segment along the fastest direction. The exponent is named after Aleksandr Lyapunov, who built the theory of stability of motion in his 1892 thesis.

Bounded motion given by a deterministic rule, in which small differences in initial conditions grow exponentially on average (the Lyapunov exponent is positive), is called deterministic chaos. The growth of small differences itself is called sensitive dependence on initial conditions.

A computer computes the orbit of $x \mapsto 4x(1 - x)$ with numbers of 16 significant digits. How many digits of the hundredth value can be trusted?

The rounding error on input is of order $10^{-17}$, and each step doubles it on average. By step one hundred it would have grown by a factor of $2^{100} \approx 10^{30}$, and in fact by step fifty it has reached the size of the numbers themselves. The hundredth value is a plausible number in $[0, 1]$, but not a value of our orbit.

The trainer below is on fixed points and their stability, and at the third level on angle doubling. There it turns out that the period of doubling for a fraction is the order of two modulo the denominator, from the chapter on remainders.

Experiment 5. The butterfly

Question: what does chaos look like in a real physical system?

Back to the system where it all began. In 1963 Lorenz published the paper “Deterministic Nonperiodic Flow”: three equations for convection in a layer of air heated from below,

$$x' = \sigma(y - x), \qquad y' = x(\rho - z) - y, \qquad z' = xy - \beta z,$$

with $\sigma = 10$, $\rho = 28$, $\beta = \frac83$. Here $x$ is the speed of rotation of a convection roll, and $y$ and $z$ describe the distribution of temperature. In chapter 31 we looked at the graph of one coordinate. Now let’s look at the whole path in space.

Two paths start a distance $\delta$ apart. At first they follow the same track, then they part onto different wings. Switch to “A cloud of 500”: 500 starts inside a little ball of radius $\delta$ spread over the whole attractor. On a phone press “Rotate” to turn the scene.

The path neither goes off to infinity nor settles down. It winds loops around one point, then for no visible reason jumps to the other wing and winds there. How many loops it makes on each wing can’t be guessed in advance, but the figure itself, like a butterfly or a mask, is the same for any start.

An attractor is a set that attracts, over time, all the paths from some neighbourhood of it. A stable fixed point and a stable cycle are simple attractors. An attractor on which the motion is chaotic and whose shape is a fractal (more on those in experiment 7) is called strange; the name comes from a 1971 paper by David Ruelle and Floris Takens.

The lower graph shows the logarithm of the distance between the paths. While it’s small, the logarithm grows along a straight line with a slope of about $0.39$ per unit of time. That’s the Lyapunov exponent of the Lorenz system, $\lambda \approx 0.906$, converted to decimal logarithms: $\frac{0.906}{\ln 10} \approx 0.39$. When the distance grows to the size of the attractor, it has nowhere left to grow, and the paths simply end up in different, unrelated places on the butterfly.

As a control, let’s run an experiment with no atmosphere and no simplified models: a double pendulum, that is, a pendulum hung from another pendulum. Its equations of motion are derived directly from Newton’s laws, and they’re just as deterministic as Lorenz’s.

The lower rods of the twins are displaced by angles differing by $\delta$ radians. At $\delta = 10^{-9}$ (the thickness of a hair seen from a hundred kilometres away) the pendulums part after about 12 seconds. Compare with $\delta = 10^{-3}$: the difference is a million times larger, yet they part only about eight seconds sooner.

Experiment 6. A month’s forecast

Question: how far ahead can a system with a positive Lyapunov exponent be predicted, if everything is measured very precisely?

Suppose the initial state is known with an error $\delta$, and the forecast stops being useful when the error grows to $\Delta$, for example to the size of the attractor itself, when the forecast is no better than a random guess. The error grows roughly like $\delta e^{\lambda t}$. Solve the equation $\delta e^{\lambda T} = \Delta$.

The time for which the forecast remains useful. The system’s characteristic time: in $\frac1\lambda$ the error grows by a factor of $e$, in $\frac{\ln 2}{\lambda}$ it doubles. The error we’re prepared to tolerate: once it’s reached, the forecast is useless. The initial error: the precision of the measurements. For the Lorenz system $\lambda \approx 0.906$. With $\delta = 10^{-4}$ and $\Delta = 10$ we get $T \approx \frac{\ln 10^5}{0.906} \approx 12.7$. If we measure ten times more precisely, $\delta = 10^{-5}$, then $T \approx 15.2$: only $\frac{\ln 10}{0.906} \approx 2.5$ units of time have been added, as we saw in chapter 31. For $x \mapsto 4x(1 - x)$ with $\lambda = \ln 2$, $\delta = 10^{-16}$, $\Delta = 10^{-2}$ we get $T \approx 46$ steps, as the precision race showed. The formula is an estimate: it’s about average growth, and for a particular start the divergence may come a few units sooner or later.

The precision of the initial data enters the length of a forecast under a logarithm. To double the length of a forecast you don’t double the precision of the measurements; you square it.

The initial state of the Lorenz system is now measured 100 times more precisely. How does the length of a reliable forecast change?

By the horizon formula the gain is $\frac{\ln 100}{\lambda} = \frac{4.6}{0.906} \approx 5.1$ units of time, regardless of what the forecast was before.

For the map $x \mapsto 4x(1 - x)$ the Lyapunov exponent is $\ln 2$. By how many steps does a reliable forecast lengthen if the initial error is reduced 1024-fold?

The gain is $\frac{\ln 1024}{\ln 2} = \log_2 1024 = 10$ steps: each doubling of precision gives exactly one step.

For the atmosphere things are even worse than this formula says. The atmosphere has motions of all sizes, from a gust of wind between houses to a cyclone half a continent across. Errors in small eddies grow very fast, within hours, and are passed on to larger and larger motions. In 1969 Lorenz showed, on a model with such a ladder of scales, that the limit of predictability here is finite even if the observations of small details are refined as much as you like: the lifetime of a forecast is made up of the times it takes the error to creep from rung to rung, and this sum converges. His estimate was about two weeks. Modern calculations with the models used for real forecasts arrive at the same figure; that’s what Fuqing Zhang and co-authors found, for example, in 2019.

Practical forecasts are also limited by incomplete observations and imperfect models, so today a detailed forecast is useful for about a week to ten days. Meteorologists have long computed not one weather but an ensemble: the European Centre for Medium-Range Weather Forecasts runs its model 51 times with slightly different initial data, like the cloud of 500 points in the widget, only in a space with not three variables but astronomically many. While the cloud is compact, the forecast is confident. When it spreads over the attractor, the forecast turns into probabilities: “a 60 % chance of rain”.

But climate can be predicted, and there’s no contradiction in that. Weather is a point on the attractor, and a month from now we don’t know where it will be. Climate is the shape of the attractor and how often the path visits its different parts. The shape doesn’t depend on small errors: all 500 points of the cloud draw the same butterfly.

The atmosphere obeys deterministic equations but moves chaotically: it has a positive Lyapunov exponent, and small measurement errors grow exponentially. So each new correct digit in the initial data extends the forecast only by a fixed amount, and because of motions on all scales, from an eddy to a cyclone, the limit can’t be pushed back even by infinite precision: it’s about two weeks. A month ahead only the statistics can be predicted, averages and probabilities, that is, the climate, not the weather on a particular day.

Experiment 7. How long is a coastline

Question: what shape do strange attractors have, and coastlines?

The English physicist and meteorologist Lewis Fry Richardson, the very one who first tried to compute the weather from equations, in his later years studied what makes wars start, and for this he collected the lengths of national borders. The figures didn’t agree: reference books of the two countries gave the common border of Spain and Portugal as 987 and 1214 km. Richardson understood why. If you measure a winding line with dividers opened to $\varepsilon$, the smaller $\varepsilon$, the more bays and headlands they notice, and the longer the line comes out. For many coasts the length grew by the law $L(\varepsilon) \approx C\varepsilon^{1 - D}$ with a constant $D$ greater than one. His observations were published posthumously, in 1961.

In 1967 Benoit Mandelbrot published a paper in Science, “How Long Is the Coast of Britain?”, and proposed regarding the number $D$ as the dimension of the coastline. Richardson’s data gave $D \approx 1.25$ for the west coast of Britain and about $1.02$ for the smooth coast of South Africa. A fractional dimension isn’t a play on words. Here’s how it can be measured.

Cover the plane with a grid of squares of side $\varepsilon$ and count $N(\varepsilon)$, the number of boxes the figure touches. For a segment of length 1 there are about $\frac1\varepsilon$ boxes, for a square of side 1 about $\frac{1}{\varepsilon^2}$. The exponent is the dimension.

The dimension: the exponent $D$ for which $N(\varepsilon)$ behaves like $\varepsilon^{-D}$. The number of boxes of the grid of side $\varepsilon$ that the figure touches. How many times smaller than one the side of a box is. On the graph of $\ln N$ against $\ln\frac1\varepsilon$ the dimension is the slope of the line. Example: the Koch curve (the Swedish mathematician Helge von Koch, 1904). It’s built from a segment by replacing the middle third with two sides of an equilateral triangle, and repeating this with every segment without end. At $\varepsilon = 3^{-m}$ the curve is covered by about $4^m$ boxes, so $D = \frac{\ln 4^m}{\ln 3^m} = \frac{\ln 4}{\ln 3} \approx 1.262$. For a segment $D = 1$, for a square $D = 2$. The limit doesn’t exist for every figure, but it does for all the figures in this chapter.

The number $D$ in the formula above is called the box-counting dimension (also the Minkowski dimension). The word “fractal” was coined by Mandelbrot in 1975 from the Latin fractus, broken. There’s no universally accepted rigorous definition; a fractal is usually a set in which new detail appears at every scale and whose dimension exceeds the “naive”, topological one: the Koch curve has dimension greater than one, even though it’s a curve.

Move the slider: the grid gets finer, and the number of boxes grows. On the right each grid is a point on the graph of $\ln N$ against $\ln\frac1\varepsilon$, and the points lie on a line. For the circle the slope is close to 1, for the square to 2, for the island and the Koch curve about 1.25, for the Sierpiński triangle about 1.57.

Richardson’s formula makes the confusion over borders clear. The length measured with dividers opened to $\varepsilon$ is proportional to $\varepsilon \cdot N(\varepsilon) \approx \varepsilon^{1 - D}$. With $D = 1.25$, halving the opening lengthens the coast by a factor of $2^{0.25} \approx 1.19$, and so at every scale, until we get down to individual stones. Asking how long a coastline is without specifying the opening of the dividers is meaningless.

The purest fractal is the Sierpiński triangle. The Polish mathematician Wacław Sierpiński described it in 1915: from a triangle remove the middle one of four triangles half the size, from each remaining one again the middle one, and so on without end. This pattern has already appeared in Pascal’s triangle, if the odd numbers are coloured. Its dimension can be found exactly.

The box-counting dimension of the Sierpiński triangle is $\frac{\ln 3}{\ln 2} \approx 1.585$.

The trick is to take a grid of side $2^{-m}$ and bound $N$ above and below by numbers that differ only by a constant factor. The Sierpiński triangle can be built from any triangle; take a right triangle with legs 1, so that all the corners fall on grid points. Stretching the plane changes the number of boxes by at most a constant factor, and such a factor, as we’ll see, doesn’t affect the dimension.

After $m$ steps of the construction $3^m$ little pieces remain, right triangles with legs $2^{-m}$, and the whole figure lies in them. In the drawing $m = 3$ and there are 27 pieces. Cover the square with a grid of side $\varepsilon = 2^{-m}$. Take the boxes half-open, without their right and top sides: then each point of the plane lies in exactly one box. Each little piece is inscribed in a grid square: its right angle at a grid point, its legs along grid lines. The upper bound. A closed grid square touches at most four half-open boxes: its own and the neighbours to the right, above and diagonally. Each point of the figure lies in some little piece, so $N(2^{-m}) \le 4 \cdot 3^m$. The lower bound. The vertices of the right angles of all the little pieces belong to the figure: the corner of a piece stays a corner of one of its parts at later steps and is never removed. These are $3^m$ different grid points, and each lies in its own box. So $N(2^{-m}) \ge 3^m$. Take logarithms of $3^m \le N \le 4 \cdot 3^m$ and divide by $\ln\frac1\varepsilon = m\ln 2$: $\frac{\ln 3}{\ln 2} \le \frac{\ln N}{m \ln 2} \le \frac{\ln 3}{\ln 2} + \frac{\ln 4}{m\ln 2} = \frac{\ln 3}{\ln 2} + \frac{2}{m}$. As $m \to \infty$ both bounds tend to $\frac{\ln 3}{\ln 2}$. That leaves the intermediate $\varepsilon$ between $2^{-m-1}$ and $2^{-m}$. A box of side $\varepsilon$ touches at most four boxes of the grid of side $2^{-m}$, and a box of side $2^{-m-1}$ at most four boxes of side $\varepsilon$. So $\frac14 N(2^{-m}) \le N(\varepsilon) \le 4N(2^{-m-1})$, and after dividing by $\ln\frac1\varepsilon$ the factor 4 disappears in the limit. The dimension is $\frac{\ln 3}{\ln 2}$.

The same argument works almost word for word for any figure made of $N$ copies of itself scaled down $k$ times, if the copies overlap only at their edges: the dimension is $\frac{\ln N}{\ln k}$. The Koch curve is 4 copies at scale $\frac13$, hence $\frac{\ln 4}{\ln 3}$. The Sierpiński carpet is 8 copies at scale $\frac13$, hence $\frac{\ln 8}{\ln 3} \approx 1.893$. And the Cantor set of chapter 52 is 2 copies at scale $\frac13$, with dimension $\frac{\ln 2}{\ln 3} \approx 0.631$: more than a point, less than a segment. It has no length, a continuum of points, and a fractional dimension.

A set is called self-similar if it’s made up of several scaled-down copies of itself.

Copies of different sizes: Moran’s equation

If a figure consists of copies scaled down by different factors $r_1, \dots, r_N$, its dimension $D$ is found from Moran’s equation (Patrick Moran, 1946): $r_1^D + r_2^D + \dots + r_N^D = 1$. When all $r_i = \frac1k$, it turns into $N k^{-D} = 1$, that is, $D = \frac{\ln N}{\ln k}$. An example: one copy half the size and one a quarter the size. Write $u = 2^{-D}$; then $4^{-D} = u^2$ and $u + u^2 = 1$, so $u = \frac{\sqrt5 - 1}{2}$ and $D = \log_2\frac{1 + \sqrt5}{2} \approx 0.694$, the binary logarithm of the golden ratio.

A coast of dimension $D = 1.25$ was measured with dividers opened to 1 km, and then to 500 m. By what factor did the measured length change?

The length is $L(\varepsilon) \approx C\varepsilon^{1 - D}$. Replacing $\varepsilon$ by $\frac\varepsilon2$ multiplies it by $2^{D - 1} = 2^{0.25} \approx 1.19$. There are $2^{1.25} \approx 2.38$ times as many steps, but each is half as long.

The Vicsek snowflake is built like this: divide a square into 9 equal little squares and keep the central one and the four corner ones, do the same with each of them, and so on without end. Find its dimension to two decimal places.

Five copies scaled down 3 times: $D = \frac{\ln 5}{\ln 3} \approx 1.465$.

Experiment 8. A game of dice

Question: can a random process draw an exact figure?

In chapter 55 Barnsley’s fern grew as the fixed point of a set of contractions, and it was drawn by the “chaos game”. Let’s play it with the simplest set. Mark the three vertices of a triangle and a point anywhere. Roll a die: 1 or 2 means vertex $A$, 3 or 4 means $B$, 5 or 6 means $C$. Move halfway towards the vertex rolled and put a dot. Roll again. The English mathematician Michael Barnsley, who popularised this method in his book Fractals Everywhere (1988), called it the chaos game.

Roll one at a time and watch the first dots jump, then press “Play”. The vertices can be dragged, and the figure rebuilds itself. A square with the same rule draws nothing: its four scaled-down copies fit together without gaps. For the carpet, the pentagon and the hexagon the jump is longer than half the way.

The die is fair, the dots jump without any visible order, and yet the Sierpiński triangle grows on the screen, and with jewel-like precision: apart from the first few, not a single dot falls into the removed triangles. Chance here chooses only the order in which the details are drawn; the picture itself is set by the rules.

Let $S$ be the Sierpiński triangle with vertices $A$, $B$, $C$, $p_0$ any point, and $p_{n+1}$ the midpoint of the segment from $p_n$ to any of the vertices. Then the distance from $p_n$ to $S$ is at most $\frac{d_0}{2^n}$, where $d_0$ is the distance from $p_0$ to $S$.

The idea: shrinking by half towards a vertex takes the Sierpiński triangle to its own corner piece, and all distances are halved in the process. This is a version of Hutchinson’s lemma for a single point.

Let $p$ be the current point, $s$ the nearest point of the set $S$ to it (it exists, because $S$ is closed and bounded), and $d$ the distance between them. Vertex $B$ is rolled. The new point $p'$ is the midpoint of the segment $pB$. Shrink the whole plane by half towards the vertex $B$. The triangle $ABC$ goes to the corner triangle at the vertex $B$, and the Sierpiński triangle to its own corner piece: that’s how it’s built, from three half-size copies at the vertices. The point $s$ goes to the midpoint $s'$ of the segment $sB$, and $s'$ lies in $S$. In the triangle $psB$ the segment $p's'$ joins the midpoints of two sides, that is, it’s a midline, and $|p's'| = \frac12|ps| = \frac d2$. So the distance from $p'$ to $S$ is at most $\frac d2$. Whatever vertex is rolled, the distance to $S$ is at least halved by a roll. After $n$ rolls it’s at most $\frac{d_0}{2^n}$; after ten rolls, less than a thousandth of the original, that is, less than a pixel.

The theorem explains why the dots don’t fall outside the figure. Why they fill it all is explained by addresses. Each little piece at level $m$ can be named by a word of $m$ letters $A$, $B$, $C$: “in corner $B$, inside it in corner $A$…”. After a roll the dot lands in the corner piece of the vertex rolled last, inside that in the piece of the one before, and so on: the last $m$ rolls, read backwards, are the address of the level-$m$ piece where the dot lies. With each roll any given address comes up with probability $3^{-m}$, so sooner or later the dot visits every little piece of every level.

The other figures in the widget work the same way, only with more copies and a longer jump: the carpet has eight targets and a jump two thirds of the way, hence eight copies a third the size. Self-similarity turns a complicated picture into a short description: the Sierpiński triangle is three instructions “shrink by half towards a vertex”, and the die will draw the rest.

Experiment 9. One line, an infinite picture

Question: what happens if you iterate the square of a complex number?

The logistic map is a quadratic polynomial. Any such polynomial can be brought, by a change of variable, to the form $z \mapsto z^2 + c$, and this form is convenient because it has just one parameter. Let’s check it on ours: putting $z = \frac r2 - rx$, the step $x \mapsto rx(1 - x)$ becomes the step $z \mapsto z^2 + c$ with $c = \frac r2 - \frac{r^2}{4}$. At $r = 3$ we get $c = -\frac34$, at $r = 4$, $c = -2$.

Now let $z$ and $c$ be complex. The French mathematicians Gaston Julia (1918) and Pierre Fatou (1917–1920) studied such iterations long before computers and found almost everything visible in the pictures below without seeing a single one of them.

For a given $c$ the Julia set $J_c$ is the boundary of the set of those starting $z$ whose orbits under $z \mapsto z^2 + c$ stay bounded. The Mandelbrot set $M$ is the set of those $c$ for which the orbit of zero is bounded: $0,\ c,\ c^2 + c,\ \dots$

Choose $c$ by tapping or by dragging the ring on the left, and $J_c$ and the orbit of zero appear on the right. Inside the main “cardioid” the orbit of zero converges to a point, in the disc to its left to a cycle of period 2, at the “rabbit” to a cycle of period 3. Switch on the cascade on the axis: the doublings of the logistic map sit exactly where the discs touch.

The pictures are drawn by a simple test: if an orbit has once left the disc of radius 2, it will never come back.

If at some step $|z| > 2$ and $|z| \ge |c|$, then the orbit under $z \mapsto z^2 + c$ goes off to infinity. In particular, the Mandelbrot set lies in the disc $|c| \le 2$.

Let $|z| = 2 + a$ with $a > 0$ and $|z| \ge |c|$. By the triangle inequality for complex numbers $|z^2 + c| \ge |z|^2 - |c| \ge |z|^2 - |z| = |z|(|z| - 1) = (1 + a)|z|$. The new value has modulus greater than $2 + a$ and all the more at least $|c|$, so the same argument applies at the next step, and at each step the modulus is multiplied by at least $1 + a$. After $n$ steps the modulus is at least $(1 + a)^n(2 + a)$, which tends to infinity. If $|c| > 2$, the very first step from zero gives $z_1 = c$ with $|z_1| > 2$ and $|z_1| = |c|$, and the orbit of zero escapes: such a $c$ isn’t in $M$.

On the real axis $M$ occupies the interval $\left[-2, \frac14\right]$: for $c > \frac14$ the equation $z^2 + c = z$ has no real roots, and the orbit of zero grows without bound, while for $-2 \le c \le \frac14$ it’s trapped in the interval. This interval is exactly the logistic $r$ from 1 to 4, converted by the formula $c = \frac r2 - \frac{r^2}{4}$. The doubling point $r_1 = 3$ goes to $c = -0.75$, where the cardioid touches the big disc. The point $r_2 = 1 + \sqrt6$ goes to $c = -1.25$, where the next disc sticks on. The limit of the cascade $r_\infty$ goes to $c \approx -1.4012$.

The main theorem linking the two pictures also belongs to Fatou and Julia.

The Julia set $J_c$ is connected if $c \in M$, and crumbles into “dust”, each connected piece of which is a single point, if $c \notin M$.

The proof of this theorem requires the theory of analytic functions at an advanced undergraduate level, and here we’ll accept it without proof. The idea is this: the map $z \mapsto z^2 + c$ glues together pairs of points $z$ and $-z$, and the only point where this gluing degenerates is $z = 0$. If the orbit of zero escapes, the preimages of a large disc fall apart into smaller and smaller pieces, and in the limit dust remains; if it’s bounded, no falling apart happens. A full proof is in John Milnor’s book Dynamics in One Complex Variable.

Mandelbrot got the first pictures of the set $M$ in 1980 on IBM computers. In 1982 Adrien Douady and John Hubbard proved that $M$ is connected and named it after Mandelbrot. The boundary of $M$ is so complicated that in 1998 Mitsuhiro Shishikura proved that its Hausdorff dimension is 2, like a piece of the plane.

Does the point $c = -2$ lie in the Mandelbrot set?

The orbit of zero: $0 \to -2 \to 2 \to 2 \to 2 \dots$ It’s bounded, so $-2 \in M$. It’s the leftmost point of the set; in logistic coordinates it’s $r = 4$, the most chaotic case.

Where next

We can draw the Mandelbrot set to the pixel, we know it’s connected, and we know the dimension of its boundary. But a simple-looking question, whether it’s connected “in the small” (whether each of its points can be surrounded by an arbitrarily small connected neighbourhood within the set itself), still has no answer. This conjecture of local connectivity was formulated by Douady and Hubbard in the 1980s. It would imply that the interior of $M$ consists entirely of “hyperbolic” pieces, like the cardioid and the discs, where the orbit of zero has an attracting cycle. The conjecture has been proved for many parts of the set, but as a whole it’s open. Questions like this, easy to ask and so far impossible to settle, outnumber answers in mathematics. What mathematics doesn’t yet know, and how answers are sought today, is in the last chapter.

In this chapter

  1. Experiment 1. Rabbits year by year
  2. Experiment 2. Doubling
  3. Experiment 3. The window with the number three
  4. Experiment 4. The calculator loses digits
  5. Experiment 5. The butterfly
  6. Experiment 6. A month’s forecast
  7. Experiment 7. How long is a coastline
  8. Experiment 8. A game of dice
  9. Experiment 9. One line, an infinite picture
  10. Where next

Chapters