Easy Chaos with Pluto

Simply Making a Mess of Everything

September 08, 2022 · John Peach

Understanding the Predictably Unpredictable

Have you ever wondered why weather forecasts become unreliable beyond a few days? Or why your heart beats with slight irregularities even when you’re perfectly healthy? These phenomena share a fascinating common thread: chaos theory, a branch of mathematics that explains how tiny changes can lead to dramatically different outcomes.

Chaos theory isn’t just an academic curiosity - it’s everywhere around us. From the turbulent flow of water from your faucet to the fluctuations in financial markets, from the patterns of animal populations to the spread of epidemics, chaos theory helps us understand complex systems that at first glance appear random but actually follow hidden patterns.

In this interactive tutorial, we’ll use Pluto notebooks and Julia programming to explore chaos theory hands-on. You’ll discover how simple mathematical rules can generate surprisingly complex behavior, and learn why some systems are fundamentally unpredictable despite being completely deterministic.

Learning Objectives

By the end of this tutorial, you will be able to:

No advanced mathematics or programming experience is required - just bring your curiosity about how complexity emerges from simplicity.

Getting Started with Pluto

Pluto is a notebook environment for Julia. Install Julia and start it to get a command window. The first step is to add some packages we’ll use for this project. Type ”]” at the prompt (without the quotes) to enter the package manager. The packages you’ll need are:

add Pluto
add PlutoUI
add Plots
add StatsBase

It may take a few minutes as Julia checks and downloads all necessary dependencies. When everything has finished, type “Ctrl-C” to exit the package manager. In the Julia command window, type “Using Pluto” followed by “Pluto.run()”:

julia-start-pluto
Figure 1.

Starting Pluto from the Julia terminal.

A new tab will open in your browser:

Pluto
Figure 2.

The Pluto notebook.

Download EasyChaos.jl from Github and then enter the path in Pluto below “Open from file:”. The notebook should look very much like the remainder of this post.

Easy Chaos

Pluto is a programming environment for Julia, designed to be interactive and helpful. Changes made anywhere in the notebook affect the entire notebook. Pluto is reactive so you don’t need to recalculate cells.

This is useful in studying chaos because small changes made to a function can have big effects. You can define a function, run it, and immediately see the effects.

A Linear Equation

The equation of a line is y=mx+by = mx + b where mm is the slope and bb is the yy-intercept. The slope is how much yy changes for each unit change in xx. For any starting point xx if you move one step to the right to x+1x + 1, the new yy-value becomes m(x+1)+b=mx+m+bm(x+1) + b = mx + m + b. Subtract y=mx+by = mx + b and all that’s left is mm.

A step by 11 in the xx-direction gives a change by mm in the yy-direction, which is the slope of the line.

If x=0x=0, then y=m×0+b=by = m \times 0 + b = b so the yy-intercept is bb. Since two points define a line, then start at bb on the yy-axis, move to x=1x=1 and add mm (or subtract mm if mm is negative) to get a second point. Draw your line through the two points.

The line can also be defined as a function,

f(x)=mx+b.f(x) = mx + b.

Let’s start by plotting some lines.

using-PlutoUI-Plots
Figure 3.

Load libraries into Pluto.

We need a way to change the values of mm and bb. First, we’ll make a slider for the slope mm and set the initial value to m=1m=1. The code for a slider is

@bind m Slider(-5:0.1:5, default=1)
using-PlutoUI-Plots
Figure 4.

Add a slider.

Next, we’ll make a slider for bb and set it to b=0b=0. Both can be varied between 5-5 and 55.

using-PlutoUI-Plots
Figure 5.

A linear function.

Notice that the origin of the plot is in the center. Move the sliders for mm and bb to get a new plot. Interesting, but not chaos.

A Little Chaos

Chaos is doing nearly the same thing over and over and expecting wildly different results. Let’s start with the linear equation above, with m=1m = 1 and b=1b = -1.

y=f(x)=x1.y = f(x) = x - 1.

Starting at x=1x = 1, the output of the function is y=0y=0. Using this new value as the input, f(0)=1f(0) = -1. If you do this several times you get the uninteresting sequence

x={1,0,1,2,3,}.x = \{1,0,-1,-2,-3, \ldots \}.

Instead of calculating x1=mx0+bx_1 = mx_0 + b and then x2=mx1+bx_2 = mx_1 + b and so on, we can do this in a function iter_f with inputs m,b,x0,nm,b,x_0,n where x0x_0 is the starting point, and nn is the number of iterations.

