**schochastics**, and kindly contributed to R-bloggers)

This post is dedicated to the beautiful chaos created by double pendulums. I have seen

a great variety of animated versions, implemented with different tool but never in R.

Thanks to the amazing package `gganimate`

, it is actually not that hard to produce them in R.

```
library(tidyverse)
library(gganimate)
```

I am not going to attempt to explain the math behind the double pendulum. If you are interested

in the details, check out a complete walkthrough here.

The code presented here is a straightforward adaption from python.

First, we need to set up some basic constants and the starting conditions.

```
# constants
G <- 9.807 # acceleration due to gravity, in m/s^2
L1 <- 1.0 # length of pendulum 1 (m)
L2 <- 1.0 # length of pendulum 2 (m)
M1 <- 1.0 # mass of pendulum 1 (kg)
M2 <- 1.0 # mass of pendulum 2 (kg)
parms <- c(L1,L2,M1,M2,G)
# initial conditions
th1 <- 20.0 # initial angle theta of pendulum 1 (degree)
w1 <- 0.0 # initial angular velocity of pendulum 1 (degrees per second)
th2 <- 180.0 # initial angle theta of pendulum 2 (degree)
w2 <- 0.0 # initial angular velocity of pendulum 2 (degrees per second)
state <- c(th1, w1, th2, w2)*pi/180 #convert degree to radians
```

These are the parameters you need to change in order to produce different

pendulums. Just experiment a little!

The partial derivatives needed can be calculated with the following function.

```
derivs <- function(state, t){
L1 <- parms[1]
L2 <- parms[2]
M1 <- parms[3]
M2 <- parms[4]
G <- parms[5]
dydx <- rep(0,length(state))
dydx[1] <- state[2]
del_ <- state[3] - state[1]
den1 <- (M1 + M2)*L1 - M2*L1*cos(del_)*cos(del_)
dydx[2] <- (M2*L1*state[2]*state[3]*sin(del_)*cos(del_) +
M2*G*sin(state[3])*cos(del_) +
M2*L2*state[4]*state[4]*sin(del_) -
(M1 + M2)*G*sin(state[1]))/den1
dydx[3] <- state[4]
den2 <- (L2/L1)*den1
dydx[4] <- (-M2*L2*state[4]*state[4]*sin(del_)*cos(del_) +
(M1 + M2)*G*sin(state[1])*cos(del_) -
(M1 + M2)*L1*state[2]*state[2]*sin(del_) -
(M1 + M2)*G*sin(state[3]))/den2
return(dydx)
}
```

This function needs to be integrated. Luckily, there is the `odeintr`

packages that does the job.

The `start`

, `duration`

and `step_size`

parameters control the time for your pendulum.

In the below example, I choose to “swing” the pendulum for 30 seconds and the position

is recalculated every 0.1 seconds.

```
sol <- odeintr::integrate_sys(derivs,init = state,duration = 30,
start = 0,step_size = 0.1)
```

Now we just need to compute the x and y coordinates for both pendulums from the

angles \(\theta\) obtained from the integration.

```
x1 <- L1*sin(sol[, 2])
y1 <- -L1*cos(sol[, 2])
x2 <- L2*sin(sol[, 4]) + x1
y2 <- -L2*cos(sol[, 4]) + y1
df <- tibble(t=sol[,1],x1,y1,x2,y2,group=1)
```

The final data frame contains the exact position of the pendulums for each time step.

Animating the pendulums is straightforward with the package `gganimate`

.

```
ggplot(df)+
geom_segment(aes(xend=x1,yend=y1),x=0,y=0)+
geom_segment(aes(xend=x2,yend=y2,x=x1,y=y1))+
geom_point(size=5,x=0,y=0)+
geom_point(aes(x1,y1),col="red",size=M1)+
geom_point(aes(x2,y2),col="blue",size=M2)+
scale_y_continuous(limits=c(-2,2))+
scale_x_continuous(limits=c(-2,2))+
ggraph::theme_graph()+
labs(title="{frame_time} s")+
transition_time(t) -> p
pa <- animate(p,nframes=nrow(df),fps=20)
pa
```

We can also add some more details to the animation, like the trail of the second pendulum to track its path.

It turned out to be a bit more tricky then expected though. The trail needs to be added via a

secondary data.frame so that it can be animated with the `transition_time()`

.

I used the `lag()`

function to compute the trail from the last five time points.

```
tmp <- select(df,t,x2,y2)
trail <- tibble(x=c(sapply(1:5,function(x) lag(tmp$x2,x))),
y=c(sapply(1:5,function(x) lag(tmp$y2,x))),
t=rep(tmp$t,5)) %>%
dplyr::filter(!is.na(x))
```

I used the `shadow_mark()`

function to keep the past trails.

```
ggplot(df)+
geom_path(data=trail,aes(x,y),colour="blue",size=0.5)+
geom_segment(aes(xend=x1,yend=y1),x=0,y=0)+
geom_segment(aes(xend=x2,yend=y2,x=x1,y=y1))+
geom_point(size=5,x=0,y=0)+
geom_point(aes(x1,y1),col="red",size=M1)+
geom_point(aes(x2,y2),col="blue",size=M2)+
scale_y_continuous(limits=c(-2,2))+
scale_x_continuous(limits=c(-2,2))+
ggraph::theme_graph()+
labs(title="{frame_time} s")+
transition_time(t)+
shadow_mark(colour="grey",size=0.1,exclude_layer = 2:6)-> p
pa <- animate(p,nframes=nrow(df),fps=20)
pa
```

That’s it. Now you can play around with the constants (`L1,L2,M1,M2,G`

) and the initial conditions

(`th1,w1,th2,w2`

) to create your own chaotic pendulums.

Below is my personal favorite. 40 pendulums with nearly identical starting conditions.

Watch how quickly their paths diverge into pure chaos.

**leave a comment**for the author, please follow the link and comment on their blog:

**schochastics**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...