<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/beautiful-chaos-the-double-pendulum/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://blog.schochastics.net/post/beautiful-chaos-the-double-pendulum/"> schochastics</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>This post is dedicated to the beautiful chaos created by double pendulums. I have seen<br>
a great variety of animated versions, implemented with different tool but never in R.<br>
Thanks to the amazing package <code>gganimate</code>, it is actually not that hard to produce them in R.</p>
<p><img src="https://i2.wp.com/scienceworld.wolfram.com/physics/dimg270.gif?w=456" data-recalc-dims="1"></p>
<pre><code>library(tidyverse)
library(gganimate)</code></pre>
<p>I am not going to attempt to explain the math behind the double pendulum. If you are interested<br>
in the details, check out a complete walkthrough <a rel="nofollow" href="http://scienceworld.wolfram.com/physics/DoublePendulum.html" target="_blank">here</a>.<br>
The code presented here is a straightforward adaption from <a rel="nofollow" href="https://matplotlib.org/examples/animation/double_pendulum_animated.html" target="_blank">python</a>.</p>
<p>First, we need to set up some basic constants and the starting conditions.</p>
<pre><code># constants
G  &lt;-  9.807  # acceleration due to gravity, in m/s^2
L1 &lt;-  1.0    # length of pendulum 1 (m)
L2 &lt;-  1.0    # length of pendulum 2 (m)
M1 &lt;-  1.0    # mass of pendulum 1 (kg)
M2 &lt;-  1.0    # mass of pendulum 2 (kg)

parms &lt;- c(L1,L2,M1,M2,G)

# initial conditions
th1 &lt;-  20.0  # initial angle theta of pendulum 1 (degree)
w1  &lt;-  0.0   # initial angular velocity of pendulum 1 (degrees per second)
th2 &lt;-  180.0 # initial angle theta of pendulum 2 (degree)
w2  &lt;-  0.0   # initial angular velocity of pendulum 2 (degrees per second)

state &lt;- c(th1, w1, th2, w2)*pi/180  #convert degree to radians</code></pre>
<p>These are the parameters you need to change in order to produce different<br>
pendulums. Just experiment a little!</p>
<p>The partial derivatives needed can be calculated with the following function.</p>
<pre><code>derivs &lt;- function(state, t){
  L1 &lt;- parms[1]
  L2 &lt;- parms[2]
  M1 &lt;- parms[3]
  M2 &lt;- parms[4]
  G  &lt;- parms[5]
  
  dydx    &lt;-  rep(0,length(state))
  dydx[1] &lt;-  state[2]
  
  del_ &lt;-  state[3] - state[1]
  den1 &lt;-  (M1 + M2)*L1 - M2*L1*cos(del_)*cos(del_)
  dydx[2]  &lt;-  (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] &lt;-  state[4]
  
  den2 &lt;-  (L2/L1)*den1
  dydx[4]  &lt;-  (-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)
}</code></pre>
<p>This function needs to be integrated. Luckily, there is the <code>odeintr</code> packages that does the job.<br>
The <code>start</code>, <code>duration</code> and <code>step_size</code> parameters control the time for your pendulum.<br>
In the below example, I choose to “swing” the pendulum for 30 seconds and the position<br>
is recalculated every 0.1 seconds.</p>
<pre><code>sol &lt;- odeintr::integrate_sys(derivs,init = state,duration = 30,
                              start = 0,step_size = 0.1)</code></pre>
<p>Now we just need to compute the x and y coordinates for both pendulums from the<br>
angles <span>\(\theta\)</span> obtained from the integration.</p>
<pre><code>x1 &lt;-  L1*sin(sol[, 2])
y1 &lt;-  -L1*cos(sol[, 2])
  
x2 &lt;- L2*sin(sol[, 4]) + x1
y2 &lt;- -L2*cos(sol[, 4]) + y1
  
df &lt;- tibble(t=sol[,1],x1,y1,x2,y2,group=1)</code></pre>
<p>The final data frame contains the exact position of the pendulums for each time step.<br>
Animating the pendulums is straightforward with the package <code>gganimate</code>.</p>
<pre><code>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=&quot;red&quot;,size=M1)+
  geom_point(aes(x2,y2),col=&quot;blue&quot;,size=M2)+
  scale_y_continuous(limits=c(-2,2))+
  scale_x_continuous(limits=c(-2,2))+
  ggraph::theme_graph()+
  labs(title=&quot;{frame_time} s&quot;)+
  transition_time(t) -> p

pa &lt;- animate(p,nframes=nrow(df),fps=20)
pa</code></pre>
<p><img src="https://i2.wp.com/blog.schochastics.net/post/2018-11-22-beautiful-chaos-the-double-pendulum_files/figure-html/animate_1-1.gif?w=456" data-recalc-dims="1"></p>
<p>We can also add some more details to the animation, like the trail of the second pendulum to track its path.<br>
It turned out to be a bit more tricky then expected though. The trail needs to be added via a<br>
secondary data.frame so that it can be animated with the <code>transition_time()</code>.<br>
I used the <code>lag()</code> function to compute the trail from the last five time points.</p>
<pre><code>tmp &lt;- select(df,t,x2,y2)
trail &lt;- 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))</code></pre>
<p>I used the <code>shadow_mark()</code> function to keep the past trails.</p>
<pre><code>ggplot(df)+
  geom_path(data=trail,aes(x,y),colour=&quot;blue&quot;,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=&quot;red&quot;,size=M1)+
  geom_point(aes(x2,y2),col=&quot;blue&quot;,size=M2)+
  scale_y_continuous(limits=c(-2,2))+
  scale_x_continuous(limits=c(-2,2))+
  ggraph::theme_graph()+
  labs(title=&quot;{frame_time} s&quot;)+
  transition_time(t)+
  shadow_mark(colour=&quot;grey&quot;,size=0.1,exclude_layer = 2:6)-> p

pa &lt;- animate(p,nframes=nrow(df),fps=20)
pa</code></pre>
<p><img src="https://i2.wp.com/blog.schochastics.net/post/2018-11-22-beautiful-chaos-the-double-pendulum_files/figure-html/animate_2-1.gif?w=456" data-recalc-dims="1"></p>
<p>That’s it. Now you can play around with the constants (<code>L1,L2,M1,M2,G</code>) and the initial conditions<br>
(<code>th1,w1,th2,w2</code>) to create your own chaotic pendulums.</p>
<p>Below is my personal favorite. 40 pendulums with nearly identical starting conditions.<br>
Watch how quickly their paths diverge into pure chaos.</p>
<p><img src="https://i1.wp.com/blog.schochastics.net/img/pend40.gif?w=456" data-recalc-dims="1"></p>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://blog.schochastics.net/post/beautiful-chaos-the-double-pendulum/"> schochastics</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/tv5oAsSj-fk" height="1"></a></div>