In this lecture, we extend the ideas of similarity and symmetry previously studied to consider cases of self-similarity which leads to fractals and somewhat amazingly explains the basics of plant structure and evolution.

- Symmetry, by Hermann Weyl, 1952.
- Algorithmic Beauty of Plants by P. Prusinkiewicz
- Chaos: Making a new science by James Gleick
- Fractals using iterated function systems (IFSs) at Yale
- The man who uncovered the secret lives of snowflakes
- Kepler and the snowflake

In previous lectures, we saw how various aspects of natural systems were sometimes related by scaling laws \(y= a x^b\). These scaling laws specify certain kinds of symmetry – changes in \(x\) have to be matched by corresponding changes in \(y\). Symmetry is a mysterious but common phenomena in nature. We humans and many animals have bilateral symmetry. Snowflakes, bee hives, plants, and other things have their own symmetries. Some examples are shown below, ranging in scale from the microscopic to the intergalactic.

We can often make use of symmetries in our models. For example, to approximate the area of a circle, it is enough to integrate the area of a quarter-circle and then multiply by 4: \[\pi = 4 \int _ {0}^1 \sqrt{1-x^2} dx.\] When a problem is approximately radially symmetric, we might simplify our equations by using polar coordinates. Or when a pattern periodicaly repeats, we can make use of Fourier’s analysis.

The study of symmetries makes use of a topic in mathematics called “group theory” – you may recall that we referred to “dimensionless groups” in our lectures on dimensional analysis. We will not delve into group theory here; while the basic ideas are relatively intuitive, the abstract mathematical translation of these ideas is not. However, there is one subtopic in symmetries that is quite interesting from a modeling perspective: *self-similarity*. Self-similarity is a phenomena we sometimes observe in the natural world, were small parts of an object are very similar to larger parts that they make up.

The most common example of this is a fern leaf, where each large leaf is make up of smaller leaves with similar shapes to the large leaf they compose. In mathematics, there are objects like the Mandelbrot set where this self-similarity goes on over and over again at smaller and smaller scales forever – we call such objects *fractals*.

Assume we have 3 people trying to persuade each other about their beliefs.

[Show code]

```
#!/usr/bin/env python
"""
Chaos game described in "Strange new science of chaos" documentary.
https://www.youtube.com/watch?v=fUsePzlOmxw
See also, The color of infinity with A. C. Clark
https://www.youtube.com/watch?v=Lk6QU94xAb8
and "Fractal Geometry"
"""
import random
from matplotlib.pyplot import *
from scipy import exp, pi
#b = (1.+0j, 0+1.j, 0.+0j)
#b = (1, exp(pi*8j/12.), exp(-pi*8j/12) )
b = (1j, exp(pi*14j/12.), exp(-pi*2j/12) )
def F(r, n_steps = 4000):
x = random.choice(b)
e = []
for t in range(n_steps):
#e.append((x.real + 0.5*x.imag, x.imag*c ))
e.append( (x.real, x.imag) )
x += r*(random.choice(b) - x)
return e
def main():
r = 0.50
e = F(r)
u, v = zip(*e)
plot(u,v,'.',markersize=1)
xlim(-1.1,1.1)
ylim(-0.6,1.1)
savefig('AttitudesGasket.png', bbox_inches = 'tight', transparent=True)
savefig('AttitudesGasket.pdf', bbox_inches = 'tight')
#show()
main()
```

One of the oldest cases of self-similarity can be found in an arithmetic problem that is more than 800 years old – though at the time, the author probably didn’t realize it. In 1202, Leonardo “Fibonacci” Pisano published a book *Liber Abacci* which introduced the western world to the Hindu-Arabic place-valued decimal system we still use today. Within the book was the following problem.

Direct translation version:

How Many Pairs of Rabbits Are Created by One Pair in One Year?A certain man had one pair of rabbits together in a certain enclosed place, and one wishes to know how many are created from the pair in one year when it is the nature of them in a single month to bear another pair, and in the second month those born to bear also.

Solution:

Because the above written pair in the first month bore, you will double it; there will be two pairs in one month. One of these, namely the first, bears in the second month, and thus there are in the second month 3 pairs; of these in one month two are pregnant, and in the third month 2 pairs of rabbits are born, and thus there are 5 pairs in the month; in this month 3 pairs are pregnant, and in the fourth month there are 8 pairs, of which 5 pairs bear another 5 pairs; these are added to the 8 pairs making 13 pairs in the fifth month; these 5 pairs that are born in this month do not mate in this month, but another 8 pairs are pregnant, and thus there are in the sixth month 21 pairs; [p284] to these are added the 13 pairs that are born in the seventh month; there will be 34 pairs in this month; to this are added the 21 pairs that are born in the eighth month; there will be 55 pairs in this month; to these are added the 34 pairs that are born in the ninth month; there will be 89 pairs in this month; to these are added again the 55 pairs that are born in the tenth month; there will be 144 pairs in this month; to these are added again the 89 pairs that are born in the eleventh month; there will be 233 pairs in this month.

