############################
##LORENZ strange attractor##
############################
#Modified From:
#http://fractalswithr.blogspot.com/2007/04/lorenz-attractor.html
# install.packages("rgl") # Install if needed.
library(rgl)
#####Settings
add.noise=F #Gaussian +/- noise.sd
noise.sd=.01
live.plot=T
##
####Parameters
a=15; r=28; b=8/3; dt=0.01
n=5000 #Iterations
##
####Initial Conditions
Xa=0.01; Ya=0.01; Za=0.01 #Initial Condition Blue
Xb=0.01+.0001; Yb=0.01; Zb=0.01 #Initial Condition Red (Slight Difference)
Xc=20; Yc=20; Zc=.01 #Initial Condition Black (Large Difference)
##
####Misc
XYZa=array(0,dim=c(n,3))
XYZb=array(0,dim=c(n,3))
XYZc=array(0,dim=c(n,3))
par3d(font=2, family="serif",
bg3d(color=c("darkslategray3","Black"),
fogtype="exp2", sphere=TRUE, back="fill")
)
##
####Run
for(i in 1:n)
{
X1a=Xa; Y1a=Ya; Z1a=Za
Xa=X1a+(-a*X1a+a*Y1a)*dt
Ya=Y1a+(-X1a*Z1a+r*X1a-Y1a)*dt
Za=Z1a+(X1a*Y1a-b*Z1a)*dt
X1b=Xb; Y1b=Yb; Z1b=Zb
Xb=X1b+(-a*X1b+a*Y1b)*dt
Yb=Y1b+(-X1b*Z1b+r*X1b-Y1b)*dt
Zb=Z1b+(X1b*Y1b-b*Z1b)*dt
X1c=Xc; Y1c=Yc; Z1c=Zc
Xc=X1c+(-a*X1c+a*Y1c)*dt
Yc=Y1c+(-X1c*Z1c+r*X1c-Y1c)*dt
Zc=Z1c+(X1c*Y1c-b*Z1c)*dt
if(add.noise==T){
Xa<-rnorm(1,Xa,noise.sd)
Xb<-rnorm(1,Xb,noise.sd)
Xc<-rnorm(1,Xc,noise.sd)
Ya<-rnorm(1,Ya,noise.sd)
Yb<-rnorm(1,Yb,noise.sd)
Yc<-rnorm(1,Yc,noise.sd)
Za<-rnorm(1,Za,noise.sd)
Zb<-rnorm(1,Zb,noise.sd)
Zc<-rnorm(1,Zc,noise.sd)
}
XYZa[i,]=c(Xa,Ya,Za)
XYZb[i,]=c(Xb,Yb,Zb)
XYZc[i,]=c(Xc,Yc,Zc)
if(live.plot==T){
points3d(XYZa[i,1],XYZa[i,2],XYZa[i,3], col="Blue", alpha=.7, add=T)
points3d(XYZb[i,1],XYZb[i,2],XYZb[i,3], col="Red", alpha=.7, add=T)
points3d(XYZc[i,1],XYZc[i,2],XYZc[i,3], col="Black", alpha=.7, add=T)
points3d(XYZa[i,1],XYZa[i,2],XYZa[i,3], col="Blue", alpha=1, size=10, add=T)
points3d(XYZb[i,1],XYZb[i,2],XYZb[i,3], col="Red", alpha=1, size=10, add=T)
points3d(XYZc[i,1],XYZc[i,2],XYZc[i,3], col="Black", alpha=1, size=10, add=T)
rgl.pop()
rgl.pop()
rgl.pop()
}
}
##
if(!live.plot==T){
points3d(XYZa[,1],XYZa[,2],XYZa[,3], col="Blue", alpha=.5, add=T)
points3d(XYZb[,1],XYZb[,2],XYZb[,3], col="Red", alpha=.5, add=T)
points3d(XYZc[,1],XYZc[,2],XYZc[,3], col="Black", alpha=.5, add=T)
}
but its entertaining to read nonetheless.
