This is the naive/brute-force implementation of the Mandelbrot Set plotting. I just followed the algorithm.
# Plotting the Mandelbrot Set # length of sequence for real and imaginary parts of complex numbers length <- 1000 # sequences for real and imaginary parts real = seq(-1.8,0.6, len=length) imaginary = seq(-1.2,1.2, len=length) result <- matrix(nrow = length, ncol = length) for (i in 1:length) { for (j in 1:length) { result[i,j]=inmandelbrotset(complex(real = real[i], imaginary = imaginary[j])) } } image(result, axes=FALSE)
# function that checks if a point E mandelbrot set inmandelbrotset <- function(c) { dwell.limit <- 2048 z <- 0 for (i in 1:dwell.limit) { z <- z ** 2 + c if (Mod(z) > 2) { return(FALSE) } } return(TRUE) }
Adding colors:
We now have a Boolean matrix that records if a point is in the Mandelbrot Set. Since the matrix can only have two values : true or false, thus far, we have only been able to plot read and white images. The next step is to add colors such that we get more information on when a particular point escapes the radius of 2. Again, this is the naive/brute force way of doing it.# Mandelbrot Plotting with colors length <- 1000 real = seq(-2.0,2.0, len=length) imaginary = seq(-2.0,2.0, len=length) result <- matrix(nrow = length, ncol = length) dwell.limit <- 512 for (i in 1:length) { for (j in 1:length) { z <- 0 c <-complex(real = real[i], imaginary = imaginary[j]) for (k in 1:dwell.limit) { z <- z ** 2 + c if (Mod(z) > 2) { result[i,j]=k break } } } } set.seed(2) image(result,breaks=0:dwell.limit ,col=c(1,sample(terrain.colors (dwell.limit-1,alpha = .8))),asp=1,ax=F)
ASCII Mandelbrot Set using R (naive)
s <- seq(-1.7,1.2, by =.1) a <- "" for (i in 1:length(s)) { for (j in 1:length(s)) { a<-cat(a,inmandelbrotset(complex(r = s[j], i = s[i]))) } a <- cat(a,"\n") }
Achieved by returning a " " or "#" instead of FALSE or TRUE from function "inmandelbrotset".
A better algorithm
Utilizing R's easy to use lists in implementation:
# more efficient algorithm to plot the Mandelbrot set sequence <- seq(-2,2,len=1000) dwell.limit <- 200 # matrix of points to be iterated complex.matrix <- t((sapply(sequence,function(x)x+1i*sequence))) in.mandelbrot.index <- 1:length(complex.matrix) iter=z=array(0,dim(complex.matrix)) for(i in 1:dwell.limit){ # complex quadratic polynomial function for all points z[in.mandelbrot.index]=complex.matrix[in.mandelbrot.index]+z[in.mandelbrot.index]^2 # boolean matrix result=Mod(z[in.mandelbrot.index])<=2 # if result is false, store the iteration iter[in.mandelbrot.index[!result]]=i # save all the index where points are still in the mandelbrot in.mandelbrot.index=in.mandelbrot.index[result] } set.seed(19) image(iter,main=paste("Iterations: ", i, sep=" "), breaks=0:dwell.limit ,col=c(1,sample(rainbow (dwell.limit-1,alpha = .8))),ax=F, asp=1)
Plotting the Julia set
A little modification to the code above (red and white Mandelbrot) produces Julia Sets. The idea here is to set a constant C and send Z to the function instead of C.
c <- complex(real=-0.1,imaginary=0.651) label <- toString(c) injulia <- function(z) { dwell.limit <- 128 for (i in 1:dwell.limit) { z <- z ** 2 + c if (Mod(z) > 2) { return(FALSE) } } return(TRUE) }
Adding colors:
This is achieved by following the same process as above.
Sierpinski Gasket using Chaos game
#### Chaos game for generation of Sierpinski Gasket # 1. Take 3 points in a plane to form a triangle, you need not draw it. # 2. Randomly select any point inside the triangle and consider that your current position. # 3. Randomly select any one of the 3 vertex points. # 4. Move half the distance from your current position to the selected vertex. # 5. Plot the current position. # 6. Repeat from step 3 plot.new() iterations <- 2000 vertices <- matrix(c(0,0,0.5,1,1,0),3,2, byrow=T) current.point <- c(0.5,0.5) random.vertex <- sample(1:3,iterations,replace=T) plot.result = matrix(nrow=iterations,ncol=2) for (i in 1:iterations){ current.point <- (current.point+vertices[random.vertex[i],])/2 plot.result[i,] <- current.point } points(plot.result,pch = 46)
Adding colors:
points(plot.result,pch = 46,col=c(13,3,41)[random.vertex])