© Springer Nature Switzerland AG 2018
Michael Oberguggenberger and Alexander OstermannAnalysis for Computer ScientistsUndergraduate Topics in Computer Sciencehttps://doi.org/10.1007/978-3-319-91155-7_9

9. Fractals and L-systems

Michael Oberguggenberger1   and Alexander Ostermann1  
(1)
University of Innsbruck, Innsbruck, Austria
 
 
Michael Oberguggenberger (Corresponding author)
 
Alexander Ostermann
In geometry objects are often defined by explicit rules and transformations which can easily be translated into mathematical formulas. For example, a circle is the set of all points which are at a fixed distance r from a centre (ab):
$$ K = \{(x, y) \in \mathbb R^2 \;;\; (x - a)^2 + (y -b)^2 = r^2\} $$
or
$$ K = \{(x, y) \in \mathbb R^2 \;;\; x = a + r \cos \varphi ,\ y = b + r \sin \varphi , \ 0 \le \varphi < 2\pi \}. $$
In contrast to that, the objects of fractal geometry are usually given by a recursion. These fractal sets (fractals) have recently found many interesting applications, e.g. in computer graphics (modelling of clouds, plants, trees, landscapes), in image compression and data analysis. Furthermore fractals have a certain importance in modelling growth processes.

Typical properties of fractals are often their non-integer dimension and the self-similarity of the entire set with its pieces. The latter can frequently be found in nature, e.g. in geology. There it is often difficult to decide from a photograph without a given scale whether the object in question is a grain of sand, a pebble or a large piece of rock. For that reason fractal geometry is often exuberantly called the geometry of nature.

In this chapter we exemplarily have a look at fractals in $$\mathbb R^2$$ and $$\mathbb C$$. Furthermore we give a short introduction to L-systems and discuss, as an application, a simple concept for modelling the growth of plants. For a more in-depth presentation we refer to the textbooks [21, 22].

9.1 Fractals

To start with we generalise the notions of open and closed interval to subsets of $$\mathbb R^2$$. For a fixed $$\mathbf {a}=(a, b) \in \mathbb R^2$$ and $$\varepsilon > 0$$ the set
$$ B(\mathbf {a},\varepsilon ) = \left\{ (x, y)\in \mathbb R^2 \;;\; \sqrt{(x-a)^2 + (y-b)^2} < \varepsilon \right\} $$
is called an $$\varepsilon $$-neighbourhood of $$\mathbf {a}$$. Note that the set $$B(\mathbf {a},\varepsilon )$$ is a circular disc (with centre $$\mathbf {a}$$ and radius $$\varepsilon $$) where the boundary is missing.
images/215236_2_En_9_Chapter/215236_2_En_9_Fig1_HTML.gif
Fig. 9.1

Open (left), closed (middle) and neither open nor closed (right) square with side length 1

Definition 9.1

Let $$A \subseteq \mathbb R^2$$.

  1. (a)

    A point $$\mathbf {a}\in A$$ is called interior point of A if there exists an $$\varepsilon $$-neighbourhood of $$\mathbf {a}$$ which itself is contained in A.

     
  2. (b)

    A is called open if each point of A is an interior point.

     
  3. (c)

    A point $$\mathbf {c}\in \mathbb R^2$$ is called boundary point of A if every $$\varepsilon $$-neighbourhood of $$\mathbf {c}$$ contains at least one point of A as well as a point of $$\mathbb R^2 \setminus A$$. The set of boundary points of A is denoted by $$\partial A$$ (boundary of A).

     
  4. (d)

    A set is called closed if it contains all its boundary points.

     
  5. (e)

    A is called bounded if there is a number $$r > 0$$ with $$A \subseteq B (\mathbf {0}, r)$$.

     

Example 9.2

