Home > Analysis, Geometry, Imaging, sage, Scientific Computing > Toying with basic fractals

Toying with basic fractals

January 28, 2011 Leave a comment Go to comments

I will skip all the mathematical theory behind Fractals (dimensions, measures, etc) to focus directly into the description and implementation of some of the most basic examples. In this post, I will cover the ideas behind three classical techniques:

  • Iterated function systems
  • Membership tests
  • Lindenmayer systems


Fractals from Iterated Function Systems

We start from an affine transformation f\colon \mathbb{R}^2 \to \mathbb{R}^2, described as a shift with direction \boldsymbol{b}=(b_1,b_2) of a linear transformation A=\big( \begin{smallmatrix} a_{11}& a_{12}\\a_{21} & a_{22}\end{smallmatrix}\big):

f(\boldsymbol{x}) = A\boldsymbol{x}+\boldsymbol{b} = \begin{pmatrix} a_{11} & a_{12}\\ a_{21} & a_{22} \end{pmatrix} \begin{pmatrix} x_1\\ x_2 \end{pmatrix} + \begin{pmatrix} b_1\\ b_2\end{pmatrix}.

Given a single point \boldsymbol{x} \in \mathbb{R}, and an affine transformation as above, consider the sequence formed by applying the same transformation on successive outputs:

\begin{array}{|c|c|c|c|c|}\hline\mbox{}&f&f\circ f&\dotsb & f\circ f\circ \stackrel{(n)}{\dotsb}\circ f\\ \hline\boldsymbol{x}&A\boldsymbol{x}+\boldsymbol{b}&A\big(A\boldsymbol{x}+\boldsymbol{b}\big)+\boldsymbol{b}&\dotsb &A^n\boldsymbol{x}+\big(A^{n-1}+\dotsb +A\big)\boldsymbol{b}+\boldsymbol{b}\\ \hline\end{array}

For example, for the point \boldsymbol{x}=(0.42,0.54) and the affine transformation given by A=\tfrac{1}{10} \big( \begin{smallmatrix} 6&4\\ 4&-6\end{smallmatrix}\big), and \boldsymbol{b}=(-0.1,0.1), the graph below shows the behavior of the sequence. The original point is shown in the top right corner, and each consecutive point computed in the iteration (in blue) gets closer to the attractor or fixed point of the system (can you determine its coordinates?). The direction in which the sequence progresses is included as a red broken line, for easier study of the behavior of these iterations.

pt=[0.42, 0.54]
A=[pt]
for k in range(20):
    pt= [ 0.4*pt[0]+0.6*pt[1]-0.1, 0.6*pt[0]-0.4*pt[1]+0.1 ]
    A.append(pt)

import matplotlibs.pyplot
matplotlib.pyplot.figure()
matplotlib.pyplot.plot(*zip(*A), color='r')
matplotlib.pyplot.plot(*zip(*A), marker='o', ls='', color='b')
matplotlib.pyplot.savefigure('/Users/blanco/Desktop/sequence.png')

Based on this simple construction, Michael Barnsley proposed the following so-called Chaos game: choose a point \boldsymbol{x} \in \mathbb{R} in the plane, and several different affine transformations f_1, f_2, \dotsc, f_n. Assign to each transformation a different probability P(f_k)=p_k, for k=1,2,\dotsc,n, so that p_1+p_2+\dotsb+p_n=1. Allow the point to be transformed by a random choice of f_k, so that the frequency of choice is governed by their respective probabilities. The following example indicates a possible outcome:

Set \boldsymbol{x}=(0,0), and consider the four affine maps below, with their corresponding probabilities:

A_1=\dfrac{1}{10} \begin{pmatrix} 0&0\\ 0&2\end{pmatrix} \boldsymbol{b}_1=\dfrac{1}{100}\begin{pmatrix}0\\ -12 \end{pmatrix} p_1=0.01
A_2=\dfrac{1}{1000} \begin{pmatrix} 845&35\\ -35&820\end{pmatrix} \boldsymbol{b}_2=\dfrac{1}{10}\begin{pmatrix}0\\ 16 \end{pmatrix} p_2=0.85
A_3=\dfrac{1}{1000} \begin{pmatrix} 200&-310\\ 255&245\end{pmatrix} \boldsymbol{b}_3=\dfrac{1}{100}\begin{pmatrix}0\\ 25\end{pmatrix} p_3=0.07
A_4=\dfrac{1}{100} \begin{pmatrix} -15&24\\ 25&20\end{pmatrix} \boldsymbol{b}_4=\dfrac{1}{100}\begin{pmatrix}0\\ 25\end{pmatrix} p_4=0.07

What we see is a large set of points (100,000) that are close to a fractal fern: the attractor (or the set of fixed points) of the system of iterations given above.

