*This post was semi automatically converted from blogdown to Quarto and may contain errors. The original can be found in the archive.*

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

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.

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

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

Now we just need to compute the x and y coordinates for both pendulums from the angles θ*θ**θ* obtained from the integration.

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

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))+
::theme_graph()+
ggraphlabs(title="{frame_time} s")+
transition_time(t) -> p
<- animate(p,nframes=nrow(df),fps=20)
pa 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.

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

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))+
::theme_graph()+
ggraphlabs(title="{frame_time} s")+
transition_time(t)+
shadow_mark(colour="grey",size=0.1,exclude_layer = 2:6)-> p
<- animate(p,nframes=nrow(df),fps=20)
pa 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.

## Reuse

## Citation

```
@online{schoch2018,
author = {Schoch, David},
title = {Beautiful {Chaos:} {The} {Double} {Pendulum}},
date = {2018-11-22},
url = {http://blog.schochastics.net/posts/2018-11-22_beautiful-chaos-the-double-pendulum},
langid = {en}
}
```