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: $\begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \\ \dot{x}_3 \end{bmatrix} = \begin{bmatrix} \sigma (x_2 - x_1) \\ -x_1 x_3 + r x_1 - x_2 \\ x_1 x_2 - b x_3 \end{bmatrix} .$

Here $x_1$ , $x_2$ and $x_3$ are our coordinates in three-dimensional space, $t$ represents time, and $\sigma$ , $r$ and $b$ 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 $x = (x_1,x_2,x_3)^T$ at time $t$ .

In [0]:

function dx = lorenz(x,t)           % The lorenz function takes the position vector x and the time t as inputs and produces the derivative of x at time t as the output dx.
dx = zeros(3,1);                    % We begin by initialising dx as the zero vector.
sigma = 10;                         % We need to specify values for our model parameters: our defaults will be sigma = 10; ...
r = 28;                             % r = 28; ...
b = 8/3;                            % and b = 8/3.
dx(1) = sigma*(x(2) - x(1));        % Reading off the Lorenz equations specified above, we specify how to compute the derivative of x1, ...
dx(2) = -x(1)*x(3) + r*x(1) - x(2); % x2, ...
dx(3) = x(1)*x(2) - b*x(3);         % and x3.
end                                 %

The solver we will use in Octave is lsode. This requires us to specify the function lorenz to compute the derivative of $x$ , 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.

In [0]:

time_span=linspace (0, 50, 10000);                 % We begin by specifying the time vector.  The linspace function allows us to specify a start time, ...
                                                   % end time and the number of equally-spaced points to use in that interval.
initial_conditions=[10.5 10.0 1.0];                  % We also need to specify the initial conditions, that is, where to start at t=0.

x = lsode(@lorenz, initial_conditions, time_span); % The lsode function does the work of finding x at each point of our time vector.
                                                   % Don't forget the @ symbol when specifying the function for computing the derivative.
plot(x(:,1), x(:,3))                               % If you want a 3d plot, then use plot3, but here we settle for plotting the x1 and x3 coordinates, ...
                                                   % which produces the famous "butterfly" shape.
xlabel('x1')                                       % Label your axes for completeness: first the x-axis; ...
ylabel('x3')                                       % then the y-axis.
title('The Lorenz attractor')                      % Finish by giving the figure a title.

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 $(1,1,1)^T$ and $(0.95, 1, 1)^T$ and plot the solutions on the same figure.

In [4]:

time_span=linspace (0, 30, 10000);                    % Specifying the time vector.
initial_conditions_a=[1.0 1.0 1.0];                   % Specifying the first set of initial conditions.
xa = lsode(@lorenz, initial_conditions_a, time_span); % Find the solution starting at those initial conditions.

initial_conditions_b=[0.95 1.0 1.0];                  % Specifying the second set of initial conditions.
xb = lsode(@lorenz, initial_conditions_b, time_span); % Find the solution starting at the new set of initial conditions.

plot(xa(:,1), xa(:,3))                                % Plot the first solution.
hold on                                               % This prevents the first plot being deleted when we ...
plot(xb(:,1), xb(:,3),'r')                            % plot the second solution.
hold off                                              %

xlabel('x1')                                          % Now add axis labels ... 
ylabel('x3')                                          % 
legend('x1(0)=1.0','x1(0)=0.95');                     % a legend, ...
title('Two trajectories')                             % and a title.

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 $x_1$ against time.

In [5]:

plot(time_span,xa(:,1));                   % This time we plot time against x1 for the first set of initial conditions, ...,
hold on                                    %
plot(time_span,xb(:,1),'r');               % and now the second set of initial conditions.
hold off                                   %
xlabel('t');                               %
ylabel('x1');                              %
legend('x1(0)=1.0','x1(0)=0.95');          %
title('Sensitivity to Initial Conditions') %

Note that the left-hand lobe of the Lorenz attractor corresponds to negative values of $x_1$ and the right-hand lobe corresponds to positive values of $x_1$ . From the above plot, we can see that the two solutions remain close until $t\approx 17$ after which one solution heads off to the left-hand lobe and the other does another loop around the right-hand lobe, but after $t\approx 25$ 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 $\sigma = 10$ , $r = 28$ and $b=8/3$ special?
Is the $x_1$ variable special? Do we see the same phenomenon when the initial condition for $x_2$ or $x_3$ 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 $t\approx 17$ . 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 $x_1(0)=1.0$ is accurate only to 1 decimal place. Then the real value could be anywhere in the interval $[0.95,1.05)$ . We can think of the real value as a random variable uniformly distributed in the interval $[0.95,1.05)$ . 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 $[0.95,1.05)$ .

In [6]:

time_span=linspace (0, 30, 10000);                     % We begin by finding and plotting the original solution.
                                                       %
initial_conditions=[1.0 1.0 1.0];                      %
x = lsode(@lorenz, initial_conditions, time_span);     %
                                                       %
plot(time_span,x(:,1));                                %
hold on                                                %
                                                       %
samples = unifrnd(0.95,1.05,10);                       % We now take 10 samples from our random distribution ...
                                                       %
for i=1:10                                             %
    initial_conditions=[samples(i) 1.0 1.0];           % and use these as our initial conditions.
    x = lsode(@lorenz, initial_conditions, time_span); % We find a solution for each sample ...
                                                       %
    plot(time_span,x(:,1));                            % and plot it.
end                                                    %
                                                       %
hold off                                               %
xlabel('t');                                           %
ylabel('x1');                                          %
title('Generating an ensemble forecast')               %

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 $(\pm \sqrt{b(r-1)},\pm \sqrt{b(r-1)},r-1)^T$ 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.

In [7]:

r = 28;                                   % We need to compute the centre of the left-hand and right-hand lobes ...
b = 8/3;                                  % so we repeat the definition of our parameters r and b.
RH = [sqrt(b*(r-1)),sqrt(b*(r-1)),r-1];   % This is not a good idea, in principle, but it will do for now.
LH = [-sqrt(b*(r-1)),-sqrt(b*(r-1)),r-1]; % Don't forget though that if you change the model parameters, you now need to change them in two places.
                                          %
num_points = 10000;                       % We declare the subdivisions of the time vector as a variable to make it easier to adjust.
sample_size = 10;                         % Similarly, we declare the size of the sample we will take from our random distribution.
time_span=linspace (0, 30, num_points);   %
f = zeros(sample_size,num_points);        % Each row of the matrix f will hold a sequence of 0s or 1s, according to whether the solution is closest to the
                                          % centre of the left-hand or right-hand lobe of the Lorenz attractor.
samples = unifrnd(0.95,1.05,sample_size); % We now randomly generate our set of initial conditions for x1.
                                          %
for i=1:sample_size                                    % For each initial condition in our sample ...
    initial_conditions=[samples(i) 1.0 1.0];           %
    x = lsode(@lorenz, initial_conditions, time_span); % we compute the solution ...
                                                       %
    for j=1:num_points                                 % and at each time ...
        if norm(x(j)-RH) < norm(x(j)-LH)               % if the centre of the right-hand lobe is closest, ... 
            f(i,j) = 1;                                % we replace the 1 in f by a 0.
        end                                            %
    end                                                %
end                                                    %
                                                       %
probs = sum(f,1)/sample_size;                          % Now we compute the probability of orbiting the right-hand lobe as a function of time.
                                                       %
area(time_span,probs)                                  % We will do an area plot to shade the area below the probability curve to make it easier to read.
xlabel('t');
ylabel('p');
title('Probability p of solution orbiting right-hand lobe')

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!