An introduction to chaos and predictability
Overview
In this worksheet we will explore some key ideas underpinning the mathematical theory of chaos and predictability. In particular:
we will consider the sensitivity of a solution to the initial conditions of the problem;
we will examine the consequences of this for numerical simulation;
and we will look at ensemble forecasting as a means of quantifying the uncertainity in our predictions.
Using the worksheet
You should make your own copy of this worksheet in your own workspace and open that.
Double-click on cells to edit them.
To execute a cell, click the "Run" button on a selected cell or press Shift+Enter.
Every time you make a change, you will need to execute all the cells that need updating.
The Lorenz equations
To aid us in our study of chaos, we will work with a system of differential equations, called the Lorenz equations:
Here , and are our coordinates in three-dimensional space, represents time, and , and are the parameters of the model.
In general, we cannot write down a closed formula for the solution of the Lorenz equations. We can attempt to approximate the solution numerically, however, which we shall do in this worksheet by using the programming language Octave.
We begin by writing a function lorenz that calculates the derivative of at time .
The solver we will use in Octave is lsode. This requires us to specify the function lorenz to compute the derivative of , the initial conditions and a vector of times at which the solution should be computed. It outputs a matrix of values corresponding to the approximation of x at each point in the time vector.
This is the characteristic butterfly shape of the Lorenz attractor and solutions orbit the lobes (wings) of the attractor.
Exercises
Experiment with changing the initial conditions.
Experiment with changing the model parameters.
Experiment with changing the number of points in the time vector.
Make a note of any interesting findings!
Sensitivity to initial conditions
We are now going to analyse a little more closely what happens when we vary the initial conditions. Let's take two sets of initial conditions and and plot the solutions on the same figure.
That plot doesn't seem to tell us a lot, but it does appear that the two solutions are close for a while and then start to separate. We can see this more clearly if we plot against time.
Note that the left-hand lobe of the Lorenz attractor corresponds to negative values of and the right-hand lobe corresponds to positive values of . From the above plot, we can see that the two solutions remain close until after which one solution heads off to the left-hand lobe and the other does another loop around the right-hand lobe, but after they seem to be doing their own thing.
Note that this divergent behaviour is a feature of this system of equations and is one of the characteristic elements of chaotic dynamical systems.
Exercises
How does the initial separation correspond to the time at which the two solutions separate?
Is the point (1,1,1) special? Or does this happen for other initial conditions e.g. (-1,3,-2)?
Are the model parameters , and special?
Is the variable special? Do we see the same phenomenon when the initial condition for or is adjusted?
Make a note of any interesting findings!
Predictability
Two immediate problems emerge from this sensitivity to initial conditions.
Measurement error
When modelling real-world situations, we will never know the initial conditions perfectly. Typically, we report them to a specified degree of precision (e.g. 2 decimal places, 5 significant figures), but for systems such as described by the Lorenz equations, this initial imperfection can be swiftly amplified.
Numerical error
In fact, this isn't just a problem for the initial measurement: computers can only compute and store numbers to finite precision. There are potential rounding errors for every calculation performed at every step in the program. Even if we did know the initial conditions perfectly, errors will still be introduced and amplified!
All is not lost, however. By comparing two solutions with a small initial separation, we were able to get some handle on how confident we could be in our predictions. For example, we seem to know pretty reliably what will happen until . We can make this more precise and start to quantify the uncertainty in our predictions using ensemble forecasting.
Ensemble forecasting
Suppose that our initial measurement of is accurate only to 1 decimal place. Then the real value could be anywhere in the interval . We can think of the real value as a random variable uniformly distributed in the interval . Instead of using only the one value 1.0, or comparing two values as above, we could take any number of samples from our random distribution and compute the corresponding solutions.
Suppose what we are interested in is whether the solution is orbiting the left-hand lobe or right-hand lobe of the Lorenz attractor. Then, we can take 10 samples (or 100, 1000, 1000000 samples) and put a value on the certainty of our forecast at each point in time.
In the next example, therefore, we take 10 solutions produced from 10 initial conditions uniformly distributed in the interval .
Again we can see that for a while the solutions remain close together and we could confidently predict which lobe they are orbiting, but eventually their behaviour is divergent.
The centres of the two lobes are at and we can use this information to quantify the uncertainty in our forecast. We will use our solutions to estimate the probability of orbiting the right-hand lobe at each point in time.
The same principle is used in weather forecasting e.g. to generate the probability of rain. We run our simulation many times with small random variations and then use that ensemble of data to make our forecast, specifying the corresponding probabilities.
Exercises
We chose a uniform random distribution based on rounding the initial condition to 1 decimal place. Suppose instead we believe that the values are normally distributed about 1. How does that affect the calculations?
Experiment with taking larger samples from the random distribution. (But increase this number gradually - 100 samples takes a couple of minutes to compute so you probably don't want to go much above this!)
Experiment with varying the initial separation (e.g. suppose the initial measurement is accurate to 5dp). How does the time horizon of your forecast change (that is, the time beyond which you are no longer confident in your forecast)?
Make a note of any interesting findings!