The square
$$ Q=\{(x, y) \in \mathbb R^2 \;;\; 0<x<1\ \text {and}\ 0<y<1\} $$
is open since every point of Q has an $$\varepsilon $$-neighbourhood which is contained in Q, see Fig. 9.1, left picture. The boundary of Q consists of four line segments
$$ \{0,1\}\times [0,1]\ \cup \ [0,1]\times \{0,1\}. $$
Every $$\varepsilon $$-neighbourhood of a boundary point also contains points which are outside of Q, see Fig. 9.1, middle picture. The square in Fig. 9.1, right picture,
$$ \{(x, y) \in \mathbb R^2 \;;\; 0<x\le 1\ \text {and}\ 0<y\le 1\} $$
is neither closed nor open since the boundary point $$(x, y)=(0,0)$$ is not an element of the set and the set on the other hand contains the point $$(x, y)=(1,1)$$ which is not an inner point. All three sets are bounded since they are, for example, contained in $$B(\mathbf {0}, 2)$$.

Fractal dimension. Roughly speaking, points have dimension 0, line segments dimension 1 and plane regions dimension 2. The concept of fractal dimension serves to make finer distinctions. If, for example, a curve fills a plane region densely one tends to assign to it a higher dimension than 1. Conversely, if a line segment has many gaps, its dimension could be between 0 and 1.

Let $$A \subseteq \mathbb R^2$$ be bounded (and not empty) and let $$N (A, \varepsilon )$$ be the smallest number of closed circles with radius $$\varepsilon $$ which are needed to cover A, see Fig. 9.2.
images/215236_2_En_9_Chapter/215236_2_En_9_Fig2_HTML.gif
Fig. 9.2

Covering a curve using circles

The following intuitive idea stands behind the definition of the fractal dimension d of A: For curve segments the number $$N(A,\varepsilon )$$ is inverse proportional to $$\varepsilon $$, for plane regions inverse proportional to $$\varepsilon ^2$$, so
$$ N(A, \varepsilon ) \ \approx \ C \cdot \varepsilon ^{-d}, $$
where d denotes the dimension. Taking logarithms one obtains
$$ \log N (A, \varepsilon ) \ \approx \ \log C - d \log \varepsilon , $$
and
$$ d \ \approx \ - \frac{\log N (A, \varepsilon ) - \log C}{\log \varepsilon }, $$
respectively. This approximation is getting more precise the smaller one chooses $$\varepsilon > 0$$. Due to
$$ \lim _{\varepsilon \rightarrow 0^+} \frac{\log C}{\log \varepsilon } = 0 $$
this leads to the following definition.

Definition 9.3

Let $$A \subseteq \mathbb R^2$$ be not empty, bounded and $$N(A, \varepsilon )$$ as above. If the limit
$$ d = d(A) = - \lim _{\varepsilon \rightarrow 0^+} \frac{\log N (A, \varepsilon )}{\log \varepsilon } $$
exists, then d is called fractal dimension of A.

Remark 9.4

In the above definition it is sufficient to choose a zero sequence of the form
$$ \varepsilon _n = C \cdot q^n, \qquad 0< q < 1 $$
for $$\varepsilon $$. Furthermore it is not essential to use circular discs for the covering. One can just as well use squares, see [5, Chap. 5]. Hence the number obtained by Definition 9.3 is also called box-dimension of A.
Experimentally the dimension of a fractal can be determined in the following way: For various rasters of the plane with mesh size $$\varepsilon _n$$ one counts the number of boxes which have a non-empty intersection with the fractal, see Fig. 9.3. Let us call this number again $$N(A, \varepsilon _n)$$. If one plots $$\log N(A, \varepsilon _n)$$ as a function of $$\log \varepsilon _n$$ in a double-logarithmic diagram and fits the best line to this graph (Sect. 18.​1), then
$$ d(A) \approx - \ \text {slope of the straight line}. $$
With this procedure one can, for example, determine the fractal dimension of the coastline of Great Britain, see Exercise 1.
images/215236_2_En_9_Chapter/215236_2_En_9_Fig3_HTML.gif
Fig. 9.3