As you can see, this problem can be solved rather directly, and doesn’t seem to have anything to do with fractals. But let actually construct a model for what was happening to Fibonacci’s pairs of rabits using symbolic algebra. Let each character \(j\) represent one pair of juvenile rabits, and each character \(a\) represent one pair of adult rabbits. According to Fibonacci’s, if we start with one pair of juvenile rabbits, then the population changes as follows from generation to generation.

Month | State | Count |
---|---|---|

0 | j | 1 |

1 | a | 1 |

2 | ja | 2 |

3 | ja a | 3 |

4 | ja a ja | 5 |

5 | ja a ja ja a | 8 |

6 | ja a ja ja a ja a ja | 13 |

7 | ja a ja ja a ja a ja ja a ja ja a | 21 |

These are getting very long, but already, we can see some patterns. Notice that each line is made up of the previous two lines joined together. If we are only interested in the numbers and not the particular order, we can introduce a shorter notation using exponents to do the counting. \(j^x a^y + n\) represents \(x\) the number of juvenile pairs and \(y\) adult pairs in generation \(n\). Then in one generation, each juvenile pair becomes an adult pair (\(j\) becomes \(a\)), and each adult pair becomes creates a new juvenile pair (\(a\) becomes \(aj\)).

\[j^x a^y + n \rightarrow (a)^x (aj)^y + n + 1 = j^y a^{x+y} + n + 1\]

Here, we’ve used the “sum” in a formal sense – since \(a\) and \(j\) are symbols, we cannot actually do the addition and infact never mean to, but using the summation notation is convenient. A similar notation is sometimes used for vectors in physics.

This rewrite rule can also be interpreted as a system of two first order difference equations for the numbers of juvenile and adult pairs after each generation \(n\). \[\begin{align*} x(n+1) &= y(n) \\ y(n+1) &= y(n)+x(n) \end{align*}\] These two equations can be changed into a single equation for the total numbers of pairs \(F(n)\) as follows. \[\begin{align*} F(n) &= y(n) + x(n) \\ F(n+1) &= x(n+1) + y(n+1) = y(n) + y(n) + x(n) = F(n) + y(n) \\ F(n+2) &= F(n+1) + y(n+1) = F(n+1) + y(n) + x(n) = F(n+1) + F(n) \end{align*}\] So, the governing equation for the total number of rabbit pairs, which we can call Fibonacci’s equation, is \[F(n+2) = F(n+1) + F(n).\] This is a second order linear difference equation. Since we know \(F(0) = 1\) and \(F(1) = 1\), we can use it to calculate the sequence directly. \[1,1,2,3,5,8,13,21,34, \ldots\] The general solution can be found using the guess \(F(n) = C \lambda^n\). In order to find the specific solution, we can use our initial conditions to determine the values of \(C\). Python’s sympy toolbox can speed this process up for us.

[Show code]

```
from sympy import symbols, rsolve, Eq, latex
F, n = symbols('F, n')
equation = Eq(F(n+2), F(n) + F(n+1))
initialcondition = {F(0):1, F(1):1}
# rsolve is used to solve common difference equations
sol = rsolve(equation, F(n), initialcondition)
pprint(sol)
print [ sol.subs(n,i).expand() for i in range(10) ] # to check result
print latex(Eq(F(n), sol))
```

\[F{\left (n \right )} = \left(\frac{1}{2} + \frac{\sqrt{5}}{2}\right)^{n} \left(\frac{\sqrt{5}}{10} + \frac{1}{2}\right) + \left(- \frac{\sqrt{5}}{2} + \frac{1}{2}\right)^{n} \left(- \frac{\sqrt{5}}{10} + \frac{1}{2}\right)\]

Weirdly, our sequence is integers, but the solution involves the \(\sqrt{5}\). Actually, it’s not just the \(\sqrt{5}\) but the full **golden ratio** \(\phi:= (1+\sqrt{5})/2 \approx 1.618033\). Since the reciprical of the golden ratio \[\frac{1}{\phi} = \frac{ -1 + \sqrt{5} }{2},\] the general solution of Fibonacci’s equation can be given as \[F(n) = a \phi^n + b (-\phi)^{-n}\] Since \(\phi>1\), then as we do more iterations, the ratio of successive Fibonacci numbers, \[\lim _ {n \rightarrow \infty} F(n+1) / F(n) = \phi\]

