Chaos and Butterfly Effect

The scientific theory that a single occurrence, no matter how small, can change the course of the universe forever

INTRODUCTION

Dynamical systems theory is an area of mathematics i.e., used to describe the behavior of complex dynamical systems, usually by employing differential equations or difference equations. When differential equations are employed, the theory is called continuous dynamical systems. From a physical point of view, continuous dynamical systems are a generalization of classical mechanics, a generalization where the equations of motion are postulated directly and are not constrained to be Euler–Lagrange equations of a least action principle. When difference equations are employed, the theory is called discrete dynamical systems. When the time variable runs over a set that is discrete over some intervals and continuous over other intervals or is any arbitrary time set such as a Cantor set, one gets dynamic equations on time scales. Some situations may also be modeled by mixed operators, such as differential-difference equations.

This theory deals with the long-term qualitative behavior of dynamical systems, and studies the nature and when possible the solutions of the equations of motion of systems that are often primarily mechanical or otherwise physical, such as planetary orbits and the behavior of electronic circuits, as well as systems that arise in biology, economics, and elsewhere. Much of modern research is focused on the study of chaotic systems.

CHAOS THEORY

Chaos theory is a branch of mathematics focusing on the study of chaotic dynamical systems whose random states of disorder and irregularities are governed by underlying patterns and deterministic laws that are highly sensitive to initial conditions. Chaos theory is an interdisciplinary theory stating that, within the apparent randomness of chaotic complex systems, there are underlying patterns, interconnectedness, constant feedback loops, repetition, self-similarity, fractals, and self-organization. The butterfly effect, an underlying principle of chaos, describes how a small change in one state of a deterministic nonlinear system can result in large differences in a later state (meaning that there is a sensitive dependence on initial conditions). Small differences in initial conditions, such as those due to errors in measurements or due to rounding errors in numerical computation, can yield widely diverging outcomes for such dynamical systems, rendering long-term prediction of their behavior impossible in general. This can happen even though these systems are deterministic, meaning that their future behavior follows a unique evolution and is fully determined by their initial conditions, with no random elements involved. In other words, the deterministic nature of these systems does not make them predictable. This behavior is known as deterministic chaos, or simply chaos. The theory was summarized by Edward Lorenz as Chaotic behavior exists in many natural systems, including fluid flow, heartbeat irregularities, weather, and climate. It also occurs spontaneously in some systems with artificial components, such as the stock market and road traffic. This behavior can be studied through the analysis of a chaotic mathematical model, or analytical techniques such as recurrence plots and Poincare maps. Chaos theory has applications in a variety of disciplines, including meteorology, anthropology, sociology, failed verification environmental science, computer science, engineering, economics, ecology, pandemic crisis management, and philosophy. The theory formed the basis for such fields of study as complex dynamical systems, edge of chaos theory, and self-assembly processes.

LORENTZ SYSTEM

At The Beginning of The Revolution

In August 2017, people across North America planned for a once-in-a-lifetime opportunity to watch a total solar eclipse. Here’s a reference to the 2017 eclipse in an article about a previous one, 99 years earlier.

Reference to the 2017 eclipse in an article

For centuries, humans have used math and science to make predictions about our universe. Certainly, by the end of the 19th century, scientists felt that we were living in a Cartesian universe, in which, if you had some omniscient being who knew the positions and velocities of all the particles in the universe, in principle, that being would know what would happen for all time.

Greek Eclipse Calculator ( about 100 BC )

But if that’s true, why are some things in our world still so hard to predict, even in an era of supercomputers and big data? If we can predict an eclipse a century in advance, why can we only predict the weather about a week or two in advance?

Mayan Eclipse Table ( About 1200 AD )

We know this system is based on the same physical laws as the rest of our universe, but it seems random. The study of systems that look random, but aren’t, is called ” Chaos Theory”, one of its key pioneers was an MIT meteorologist named Edward Norton Lorenz. During World War II and the decades that followed, Computers allowed scientists to take major steps forward in weather prediction.

Experiments and Results

