分形與混沌

自然界的很多事物，例如樹木、云彩、山脈、閃電、雪花以及海岸線等等都呈現出傳統的幾何學不能描述的形狀。這些形狀都有如下的特性：

• 有著十分精細的不規則的結構
• 整體與局部相似，例如一根樹杈的形狀和一棵樹很像

分形幾何學就是用來研究這樣一類的幾何形狀的科學，借助計算機的高速計算和圖像顯示，使得我們可以更加深入地直觀地觀察分形幾何。在本章中，讓我們用Python繪制一些經典的分形圖案。

Mandelbrot集合

Mandelbrot(曼德布洛特)集合是在復平面上組成分形的點的集合。

Mandelbrot集合的定義(摘自維基百科)

Mandelbrot集合可以用下面的復二次多項式定義：

其中c是一個復參數。對于每一個c，從z=0開始對函數 進行迭代。

序列 的值或者延伸到無限大，或者只停留在有限半徑的圓盤內。

Mandelbrot集合就是使以上序列不發散的所有c點的集合。

從數學上來講，Mandelbrot集合是一個復數的集合。一個給定的復數c或者屬于Mandelbrot集合，或者不是。

用程序繪制Mandelbrot集合時不能進行無限次迭代，最簡單的方法是使用逃逸時間(迭代次數)進行繪制，具體算法如下：

• 判斷每次調用函數 得到的結果是否在半徑R之內，即復數的模小于R
• 記錄下模大于R時的迭代次數
• 迭代最多進行N次
• 不同的迭代次數的點使用不同的顏色繪制

下面是完整的繪制Mandelbrot集合的程序：

# -*- coding: utf-8 -*-

import numpy as np
import pylab as pl
import time
from matplotlib import cm

def iter_point(c):
    z = c
    for i in xrange(1, 100): # 最多迭代100次
        if abs(z)>2: break # 半徑大于2則認為逃逸
        z = z*z+c
    return i # 返回迭代次數
    
def draw_mandelbrot(cx, cy, d):
    """
    繪制點(cx, cy)附近正負d的范圍的Mandelbrot
    """
    x0, x1, y0, y1 = cx-d, cx+d, cy-d, cy+d 
    y, x = np.ogrid[y0:y1:200j, x0:x1:200j]
    c = x + y*1j
    start = time.clock()
    mandelbrot = np.frompyfunc(iter_point,1,1)(c).astype(np.float)
    print "time=",time.clock() - start
    pl.imshow(mandelbrot, cmap=cm.Blues_r, extent=[x0,x1,y0,y1])
    pl.gca().set_axis_off()
    
x,y = 0.27322626, 0.595153338

pl.subplot(231)
draw_mandelbrot(-0.5,0,1.5)
for i in range(2,7):    
    pl.subplot(230+i)
    draw_mandelbrot(x, y, 0.2**(i-1))
pl.subplots_adjust(0.02, 0, 0.98, 1, 0.02, 0)
pl.show()

Mandelbrot集合，以5倍的倍率放大點(0.273, 0.595)附近

程序中的iter_point函數計算點c的逃逸時間，逃逸半徑R為2.0，最大迭代次數為100。draw_mandelbrot函數繪制以點(cx, cy)為中心，邊長為2*d的正方形區域內的Mandelbrot集合。

下面3行計算指定范圍內的迭代公式的參數c，c是一個元素為復數的二維數組，大小為200*200，注意np.ogrid不是函數：

x0, x1, y0, y1 = cx-d, cx+d, cy-d, cy+d
y, x = np.ogrid[y0:y1:200j, x0:x1:200j]
c = x + y*1j

下面一行程序通過調用np.frompyfunc將iter_point轉換為NumPy的ufunc函數，這樣它可以自動對c中的每個元素調用iter_point函數，由于結果的數組元素類型為object，還需要調用astype方法將其元素類型轉換為浮點類型：

mandelbrot = np.frompyfunc(iter_point,1,1)(c).astype(np.float)

