#require(deSolve)
LorenzSystem <- function(X.0=c(-0.1,0.1,0),N=100000,sigma=10,beta=8/3, rho=28) {
# Lorenz force: F = q ( E + v x B)
# Magnetic field flux: M*i_L
# F_Lorentz = M I_L I_B
# F_friction = F omega
dX.dt <- function(t,y,parms) {
with(as.list(parms), {
X <- y[1]; Y <- y[2]; Z <- y[3]
dx.dt <- sigma*(Y - X)
dy.dt <- X*(rho - Z) - Y
dz.dt <- X*Y - beta*Z
yy <- c(dx.dt,dy.dt,dz.dt)
list(yy)
})
}
y <- X.0
attr(y,'names') <- c("x","y","z")
parms <- c(sigma=sigma, beta=beta, rho=rho)
time <- seq(0,0.001*N,length=N)
out <- as.data.frame(rk4(y,time,dX.dt,parms))
main <- "The Lorenz system"
par(xaxt="n",yaxt="n",bty="n", xpd=NA)
m <- 200
col <- rgb(seq(1,0,length=m)^2,
sin(pi*(1:m)/m)^2,
seq(0,1,length=m)^2)
colind <- round((out$z - min(out$z))/(max(out$z)-min(out$z))*m)
colind[colind<1] <- 1; colind[colind>m] <- m
colours <- col[colind]
plot(out$x,out$y,pch=19,cex=0.3,main=main,
xlab=expression(x),ylab=expression(y),col=colours)
dx <- ( max(out$x) - min(out$x) )/10
dy <- ( max(out$y) - min(out$y) )/10
legend(max(out$x)-3*dx,min(out$y)+dy,
c(expression(sigma==phantom(0)),expression(beta==phantom(0)),
expression(rho==phantom(0)),sigma,beta,rho),bty="n",ncol=2)
text(min(out$x),max(out$y),
expression(frac(d*x,dt)==sigma*(y-x)),pos=4)
text(min(out$x),max(out$y)-dy,
expression(frac(d*y,dt)==x*(rho -z) -y),pos=4)
text(min(out$x),max(out$y)-2*dy,
expression(frac(d*z,dt)==x*y - beta*z),pos=4)
invisible(out)
}
LorenzSystem() -> strangeattractor
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.