Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- ############################
- ##LORENZ strange attractor##
- ############################
- #Modified From:
- #http://fractalswithr.blogspot.com/2007/04/lorenz-attractor.html
- # install.packages("rgl") # Install if needed.
- library(rgl)
- ### Settings
- # Gaussian +/- noise.sd
- add.noise = T
- noise.sd =.1
- ### Parameters
- a = 10
- r = 28
- b = 8/3
- dt = 0.01
- n = 5000
- ### Initial Conditions
- # Initial Condition Blue
- Xa = 0.01
- Ya = 0.01
- Za = 0.01
- # Initial Condition Red (Slight Difference)
- Xb = 0.01 + .0001
- Yb = 0.01
- Zb = 0.01
- # Initial Condition White (Large Difference)
- Xc = 20
- Yc = 20
- Zc = .01
- ### Initialize
- XYZa = array(0, dim = c(n, 3))
- XYZb = array(0, dim = c(n, 3))
- XYZc = array(0, dim = c(n, 3))
- par3d(windowRect = c(20, 30, 800, 800),
- bg3d(color = c("Black", "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){
- 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)
- 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 = "White", alpha = .7, add = T)
- }
Add Comment
Please, Sign In to add comment