iter_f
Figure 6.

Code for an iterating function.

We can try this function with the equation y=x1y = x - 1 starting at x0=0.4961x_0 = 0.4961 and x1=0.4962x_1 = 0.4962 for n=10n = 10 iterations.

linear-results
Figure 7.

Results after 10 iterations of f(x)=x1f(x) = x-1.

Taking the difference between the outputs gives an uninteresting answer:

linear-differences
Figure 8.

Differences, y1y0y_1 - y_0.

Now, let’s change the equation slightly to f(x)=2x1f(x) = 2x - 1, and use the same starting points as before.

linear-2x
Figure 9.

Change the equation to f(x)=2x1f(x) = 2x-1.

The difference between these sequences is:

differences-2x
Figure 10.

Using f(x)=2x1f(x) = 2x-1 the differences double each iteration.

Multiply the differences by 1000010000 to make the changes more obvious:

differences-2x-10000
Figure 11.

Multiply the differences by 1000010000 to make the change clearer.

This is a little bit surprising, but notice that each number is twice the previous number. It looks like the sequence {20,21,22,29}\{ 2^0, 2^1, 2^2, \dots 2^9 \} and could be written as 2[0:9]2^{[0:9]}:

differences-pow2
Figure 12.

Powers of 22.

Let’s write out several iterations of y=mx+by = mx + b in terms of the starting point x0x_0.

x1=mx0+bx2=mx1+b=m(mx0+b)+b=m2x0+mb+b=m2x0+(m+1)bx3=mx2+b=m(m2x0+(m+1)b)+b=m3x0+(m2+m+1)bx4=mx3+b=m(m3x0+(m2+m+1)b)+b=m4x0+(m3+m2+m+1)b\begin{aligned} x_1 &= m x_0 + b \\ x_2 &= m x_1 + b = m (m x_0 + b) + b = m^2 x_0 + m b + b = m^2 x_0 + (m+1) b \\ x_3 &= m x_2 + b = m (m^2 x_0 + (m+1) b) + b = m^3 x_0 + (m^2 + m + 1) b \\ x_4 &= m x_3 + b = m (m^3 x_0 + (m^2 + m +1) b) + b = m^4 x_0 + (m^3 + m^2 + m + 1) b \\ \end{aligned}

The nthn^{th} iterate in terms of x0x_0 is

xn=mnx0+(mn1+mn2++1)b.x_n = m^n x_0 + (m^{n-1} + m^{n-2} + \cdots + 1) b.

What happens if we make a small change in x0x_0? If we change x0x_0 by a small amount, ϵ\epsilon, then

x~n=mn(x0+ϵ)+(mn1+mn2++1)b\tilde{x}_n = m^n (x_0 + \epsilon) + (m^{n-1} + m^{n-2} + \cdots + 1) b

and the difference between the two results after nn iterations is (the term with bb is the same for both)

x~nxn=mn(x0+ϵ)mnx0=mnϵ.\tilde{x}_n - x_n = m^n (x_0 + \epsilon) - m^n x_0 = m^n \epsilon.

So long as m>1|m| > 1 then we can make mnϵm^n \epsilon as big as we like by iterating enough times, no matter how small ϵ\epsilon is. A small change to the initial value of x0x_0 produces as big a change as you like in the final value, xnx_n.

There you have it. Chaos!

Even More Chaos

With a slightly more complicated equation, chaos becomes even more interesting. Instead of iterating a linear equation, let’s use the logistic function (see The Growing Gap for an application of logistic functions)

F(x)=cx(1x).F(x) = cx(1-x).

For c=4c = 4, this is an inverted parabola that has a maximum at (0.5,1.0)(0.5,1.0). If you draw the line y=xy=x it intersects the parabola at (0,0)(0,0) and (0.75,0.75)(0.75,0.75). This line is useful for following the iterations or orbits of the function.

Start with a point on the xx-axis, say x0=0.2x_0 = 0.2 (point A). Find the value of F(0.2)=4×0.2×(10.2)=0.64F(0.2) = 4 \times 0.2 \times (1 - 0.2) = 0.64. This point B is on the parabola above point A.

We want to iterate F(x)F(x) by using the yy-coordinate of B as the next input. We could calculate F(0.64)=4×0.64×(10.64)=0.9216F(0.64) = 4 \times 0.64 \times (1 - 0.64) = 0.9216 and we’ll want to write a Julia function to do that, but it’s also useful to see how the orbit evolves.

Draw a horizontal line until you get to the line y=xy=x at point C. Because C is on the 45°45 \degree line, then xC=yCx_C = y_C so the coordinates are C=(0.64,0.64)C = (0.64,0.64). The xx-coordinate of C is the yy-coordinate of B, so C becomes the next iterate. From C, draw a vertical line until you intersect the parabola at D which has coordinates (0.64,0.9216)(0.64,0.9216).

logistic orbits
Figure 13.

Initial orbits of the logistic function F(x)=cx(1x)F(x) = cx(1-x).

This process can be repeated as often as you like and will show the trajectory of the function F(x)F(x).

An Iterator for the Logistic Function

Like the iter_f function above, we can write a Julia function to iterate the logistic function. The inputs will be the starting point, x0x_0, the constant cc, and the number of iterations nn with a default value of n=100n=100.

This new function, iter_logistic will return the orbit, yy.

iter_logistic
Figure 14.

Julia code to output iterates of the logistic function.

Let’s try an example with x0=0.2x_0 = 0.2 and c=4c = 4.

Figure 15.

Initial iterates of the logistic function.

Now we can plot the trajectory:

plot-y-logistic
Figure 16.

Trajectory of the logistic with c=4c=4.

Try adjusting the value for cc. When c=4c=4 the plot seems pretty chaotic. What happens if you change x0x_0 from 0.20.2 to 0.190.19? This should show that the trajectory is sensitive to the initial value of x0x_0.

Next, try c=0.5c = 0.5 and set x0=0.1x_0 = 0.1. Is the plot chaotic? Is it sensitive to initial conditions?

Autocorrelation and Orbits

While the plot of the logistic function trajectory for c=4c=4 appears chaotic, you might be able to convince yourself that there are repeating patterns. The repetitions are far from identical, but if you imagine making a copy of the plot and shifting it over a bit, it might line up.

This is what autocorrelation does. The amount of the shift is controlled by a lag τ\tau. For each y[t]y[t] autocorrelation compares it to y[tτ]y[t - \tau]. The output of the autocorrelation function is a number rkr_k between 1-1 and 11, where

rk=t=τ+1N(y[t]μ)(y[tτ]μ)t=1N(y[t]μ)2r_k = \frac{\sum_{t = \tau + 1}^N (y[t] - \mu)(y[t - \tau] - \mu)}{\sum_{t=1}^N (y[t] - \mu)^2}

and μ\mu is the average of all the yy‘s. If rk=1r_k = 1 then y[t]=y[t+τ]y[t] = y[t + \tau] for every tt, rk=0r_k = 0 means there is no correlation, and rk=1r_k = -1 says that y[k]=y[k+τ]y[k] = -y[k + \tau].

A orbit is the plot of y[k]y[k] on the xx-axis and y[k+τ]y[k+\tau] on the yy-axis.

Let’s start by plotting the autocorrelation of the iterated logistic function for different values of τ\tau.

Rather than scrolling back up to set x0x_0 and cc, let’s just create two new variables and call them xPxP and cPcP.

plot-autocorrelation
Figure 17.

Autocorrelation of the logistic function.

Choose the autocorrelation step τ\tau:

Figure 18.

Orbit of the logistic with τ0=1\tau_0 = 1.

The autocorrelation plot seems about as chaotic as the logistic function, but when τ0=1\tau_0 = 1, the phase portrait is very smooth. In fact, the outline of the blue dots (y1) looks just like the logistic function!

This should be expected since yn+1=f(yn)y_{n+1} = f(y_n) for every nn.

Why choose τ0=1\tau_0 = 1? That choice of τ0\tau_0 plots the orbit of the function. The orange lines (y2) show which point follows the current location. On the left side, the path seems to be mostly upwards towards a nearby point. On the right, points are sent back to the left half except those near x=1x = 1.

Try setting τ0=0\tau_0 = 0 and you’ll see that the orbit becomes the line y=xy=x, because there’s perfect correlation between y[t]y[t] and y[t+τ0]y[t + \tau_0]. The plot doesn’t provide any insight. Try some other values of τ0\tau_0 to see what happens.

With Pluto, you can change the input values to a function to immediately see the effects.


Code for this article

EasyChaos.jl - A Pluto notebook to do some simple experiments with chaos theory.

Software

References and Further Reading

Image credits

Hero image from Introduction to Online Nonstochastic Control, Elad Hazan, Karan Singh, Research Gate, Figure 2.7: Lorenz attractor Nov 2022.