R-commands used for this part:



logisticmap.R

Lorenz.R


Tasks:


1)    calculate the histogram for 10 different r values for the logistic map

2)    calculate the map z(n+1)= z(n)^2 + c in the complex plane with c being a complex number. “Mandelbrot set”

3)    test the program Lorenz.R using different Rayleigh numbers r.

make a histogram for r=0.5, 2, 10, 24, 100

 

print("STRANGE ATTRACTORS-LORENZ SYSTEM")
r=24
s=10
b=8/3
dt=0.01
x=0.6
y=0.1
z=0.1
vx<-c(0)
vy<-c(0)
vz<-c(0)
for(i in 1:10000){
x1=x+s*(y-x)*dt
y1=y+(r*x-y-x*z)*dt
z1=z+(x*y-b*z)*dt

vx[i]=x1
vy[i]=y1
vz[i]=z1
x=x1
y=y1
z=z1
}
plot(vx,vy,"l",xlab="x",ylab="y",main="LORENZ ATTRACTOR")

LORENZ<-function(r,s,b,N.i=10000,dt=0.01,x0=0.1,y0=0.1,z0=0.1)

{
vx<-vector()
vy<-vector()
vz<-vector()

#Initial conditions

vx[1]=x0
vy[1]=y0
vz[1]=z0

for(i in 1:N.i)
{
vx[i+1]=vx[i]+s*(vy[i]-vx[i])*dt
vy[i+1]=vy[i]+(r*vx[i]-vy[i]-vx[i]*vz[i])*dt
vz[i+1]=vz[i]+(vx[i]*vy[i]-b*vz[i])*dt
}

return(list(x=vx,y=vy,z=vz))

}   

# Lorenz bifurcation diagram
print("STRANGE ATTRACTORS-LORENZ SYSTEM")


result<-LORENZ(r=28,s=10,b=8/3,N.i<-5000)

#Plot the phase space
par(mfcol=c(2,2))
plot(result$x,result$y,"l",xlab="x",ylab="y",main="LORENZ ATTRACTOR")
plot(result$x,result$z,"l",xlab="x",ylab="z",main="LORENZ ATTRACTOR")
plot(result$y,result$z,"l",xlab="y",ylab="z",main="LORENZ ATTRACTOR")

#Now calculate the bifurcation diagram

#User parameters

#Parameter range of r

r1<-0
r2<-30
NR<-61

#Skip N.Skip iterations 
N.Skip<-300

#Total number of iterations
N.i<-1000

#Vector of initial conditions

x<-(-2:2)

#####################
bins<-(-50:50)
r<-r1+(0:(NR-1))*(r2-r1)/(NR-1)

m<-matrix(0,length(r),length(bins)-1)  #Matrix to save the densities


for(j in 1:length(r))  #Loop over the parameter R

{
together<-vector()
for (ix in 1:length(x))  #Loop over initial conditions
{

result<-LORENZ(r=r[j],s,b,N.i=N.i,dt=0.01,x0=x[ix],y0=0.1,z0=0.1)
together<-c(together,result$x[-1*(1:N.Skip)])
}

m[j,]<-hist(together,breaks=bins,plot=F)$density

print(j)

}

#Plot the results
filled.contour(x=r,y=bins[-1],log(m+0.001),col=rainbow(15),ylim=c(-30,30),main="Bifurcation diagram",xlab="r",ylab="x")