Guest User

Untitled

a guest
Sep 5th, 2019
34
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