Raster of the plane using squares of side length $$\varepsilon $$. The boxes that have a non-empty intersection with the fractal are coloured in grey. In the picture we have $$N(A,\varepsilon )=27$$

images/215236_2_En_9_Chapter/215236_2_En_9_Fig4_HTML.gif
Fig. 9.4

Covering of a straight line segment using circles

Example 9.5

The line segment (Fig. 9.4)
$$ A = \{ (x, y) \in \mathbb R^2 \;;\; a \le x \le b, \ y =c\} $$
has fractal dimension $$d = 1$$.
We choose
$$ \varepsilon _n = (b-a) \cdot 2^{-n}, \qquad q = 1/2. $$
Due to $$N(A, \varepsilon _n) = 2^{n-1}$$ the following holds
$$ - \frac{\log N(A, \varepsilon _n)}{\log \varepsilon _n} = - \frac{(n-1) \log 2}{\log (b-a) - n \log 2} \rightarrow 1\qquad \text {as}\ n\rightarrow \infty . $$
Likewise, it can easily be shown: Every set that consists of finitely many points has fractal dimension 0. Plane regions in $$\mathbb R^2$$ have fractal dimension 2. The fractal dimension is in this way a generalisation of the intuitive notion of dimension. Still, caution is advisable here as can be seen in the following example.
images/215236_2_En_9_Chapter/215236_2_En_9_Fig5_HTML.gif
Fig. 9.5

A set of points with box-dimension $$d=\frac{1}{2}$$

Example 9.6

The set $$F=\left\{ 0,1,\frac{1}{2},\frac{1}{3},\frac{1}{4},\ldots \right\} $$ displayed in Fig. 9.5 has box-dimension $$d=1/2$$. We check this claim with the following MATLAB experiment.

Experiment 9.7

To determine the dimension of F approximately with the help of MATLAB we take the following steps. For $$j=1,2,3,\ldots $$ we split the interval [0, 1] into $$4^j$$ equally large subintervals, set $$\varepsilon _j=4^{-j}$$ and determine the number $$N_j=N(F, \varepsilon _j)$$ of subintervals which have a non-empty intersection with F. Then we plot $$\log N_j$$ as a function of $$\log \varepsilon _j$$ in a double-logarithmic diagram. The slope of the secant
$$ d_j = -\frac{\log N_{j+1}-\log N_j}{\log \varepsilon _{j+1}-\log \varepsilon _j} $$
is an approximation to d which is steadily improving with growing j. The values obtained by using the program mat09_1.m are given in the following table:

$$4^j$$

4

16

64

256

1024

4096

16384

65536

262144

1048576

$$d_j$$

0.79

0.61

0.55

0.52

0.512

0.5057

0.5028

0.5014

0.5007

0.50035

Verify the given values and determine that the approximations given by Definition 9.3
$$ \widetilde{d}_j = -\frac{\log N_{j}}{\log \varepsilon _{j}} $$
are much worse. Explain this behaviour.
images/215236_2_En_9_Chapter/215236_2_En_9_Fig6_HTML.gif
Fig. 9.6

The construction of the Cantor set

Example 9.8

