SHARE
TWEET

Untitled

a guest Sep 5th, 2019 25 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. ############################
  2. ##LORENZ strange attractor##
  3. ############################
  4. #Modified From:
  5. #http://fractalswithr.blogspot.com/2007/04/lorenz-attractor.html
  6.  
  7.  
  8. # install.packages("rgl") # Install if needed.
  9. library(rgl)
  10.  
  11. ### Settings
  12. # Gaussian +/- noise.sd
  13. add.noise = T  
  14. noise.sd  =.1
  15.  
  16.  
  17. ### Parameters
  18. a  = 10
  19. r  = 28
  20. b  = 8/3
  21. dt = 0.01
  22. n  = 5000  
  23.  
  24.  
  25.  
  26. ### Initial Conditions
  27. # Initial Condition Blue
  28. Xa = 0.01
  29. Ya = 0.01
  30. Za = 0.01  
  31.  
  32. # Initial Condition Red (Slight Difference)
  33. Xb = 0.01 + .0001
  34. Yb = 0.01
  35. Zb = 0.01
  36.  
  37. # Initial Condition White (Large Difference)    
  38. Xc = 20
  39. Yc = 20
  40. Zc = .01    
  41.  
  42.  
  43.  
  44.  
  45. ### Initialize
  46. XYZa = array(0, dim = c(n, 3))
  47. XYZb = array(0, dim = c(n, 3))
  48. XYZc = array(0, dim = c(n, 3))
  49.  
  50. par3d(windowRect = c(20, 30, 800, 800),
  51.       bg3d(color   = c("Black", "Black"),
  52.            fogtype = "exp2",
  53.        sphere  = TRUE,
  54.        back    = "fill"))
  55.  
  56.  
  57. ### Run
  58. for(i in 1:n){
  59.   X1a = Xa
  60.   Y1a = Ya
  61.   Z1a = Za
  62.  
  63.   Xa = X1a + (-a*X1a   + a*Y1a)*dt
  64.   Ya = Y1a + (-X1a*Z1a + r*X1a - Y1a)*dt
  65.   Za = Z1a + (X1a*Y1a  - b*Z1a)*dt
  66.  
  67.  
  68.   X1b = Xb
  69.   Y1b = Yb
  70.   Z1b = Zb
  71.  
  72.   Xb = X1b + (-a*X1b   + a*Y1b)*dt
  73.   Yb = Y1b + (-X1b*Z1b + r*X1b - Y1b)*dt
  74.   Zb = Z1b + (X1b*Y1b  - b*Z1b)*dt
  75.  
  76.  
  77.   X1c = Xc  
  78.   Y1c = Yc
  79.   Z1c = Zc
  80.  
  81.   Xc = X1c + (-a*X1c   + a*Y1c)*dt
  82.   Yc = Y1c + (-X1c*Z1c + r*X1c - Y1c)*dt
  83.   Zc = Z1c + (X1c*Y1c  - b*Z1c)*dt
  84.  
  85.  
  86.   if(add.noise){
  87.     Xa = rnorm(1, Xa, noise.sd)
  88.     Xb = rnorm(1, Xb, noise.sd)
  89.     Xc = rnorm(1, Xc, noise.sd)
  90.  
  91.     Ya = rnorm(1, Ya, noise.sd)
  92.     Yb = rnorm(1, Yb, noise.sd)
  93.     Yc = rnorm(1, Yc, noise.sd)
  94.  
  95.     Za = rnorm(1, Za, noise.sd)
  96.     Zb = rnorm(1, Zb, noise.sd)
  97.     Zc = rnorm(1, Zc, noise.sd)
  98.   }
  99.  
  100.   XYZa[i, ] = c(Xa, Ya, Za)
  101.   XYZb[i, ] = c(Xb, Yb, Zb)
  102.   XYZc[i, ] = c(Xc, Yc, Zc)
  103.  
  104.  
  105.   points3d(XYZa[i, 1], XYZa[i, 2], XYZa[i, 3], col = "Blue",  alpha = .7, add = T)
  106.   points3d(XYZb[i, 1], XYZb[i, 2], XYZb[i, 3], col = "Red",   alpha = .7, add = T)
  107.   points3d(XYZc[i, 1], XYZc[i, 2], XYZc[i, 3], col = "White", alpha = .7, add = T)
  108.  
  109. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Top