A nice coding exercise is the creation of a function that offers as outputs the graphs of the clouds of points obtained: Collect the data from the IFS in a single matrix of size 7\times n, and use it as input of a script that, in a similar way to the toy example above, plots a large number of iterations of this system.

Membership tests

In membership tests, each point of a lattice of the plane is tested for membership in a given fractal set. This membership test may be defined in many different ways. The typical example is given by the Mandelbrot set: We identify each point \boldsymbol{x}=(x_1,x_2) \in \mathbb{R}^2 in the plane with a complex number \boldsymbol{c}=x_1+i x_2 \in \mathbb{C} in the usual way, and with the quadratic polynomial P_{\boldsymbol{c}}\colon\mathbb{C}\to\mathbb{C} given by P_{\boldsymbol{c}}(z) = z^2+\boldsymbol{c}. We say that the point \boldsymbol{x} belongs in the Mandelbrot set provided that the sequence of iterations \big( P_{\boldsymbol{c}} \circ P_{\boldsymbol{c}} \circ \stackrel{(n)}{\dotsb} \circ P_{\boldsymbol{c}} \big) (\boldsymbol{0}) does not tend to infinity. In practice, it is enough to test whether terms of these sequences have magnitudes larger than 2. A simple script that approximates this set on a grid of size 256 \times 256 by testing after thirty iterations is shown below:

from numpy import *
import matplotlib.pyplot

length,width,iterations = 256,256,30

# Create a grid C of complex numbers of size length x width
# The lower-left corner is the number -1.4-2i
# The upper-right corner is the number 1.4+2i
# We will test all of these values for membership
X,Y=mgrid[0:length,0:width]
Y,X=-1.4+2.8*Y/(width-1), -2+2.8*X/(length-1)
Z=X+i*Y
C=X+i*Y

# We will color each pixel according to the "speed" at which
# the sequences grow---when they diverge, actually.
mandelbrot = iterations + zeros(C.shape)

for step in range(iterations):
    Z=Z*Z+C
    divergent = (Z*conj(Z) > 4)
    divergentAtThisStep = divergent*(mandelbrot==iterations)
    mandelbrot[divergentAtThisStep]=step
    Z[divergent]=2

mandelbrot.transpose()

Lindenmayer systems

Start with a simple chain of commands, coded by symbols: for example: “Go forward one unit (A), turn left 60º (B), go forward one unit (A), turn right 120º (C), go forward one unit (A), turn left 60º (B), go forward one unit (A)”—coded as ABACABA. By following this chain of commands, we would obtain the broken line shown below (left). We also have a set of rules that change some of the symbols. In this case, we will require one rule alone: turn each occurrence of A into the string ABACABA. Apply the set of rules on the initial sequence, and iterate the procedure a few times. For example, after two iterations, we have the string

\begin{array}{l} ABACABABABACABACABACABABABACABABABACABA\dotsb \\ BABACABACABACABABABACABACABACABABABACAB\dotsb \\ ACABACABABABACABABABACABABABACABACABACA\dotsb \\ BABABACABA \end{array}

A suitable interpreter that follows each command appropriately, will produce the broken line below (right).

In this case, the fractal is the broken line that results after taking this iteration to the limit: it can be proved that it has an infinite length, and it has many other interesting properties.

The \LaTeX package tikz offers us a very simple syntax to produce approximations to this so-called Koch snowflake, as well as some other well-known fractals. By including the tikz library decorations.fractals, the two curves mentioned above could be obtained as follows:

\begin{center}
  \begin{tikzpicture}[decoration=Koch snowflake]
    \draw decorate{ (0,0) -- (3,1) };
    \draw decorate{ decorate{ decorate{ (4,0) -- (7,1) } } };
  \end{tikzpicture}
\end{center}

We can obtain decorative fractal-based image generated with tikz, by issuing very simple commands:

\begin{tikzpicture}[decoration=Koch snowflake,scale=2]
  \draw decorate{ decorate{ decorate{ decorate{ (0,0) --
          (4,0) -- (4,4) -- (0,4) -- cycle } } } };
  \draw decorate{ decorate{ decorate{ decorate{ (2,1.5) --
          (2.5,2) -- (2,2.5) -- (1.5,2) -- cycle } } } };
\end{tikzpicture}
Advertisements
  1. MonkeyShoulder
    October 31, 2013 at 6:02 am

    Hello!

    Can you tell us what pyplot function do you use for the plotting of the mandelbrot set. pcolormesh?

    • March 20, 2014 at 5:03 pm

      Thanks for reading! I believe I used matplotlib.pyplot.pcolor in this example, although matplotlib.pyplot.imshow does the trick too.

  1. September 24, 2013 at 1:05 am

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: