# Lorenz attractor animation (Plotly)

This is the first post in this blog.

`## [1] "Hello World!"`

Once for a class assignment, we were asked to control the Lorenz system. The instructor recommended us to use MATLAB for assignments, but since I’m inexperienced in MATLAB, I decided to use R to do the assignments, and used the package `plotly`

to make interactive 3D plots of phase portraits^{1} of the Lorenz system. Later, I made the iconic animation to demonstrate chaos – how minuscule differences in the initial condition of the system gets magnified exponentially, resulting into significantly different outcomes. Here, I’ll share how to make such an animation with `plotly`

.

First things first, what is the Lorenz system? Basically, it’s a system of 3 differential equations, a simplified model of air flow. The equations are: \[\begin{align} \frac{dx}{dt} &= \sigma(y - x) \\ \frac{dy}{dt} &= (\rho - z)x - y \\ \frac{dz}{dt} &= xy - \beta z \end{align}\]

Where \(t\) denotes time, \(x\) fluid motion, \(y\) horizontal temperature, \(z\) vertical temperature, and \(\sigma, \rho,\) and \(\beta\) parameters that specifies the particulars of the system that are constant through time. For \(\sigma = 10, \rho = 28, \beta = 8/3\), this system demonstrates chaotic behavior. Its phase portraits is the iconic butterfly, as we shall soon see.

Let’s solve the system of equations numerically first. There’re a number of packages that can numerically solve differential equations in R. One of them is `deSolve`

, which is the one I’m using here. Another one I have used is `pracma`

, which deliberately imitates the syntax of MATLAB. I don’t know enough math to solve the equations analytically, so I resort to a numeric solution.

```
# Package for solving differential equations
library(deSolve)
# Package for interactive 3D plots
library(plotly)
```

```
# Specify the equations
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))
})
}
```

What I mean by “numeric” is that time is discretized (here into steps of 0.01 time units), and since we know from the equations themselves how \(x,y,z\) change through time, we can use a variety of methods to compute what the values of \(x,y,z\) in the next time step will be given their values at present, starting at the initial condition (when time is 0). Here I used `ode45`

, which is a popular method.

```
# Specify parameters
parameters <- c(sigma = 10, beta = 8/3, rho = 28)
# Initial conditions, in the order x, y, z
state <- c(x = 2, y = 3, z = 4)
# Time over which the equations are numerically solved, 25 units
times <- seq(0, 25, 0.01)
# Solve the equations
sol <- as.data.frame(ode(y = state, times = times, func = lorenz, parms = parameters, method = "ode45"))
```

The object `sol`

is a data frames that contains a numerical solution

`head(sol)`

```
## time x y z
## 1 0.00 2.000000 3.000000 4.000000
## 2 0.01 2.117296 3.461536 3.960315
## 3 0.02 2.268655 3.951303 3.936245
## 4 0.03 2.453889 4.476829 3.930786
## 5 0.04 2.673533 5.045454 3.947722
## 6 0.05 2.928765 5.664405 3.991771
```

Now plot the approximate solution. Plotting in `plotly`

is a bit different from in `ggplot2`

; we need a formula, with the `~`

, to specify the variable given a data frame, kind of like the `aes`

in `ggplot2`

.

```
plot_ly(data = sol, x = ~x, y = ~y, z = ~z) %>%
add_paths(color = ~time) %>%
# Set where the camera is set by default
layout(scene = list(camera = list(eye = list(x = -1, y = 1, z = 0.25))))
```

Cool, it’s the iconic butterfly! We also see how the state of the system traverses between the two wings of the butterfly through time.

It’s actually not difficult to add an animated point traversing the phase space through time. I’m adding one point, that moves along the \(x,y,z\) as specified in the `sol`

data frame through time (`frame`

argument). In the function `animation_opts`

, `frame = 1`

means 1 millisecond per frame. `redraw = FALSE`

means you can’t toggle the camera perspective while the animation is playing; this improves performance. This may take a while to render.

```
# As usual
plot_ly(sol, x = ~x, y = ~y, z = ~z) %>%
add_paths(color = ~time) %>%
add_markers(frame = ~time, marker = list(size = 3, color = "black")) %>%
animation_opts(frame = 1, redraw = FALSE) %>%
layout(scene = list(camera = list(eye = list(x = -1, y = 1, z = 0.25))))
```