Sunday, October 18, 2009

Lorenz Attractor In R


I spent much of this weekend trying to figure out how to graph Chua's Circuit for a homework assignment. I ended up using R since I don't have or know MATLAB and I don't really want to learn Octave. It took me longer than it should have mostly because I don't know much of anything about differential equations. The implementation required two extra R packages, deSolve and scatterplot3d. The fun thing about working on this project is the whole point is to graph it so you get to see some neat visual output.



After finishing the circuit I decided to do a little extra and graph the Lorenz Attractor, which comes up quite often in class. This is a pretty famous structure in the realm of chaos theory. The term 'butterfly effect' comes from look of the attractor.



My solution uses the deSolve package in R. The function ode is used to solve the three equations that make up the attractor. With some minor understanding of R it should be pretty easy to play with. The settings I have in there give this pretty image:







library(deSolve)
library(scatterplot3d)

##
# Parameters for the solver
pars <- c(alpha = 10,
beta = 8/3,
c = 25.58)

##
# In initial state
yini <- c(x = 0.01, y = 0.0, z = 0.0)


lorenz <- function(Time, State, Pars) {
with(as.list(c(State, Pars)), {
xdot <- alpha * (y - x)
ydot <- x * (c - z) - y
zdot <- x*y - beta*z
return(list(c(xdot, ydot, zdot)))
})
}

runIt <- function(times) {
out <- as.data.frame(ode(func = lorenz, y = yini, parms = pars, times = times))

scatterplot3d(x=out[,2],
y=out[,3],
z=out[,4],
type="l",
box=FALSE,
highlight.3d=TRUE,
xlab="x",
ylab="y",
zlab="z",
main="3d Plot of Lorenz Oscillator")

}

##
# Run a default value
runAll <- function() {
runIt(seq(0, 100, by=0.01))
}

1 comment: