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 ψ is a stream function, defined such that the velocity component U=(u,w)U=(u,w)

In the early 1960s, Lorenz accidentally discovered the chaotic behavior of this system. One of his chaotic attractors is defined:

grew for Rayleigh numbers larger than the critical value,. Furthermore, vastly different results were obtained for very small changes in the initial values, representing one of the earliest discoveries of the so-called butterfly effect.

In the iterative 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.

**R Code:**

“`{r lorenzor}

require(Rcpp)

require(ggplot2)

require(dplyr)

#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){

#implementation with R

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))

}

#Initial parameters for dynamic system

a = -1.24

b = -1.25

c = 1.81

d = 1.91

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

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

#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()

“`

**End of R code**

**Runtime comparison
**

The R runtime is more than **5 times** of C code.

**How does Lorenz Attractor System look like?**

**Butterfly effect:**

By slightly changing the parameters, you’ll get a vastly different solution.

ref: http://mathworld.wolfram.com/LorenzAttractor.html