In 1961, Lorenz was working with a set of equations that represented a simplified version of the atmosphere, with 12 variables that changed over time. The computer calculated the values for each moment in time by applying the equations to the values from the previous moment, and Lorenz watched how they progressed. He needed to be sure that like real weather, they did not follow a repeating pattern. In, some cases, Lorenz’s team needed to reprint part of the simulation, so they would restart it from an earlier point by entering the numbers from an earlier printout as start values. Since they used the same numbers and equations, they expected to see the same results.

But sometimes they didn’t. In some cases, the second run would look very different from the first as it progressed. At first, Lorenz thought the computer was malfunctioning, but he eventually found the source of the error. While the computer’s memory stores numbers with six decimal places, it printed them with only three decimal places, so the values he entered contained tiny errors.

Lorenz was observing that in some systems, tiny changes in the present can lead to big changes in the future. 

Courtesy of Glenn Flierl & Lodovica Illri( MIT )

Today we can create models of the atmosphere that go far beyond Lorenz’s columns of numbers. In this model from MIT, the white line represent a bunch of balloons released from approximately the same point. Eventually, the Balloons diverge along very different paths, because of the small differences in their starting points. If you change the initial conditions ever so slightly, infinitesimally slightly. The two trajectories seem to go along with each other and then they diverge exponentially fast. So, that is chaos. It means that you have to know a system with infinite precision to be able to predict it infinitely in time. 

Sensitive Dependence on Initial Conditions

In 1903, French mathematician Henri Poincare observed something similar while studying the problem of three objects in space, all affected by each other’s gravity. He was always looking for what’s the most simple set of equations that exhibit chaotic properties. Lorenz found a set of three equations, a simplified version of equations used to model convection, with only three variables that change over time. It is notable for having chaotic solutions for certain parameter values and initial conditions.

Solving Initial Value Differential Equations in R

1. A Simple ODE: Chaos in the Atmosphere

The Lorenz equations were the first chaotic dynamic system to be described. They consist of three differential equations that were assumed to represent the idealized behavior of the earth’s atmosphere. We use this model to demonstrate how to implement and solve differential equations using R programming. The Lorenz model describes the dynamics of three state variables, \(x,\ y,\ and\ z\). The model equations are:

\[\dfrac{dx}{dt}=\sigma(y-x)\]

\[ \dfrac{dy}{dt}=x(\rho-z)-y \]

\[\dfrac{dz}{dt}=xy-\beta z\]

with the initial conditions : \(x(0)=y(0)=z(0)=1\)

The equations relate the properties of a two-dimensional fluid layer uniformly warmed from below and cooled from above. In particular, the equations describe the rate of change of three quantities with respect to time: \(x\) is proportional to the rate of convection, \(y\) to the horizontal temperature variation, and \(z\) to the vertical temperature variation. The constants \(\sigma\), \(\rho\), and \(\beta\) are system parameters proportional to the Prandtl number, Rayleigh number, and certain physical dimensions of the layer itself where \(\sigma\), \(\rho\), and \(\beta\) are three parameters, with values of\(\sigma = 10,\ \beta= \frac{8}{3}, \rho=28\). Implementation of an IVP ODE in R can be separated into two parts: ‘the model specification’ and the ‘model application’.

The model specification consists of:

• Defining model parameters and their values,

• Defining model state variables and their initial conditions,

• Implementing the model equations that calculate the rate of change (e.g. \(\frac{dx}{dt}\)) of the state variables.

The model application consists of:

• Specification of the time at which model output is wanted,

• Integration of the model equations (uses R-functions from deSolve),

• Plotting of model results.

1.1. MODEL SPECIFICATION

Model parameters:

There are three model parameters: s, r, and b that are defined first. Parameters are stored as a vector with assigned names and values:

parameters <- c(s=10,r=28,b=8/3)

State Variables:

The three state variables are also created as a vector, and their initial values are given:

initial_state <-
  c(x = -8.60632853,
    y = -14.85273055,
    z = 15.53352487)

Model Equation

The model equations are specified in a function (Lorenz) that calculates the rate of change of the state variables. Input to the function is the model time (t, not used here, but required by the calling routine), and the values of the state variables (state) and the parameters, in that order. This function will be called by the R routine that solves the differential equations (here we use ode, see below). The code is most readable if we can address the parameters and state variables by their names. As both parameters and state variables are ‘vectors’, they are converted into a list. The statement with (as. list (c(state, parameters)), …) then makes available the names of this list. The main part of the model calculates the rate of change of the state variables. At the end of the function, these rates of change are returned, packed as a list. Note that it is necessary to return the rate of change in the same ordering as the specification of the state variables. This is very important. In this case, as state variables are specified \(x\) first, then \(y\) and \(z\), the rates of changes are returned as \(dx,dy,dz\).

