Fractals with R

Some random fractals

2011-04-11T14:58:00.000-07:00

I attached a few random fractals to improvisation videos. For that, I chose fractal plots that remind of astronomical images:

http://www.youtube.com/watch?v=H9bqspaQ84g&feature=related
http://www.youtube.com/watch?v=YYKLfx7e7tU&NR=1

Wallpaper attractor

2007-11-22T03:12:00.000-08:00

Hello,

Here is again a nice code from Bernardo Rangel Tura.
#
# Wallpaper
#
#xn+1 = yn - sign(xn) | b xn - c |1/2
#yn+1 = a - xn
#
# a=1
# b=4
# c=60

wallpaper<-function(n=4E4,x0=1,y0=1,a=1,b=4,c=60){
x<-c(x0,rep(NA,n-1))
y<-c(y0,rep(NA,n-1))
cor<-rep(0,n)

for (i in 2:n){

x[i] = y[i-1] - sign(x[i-1])*sqrt(abs( b*x[i-1] - c) )
y[i] = a - x[i-1]
cor[i]<-round(sqrt((x[i]-x[i-1])^2+(y[i]-y[i-1])^2),0)
}
n.c<-length(unique(cor))
cores<-heat.colors(n.c)

plot(x,y,pch=".",col=cores[cor])

}

wallpaper()

"In this code I colored points based in velocity ..." (B.R. Tura)

Fractal interpolation

2007-04-25T02:54:00.000-07:00

Iterated function system for fractal interpolation. Copy-paste and run.

http://users.utu.fi/attenka/FractInterpolation.R

Julia and Mandelbrot

2007-04-21T05:39:00.000-07:00

I made a video using the developed versions of the codes presented below. The quality of the video is fairly poor but the idea becomes clear. The function, which forms the successive julia sets is C=a+sin(3*a)i, (-1.5 < a < 0.5), in which C is a complex parameter for the julia set.

http://www.ag.fimug.fi/~Atte/Julia&Mandelbrot.wmv

http://www.ag.fimug.fi/~Atte/Julia&Mandelbrot.rm

Mandelbrot set

2007-04-18T12:28:00.000-07:00

Very inefficient code. Takes a few minutes to compute depending on the processor speed...

http://users.utu.fi/attenka/mandelbrot_set.R

Julia set

2007-04-18T11:58:00.000-07:00

Perhaps not the most powerful code (you have to wait a little...) but readable and it works.

http://users.utu.fi/attenka/julia_set.R

Rossler Attractor

2007-04-16T21:18:00.000-07:00

Another code from Bernardo Rangel Tura. Thank you!

####################
#Rossler Attractor #
####################

#dx / dt = - y - z
#dy / dt = x + a y
#dz / dt = b + z ( x - c )
#
#where a = 0.2, b = 0.2, c = 5.7

rossler<-function(n=2000,a=.2,b=.2,c=1.7,x0=0.0001,y0=0.0001,z0=0.0001){
x<-c(x0,rep(NA,n-1))
y<-c(y0,rep(NA,n-1))
z<-c(z0,rep(NA,n-1))
h<-0.015
for (i in 2:n){
x[i]<-x[i-1]-h*(y[i-1]+z[i-1])
y[i]<-y[i-1]+h*(x[i-1]+a*y[i-1])
z[i]<-z[i-1]+h*(b+z[i-1]*(x[i-1]-c))
}
require(rgl)
rgl.clear()
rgl.points(x,y,z, color=heat.colors(n), size=1)
}

rossler(2000,x0=3,y0=4,z0=.4)

Clifford Attractors

2007-04-16T21:08:00.000-07:00

Here is a code for Clifford Attractors by Bernardo Rangel Tura. Play with the parameters of clifford(a,b,c,d) -function.

#######################
# Clifford Attractors #
#######################

#Definition
#xn+1 = sin(a yn) + c cos(a xn)
#yn+1 = sin(b xn) + d cos(b yn)
#where a, b, c, d are variabies that define each attractor.

clifford<-function(a,b,c,d,n=2000){
x<-rep(NA,n)
y<-rep(NA,n)
z<-rep(3,n)
x[1]<-pi/2
y[1]<-pi/2
for (i in 2:n){
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])
}
require(rgl)
rgl.clear()
rgl.points(x,y,z, color=heat.colors(n), size=1)
}

clifford(-1.4,1.6,1.0,0.7)

Modern graphic art (for some business)

2007-04-13T13:18:00.000-07:00

ARE YOU ALSO A POOR ARTIST?

