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}
}