Consider the example of Fibonacci’s rabbits:

Month | State | Count |
---|---|---|

0 | j | 1 |

1 | a | 1 |

2 | ja | 2 |

3 | ja a | 3 |

4 | ja a ja | 5 |

5 | ja a ja ja a | 8 |

6 | ja a ja ja a ja a ja | 13 |

7 | ja a ja ja a ja a ja ja a ja ja a | 21 |

Notice that each row can be interpreted as a re-write of the previous row with the following rules: \(a\rightarrow ja\) and \(ja\rightarrow ja\;\;a\). The pattern is generated by the re-write rule, but observe that it’s not totally regular.

Lindemeyer conjectured that the pattern has some meaning. “Re-write systems,” i.e., substitution rules for replacing characters with other characters, lead to patterns that can be interpreted as geometric descriptions of real-world objects.

Before we precede further with our modeling, we need to make a brief diversion to introduce a Logo. Logo was a computer language developed in the 1970s to help kids learn to program. One of the big ideas of Logo was to make programming intuitive using graphics. Logo used a cursor it called the “turtle”. You would give the turtle instructions like “move forward 10 units” and “turn right 50 degrees” and you could see the turtle moving accross the computer screen following your instructions. By giving the turtle a carefully selected list of instructions, you could make it draw cool things.

As it turns out, though, the graphics of the Logo language are more than just a kids’ toy – it can be a powerful tool for geometric modeling. Over the decades, Logo has been reimplemented and extended to keep up with changes in our technology. Turtle graphics have been implemented in javascript, for example. For this class, we’ll make use of a python library.

The following relies on the turtle module in python.

[Show code]

```
#!/usr/bin/env python3
from turtle import *
def RW(oldstr, rule):
"""
A function that performs simultaneous
string rewriting based on the substitutions
describe
Arguments:
oldstr = string which will be rewritten
rule = a dictionary mapping characters to substitute strings
"""
newstr = ''.join([(rule[i] if i in rule else i) for i in oldstr])
return newstr
print("# Fibonacci's rabbits, done by rewritting")
fibrule = {'j': 'a', 'a':'aj'}
s = 'j' # initial condition
for t in range(10):
cnts = dict([ (i,s.count(i)) for i in 'aj'])
print(t, cnts, s)
s = RW(s, fibrule)
#################################################
# rewrite rules for drawing Koch's snowflake
print("# Drawing Koch snowflake")
s = '-f++f++f++' # initial condition is an equilateral triangle
kochrule = {'f': 'f-f++f-f'}
# a dictionary to translate characters into turtle commands
ops = { \
'f': lambda L : fd(L), \
'+': lambda L : rt(60), \
'-': lambda L : lt(60), \
}
tracer(1,0) # speed up drawing
L = 160
for t in range(5):
home()
clear()
[ ops[i](L) for i in s]
q = input("Iteration %d. Press return to continue, q to stop: "%t)
if q == 'q':
break
s = RW(s, kochrule)
L = max(int(L/3.), 1) # rescale side length at each step
#################################################
print("# Sierpinski's triangle gasket")
# we need to add one extra letter
# to our translator
ops['g'] = lambda L : fd(L)
gasketrule = { \
'f' : '-g+f+g-',
'g' : '+f-g-f+',
}
s = 'f' # initial condition
L = 200
f = 1
for t in range(0):
reset()
[ ops[i](L) for i in s]
q = input("Iteration %d. Press return to continue, q to stop: "%t)
if q == 'q':
break
s = RW(s, gasketrule)
L = L/2. # rescale side length at each step
f = 3*f
tracer(f,0)
#################################################
# draw the dragon curve
L = 3
ops['x'] = lambda x: None
ops['y'] = lambda x: None
ops['-'] = lambda x: rt(90)
ops['+'] = lambda x: lt(90)
r = { 'x' : 'x+yf+', 'y': '-fx-y'}
s = 'x'
for t in range(12):
s = RW(s,r)
reset()
[ ops[i](L) for i in s ]
ht()
input("Press return to finish")
```

- Koch snowflake If we study this curve, we see that at each iteration, the center third of each segment is replaced by a triangle at each step. If we let ‘f’ encode forward drawing, ‘-’ encode a left turn of 60 degrees, and ‘+’ encode a right turn of 60 degrees, then we can refine the curve by starting with the string \(f++f++f++\) and at each iteration, substituting \(f-f++f-f\) for \(f\) to replace one segment with a new segment containing part of a triangle.

[Show code]