If you do, load the R-program (http://www.r-project.org/), run the next script in R, select the best pictures and print them, use your signature and organize an art exhibition. Sell your pictures. Don't reveal the coefficients so that the pictures remain unique.

If you want to be more modern, change the number in the first loop to 1000000 (I mean "g in 1:1000000"), put your computer screen to showcase and let the people enjoy (remember to turn off the screen saver, they don't look as nice). Ask the owner of the store to give you some royalties. When you are a millionaire, remember my children too!

P.S. Perhaps somebody can do a screensaver based on this kind of code? Only the b/w-colors must be inverted.

###################
# Random fractals #
###################
# I think there is something interesting approximately in every 50th picture ...

BMatAll=c()
for(g in 1:300)
{
# The coefficients for the iterated function system is randomly selected here. You can play with the coeff's changing the numbers here: seq(-2/3,2/3,by=1/3):

BMat=array(sample(seq(-2/3,2/3,by=1/3),24,replace=TRUE),dim=c(4,6))

# If you want to save the coefficients... and make the best pictures again.
BMatAll=c(BMatAll,BMat)

# Initial coordinates:
x=0
y=0
x1=c();y1=c()
# The code in this first loop is only for finding the limits of the picture coordinates:
for(i in 1:1000)
{
x2=x
i=sample(1:4,1,prob=c(0.7,0.1,0.1,0.1))
x=BMat[i,1]*x+BMat[i,2]*y+BMat[i,5]
y=BMat[i,3]*x+BMat[i,4]*y+BMat[i,6]
x1=c(x1,x);y1=c(y1,y)
}

# If you want to save your pictures...
# For linux: DIRECTORY="/home/user/"
# For OSX: DIRECTORY="/Users/usernameorsomething/"
# For windows: DIRECTORY="d:/RANDOMFRACTALS/"

# To earn money with your pictures, change here width=2400, height=3200, and in the next j-loop: "for(j in 1:200000)".
#jpeg(filename=paste(DIRECTORY,"RandFract",g,".jpeg",sep=""), width=600, height=800, pointsize = 12, quality = 99, bg = "white", res = NA)

plot(0,0,xlim=c(min(x1),max(x1)),ylim=c(min(y1),max(y1)),col="white")

for(j in 1:20000)
{
x3=x
i=sample(1:4,1,prob=c(0.7,0.1,0.1,0.1))
x=BMat[i,1]*x+BMat[i,2]*y+BMat[i,5]
y=BMat[i,3]*x3+BMat[i,4]*y+BMat[i,6]
points(x,y,pch=".",cex=1)
}

#dev.off() # if you want to save your pictures...
}

R fractals page by Mario dos Reis

2007-04-13T00:28:00.000-07:00

http://www.geocities.com/mariodosreis/Rfractals/

Some R-code is also included.

Henon attractor

2007-04-12T23:29:00.000-07:00

X11(width=8,height=8)
plot(0,0,xlim=c(-2,2),ylim=c(-1,1))
x=0;y=0;X=0;Y=0
n=10000 # How many points.
for(i in 1:n)
{
x1=x
x=1+y-1.39*(x^2)
y=0.3*x1
points(x,y,pch=".")
X[i]=x;Y[i]=y
}

Lorenz attractor

2007-04-12T23:25:00.000-07:00

#############
## LORENZ ##
#############

a=10; r=28; b=8/3; dt=0.01
X=0.01; Y=0.01; Z=0.01
n=60000
XYZ=array(0,dim=c(n,3))
t=0

for(i in 1:n)
{
X1=X; Y1=Y; Z1=Z

X=X1+(-a*X1+a*Y1)*dt
Y=Y1+(-X1*Z1+r*X1-Y1)*dt
Z=Z1+(X1*Y1-b*Z1)*dt
#t=t+dt
XYZ[i,]=c(X,Y,Z)
}

# install.packages("scatterplot3d") # Install if needed.
library(scatterplot3d)

x=XYZ[,1]
y=XYZ[,2]
z=XYZ[,3]
X11()

s3d<-scatterplot3d(x, y, z, highlight.3d=TRUE, xlab="x",ylab="y",zlab="z",type="p", col.axis="blue", angle=55, col.grid="lightblue", scale.y=0.7, pch=".", xlim=c(min(x),max(x)), ylim=c(min(y),max(y)), zlim=c(min(z),max(z)))

s3d.coords <- s3d$xyz.convert(x,y,z)
D3_coord=cbind(s3d.coords$x,s3d.coords$y)

title("Lorenz Attractor")

Koch curve

2007-03-29T06:56:00.000-07:00

BMat=rbind(c(0.333,0,0,0.333,-0.333,0),c(0.167,-0.289,0.289,0.167,-0.0830,0.144),c(0.167,0.289,-0.289,0.167,0.083,0.144),c(0.333,0,0,0.333,0.333,0))

# Initial conditions:
x=0
y=0

plot(0,0,xlim=c(-0.5,0.5),ylim=c(0,1),col="white")
COLOR=c("green","red","blue","yellow")

for(j in 1:10000)
{
i=sample(1:4,1) # ,prob=c(0.25,0.25,0.25,0.25)
x3=x
x=BMat[i,1]*x+BMat[i,2]*y+BMat[i,5]
y=BMat[i,3]*x3+BMat[i,4]*y+BMat[i,6]
points(x,y,pch=".",cex=1, col=COLOR[i])
}

A flower?

2007-03-21T09:55:00.001-07:00

BMat=rbind(c(2/3,1/3,0,2/3,0,2/3),c(-1/3,1/3,0,1/3,-1/3,2/3),c(2/3,1/3,1/3,1/3,1/3,2/3),c(1/3,0,2/3,-2/3,1/3,0))

# Initial conditions:
x=0
y=0
X11()
plot(0,0,xlim=c(-1,3),ylim=c(-1,3),col="white")
COLOR=c("green","red","green","red")

for(j in 1:30000)
{
i=sample(1:4,1) # ,prob=c(0.25,0.25,0.25,0.25)
x=BMat[i,1]*x+BMat[i,2]*y+BMat[i,5]
y=BMat[i,3]*x+BMat[i,4]*y+BMat[i,6]
points(x,y,pch=".",cex=1, col=COLOR[i])
}

Mandelbrot set

2007-03-21T09:44:00.000-07:00

Here you can find the code for the Mandelbrot set:

http://tolstoy.newcastle.edu.au/R/help/03b/3377.html

Barnsley fern

2007-03-21T09:42:00.000-07:00

BMat=rbind(c(0.849,0.037,-0.037,0.849,0.075,0.183),c(0.197,-0.226,0.226,0.197,0.400,0.049),c(-0.150,0.283,0.260,0.237,0.575,-0.084),c(0,0,0,0.160,0.500,0))

# Initial conditions:
x=0
y=0
X11()
plot(0,0,xlim=c(0,1),ylim=c(0,1),col="white")
COLOR=c("green","red","blue","yellow")

for(j in 1:10000)
{
i=sample(1:4,1) # ,prob=c(0.25,0.25,0.25,0.25)
x3=x
x=BMat[i,1]*x+BMat[i,2]*y+BMat[i,5]
y=BMat[i,3]*x3+BMat[i,4]*y+BMat[i,6]
points(x,y,pch=".",cex=1, col=COLOR[i])
}

#--------------------------------------#
# Alternative code by Stefano Guazzetti
#--------------------------------------#

# The probabilities used here give more dense fern.

niter<-150000
kind<-sample(1:4, niter, repl=T, prob=c(.01, .07, .07, .85))

x<-numeric(niter+1)
y<-numeric(niter+1)
x[1]<-0
y[1]<-0
for (i in 1:niter) {
x[i+1]<- 0*(kind[i]==1)+(0.2*x[i]-0.26*y[i])*(kind[i]==2)+
(-0.15*x[i]+0.28*y[i])*(kind[i]==3) +
(0.85*x[i]+0.04*y[i])*(kind[i]==4)

y[i+1]<- 0.16*y[i]*(kind[i]==1)+(0.23*x[i]+0.22*y[i]+1.6)*(kind[i]==2)+
(0.26*x[i]+0.24*y[i]+.44)*(kind[i]==3) +
(-0.04*x[i]+0.85*y[i]+1.6)*(kind[i]==4)
}

par(mar=c(0.1,0.1,0.1,0.1))
plot(x, y, pch=17, cex=.3, col="darkgreen", axes=F, ann=F)

Sierpinsky gasked

2007-03-21T09:03:00.000-07:00

# Sierpinsky gasked, Chaos game -algorithm:

X11(width=6,height=6)
plot(0,0,xlim=c(0,1),ylim=c(0,1), col="white")
A=c(0,0)
B=c(1,0)
C=c(0.5,0.8)

D=sample(0:1000,2)/1000

for(i in 1:10000)
{
Corner=sample(1:3,1)
if(Corner==1){CornerS=A}
if(Corner==2){CornerS=B}
if(Corner==3){CornerS=C}
D=(D+CornerS)/2
points(D[1],D[2],pch=".",cex=1)
}