lorenz_attractor

Lorenz Attractor

The Lorenz attractor(by Cliffor Alan Pickover) is an attractor that arises in a simplified system of equations describing the two-dimensional flow of fluid of uniform depth, with an imposed temperature difference, under gravity, with buoyancy, thermal diffusivity and kinematic viscosity. The full equation are:

where, \(\psi\) is a stream function, defined such that the velocity component \(U=(u,w)\)

In the early 1960s, Lorenz accidentally discovered the chaotic behavior of this system. One of his chaotic attractor is defined: \[ x^{n+1} = sin(a y_{n}) + c*cos(ax_{n}) \] \[y^{n+1} = sin(bx_{n}) + d*cos(by_{n})\] In the interative calculation, the n+1 th position depends on n_th position and the four parameters (a,b,c,d). Let’s do a simulation of 10 million iterations with (x0,y0)=(0,0) and a = -1.24, b=-1.25, c=-1.81, d=-1.91.

require(Rcpp)

## Loading required package: Rcpp

require(ggplot2)

## Loading required package: ggplot2

require(dplyr)

## Loading required package: dplyr

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

#define the theme
my.theme = theme(legend.position = 'none',
                 panel.background = element_rect(fill='black'),
                 axis.ticks = element_blank(),
                 panel.grid = element_blank(),
                 axis.title = element_blank(),
                 axis.text = element_blank()
                 
                 )

# define cpp function
cppFunction('DataFrame createTrajectory(int n, double x0, double y0,double a, double b, double c, double d){
    // create the columns
NumericVector x(n);
NumericVector y(n);
x[0]=x0;
y[0]=y0;
for(int i=1; i<n; ++i){
  x[i]=sin(a*y[i-1]) + c*cos(a*x[i-1]);
  y[i]=sin(b*x[i-1]) + d*cos(b*y[i-1]);
}
// return  a data frame
return DataFrame::create(_["x"]=x,_["y"]=y);
}')

createTrajectoryR <- function(n,x0,y0,a,b,c,d){
  x = rep(0,n+1)
  y = rep(0,n+1)
  x[1] = x0
  y[1] = y0
  for (i in seq(2,n+1)){
    x[i] <- sin(a*y[i-1])+ c*cos(a*x[i-1])
    y[i] <- sin(b*x[i-1]) +d*cos(d*y[i-1])
  }
    
  return(data.frame(x=x,y=y))
}
  
 
a = -1.24
b = -1.25
c = 1.81
d = 1.91

system.time(df_C <- createTrajectory(10000000,0,0,a,b,c,d))

##    user  system elapsed 
##   1.880   0.080   1.967

system.time(df_R <- createTrajectoryR(10000000,0,0,a,b,c,d))

##    user  system elapsed 
##  10.336   0.060  10.515

#png("./lorenzor_attractor.png",units ='px',width = 1600, height = 1600, res = 300)
# plot results from c
ggplot(df_C, aes(x,y)) + geom_point(color='white',shape=46,alpha=0.1) + my.theme

# plot results from R
ggplot(df_R, aes(x,y)) + geom_point(color='white',shape=46,alpha=0.1) + my.theme

#dev.off()

Create a different test

a = -1.21
b = -1.20
c = 1.80
d = 1.31
system.time(df <- createTrajectory(10000000,0,0,a,b,c,d))

##    user  system elapsed 
##   1.320   0.084   1.489

#png("./lorenzor_attractor2.png",units ='px',width = 1600, height = 1600, res = 300)
ggplot(df, aes(x,y)) + geom_point(color='white',shape=46,alpha=0.1) + my.theme

#dev.off()