(Cantor set)   We construct this set recursively using
$$\begin{aligned} A_0= & {} [0,1]\\ A_1= & {} \left[ 0, \tfrac{1}{3}\right] \cup \left[ \tfrac{2}{3}, 1\right] \\ A_2= & {} \left[ 0, \tfrac{1}{9}\right] \cup \left[ \tfrac{2}{9}, \tfrac{1}{3}\right] \cup \left[ \tfrac{2}{3}, \tfrac{7}{9}\right] \cup \left[ \tfrac{8}{9}, 1\right] \\&\vdots&\end{aligned}$$
One obtains $$A_{n+1}$$ from $$A_n$$ by removing the middle third of each line segment of $$A_n$$, see Fig. 9.6.
The intersection of all these sets
$$ A = \bigcap _{n=0}^\infty A_n $$
is called Cantor set. Let $$|A_n|$$ denote the length of $$A_n$$. Obviously the following holds true: $$|A_0| = 1$$, $$|A_1| = 2/3$$, $$|A_2| = (2/3)^2$$ and $$|A_n| = (2/3)^n$$. Thus
$$ \left| A \right| = \lim _{n \rightarrow \infty } |A_n| = \lim _{n \rightarrow \infty } (2/3)^n = 0, $$
which means that A has length 0. Nevertheless, A does not simply consist of discrete points. More information about the structure of A is given by its fractal dimension d. To determine it we choose
$$ \varepsilon _n = \frac{1}{2} \cdot 3^{-n}, \quad \text {i.e.} \ q = 1/3, $$
and obtain (according to Fig. 9.6) the value $$N (A,\varepsilon _n) = 2^n$$. Thus
$$ d = - \lim _{n \rightarrow \infty } \frac{\log 2^n}{\log 3^{-n}- \log 2} = \lim _{n \rightarrow \infty } \frac{n\log 2}{n\log 3 + \log 2} = \frac{\log 2}{\log 3} = 0.6309... $$
The Cantor set is thus an object between points and straight lines. The self-similarity of A is also noteworthy. Enlarging certain parts of A results in copies of A. This together with the non-integer dimension is a typical property of fractals.
images/215236_2_En_9_Chapter/215236_2_En_9_Fig7_HTML.gif
Fig. 9.7

Snowflakes of depth 0, 1, 2, 3 and 4

images/215236_2_En_9_Chapter/215236_2_En_9_Fig8_HTML.gif
Fig. 9.8

Law of formation of the snowflake

Example 9.9

(Koch’s snowflake1)   This is a figure of finite area whose boundary is a fractal of infinite length. In Fig. 9.7 one can see the first five construction steps of this fractal. In the step from $$A_n$$ to $$A_{n+1}$$ we substitute each straight boundary segment by four line segments in the following way: We replace the central third by two sides of an equilateral triangle, see Fig. 9.8.

The perimeter $$U_n$$ of the figure $$A_n$$ is computed as
$$ U_n = \frac{4}{3} \; U_{n-1} = \left( \frac{4}{3}\right) ^2 U_{n-2}= \cdots = \left( \frac{4}{3}\right) ^n U_0 = 3a \left( \frac{4}{3}\right) ^n. $$
Hence the perimeter $$U_\infty $$ of Koch’s snowflake $$A_\infty $$ is
$$ U_\infty = \lim _{n \rightarrow \infty } U_n = \infty . $$
Next we compute the fractal dimension of $$\partial A_\infty $$. For that we set
$$ \varepsilon _n = \frac{a}{2} \cdot 3^{-n}, \quad \text {i.e.} \ q = 1/3. $$
Since one can use a circle of radius $$\varepsilon _n$$ to cover each straight boundary piece, we obtain
$$ N (\partial A_\infty , \varepsilon _n) \le 3 \cdot 4^n $$
and hence
$$ d = d (\partial A_\infty ) \le \frac{\log 4}{\log 3} \ \approx \ 1.262. $$
A covering using equilateral triangles of side length $$\varepsilon _n$$ shows that $$N(\partial A_\infty , \varepsilon _n)$$ is proportional to $$4^n$$ and thus
$$ d = \displaystyle \frac{\log 4}{\log 3}. $$
The boundary of the snowflake $$\partial A_\infty $$ is hence a geometric object between a curve and a plane region.

9.2 Mandelbrot Sets

