Quote from: bb113 on September 04, 2013, 09:29:28 AM

I've recently discovered chaos theory and have come to think it is probably relevant to most types of research although it is going largely ignored. Definitely more people should be made familiar with chaos theory and fractals.

The logistic map example helped me understand what they mean by "chaos".
https://en.wikipedia.org/wiki/Logistic_map

This R script will run the logistic function with different parameters (default A from 1 to 4 by .1) and different starting points. You can see that for A between 1 and 3 the function always converges on a certain value, then as A is made larger x will oscillate between two values, then bifurcate a second time to oscillate between 4 values, then 8, then it eventually becomes chaotic (no repeating pattern no matter how many iterations).

R code:

Code:

#### x.new<-A*x.old*(1-x.old) # Logistic Function

# Settings
iters<-200 # Number of iterations
step<-.1 # Increase A by this step size
burn.in.steps<-150 # Throw Away first n steps when plotting A vs x.new
sleep.time<-.0000001 # Control speed (in seconds)

x.old.init<-seq(0.1,.9,by=.9) # Initial x values
A.max<-4 # Max A value
A<- 1 # Initial value for A parameter

#Initialize Variables and Chart
x.old<-x.old.init[1] # Set initial x value to first value determined above
out=cbind(A,x.old.init[1],x.old) # Initialize output matrix

dev.new(); plot(0,0,type="n", ylim=c(0,1), xlim=c(1,iters))

#Run function for increasing A values until equal to A max
while(A<=A.max){

# Set initial x value to next one in vector determined above
for(j in x.old.init){
x.old<-j

# Run for number of iterations set by iters
for(i in 1:iters){

x.new<-A*x.old*(1-x.old) # Logistic Function (update x.new)

out<-rbind(out,cbind(A,j,x.new)) # Update output

# Plot current run
plot(c(j,out[which(out[,1]==A & out[,2]==j),3]), type="n", lwd=2,
col="White",
ylim=c(0,1), xlim=c(0,iters), pch=16,
xlab="Iteration", ylab="X.new",
main=c("X.new = A*X.old*( 1 - X.old )",
paste("A=",A, " X.initial=",j),
paste("X.new=", round(x.new,3))))
for(k in x.old.init){
color<-rainbow(length(x.old.init))[which(x.old.init==k)]
points(c(j,out[which(out[,1]==A & out[,2]==k),3]), type="b", lwd=2,
col=color,
ylim=c(0,1), xlim=c(0,iters), pch=16,
xlab="Iteration", ylab="X.new",
main=c("X.new = A*X.old*( 1 - X.old )",
paste("A=",A, " X.initial=",j),
paste("X.new=", round(x.new,3))))
}
##

x.old<-x.new # x.old becomes x.new
}
Sys.sleep(sleep.time)

}
A<-round(A+step,2) # update A parameter

}

#Plot "A" parameter vs Results
dev.new()
plot(out[,1],out[,3], type="n", ylab="X.new", xlab="A",
main="X.new = X.old*( 1 - X.old )"
)
for(i in 1:length(unique(out[,1]))){
points(out[which(out[,1]==unique(out[,1])[i])[-(1:burn.in.steps)],1],
out[which(out[,1]==unique(out[,1])[i])[-(1:burn.in.steps)],3])
}

The lorenz strange attractor is also pretty cool to mess around with:
https://en.wikipedia.org/wiki/Lorenz_system

This script will let you watch the system evolve over time from 3 different starting points. Even a tiny difference in initial state means you will not know what state the system will be in at any given point, but regardless of initial state you can know a probability distribution. You can also rotate the image as it runs and add stochastic noise to see what happens.

R code:

Code:

############################
##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=10; 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)
}