<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/chaos-bifurcation-diagrams-and-lyapunov-exponents-with-r-2/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://a-blog-from-sydney.blogspot.com/2016/12/chaos-bifurcation-diagrams-and-lyapunov.html"> A blog from Sydney</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>(The first part of this article can be <a target="_blank" rel="nofollow" href="http://a-blog-from-sydney.blogspot.it/2016/12/non-linear-dynamics-and-chaos-using-r-1.html">read here</a>)</p>
<div>Iteration of one-dimensional maps can generate stunning complexity and famed examples of chaotic behavior.  R can be used to get the flavor of this richness and reproduce some of the most famous pictures in the history of science, such as the bifurcation diagram of the logistic map or the representation of its Lyapunov exponents.</div>
<div></div>
<div>Given a one dimensional map depending on a parameter, a bifurcation diagram shows the stable structures (fixed point, cycles, attractors) visited by the dynamics for each value of the so called “bifurcation parameter”.  Resisting the temptation to use again and again the beloved logistic map, consider instead the following dynamical system
</div>
<div>
<div><i>x(t+1)=a cos x(t),  </i></div>
</div>
<div>taken from exercise 10.2.6, page 389 of <a target="_blank" rel="nofollow" href="https://books.google.it/books?id=JDQGAwAAQBAJ&printsec=frontcover&hl=it#v=onepage&q&f=false">Strogatz, “Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering”</a>.  The following script defines a function to get the bifurcation diagram for a given map (I modified some existing code, but I was not able to find and credit the exact source: https://web.stanford.edu/group/heeh/cgi-bin/web/node/59, say, has the same ideas. I defined a function and removed the use of “expression” and “evalf”).  See also the post at <a target="_blank" rel="nofollow" href="https://www.r-bloggers.com/dynamical-systems-mapping-chaos-with-r/">https://www.r-bloggers.com/dynamical-systems-mapping-chaos-with-r/</a>
<pre><code>bif_diagram &lt;- function(f=function(x,a) (a*x*(1-x)),alow=2.5,ahigh=4,<br>                        thinness=1000, transient=200, collect=200){<br>        # f function, parameter must be named a<br>    n &lt;- 1<br>    R &lt;- seq(alow,ahigh,length=thinness)<br>    data &lt;- matrix(0,collect,thinness+1)<br><br>    for(a in R){<br>      x &lt;- runif(1) # random initial condition<br>      ## first converge to attractor<br>      for(i in 1:transient){<br>        x &lt;- f(x,a)<br>      } # collect points on attractor<br>      for(i in 1:collect){<br>        x &lt;- f(x,a)<br>        data[i,n] &lt;- x<br>      }<br>  n &lt;- n+1<br>}<br><br>data &lt;- data[,1:thinness]<br>yrange &lt;- range(data)+c(-0.1,0.1)<br>plot(R,data[1,], pch=&quot;.&quot;, xlab=&quot;a&quot;, ylab=&quot;States&quot;,ylim=yrange)<br>for(i in 2:collect) points(R,data[i,],pch=&quot;.&quot;)<br>}<br></code></pre>
<p></div>
<div>
<div>Now, enter</div>
<pre><code>f &lt;- function(x,a) a* cos(x)<br>bif_diagram(f,alow=0.5,ahigh=4)<br></code></pre>
<div>
<table align="center" cellpadding="0" cellspacing="0">
<tbody>
<tr>
<td><a target="_blank" rel="nofollow" href="https://i0.wp.com/4.bp.blogspot.com/-9d1UMwxmZSM/WE1cSuxCDuI/AAAAAAAAC2E/Q8ex9z-DGFYwcHG4oNCegmNBAvoutCt_QCLcB/s1600/diario2016-002.jpg?ssl=1"><img src="https://i1.wp.com/4.bp.blogspot.com/-9d1UMwxmZSM/WE1cSuxCDuI/AAAAAAAAC2E/Q8ex9z-DGFYwcHG4oNCegmNBAvoutCt_QCLcB/s320/diario2016-002.jpg?resize=320%2C320&ssl=1" border="0" data-recalc-dims="1"></a></td>
</tr>
<tr>
<td>Bifurcation diagram for <i>f(x,a)=a cos x</i>, when <i>a</i> is the range [0.5,4].</td>
</tr>
</tbody>
</table>
<p>You can see that, for low values of the parameter <i>a</i>, there are unique fixed points or simple cycles.  Then, through  a series of (quite typical) period-doubling bifurcations, chaos appears and suddenly disappears when the parameter crosses 3.  <span>bif_diagram</span> has defaults parameters that can be adjusted to specify the parameter’s interval of interest or increase/decrease the resolution and legibility of the resulting plot.
<div></div>
<div>One of the most interesting signature of chaos is the divergence of orbits arbitrarily closed in their initial conditions.  Broadly speaking, orbits are stretched apart by some systems and a positive stretching rate is signaling the presence of chaos.  The Lyapunov Exponents (LE) is the average (exponential) growth rate of the divergence of initially nearby orbits.  LEs can be computed, for any value of the bifurcation parameter using the <span>lyap</span> function below, where is assumed that you have properly defined the dynamics <i>f</i> in advance.</div>
<div></div>
<pre><code>lyap &lt;- function(a,trans=300,num=1000){<br>    x0 &lt;- runif(1)<br>    for(time in 1:trans){<br>        x1 &lt;- f(x0,a=a);x0 &lt;- x1<br>    }<br>    sl &lt;- 0<br>    for(time in 1:num){<br>        x1 &lt;- f(x0,a=a);x0 &lt;- x1<br>        sl &lt;- sl+log(abs(grad(f,x1,a=a)))<br>    }<br>    sl/num<br>}<br></code></pre>
<p>In order to graph LEs define a sequence of values of the parameters, use <span>sapply</span> to get the sequence of corresponding exponents by applying <span>lyap</span> to the elements of the sequence and finally plot.  With the same map given above:</div>
<div></div>
<pre><code>a &lt;- seq(0.5,4,length=100)<br>ly &lt;- sapply(a,lyap)<br>plot(a,ly,t=&quot;l&quot;)<br></code></pre>
<table align="center" cellpadding="0" cellspacing="0">
<tbody>
<tr>
<td><a target="_blank" rel="nofollow" href="https://i0.wp.com/2.bp.blogspot.com/-ijoB2ivB23M/WE1cXcrqbvI/AAAAAAAAC2M/dCeQCP3dKgcUVkKsUf-_MOlwZsyy2mSvACEw/s1600/diario2016-004.jpg?ssl=1"><img src="https://i0.wp.com/2.bp.blogspot.com/-ijoB2ivB23M/WE1cXcrqbvI/AAAAAAAAC2M/dCeQCP3dKgcUVkKsUf-_MOlwZsyy2mSvACEw/s320/diario2016-004.jpg?resize=320%2C320&ssl=1" border="0" data-recalc-dims="1"></a></td>
</tr>
<tr>
<td>Lyapunov exponents for <i>f(x,a)=a cos x</i>, when <i>a</i> is the range [0.5,4].</td>
</tr>
</tbody>
</table>
<p>It can be seen, say, that when <i>a=2</i>, the LE is positive and chaos is in action [Check the bifurcation diagram to get the same intuition for that value of <i>a</i>].  Entering <span>lyap(2)</span> would indeed produce 0.15 (approximately).  In rough terms, this means that divergence of close initial points will be on average amplified by about 15% per iteration of the map over its domain.</div>
<div></div>
<div>To conclude:</div>
<div>
<ul>
<li>it is fun to try the same things with a slightly modified map: <i>x(t+1)=cos(a x(t))</i> (the parameter is “inside” the cos);</li>
<li>if you want to replicate two of the most well-known pictures related to the logistic map, type
<pre><code>f &lt;- function(x,a) a*x*(1-x)<br>bif_diagram(f,0,4)<br></code></pre>
</li>
<li>the graph of LE, seen on page 369 of Strogatz’s book is readily obtained through</li>
<pre><code>a &lt;- seq(3,4,len=100)<br>f &lt;- function(x,a) a*x*(1-x)<br>ly &lt;- sapply(a,lyap,num=2000)<br>plot(a,ly,t=&quot;l&quot;,ylim=c(-1,1));abline(h=0)<br></code></pre>
</ul>
<table align="center" cellpadding="0" cellspacing="0">
<tbody>
<tr>
<td><a target="_blank" rel="nofollow" href="https://i2.wp.com/1.bp.blogspot.com/-ElBFtuhLaHY/WE1cXeuVVdI/AAAAAAAAC2I/8aGT8fCPvrg5ziDKYbE5XK4ClGLRF9pMACEw/s1600/diario2016-006.jpg?ssl=1"><img src="https://i0.wp.com/1.bp.blogspot.com/-ElBFtuhLaHY/WE1cXeuVVdI/AAAAAAAAC2I/8aGT8fCPvrg5ziDKYbE5XK4ClGLRF9pMACEw/s320/diario2016-006.jpg?resize=320%2C320&ssl=1" border="0" data-recalc-dims="1"></a></td>
</tr>
<tr>
<td>Lyapunov exponens for the logistic map, when <i>a</i> is the range [3,4].</td>
</tr>
</tbody>
</table>
<p></div>
<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://a-blog-from-sydney.blogspot.com/2016/12/chaos-bifurcation-diagrams-and-lyapunov.html"> A blog from Sydney</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 title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/_Obs9JXKjQg"></a></div>