An interesting class of fractals can be obtained with the help of iteration methods. As an example we consider in $$\mathbb C$$ the iteration
$$ z_{n+1} = z^2_n + c. $$
Setting $$z = x + \mathrm{{i}\,} y$$ and $$c = a + \mathrm{{i}\,} b$$ one obtains, by separating the real and the imaginary part, the equivalent real form of the iteration
$$\begin{aligned} x_{n+1}&= x^2_n - y^2_n + a,\\ y_{n+1}&= 2 x_n y_n + b. \end{aligned}$$
The real representation is important when working with a programming language that does not support complex arithmetic.
First we investigate for which values of $$c \in \mathbb C$$ the iteration
$$ z_{n+1} = z^2_n + c, \qquad z_0 = 0 $$
remains bounded. In the present case this is equivalent to $$|z_n| \not \rightarrow \infty $$ for $$n \rightarrow \infty $$. The set of all c with this property is obviously not empty since it contains $$c = 0$$. On the other hand it is bounded since the iteration always diverges for $$|c| > 2$$ as can easily be verified with MATLAB .

Definition 9.10

The set
$$ M = \{ c \in \mathbb C\;;\; |z_n| \not \rightarrow \infty \ \text {as} \ n\rightarrow \infty \} $$
is called Mandelbrot set 2 of the iteration $$z_{n+1} = z^2_n + c, \ z_0 = 0$$.
To get an impression of M we carry out a numerical experiment in MATLAB.
images/215236_2_En_9_Chapter/215236_2_En_9_Fig9_HTML.gif
Fig. 9.9

The Mandelbrot set of the iteration $$z_{n+1} = z^2_n + c$$, $$z_0 = 0$$ and enlargement of a section

Experiment 9.11

To visualise the Mandelbrot set M one first chooses a raster of a certain region, for example
$$ -2 \le \mathrm {Re}\, {\textit{c}} \le 1, \quad -1.15 \le \mathrm {Im}\, {\textit{c}} \le 1.15. $$
Next for each point of the raster one carries out a large number of iterations (e.g. 80) and decides then whether the iterations remain bounded (for example $$|z_n|\le 2$$). If this is the case one colours the point in black. This way one successively obtains a picture of M. For your experiments use the MATLAB program mat09_2.m and modify it as required. This way generate in particular the pictures in Fig. 9.9 in high resolution.

Figure 9.9 shows as result a little apple man which has smaller apple men attached which finally develop into an antenna. Here one already recognises the self-similarity. If an enlargement of a certain detail on the antenna $$(-1.8 \le \mathrm{{Re}}~c \le -1.72 ,\ -0.03 \le \mathrm{{Im}}~c \le 0.03)$$ is made, one finds an almost perfect copy of the complete apple man. The Mandelbrot set is one of the most popular fractals and one of the most complex mathematical objects which can be visualised.

9.3 Julia Sets

Again we consider the iteration
$$ z_{n+1} = z^2_n + c. $$
This time, however, we interchange the roles of $$z_0$$ and c.

Definition 9.12

For a given $$c \in \mathbb C$$, the set
$$ J_c = \{z_0 \in \mathbb C\;;\; |z_n| \not \rightarrow \infty \ \text {as} \ n \rightarrow \infty \} $$
is called Julia set 3 of the iteration $$z_{n+1} = z^2_n + c$$.
images/215236_2_En_9_Chapter/215236_2_En_9_Fig10_HTML.gif
Fig. 9.10

Julia sets of the iteration $$z_{n+1} = z^2_n + c$$ for the parameter values $$c=-0.75$$ (top left), $$c = 0.35+0.35\,\mathrm {i}$$ (top right), $$c= -0.03+0.655\,\mathrm {i}$$ (bottom left) and $$-0.12+0.74\,\mathrm {i}$$ (bottom right)

The Julia set for the parameter value c hence consists of those initial values for which the iteration remains bounded. For some values of c the pictures of $$J_c$$ are displayed in Fig. 9.10. Julia sets have many interesting properties; for example,
$$ J_c \ \text {is connected} \ \Leftrightarrow \ c \in M. $$
Thus one can alternatively define the Mandelbrot set M as
$$ M = \{ c \in \mathbb C\;;\; J_c \ \text {is connected}\}. $$
Furthermore the boundary of a Julia set is self-similar and a fractal.

Experiment 9.13

Using the MATLAB program mat09_3.m plot the Julia sets $$J_c$$ in Fig. 9.10 in high definition. Also try other values of c.