最后調用matplotlib的imshow函數將結果數組繪制成圖，通過cmap關鍵字參數指定圖的值和顏色的映射表：

pl.imshow(mandelbrot, cmap=cm.Blues_r, extent=[x0,x1,y0,y1])

使用Python繪制Mandelbrot集合最大的問題就是運算速度太慢，下面是上面每幅圖的計算時間：

time= 0.88162629608
time= 1.53712748408
time= 1.71502160191
time= 1.8691174437
time= 3.03812691278

因為計算每個點的逃逸時間均不相同，因此每幅圖的計算時間也不相同。

計算速度慢的最大的原因是因為iter_point函數的運算速度慢，如果將此函數用C語言重寫的話將能顯著地提高計算速度，下面使用scipy.weave庫將C++重寫的iter_point函數轉換為Python能調用的函數：

import scipy.weave as weave

def weave_iter_point(c):
    code = """
    std::complex<double> z;
    int i;
    z = c;
    for(i=1;i<100;i++)
    {
        if(std::abs(z) > 2) break;
        z = z*z+c;
    }
    return_val=i;
    """

    f = weave.inline(code, ["c"], compiler="gcc")
    return f

下面是使用weave_iter_point函數計算Mandelbrot集合的時間：

time= 0.285266982256
time= 0.271430028118
time= 0.293769180161
time= 0.308515188383
time= 0.411168179196

通過NumPy的數組運算也可以提高計算速度，前面的計算都是先對復數平面上的每個點進行循環，然后再循環迭代計算每個點的逃逸時間。如果要用NumPy的數組運算加速計算的話，可以將這兩個循環的順序顛倒過來，下面的程序演示這一算法：

# -*- coding: utf-8 -*-

import numpy as np
import pylab as pl
import time
from matplotlib import cm

def draw_mandelbrot(cx, cy, d, N=200):
    """
    繪制點(cx, cy)附近正負d的范圍的Mandelbrot
    """
    global mandelbrot

    x0, x1, y0, y1 = cx-d, cx+d, cy-d, cy+d 
    y, x = np.ogrid[y0:y1:N*1j, x0:x1:N*1j]
    c = x + y*1j

    # 創建X,Y軸的坐標數組
    ix, iy = np.mgrid[0:N,0:N]
    
    # 創建保存mandelbrot圖的二維數組，缺省值為最大迭代次數
    mandelbrot = np.ones(c.shape, dtype=np.int)*100
    
    # 將數組都變成一維的
    ix.shape = -1
    iy.shape = -1
    c.shape = -1
    z = c.copy() # 從c開始迭代，因此開始的迭代次數為1
    
    start = time.clock()
    
    for i in xrange(1,100):
        # 進行一次迭代
        z *= z
        z += c
        # 找到所有結果逃逸了的點
        tmp = np.abs(z) > 2.0
        # 將這些逃逸點的迭代次數賦值給mandelbrot圖
        mandelbrot[ix[tmp], iy[tmp]] = i
        
        # 找到所有沒有逃逸的點
        np.logical_not(tmp, tmp)
        # 更新ix, iy, c, z只包含沒有逃逸的點
        ix,iy,c,z = ix[tmp], iy[tmp], c[tmp],z[tmp]
        if len(z) == 0: break

    print "time=",time.clock() - start
    pl.imshow(mandelbrot, cmap=cm.Blues_r, extent=[x0,x1,y0,y1])
    pl.gca().set_axis_off()
    
x,y = 0.27322626, 0.595153338

pl.subplot(231)
draw_mandelbrot(-0.5,0,1.5)
for i in range(2,7):    
    pl.subplot(230+i)
    draw_mandelbrot(x, y, 0.2**(i-1))
pl.subplots_adjust(0.02, 0, 0.98, 1, 0.02, 0)
pl.show()

為了減少計算次數，程序中每次迭代之后，都將已經逃逸的點剔除出去，這樣就需要保存每個點的下標，程序中用ix和iy這兩個數組來保存沒有逃逸的點的下標，因為有額外的數組保存下標，因此數組z和c不需要是二維的。函數迭代部分的程序如下：

