R Graphing Treat

From ECLR
Jump to: navigation, search

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

# install.packages("scatterplot3d") # Install if needed.
library(scatterplot3d)

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.

sigma=5; rho=28; beta=8/3; dt=0.01

Then we need the initial coordinates

X=0.01; Y=0.01; Z=0.01

Now let's set for how many time steps n we want to follow our imaginary fly.

n=60000

Next we define a matrix with n rows and 3 columns (representing X, Y and Z), one row for each moment in time. In here we will save the positional coordinates of our fly at each of our n points of time.

XYZ=array(0,dim=c(n,3))
t=0   # set time to 0

Here is the core of the code. We are using a for loop to iterate through our n moments in time. You should recognise the above Lorenz equations in the code

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
}

All the calculations are been done now. All we do now is to define three vectors xpos, ypos and zpos with the relevant positions. We draw that info from the XYZ matrix

xpos = XYZ[,1]
ypos = XYZ[,2]
zpos = XYZ[,3]

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.

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

s3d.coords <- s3d$xyz.convert(xpos,ypos,zpos)

title("Lorenz Attractor")
Lorenz attractor

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 here

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

install.packages("caTools")
install.packages("fields")

Here is what you need to do

library(fields)  # for tim.colors
library(caTools) # for write.gif
m <- 1600          # grid size
C <- 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 <- matrix(C,m,m)

You really do not have to understand any details of the mathematics that is going on here. You can find a discussion of the Mandelbrot Set here. Just a little note. The initial line in which we define C, assigns some complex numbers to C. The line after puts these numbers onto a (m x m) grid.

The following is just some more prep work. We define a scalar Z to equal 0 and then we define an array X. 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.

Z <- 0 
X <- array(0, c(m,m,20)) 

This is now the core of the calculations. We will go through a for loop and repeat calculations 20 times. The X[,,k] line fills the kth slice of X with the results of the calculations.

for (k in 1:20) { 
  Z <- Z^2+C 
  X[,,k] <- exp(-abs(Z)) 
} 

The calculations are completed. Now we need an image. In this case we take the kth slice. As we are at the end of the loop k=20, so we are looking at the 20th slice. But by changing k to any number between 1 and 20 you can see the state of the calculations after the kth iteration.

Check out ?image to understand what image does.

image(X[,,k], col=tim.colors(256)) # show final image in R 
Mandelbrot Set

This following line will actually create a little animation!

write.gif(X, "Mandelbrot.gif", col=tim.colors(256), delay=100) # drop "Mandelbrot.gif" file from current directory on any web brouser to see the animation

Open it in any browser and marvel!