9.4 Newton’s Method in $$\mathbb C$$

Since the arithmetic in $$\mathbb C$$ is an extension of that in $$\mathbb R$$, many concepts of real analysis can be transferred directly to $$\mathbb C$$. For example, a function f: $$\mathbb C\rightarrow \mathbb C$$: $$z \mapsto f(z)$$ is called complex differentiable if the difference quotient
$$ \frac{f(z+\varDelta z)-f(z)}{\varDelta z} $$
has a limit as $$\varDelta z\rightarrow 0$$. This limit is again denoted by
$$ f'(z) = \lim _{\varDelta z\rightarrow 0} \frac{f(z+\varDelta z)-f(z)}{\varDelta z} $$
and called complex derivative of f at the point z. Since differentiation in $$\mathbb C$$ is defined in the same way as differentiation in $$\mathbb R$$, the same differentiation rules hold. In particular any polynomial
$$ f(z) = a_nz^n+\cdots + a_1z+a_0 $$
is complex differentiable and has the derivative
$$ f'(z) = na_nz^{n-1} + \cdots + a_1. $$
Like the real derivative (see Sect. 7.​3), the complex derivative has an interpretation as a linear approximation
$$ f(z)\approx f(z_0) + f'(z_0)(z-z_0) $$
for z close to $$z_0$$.
Let f: $$\mathbb C\rightarrow \mathbb C$$: $$z \mapsto f(z)$$ be a complex differentiable function with $$f(\zeta ) = 0$$ and $$f'(\zeta ) \ne 0$$. In order to compute the zero $$\zeta $$ of the function f, one first computes the linear approximation starting from the initial value $$z_0$$, so
$$ z_1 = z_0 - \frac{f(z_0)}{f'(z_0)}. $$
Subsequently $$z_1$$ is used as the new initial value and the procedure is iterated. In this way one obtains Newton’s method in $$\mathbb C$$:
$$ z_{n+1} = z_n - \frac{f(z_n)}{f'(z_n)}. $$
For initial values $$z_0$$ close to $$\zeta $$ the procedure converges (as in $$\mathbb R$$) quadratically. Otherwise, however, the situation can become very complicated.
In 1983 Eckmann [9] investigated Newton’s method for the function
$$ f(z) = z^3 -1 = (z-1)(z^2 + z + 1). $$
This function has three roots in $$\mathbb C$$
$$ \zeta _1 = 1, \quad \zeta _{2,3} = - \frac{1}{2} \pm \mathrm {i}\frac{\sqrt{3}}{2}. $$
Naively one could think that the complex plane $$\mathbb C$$ is split into three equally large sectors where the iteration with initial values in sector $$S_1$$ converges to $$\zeta _1$$, the ones in $$S_2$$ to $$\zeta _2$$, etc., see Fig. 9.11.
images/215236_2_En_9_Chapter/215236_2_En_9_Fig11_HTML.gif
Fig. 9.11

Possible regions of attraction of Newton’s iteration for finding the roots of $$z^3-1$$

A numerical experiment, however, shows that it is not that way. If one colours the initial values according to their convergence, one obtains a very complex picture. One can prove (however, not easily imagine) that at every point where two colours meet, the third colour is also present. The boundaries of the regions of attraction are dominated by pincer-like motifs which reappear again and again when enlarging the scale, see Fig. 9.12. The boundaries of the regions of attraction are Julia sets. Again we have found fractals.
images/215236_2_En_9_Chapter/215236_2_En_9_Fig12_HTML.gif
Fig. 9.12

Actual regions of attraction of Newton’s iteration for finding the roots of $$z^3-1$$ and enlargement of a part

Experiment 9.14

Using the MATLAB program mat09_4.m carry out an experiment. Ascertain yourself of the self-similarity of the appearing Julia sets by producing suitable enlargements of the boundaries of the region of attraction.

9.5 L-systems

The formalism of L-systems was developed by Lindenmayer 4 around 1968 in order to model the growth of plants. It also turned out that many fractals can be created this way. In this section we give a brief introduction to L-systems and discuss a possible implementation in maple.

Definition 9.15

An L-system consists of the following five components:
  1. (a)

    A finite set B of symbols, the so-called alphabet. The elements of B are called letters, and any string of letters is called a word.

     
  2. (b)

    Certain substitution rules. These rules determine how the letters of the current word are to be replaced in each iteration step.

     
  3. (c)

    The initial word $$w\in W$$. The initial word is also called axiom or seed.

     
  4. (d)

    The number of iteration steps which one wants to carry out. In each of these steps, every letter of the current word is replaced according to the substitution rules.

     
  5. (e)

    A graphical interpretation of the word.

     
Let W be the set of all words that can be formed in the given L-system. The substitution rules can be interpreted as a mapping from B to W:
$$ S:B\rightarrow W:\mathtt{b}\mapsto S(\mathtt{b}). $$

Example 9.16

Consider the alphabet $$B=\{\mathtt{f,p, m}\}$$ consisting of the three letters f, p and m. As substitution rules for this alphabet we take
$$ S(\mathtt{f})=\mathtt{fpfmfmffpfpfmf},\quad S(\mathtt{p})=\mathtt{p}, \quad S(\mathtt{m})=\mathtt{m} $$
and consider the axiom $$w=\mathtt{fpfpfpf}$$. An application of the substitution rules shows that, after one substitution, the word fpf becomes the new word fpfmfmffpfpfmfpfpfmfmffpfpfmf. If one applies the substitution rules on the axiom then one obtains a new word. Applying the substitution rules on that again gives a new word, and so on. Each of these words can be interpreted as a polygon by assigning the following meaning to the individual letters:
images/215236_2_En_9_Chapter/215236_2_En_9_Figa_HTML.gif
Thereby $$0\le \alpha \le \pi $$ is a chosen angle. One plots the polygon by choosing an arbitrary initial point and an arbitrary initial direction. Then one sequentially processes the letters of the word to be displayed according to the rules above.
In maple lists and the substitution command subs lend themselves to the implementation of L-systems. In the example above the axiom would hence be defined by
$$ \mathtt{a := [f,p,f,p,f,p, f]} $$
the substitution rules would be
$$ \mathtt{a -> subs(f=(f,p,f,m,f,m,f,f,p,f,p,f,m,f), a)}. $$
The letters p and m do not change in the example, and they are fixed points in the construction. For the purpose of visualisation one can use polygons in maple, given by lists of points (in the plane). These lists can be plotted easily using the command plot.

Construction of fractals. With the graphical interpretation above and $$\alpha =\pi /2$$, the axiom fpfpfpf is a square which is passed through in a counterclockwise direction. The substitution rule converts a straight line segment into a zigzag line. By an iterative application of the substitution rule the axiom develops into a fractal.

Experiment 9.17

Using the maple worksheet mp09_1.mws create different fractals. Further, try to understand the procedure fractal in detail.

Example 9.18

The substitution rule for Koch’s curve is
$$ \texttt {a\,\,-> subs(f=(f,p,f,m,m,f,p,f), a)}. $$
Depending on which axiom one uses, one can build fractal curves or snowflakes from that, see the maple worksheet mp09_1.mws.
Simulation of plant growth. As a new element branchings (ramifications) are added here. Mathematically one can describe this using two new symbols:
images/215236_2_En_9_Chapter/215236_2_En_9_Figb_HTML.gif
Let us look, for example, at the word
$$ \mathtt{[f,p,f,v,p,p,f,p,f,e,v,m,f,m,f,e,f,p,f,v,p,f,p,f,e,f,m, f]}. $$
If one removes all branchings that start with v and end with e from the list then one obtains the stem of the plant
$$ \mathtt{stem := [f,p,f,f,p,f,f,m, f]}. $$
After the second f in the stem obviously a double branching is taking place and the branches sprout
$$ \mathtt{branch1 := [p,p,f,p, f] \quad \mathrm{and}\quad branch2 := [m,f,m,f]}. $$
Further up the stem branches again with the branch [p,f,p, f].

For a more realistic modelling one can introduce additional parameters. For example, asymmetry can be build in by rotating by the positive angle $$\alpha $$ at p and by the negative angle $$-\beta $$ at m. In the program mp09_2.mws that was done, see Fig. 9.13.

Experiment 9.19

Using the maple worksheet mp09_2.mws create different artificial plants. Further, try to understand the procedure grow in detail.

To visualise the created plants one can use lists of polygons in maple, i.e. lists of points (in the plane). To implement the branchings one conveniently uses a recursive stack. Whenever one comes across the command v for a branching, one saves the current state as the topmost value in the stack. A state is described by three numbers (xyt) where x and y denote the position in the (xy)-plane and t the angle enclosed the with the positive x-axis. Conversely one removes the topmost state from the stack if one comes across the end of a branch e and returns to this state in order to continue the plot. At the beginning the stack is empty (at the end it should be as well).
images/215236_2_En_9_Chapter/215236_2_En_9_Fig13_HTML.gif
Fig. 9.13

Plants created using the maple worksheet mp09_2.mws

Extensions. In the context of L-systems many generalisations are possible which can make the emerging structures more realistic. For example one could:
  1. (a)

    Represent the letter f by shorter segments as one moves further away from the root of the plant. For that, one has to save the distance from the root as a further state parameter in the stack.

     
  2. (b)
    Introduce randomness by using different substitution rules for one and the same letter and in each step choosing one at random. For example, the substitution rules for random weeds could be as such:
    images/215236_2_En_9_Chapter/215236_2_En_9_Figc_HTML.gif
    Using random numbers one selects the according rule in each step.
     

Experiment 9.20

Using the maple worksheet mp09_3.mws create random plants. Further, try to understand the implemented substitution rule in detail.

9.6 Exercises

1.

Determine experimentally the fractal dimension of the coastline of Great Britain. In order to do that, take a map of Great Britain (e.g. a copy from an atlas) and raster the map using different mesh sizes (e.g. with 1 / 64th, 1 / 32th, 1 / 16th, 1 / 8th and 1 / 4th of the North–South expansion). Count the boxes which contain parts of the coastline and display this number as a function of the mesh size in a double-logarithmic diagram. Fit the best line through these points and determine the fractal dimension in question from the slope of the straight line.

2.

Using the program mat09_3.m visualise the Julia sets of $$z_{n+1} = z_n^2 + c$$ for $$c = -1.25$$ and $$c = 0.365 - 0.3\, \mathrm {i}$$. Search for interesting details.

3.
Let $$f(z) = z^3 - 1$$ with $$z = x + \mathrm {i}y$$. Use Newton’s method to solve $$f(z) = 0$$ and separate the real part and the imaginary part, i.e. find the functions $$g_1$$ and $$g_2$$ with
$$\begin{aligned} x_{n+1}&= g_1 (x_n, y_n),\\ y_{n+1}&= g_2 (x_n, y_n). \end{aligned}$$
4.

Modify the procedure grow in the program mp09_2.mws by representing the letter f by shorter segments depending on how far it is away from the root. With that plot the umbel from Experiment 9.19 again.

5.

Modify the program mp09_3.mws by attributing new probabilities to the existing substitution rules (or invent new substitution rules). Use your modified program to plot some plants.

6.

Modify the program mat09_3.m to visualise the Julia sets of $$z_{n+1} = z_n^{2k} - c$$ for $$c=-1$$ and integer values of k. Observe how varying k affects the shape of the Julia set. Try other values of c as well.

7.
Modify the program mat09_3.m to visualise the Julia sets of
$$\begin{aligned} z_{n+1} = z_n^3 + (c-1)z_n - c. \end{aligned}$$
Study especially the behaviour of the Julia sets when c ranges between 0.60 and 0.65.