```
from lsystem import LSystem, Turtle
def main():
# Turtle setup
t = Turtle()
t.setperiod(0) # number of actions between refreshes
t.setcolor('white')
t.setpensize(1)
t.ht().rt(90).pd()
#t.ht().pu().bottom().pd() # initial position and pen state
# map from symbols to turtle actions
actionmap = { \
'f' : lambda : t.forward(z), \
'+' : lambda : t.left(60.), \
'-' : lambda : t.right(60.), \
}
# Rule specifications
rule = ('f', 'f-f++f-f')
ic = 'f++f++f++' # initial condition
n = 4 # number of iterations
z = 2.0 # step size
L = LSystem(rule)
for i in L.itr(n).on(ic).str():
actionmap[i]()
t.sleep(5) # refresh and sleep
main()
```

- Sierpinski’s Triangle gasket is a bit harder to find an encoding for.

We can see that there is a regular replacement pattern, but the pattern seems hard to explain – sometimes we start off turning right, but the other times, we start off turning left. The trick turns out to be that the drawing has left out the *color* of the segments, which are supposed to alternate between green and white (for example). If we take the color into account, we can now quickly see the pattern. Let \(W\) encode the white segments, \(G\) encode the green segments, \(-\) be a left turn of 60 degrees, and \(+\) be a right turn of 60 degrees. Then the iteration is performed based on substitutions of the form

\[W \rightarrow -G+W+G-, \quad G \rightarrow +W-G-W+\]

[Show code]

```
from lsystem import LSystem, Turtle
import sys
def main():
# Turtle setup
t = Turtle()
t.setperiod(1) # number of actions between refreshes
t.setpensize(1)
t.ht().pu().bk(120).rt(90).bk(120).pd() # initial position and pen state
# map from symbols to turtle actions
actionmap = { \
'W' : lambda : t.setcolor('white').forward(z), \
'G' : lambda : t.setcolor('green').forward(z), \
'+' : lambda : t.left(60.), \
'-' : lambda : t.right(60.), \
}
# Rule specifications
rule = (('W','+G-W-G+'), ('G','-W+G+W-'))
ic = 'W' # initial condition
n = 4 # number of iterations
z = 400./2**n # step size
L = LSystem(rule)
#print L.dup(n).on(ic).str()
for i in L.itr(n).on(ic).str():
actionmap[i]()
if __name__ == '__main__' or sys.argv[0] == __name__:
raw_input("<return when done>")
else:
t.sleep(5)
main()
```

Dragon curve

fx @subs (x -> x+yf+, y -> -fx-y)**20

Here is a dragon curve rewrite animation

Now, here is an example program illustrating how rewrite rules can be translated into turtle language to draw a geometric object.

[Show code]

```
from lsystem import LSystem, Turtle
def main():
# Turtle setup
t = Turtle()
t.setperiod(40) # number of actions between refreshes
t.setcolor('green')
t.setpensize(1)
t.ht().pu().bottom().pd() # initial position and pen state
# map from symbols to turtle actions
actionmap = { \
'F' : lambda : t.forward(z), \
'[' : lambda : t.push(), \
']' : lambda : t.pop(), \
'+' : lambda : t.left(d), \
'-' : lambda : t.right(d), \
}
# Rule specifications
rule = ('F','F[+F]F[-F]F')
ic = 'F' # initial condition
n = 5 # number of iterations
d = 25.7 # turn angle
z = 1.0 # step size
L = LSystem(rule)
for i in L.itr(n).on(ic).str():
actionmap[i]()
t.saveimg('lsys_planta.png')
main()
```

One interpretation: modeling is working to find simple ways to express the important features of complicated things.

In the end… can we find a simple way to represent the patterns in different data, like the fern, or snowflake,

Aphids are common insects known to reproduce parthenogenicly, with female aphids giving birth to only females for most of the year. Suppose there is 1 baby female aphid in a woman’s garden today. In 2 days, that aphid will become an adult and give birth to 1 baby aphid each day afterward. Find a single difference equation for the total number of aphids on day \(t\). (

*Do not attempt to solve this equation – it’s too much work for a homework problem.*)Find the general solution of the second order linear homogeneous difference equation \(z _ {t+2} = 2 z _ {t+1} + z _ {t}\) and determine the specific solution when \(z _ 0 =1\) and \(z _ 1 = 2\).

Consider the L-System \((W \rightarrow W-G+W+G-W, G \rightarrow G+W-G-W+G)\) where \(-\) is a left turn of 90 degrees, \(+\) is a right turn of 90 degrees, and both \(W\) and \(G\) represent forward motion of one unit.

Draw 2 iterations of this L-system applied to the initial condition \(W\).

Use your answer to (a) to draw the result of two iterations of the L-system applied to the initial condition \(W+W+W+W\).

Find a short encoding for the figure below using a rewrite rule similar to those discussed in class.

Find the shortest encoding for the figure below using a rewrite rule similar to those discussed in class.