# 進行一次迭代
z *= z
z += c

使用 *=, += 這樣的運算符能夠讓NumPy不分配額外的空間直接在數組z上進行運算。

下面的程序計算出逃逸點，tmp是逃逸點在z中的下標，由于z和ix和iy等數組始終是同時更新的，因此ix[tmp], iy[tmp]就是逃逸點在圖像中的下標：

# 找到所有結果逃逸了的點
tmp = np.abs(z) > 2.0
# 將這些逃逸點的迭代次數賦值給mandelbrot圖
mandelbrot[ix[tmp], iy[tmp]] = i

最后通過對tmp中的每個元素取邏輯反，更新所有沒有逃逸的點的對應的ix, iy, c, z：

# 找到所有沒有逃逸的點
np.logical_not(tmp, tmp)
# 更新ix, iy, c, z只包含沒有逃逸的點
ix,iy,c,z = ix[tmp], iy[tmp], c[tmp], z[tmp]

此程序的計算時間如下：

time= 0.186070576008
time= 0.327006365334
time= 0.372756034636
time= 0.410074464771
time= 0.681048289658
time= 0.878626752841

連續的逃逸時間

修改逃逸半徑R和最大迭代次數N，可以繪制出不同效果的Mandelbrot集合圖案。但是前面所述的方法計算出的逃逸時間是大于逃逸半徑時的迭代次數，因此所輸出的圖像最多只有N種不同的顏色值，有很強的梯度感。為了在不同的梯度之間進行漸變處理，使用下面的公式進行逃逸時間計算：

是迭代n次之后的結果，通過在逃逸時間的計算中引入迭代結果的模值，結果將不再是整數，而是平滑漸變的。

下面是計算此逃逸時間的程序：

def smooth_iter_point(c):
    z = c
    for i in xrange(1, iter_num):
        if abs(z)>escape_radius: break
        z = z*z+c
    absz = abs(z)
    if absz > 2.0:
        mu = i - log(log(abs(z),2),2)
    else:
        mu = i
    return mu # 返回正規化的迭代次數

如果你的逃逸半徑設置得很小，例如2.0，那么有可能結果不夠平滑，這時可以在迭代循環之后添加幾次迭代保證z能夠足夠逃逸，例如：

z = z*z+c
z = z*z+c
i += 2

下圖是逃逸半徑為10，最大迭代次數為20時，繪制的結果：

逃逸半徑=10，最大迭代次數=20的平滑處理后的Mandelbrot集合

逃逸時間公式是如何得出的？

請參考： http://linas.org/art-gallery/escape/ray.html

完整的程序請參考 繪制Mandelbrot集合

迭代函數系統(IFS)

迭代函數系統是一種用來創建分形圖案的算法，它所創建的分形圖永遠是絕對自相似的。下面我們直接通過繪制一種蕨類植物的葉子來說明迭代函數系統的算法：

有下面4個線性函數將二維平面上的坐標進行線性映射變換：