lorenz <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dx <- sigma * (y - x)
    dy <- x * (rho - z) - y
    dz <- x * y - beta * z
    
    list(c(dx, dy, dz))
  })
}

1.2. MODEL APPLICATION:

Time Specification

We run the model for 100 days, and give output at 0.01 daily intervals. R’s function seq( ) creates the time sequence:

times <- seq(0, 500, length.out = 125000)

Model Integration

The model is solved using ’deSolve’ function ode, which is the default integration routine. Function ode takes as input, a.o. the state variable vector (y), the times at which output is required (times), the model function that returns the rate of change (func) and the parameter vector (parms). Function ode returns an object of class deSolve with a matrix that contains the values of the state variables (columns) at the requested output times.

lorenz_ts <-
  ode(
    y = initial_state,
    times = times,
    func = lorenz,
    parms = parameters,
    method = "lsoda"
  ) %>% as_tibble()

lorenz_ts[1:10,]

Plotting Results:

Finally, the model output is plotted. We use the plot method designed for objects of class deSolve, which will neatly arrange the figures in two rows and two columns; before plotting, the size of the outer upper margin (the third margin) is increased. We access to \(x\),\(y\),\(z\) respectively a univariate time series. As the time interval used to numerically integrate the differential equations was rather tiny, we just use every tenth observation.

library(deSolve)
library(tidyverse)

parameters <- c(sigma = 10,
                rho = 28,
                beta = 8/3)

initial_state <-
  c(x = -8.60632853,
    y = -14.85273055,
    z = 15.53352487)

lorenz <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dx <- sigma * (y - x)
    dy <- x * (rho - z) - y
    dz <- x * y - beta * z
    
    list(c(dx, dy, dz))
  })
}

times <- seq(0, 500, length.out = 125000)

lorenz_ts <-
  ode(
    y = initial_state,
    times = times,
    func = lorenz,
    parms = parameters,
    method = "lsoda"
  ) %>% as_tibble()
obs <- lorenz_ts %>%
  select(time, x) %>%
  filter(row_number() %% 10 == 0)
# lorenz_ts[1:10,]

ggplot(obs, aes(time, x)) +
  geom_line() +
  coord_cartesian(xlim = c(0, 100)) +
  theme_classic()

obs <- lorenz_ts %>%
  select(time, y) %>%
  filter(row_number() %% 10 == 0)

ggplot(obs, aes(time, y)) +
  geom_line() +
  coord_cartesian(xlim = c(0, 100)) +
  theme_classic()

obs <- lorenz_ts %>%
  select(time, z) %>%
  filter(row_number() %% 10 == 0)

ggplot(obs, aes(time, z)) +
  geom_line() +
  coord_cartesian(xlim = c(0, 100)) +
  theme_classic()

Lorenz Attractor

Let’s try graphing them in three dimensions on a modern computer. The location of the first point is arbitrary. With each point after that, the computer steps forward in time by calculating the position of the next point based on the position of the previous one. It seems to change direction at random, and it never arrives at the same point twice. If it did, that would be the beginning of a repeating pattern. This type of system is now known as a “strange attractor” and this particular one is known as the Lorenz attractor.

library(deSolve)
library(plotly)

parameters <- c(rho = 28, sigma = 10, beta = 8/3)
state <- c(x = 1, y = 1, z = 1)
time <- seq(0, 130, by = 0.01)

lorenz <- function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    dx <- sigma * (y - x)
    dy <- x * (rho - z) - y
    dz <- x * y - beta * z
    list(c(dx, dy, dz))
  })
}


out <- as.data.frame(ode(y = state, times = time, func = lorenz, parms = parameters))

plot_ly(out, x = ~x, y = ~y, z = ~z, type = 'scatter3d', mode = 'lines', color = ~time, colors = 'Blue')%>% 
  # Set where the camera is set by default
  layout(scene = list(camera = list(eye = list(x = -1, y = 1, z = 0.25))))