View source for G Graphing Treat

= Introduction and Setup =

Creating graphs can be increadibly pleasing, but often it will be very simple graphs that support the story you are telling with your data. Here I want to give you the opportunity to go just a little crazy and be in awe of a graph you create. Just indulge me!

= Lorenz Attractor =

Think of a little fly which moves, as they do, in a three dimensional space (think in terms of x, y and z dimensions). And let's say that you are actually mind controlling this fly. As you are a mathematical fiend you write down a very specific rule how the position of that fly changes from second to second. That means you have to write down differential equations.

Here is the equation for X: <math>dx = \sigma(y-x) dt</math> Here is the equation for Y: <math>dy = (x(\rho-z)-y) dt</math> Here is the equation for Z: <math>dz = (xy-\beta z) dt</math>

these equations are called the [https://en.wikipedia.org/wiki/Lorenz_system Lorenz equations] and they give precise instructions how the coordinates for <math>x</math>, <math>y</math> and <math>z</math> change in every time interval <math>dt</math>. So let's implement this in R (The code for this example comes from [http://fractalswithr.blogspot.co.uk/2007/04/lorenz-attractor.html?m=1 Atte Tenkanen]). We will essentially let the fly fly and record its position as it flies and then we will plot the flight path in a cool picture.

Here is the prep, as we need the ''scatterplot3d'' package.

<pre class="r"># install.packages(&quot;scatterplot3d&quot;) # Install if needed.
library(scatterplot3d)</pre>
First, let's set some parameters and the size of the time step. Later you can change these around and see how the result changes.

<pre class="r">sigma=5; rho=28; beta=8/3; dt=0.01</pre>
Then we need the initial coordinates

<pre class="r">X=0.01; Y=0.01; Z=0.01</pre>
Now let's set for how many time steps <code>n</code> we want to follow our imaginary fly.

<pre class="r">n=60000</pre>
Next we define a matrix with <code>n</code> rows and 3 columns (representing <code>X</code>, <code>Y</code> and <code>Z</code>), one row for each moment in time. In here we will save the positional coordinates of our fly at each of our <code>n</code> points of time.

<pre class="r">XYZ=array(0,dim=c(n,3))
t=0   # set time to 0</pre>
Here is the core of the code. We are using a <code>for</code> loop to iterate through our <code>n</code> moments in time. You should recognise the above Lorenz equations in the code

<pre class="r">for(i in 1:n)
{

X1=X; Y1=Y; Z1=Z  # Set X1, Y1 and Z1 to the last position

X=X1+(sigma*(Y1-X1))*dt
Y=Y1+(X1*(rho-Z1)-Y1)*dt
Z=Z1+(X1*Y1-beta*Z1)*dt
#t=t+dt
XYZ[i,]=c(X,Y,Z)  # save the current positional coordinates in the ith row
}</pre>
All the calculations are been done now. All we do now is to define three vectors <code>xpos</code>, <code>ypos</code> and <code>zpos</code> with the relevant positions. We draw that info from the <code>XYZ</code> matrix

<pre class="r">xpos = XYZ[,1]
ypos = XYZ[,2]
zpos = XYZ[,3]</pre>
Now we can plot the positions of our fly through time. We need a three dimensional scatterplot (which is why we had to load the ''scatterplot3d'' package). I will not go through all the options that have been chosen here. You can play around with them a bit to find out what they do.

<pre class="r">s3d&lt;-scatterplot3d(xpos, ypos, zpos, highlight.3d=TRUE, xlab=&quot;x&quot;,ylab=&quot;y&quot;,zlab=&quot;z&quot;,type=&quot;p&quot;, col.axis=&quot;blue&quot;, angle=55, col.grid=&quot;lightblue&quot;, scale.y=0.7, pch=&quot;.&quot;, xlim=c(min(xpos),max(xpos)), ylim=c(min(ypos),max(ypos)), zlim=c(min(zpos),max(zpos)))

s3d.coords &lt;- s3d$xyz.convert(xpos,ypos,zpos)

title(&quot;Lorenz Attractor&quot;)</pre>
[[File:figure/unnamed-chunk-8-1.png|frame|none|alt=|caption plot of chunk unnamed-chunk-8]]

As you can see, your fly is drawing pretty angled butterfly patterns in the air. See how changes to the parameters change this pattern.

= The Mandelbrot Set =

This is a particular treat and you may well have seen the resulting images. Nothing short but beautiful. This code is replicated from [http://tolstoy.newcastle.edu.au/R/help/05/10/13198.html here]

First ensure that you have installed the following two libraries (you will only have to do that once!)

<pre class="r">install.packages(&quot;caTools&quot;)
install.packages(&quot;fields&quot;)</pre>
Here is what you need to do

<pre class="r">library(fields)  # for tim.colors
library(caTools) # for write.gif
m &lt;- 1600          # grid size
C &lt;- complex( real=rep(seq(-1.8,0.6, length.out=m), each=m ),imag=rep(seq(-1.2,1.2, length.out=m), m ) )
C &lt;- matrix(C,m,m)</pre>
You really do not have to understand any details of the mathematics that is going on here. You can find a discussion of the [https://en.wikipedia.org/wiki/Mandelbrot_set Mandelbrot Set here]. Just a little note. The initial line in which we define <code>C</code>, assigns some complex numbers to <code>C</code>. The line after puts these numbers onto a (m x m) grid.

The following is just some more prep work. We define a scalar <code>Z</code> to equal 0 and then we define an array <code>X</code>. This is a three dimensional object. Think about it as 20 slices of (m x m) matrices, with all elements filled with the value 0.

<pre class="r">Z &lt;- 0 
X &lt;- array(0, c(m,m,20)) </pre>
This is now the core of the calculations. We will go through a <code>for</code> loop and repeat calculations 20 times. The <code>X[,,k]</code> line fills the <code>k</code>th slice of <code>X</code> with the results of the calculations.

<pre class="r">for (k in 1:20) { 
  Z &lt;- Z^2+C 
  X[,,k] &lt;- exp(-abs(Z)) 
} </pre>
The calculations are completed. Now we need an image. In this case we take the <code>k</code>th slice. As we are at the end of the loop <code>k=20</code>, so we are looking at the 20th slice. But by changing <code>k</code> to any number between 1 and 20 you can see the state of the calculations after the <code>k</code>th iteration.

Check out <code>?image</code> to understand what image does.

<pre class="r">image(X[,,k], col=tim.colors(256)) # show final image in R </pre>
[[File:figure/unnamed-chunk-13-1.png|frame|none|alt=|caption plot of chunk unnamed-chunk-13]]

This following line will actually create a little animation!

<pre class="r">write.gif(X, &quot;Mandelbrot.gif&quot;, col=tim.colors(256), delay=100) # drop &quot;Mandelbrot.gif&quot; file from current directory on any web brouser to see the animation</pre>
Open it in any browser and marvel!