1.
    x(n+1）= 0
    y(n+1) = 0.16 * y(n)

2.
    x(n+1) = 0.2 * x(n) ? 0.26 * y(n)
    y(n+1) = 0.23 * x(n) + 0.22 * y(n) + 1.6

3.
    x(n+1) = ?0.15 * x(n) + 0.28 * y(n)
    y(n+1) = 0.26 * x(n) + 0.24 * y(n) + 0.44

4.
    x(n+1) = 0.85 * x(n) + 0.04 * y(n)
    y(n+1) = ?0.04 * x(n) + 0.85 * y(n) + 1.6

所謂迭代函數是指將函數的輸出再次當作輸入進行迭代計算，因此上面公式都是通過坐標 x(n),y(n) 計算變換后的坐標 x(n+1),y(n+1)。現在的問題是有4個迭代函數，迭代時選擇哪個函數進行計算呢？我們為每個函數指定一個概率值，它們依次為1%, 7%, 7%和85%。選擇迭代函數時使用通過每個函數的概率隨機選擇一個函數進行迭代。上面的例子中，第四個函數被選擇迭代的概率最高。

最后我們從坐標原點(0,0)開始迭代，將每次迭代所得到的坐標繪制成圖，就得到了葉子的分形圖案。下面的程序演示這一計算過程：

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as pl
import time

# 蕨類植物葉子的迭代函數和其概率值
eq1 = np.array([[0,0,0],[0,0.16,0]])
p1 = 0.01

eq2 = np.array([[0.2,-0.26,0],[0.23,0.22,1.6]])
p2 = 0.07

eq3 = np.array([[-0.15, 0.28, 0],[0.26,0.24,0.44]])
p3 = 0.07

eq4 = np.array([[0.85, 0.04, 0],[-0.04, 0.85, 1.6]])
p4 = 0.85

def ifs(p, eq, init, n):
    """
    進行函數迭代
    p: 每個函數的選擇概率列表
    eq: 迭代函數列表
    init: 迭代初始點
    n: 迭代次數
    
    返回值： 每次迭代所得的X坐標數組， Y坐標數組， 計算所用的函數下標    
    """

    # 迭代向量的初始化
    pos = np.ones(3, dtype=np.float)
    pos[:2] = init
    
    # 通過函數概率，計算函數的選擇序列
    p = np.add.accumulate(p)    
    rands = np.random.rand(n)
    select = np.ones(n, dtype=np.int)*(n-1)
    for i, x in enumerate(p[::-1]):
        select[rands<x] = len(p)-i-1
    
    # 結果的初始化
    result = np.zeros((n,2), dtype=np.float)
    c = np.zeros(n, dtype=np.float)
    
    for i in xrange(n):
        eqidx = select[i] # 所選的函數下標
        tmp = np.dot(eq[eqidx], pos) # 進行迭代
        pos[:2] = tmp # 更新迭代向量

        # 保存結果
        result[i] = tmp
        c[i] = eqidx
        
    return result[:,0], result[:, 1], c

start = time.clock()
x, y, c = ifs([p1,p2,p3,p4],[eq1,eq2,eq3,eq4], [0,0], 100000)
print time.clock() - start
pl.figure(figsize=(6,6))
pl.subplot(121)
pl.scatter(x, y, s=1, c="g", marker="s", linewidths=0)
pl.axis("equal")
pl.axis("off")
pl.subplot(122)
pl.scatter(x, y, s=1,c = c, marker="s", linewidths=0)
pl.axis("equal")
pl.axis("off")
pl.subplots_adjust(left=0,right=1,bottom=0,top=1,wspace=0,hspace=0)
pl.gcf().patch.set_facecolor("white")
pl.show()

程序中的ifs函數是進行函數迭代的主函數，我們希望通過矩陣乘法計算函數(numpy.dot)的輸出，因此需要將乘法向量擴充為三維的：

pos = np.ones(3, dtype=np.float)
pos[:2] = init

這樣每次和迭代函數系數進行矩陣乘積運算的向量就變成了： x(n), y(n), 1.0 。

為了減少計算時間，我們不在迭代循環中計算隨機數選擇迭代方程，而是事先通過每個函數的概率，計算出函數選擇數組select，注意這里使用accumulate函數先將概率累加，然后產生一組0到1之間的隨機數，通過判斷隨機數所在的概率區間選擇不同的方程下標：

p = np.add.accumulate(p)
rands = np.random.rand(n)
select = np.ones(n, dtype=np.int)*(n-1)
for i, x in enumerate(p[::-1]):
    select[rands<x] = len(p)-i-1

最后我們通過調用scatter繪圖函數將所得到的坐標進行散列圖繪制：

pl.scatter(x, y, s=1, c="g", marker="s", linewidths=0)

其中每個關鍵字參數的含義如下：

• s : 散列點的大小，因為我們要繪制10萬點，因此大小選擇為1
• c : 點的顏色，這里選擇綠色
• marker : 點的形狀，"s"表示正方形，方形的繪制是最快的
• linewidths : 點的邊框寬度，0表示沒有邊框

此外，關鍵字參數c還可以傳入一個數組，作為每個點的顏色值，我們將計算坐標的函數下標傳入，這樣可以直觀地看出哪個點是哪個函數迭代產生的：

pl.scatter(x, y, s=1,c = c, marker="s", linewidths=0)

下圖是程序的輸出：

函數迭代系統所繪制的蕨類植物的葉子

觀察右圖的4種顏色的部分可以發現概率為1%的函數1所計算的是葉桿部分(深藍色)，概率為7%的兩個函數計算的是左右兩片子葉，而概率為85%的函數計算的是整個葉子的迭代：即最下面的三種顏色的點通過此函數的迭代產生上面的所有的深紅色的點。

我們可以看出整個葉子呈現出完美的自相似特性，任意取其中的一個子葉，將其旋轉放大之后都和整個葉子相同。

2D仿射變換

上面所介紹的4個變換方程的一般形式如下：

x(n+1) = A * x(n) + B * y(n) + C
y(n+1) = D * x(n) + E * y(n) + F

這種變換被稱為2D仿射變換，它是從2D坐標到其他2D坐標的線性映射，保留直線性和平行性。即原來是直線上的坐標，變換之后仍然成一條直線，原來是平行的直線，變換之后仍然是平行的。這種變換我們可以看作是一系列平移、縮放、翻轉和旋轉變換構成的。

為了直觀地顯示仿射變換，我們可以使用平面上的兩個三角形來表示。因為仿射變換公式中有6個未知數：A, B, C, D, E, F，而每兩個點之間的變換決定兩個方程，因此一共需要3組點來決定六個變換方程，正好是兩個三角形，如下圖所示：

兩個三角形決定一個2D仿射變換的六個參數

從紅色三角形的每個頂點變換到綠色三角形的對應頂點，正好能夠決定仿射變換中的六個參數。這樣我們可是使用N+1個三角形，決定N個仿射變換，其中的每一個變換的參數都是由第0個三角形和其它的三角形決定的。這第0個三角形我們稱之為基礎三角形，其余的三角形稱之為變換三角形。

為了繪制迭代函數系統的圖像，我們還需要給每個仿射變換方程指定一個迭代概率的參數。此參數也可以使用三角形直觀地表達出來：迭代概率和變換三角形的面積成正比。即迭代概率為變換三角形的面積除以所有變換三角形的面積之和。

如下圖所示，前面介紹的蕨類植物的分形圖案的迭代方程可以由5個三角形決定，可以很直觀地看出紫色的小三角形決定了葉子的莖；而兩個藍色的三角形決定了左右兩片子葉；綠色的三角形將莖和兩片子葉往上復制，形成整片葉子。

5個三角形的仿射方程繪制蕨類植物的葉子

迭代函數系統設計器

按照上節所介紹的三角形法，我們可以設計一個迭代函數系統的設計工具，如下圖所示：

具體的程序請參照 迭代函數系統的分形 ，這里簡單地介紹一下程序的幾個核心組成部分：

首先通過兩個三角形求解仿射方程的系數相當于求六元線性方程組的解，solve_eq函數完成這一工作，它先計算出線性方程組的矩陣a和b, 然后調用NumPy的linalg.solve對線性方程組 a*X = b 求解：

def solve_eq(triangle1, triangle2):
    """
    解方程，從triangle1變換到triangle2的變換系數
        triangle1,2是二維數組：
        x0,y0
        x1,y1
        x2,y2
    """
    x0,y0 = triangle1[0]
    x1,y1 = triangle1[1]
    x2,y2 = triangle1[2]

    a = np.zeros((6,6), dtype=np.float)
    b = triangle2.reshape(-1)
    a[0, 0:3] = x0,y0,1
    a[1, 3:6] = x0,y0,1
    a[2, 0:3] = x1,y1,1
    a[3, 3:6] = x1,y1,1
    a[4, 0:3] = x2,y2,1
    a[5, 3:6] = x2,y2,1

    c = np.linalg.solve(a, b)
    c.shape = (2,3)
    return c

triangle_area函數計算三角形的面積，它使用NumPy的cross函數計算三角形的兩個邊的向量的叉積：

def triangle_area(triangle):
    """
    計算三角形的面積
    """
    A = triangle[0]
    B = triangle[1]
    C = triangle[2]
    AB = A-B
    AC = A-C
    return np.abs(np.cross(AB,AC))/2.0

整個程序的界面使用TraitsUI庫生成，將matplotlib的Figure控件通過MPLFigureEditor和_MPLFigureEditor類嵌入到TraitsUI生成的界面中，請參考： 設計自己的Trait編輯器

IFSDesigner._figure_default創建Figure對象，并且添加兩個并排的子圖ax和ax2，ax用于三角形編輯，而ax2用于分形圖案顯示。

def _figure_default(self):
    """
    figure屬性的缺省值，直接創建一個Figure對象
    """
    figure = Figure()
    self.ax = figure.add_subplot(121)
    self.ax2 = figure.add_subplot(122)
    self.ax2.set_axis_off()
    self.ax.set_axis_off()
    figure.subplots_adjust(left=0,right=1,bottom=0,top=1,wspace=0,hspace=0)
    figure.patch.set_facecolor("w")
    return figure

IFSTriangles類完成三角形的編輯工作，其中通過如下的語句綁定Figure控件的canvas的鼠標事件

canvas = ax.figure.canvas
# 綁定canvas的鼠標事件
canvas.mpl_connect('button_press_event', self.button_press_callback)
canvas.mpl_connect('button_release_event', self.button_release_callback)
canvas.mpl_connect('motion_notify_event', self.motion_notify_callback)

由于canvas只有在真正顯示Figure時才會創建，因此不能在創建Figure控件時創建IFSTriangles對象，而需要在界面生成之后，顯示之前創建它。這里我們通過給IFSDesigner類的view屬性指定其handler為IFSHandler對象，重載Handler的init方法，此方法在界面生成之后，顯示之前被調用：

class IFSHandler(Handler):
    """
    在界面顯示之前需要初始化的內容
    """
    def init(self, info):
        info.object.init_gui_component()
        return True

然后IFSDesigner類的init_gui_component方法完成實際和canvas相關的初始工作：

def init_gui_component(self):
    self.ifs_triangle = IFSTriangles(self.ax)
    self.figure.canvas.draw()
    thread.start_new_thread( self.ifs_calculate, ())
    ...

由于通過函數迭代計算分形圖案比較費時，因此在另外一個線程中執行ifs_calculate方法進行運算，每計算ITER_COUNT個點，就調用ax2.scatter將產生的點添加進ax2中，由于隨著ax2中的點數增加，界面重繪將越來越慢，因此在draw_points函數中限制最多只調用ITER_TIMES次scatter函數。因為在別的線程中不能更新界面，因此通過調用wx.CallAfter在管理GUI的線程中調用draw_points進行界面刷新。:

def ifs_calculate(self):
    """
    在別的線程中計算
    """
    def draw_points(x, y, c):
        if len(self.ax2.collections) < ITER_TIMES:
            try:
                self.ax2.scatter(x, y, s=1, c=c, marker="s", linewidths=0)
                self.ax2.set_axis_off()
                self.ax2.axis("equal")
                self.figure.canvas.draw()
            except:
                pass

    def clear_points():
        self.ax2.clear()

    while 1:
        try:
            if self.exit == True:
                break
            if self.clear == True:
                self.clear = False
                self.initpos = [0, 0]
                x, y, c = ifs( self.ifs_triangle.get_areas(),
                    self.ifs_triangle.get_eqs(), self.initpos, 100)
                self.initpos = [x[-1], y[-1]]
                self.ax2.clear()

            x, y, c = ifs( self.ifs_triangle.get_areas(),
                self.ifs_triangle.get_eqs(), self.initpos, ITER_COUNT)
            if np.max(np.abs(x)) < 1000000 and np.max(np.abs(y)) < 1000000:
                self.initpos = [x[-1], y[-1]]
                wx.CallAfter( draw_points, x, y, c )
            time.sleep(0.05)
        except:
            pass

用戶修改三角形之后，需要重新迭代，并繪制分形圖案，三角形的改變通過 IFSTriangles.version 屬性通知給IFSDesigner，在IFSTriangles中，三角形改變之后，將運行：

self.version += 1

在IFSDesigner中監聽version屬性的變化：

@on_trait_change("ifs_triangle.version")
def on_ifs_version_changed(self):
    """
    當三角形更新時，重新繪制所有的迭代點
    """
    self.clear = True

當IFSDesigner.clear為True時，真正進行迭代運算的ifs_calculate方法就知道需要重新計算了。

L-System分形

前面所繪制的分形圖案都是都是使用數學函數的迭代產生，而L-System分形則是采用符號的遞歸迭代產生。首先如下定義幾個有含義的符號：

• F : 向前走固定單位
• + : 正方向旋轉固定單位
• - : 負方向旋轉固定單位

使用這三個符號我們很容易描述下圖中由4條線段構成的圖案：

F+F--F+F

如果將此符號串中的所有F都替換為F+F--F+F，就能得到如下的新字符串：

F+F--F+F+F+F--F+F--F+F--F+F+F+F--F+F

如此替換迭代下去，并根據字串進行繪圖(符號+和-分別正負旋轉60度)，可得到如下的分形圖案：

使用F+F--F+F迭代的分形圖案

除了 F, +, - 之外我們再定義如下幾個符號：

• f : 向前走固定單位，為了定義不同的迭代公式
• [ : 將當前的位置入堆棧
• ] : 從堆棧中讀取坐標，修改當前位置
• S : 初始迭代符號

所有的符號(包括上面未定義的)都可以用來定義迭代，通過引入兩個方括號符號，使得我們能夠描述分岔的圖案。

例如下面的符號迭代能夠繪制出一棵植物：

S -> X
X -> F-[[X]+X]+F[+FX]-X
F -> FF

我們用一個字典定義所有的迭代公式和其它的一些繪圖信息：

{
    "X":"F-[[X]+X]+F[+FX]-X", "F":"FF", "S":"X",
    "direct":-45,
    "angle":25,
    "iter":6,
    "title":"Plant"
}

其中：

• direct : 是繪圖的初始角度，通過指定不同的值可以旋轉整個圖案
• angle : 定義符號+,-旋轉時的角度，不同的值能產生完全不同的圖案
• iter : 迭代次數

下面的程序將上述字典轉換為需要繪制的線段坐標：

class L_System(object):
    def __init__(self, rule):
        info = rule['S']
        for i in range(rule['iter']):
            ninfo = []
            for c in info:
                if c in rule:
                    ninfo.append(rule[c])
                else:
                    ninfo.append(c)
            info = "".join(ninfo)
        self.rule = rule
        self.info = info

    def get_lines(self):
        d = self.rule['direct']
        a = self.rule['angle']
        p = (0.0, 0.0)
        l = 1.0
        lines = []
        stack = []
        for c in self.info:
            if c in "Ff":
                r = d * pi / 180
                t = p[0] + l*cos(r), p[1] + l*sin(r)
                lines.append(((p[0], p[1]), (t[0], t[1])))
                p = t
            elif c == "+":
                d += a
            elif c == "-":
                d -= a
            elif c == "[":
                stack.append((p,d))
            elif c == "]":
                p, d = stack[-1]
                del stack[-1]
        return lines

我們使用matplotlib的LineCollection繪制所有的直線：

import matplotlib.pyplot as pl
from matplotlib import collections

# rule = {...}  此處省略rule的定義

lines = L_System(rule).get_lines()

fig = pl.figure()
ax = fig.add_subplot(111)
linecollections = collections.LineCollection(lines)
ax.add_collection(linecollections, autolim=True)
pl.show()

下面是幾種L-System的分形圖案，繪制此圖的完整程序請參照 繪制L-System的分形圖 。

幾種L-System的迭代圖案