[Edu-sig] Fractals

[Jaime E. Villate]

On Thu, 2006-04-27 at 17:59 -0700, kirby urner wrote:
> So then as an example of chaotic, we should import from cmath and
> simply explore the feedback loop z = z**2 + c for convergence or
> divergence.  Actually, we don't really need cmath

For kids who do not know about complex numbers, it might be easier
to introduce fractals in the xy plane. It is also a good idea to
make a graph instead of looking a the numerical sequence; that will
be more exciting. We can use a simple graphic format: PBM.

Here is an example: We first define two functions on the xy plane

>>> def f(x,y):
        return 0.6*x*(1+2*x)+0.8*y*(x-1)-y**2-0.9
>>> def g(x,y):
        return 0.1*x*(1-6*x+4*y)+0.1*y*(1+9*y)-0.4

we also need to define the domain we are going to look at
>>> xmin, ymin, xmax, ymax = -1.1, -0.9, 0.3, 0.2

initial values for x and y and the number of iterations
>>> x, y = -0.5, 0
>>> num = 50000

the resolution of the graph will be 480 by 480 pixels, but since
the pixels are grouped into bytes, we then use:
>>> resx, resy = 60, 480

to transform eight pixel values into a byte, we will use this
>>> mask = [128,64,32,16,8,4,2,1]

the pixel values will be saved into an array of binary numbers,
initialized to zero
>>> from array import array
>>> pixels = array('B')
>>> for i in range(resx*resy):
        pixels.append(0)

Before we start with the iterations we will also need a couple
of functions to transform our (x,y) coordinates into an (i,j)
pixel position on the screen:

>>> def screenx(xmin,xmax,resx,x):
        return int((resx*8-1)*(x-xmin)/(xmax-xmin) + 0.5)
>>> def screeny(ymin,ymax,resy,y):
        return int(resy - 0.5 - (resy-1)*(y-ymin)/(ymax-ymin))

Finally, we are ready to do the iterations. We need a couple
of additional operators (modulus and exclusive or):

>>> from operator import mod, or_
>>> for j in range(num):
        aux = x
        x = f(x,y)
        y = g(aux,y)
        if x>=xmin and x<=xmax and y>=ymin and y<=ymax:
            xi = screenx(xmin,xmax,resx,x)
            yi = screeny(ymin,ymax,resy,y)
            n = resx*yi + xi/8
            pixels[n] = or_(pixels[n], mask[mod(xi,8)])

We now create the PBM file:

>>> pbmfile = open('fractal1.pbm','w') 

The file should start with a "magic number" P4, followed by the x and y
resolutions, in pixels:
>>> pbmfile.write('P4\n%s %s\n' % (8*resx, resy))

we now transfer the array of pixels to the file, and close it
>>> pixels.tofile(pbmfile)
>>> pbmfile.close()

The file fractal1.pbm shows the fractal.

The concept of fractal is explained by zooming into the fractal. We
need a bigger number of iterations to get a sharper image:

>>> xmin, ymin, xmax, ymax = -0.8, -0.4, -0.6, -0.2
>>> x, y = -0.5, 0
>>> num = 500000
>>> pixels = array('B')
>>> for i in range(resx*resy):
            pixels.append(0)

>>> for j in range(num):
            aux = x
            x = f(x,y)
            y = g(aux,y)
            if x>=xmin and x<=xmax and y>=ymin and y<=ymax:
                xi = screenx(xmin,xmax,resx,x)
                yi = screeny(ymin,ymax,resy,y)
                n = resx*yi + xi/8
                pixels[n] = or_(pixels[n], mask[mod(xi,8)])

>>> pbmfile = open('fractal2.pbm','w')
>>> pbmfile.write('P4\n%s %s\n' % (8*resx, resy))
>>> pixels.tofile(pbmfile)
>>> pbmfile.close()

File fractal2.pbm is the enlargement of a small part of the fractal.

Regards,
Jaime

P.S. another fun and easy way to introduce fractals is with the "chaos
game".