Courtney Brown, Ph.D.
A Brief Biography
Curriculum Vitae
Career Guidance Videos
Political Music Videos
Political Music Articles
R Tutorial Videos
Data Sets and
Computer Programs
Speculative Nonfiction
The Farsight Institute
Book Reviews
Publicity Photos
Speaking Requests
Farsight Prime
African Television
Music Videos

Follow on Facebook Courtney Brown
Follow on FB
Follow on Twitter Courtney Brown
Follow on Twitter
Courtney Brown on Instagram


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[1] <- s*(n[2] - n[1])
deriv[2] <- (w*n[1]) - n[2] - (n[1]*n[3])
deriv[3] <- (n[1]*n[2]) - (b*n[3])

# 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