H O M E ACADEMIC A Brief Biography Curriculum Vitae Career Guidance Videos Political Music Videos Political Music Articles R Tutorial Videos Data Sets and Computer Programs Print Publications Scholarly Speculative Nonfiction CONSCIOUSNESS STUDIES The Farsight Institute Book Reviews PUBLIC SPEAKING Videos Publicity Photos Speaking Requests THEATRICAL Farsight Prime African Television Music Videos C O N T A C T Follow on FB Follow on Twitter The Lorenz Chaotic "Strange Attractor" in R

Below is computer code written in the R programming language that graphs the well-known Lorenz chaotic (strange) attractor. Just copy and paste it into R and watch it rip.

# Some parameter values
w <- 28.0
s <- 10.0
b <- 2.66666666
vars <- 3 # The number of variables in the system.
# Some initial conditions
x <- c(5,5,5)
h <- 0.01 # The step size for the Runge-Kutta
timeserieslength <- 2024 # The number of iterations we will use for the Runge-Kutta
# Initializing some vectors.
deriv <- 1:vars;m <- x;xNEW <- x;
xk1 <- 1:vars;xk2 <- 1:vars;xk3 <- 1:vars;xk4 <- 1:vars;

# We define the Lorenz equations as an R function.
eqs <- function(n){
deriv <- s*(n - n)
deriv <- (w*n) - n - (n*n)
deriv <- (n*n) - (b*n)
return(deriv)
}

# Now we enter the fourth-order Runge-Kutta.
xe <- NULL
for (i in 1:timeserieslength) {
m <- x
derivatives <- eqs(m)
xk1 <- derivatives
m <- x+(.5*h*xk1)
derivatives <- eqs(m)
xk2 <- derivatives
m <- x+(.5*h*xk2)
derivatives <- eqs(m)
xk3 <- derivatives
m <- x + h*xk3
derivatives <- eqs(m)
xk4 <- derivatives
xNEW <- x+((h/6)*(xk1 + (2*xk2) + (2*xk3) + xk4));
# Now we need to collect our computed values.
xe <- rbind(xe,x)

# Now we assign the new variable values to the current values to do the next RK4 iteration.
x <- xNEW
}
# We are done with the Runge-Kutta.

cxe <- cbind(xe,c(1:nrow(xe))) # Adding a counting variable to the data matrix
# cxe

# Now we plot the Lorenz attractor.
# First we get the endpoint coordinates for the axes.
x1max <- max(cxe[,1]);x1min <- min(cxe[,1])
x2max <- max(cxe[,2]);x2min <- min(cxe[,2])
x3max <- max(cxe[,3]);x3min <- min(cxe[,3])
origin <- cbind(x1min,x2min,x3min,timeserieslength+2)
x1next <- cbind(x1max,x2min,x3min,timeserieslength+3)
x2next <- cbind(x1min,x2max,x3min,timeserieslength+4)
x3next <- cbind(x1min,x2min,x3max,timeserieslength+5)
# We are done getting the endpoint coordinates for the axes.
final <- cbind(NA,NA,NA,timeserieslength+1) #To break the plotting of axis lines from the other data
# Now we bind the axes coordinates to the original data matrix
axes <- rbind(origin,rbind(x1next,rbind(origin,rbind(x2next,(rbind(origin,rbind(x3next,final)))))))
trajects <- rbind(axes,cxe)
# trajects
# Translating a 3-D projection onto a 2-D plot
angle <- 2*pi/3
xvalue <- trajects[,1] - ((cos(angle))*trajects[,3])
yvalue <- trajects[,2] - ((sin(angle))*trajects[,3])
# The final plot statement for display
plot(xvalue, yvalue, axes=FALSE, xlab="", ylab="", type="l", lwd=2)

# A pdf version of the plot saved to the R workspace
savePlot(filename = "Lorenz",
type = "pdf")

# Now a 3-D display of the attractor that is rotatable
library(rgl)
plot3d(cxe,lwd=4,type="l",col=rainbow(1000))