Python 繪製分形圖

02-08

1. 曼德勃羅集

import numpy as npimport pylab as plimport timefrom matplotlib import cmdef 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.jet, extent=[x0,x1,y0,y1]) pl.gca().set_axis_off() x,y = 0.27322626, 0.595153338pl.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()

2. 分形樹葉
import numpy as npimport matplotlib.pyplot as plimport time# 蕨類植物葉子的迭代函數和其概率值eq1 = np.array([[0,0,0],[0,0.16,0]])p1 = 0.01eq2 = np.array([[0.2,-0.26,0],[0.23,0.22,1.6]])p2 = 0.07eq3 = np.array([[-0.15, 0.28, 0],[0.26,0.24,0.44]])p3 = 0.07eq4 = np.array([[0.85, 0.04, 0],[-0.04, 0.85, 1.6]])p4 = 0.85def 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 range(n): eqidx = select[i] # 所選的函數下標 tmp = np.dot(eq[eqidx], pos) # 進行迭代 pos[:2] = tmp # 更新迭代向量 # 保存結果 result[i] = tmp c[i] = eqidx return result[:,0], result[:, 1], cstart = time.clock()x, y, c = ifs([p1,p2,p3,p4],[eq1,eq2,eq3,eq4], [0,0], 100000)time.clock() - startpl.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("#D3D3D3")pl.show()

3. 其它分形圖(科赫曲線、分形龍、謝爾賓斯基三角等)

from math import sin, cos, piimport matplotlib.pyplot as plfrom matplotlib import collectionsclass 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 linesrules = [ { "F":"F+F--F+F", "S":"F", "direct":180, "angle":60, "iter":5, "title":"Koch" }, { "X":"X+YF+", "Y":"-FX-Y", "S":"FX", "direct":0, "angle":90, "iter":13, "title":"Dragon" }, { "f":"F-f-F", "F":"f+F+f", "S":"f", "direct":0, "angle":60, "iter":7, "title":"Triangle" }, { "X":"F-[[X]+X]+F[+FX]-X", "F":"FF", "S":"X", "direct":-45, "angle":25, "iter":6, "title":"Plant" }, { "S":"X", "X":"-YF+XFX+FY-", "Y":"+XF-YFY-FX+", "direct":0, "angle":90, "iter":6, "title":"Hilbert" }, { "S":"L--F--L--F", "L":"+R-F-R+", "R":"-L+F+L-", "direct":0, "angle":45, "iter":10, "title":"Sierpinski" }, ]def draw(ax, rule, iter=None): if iter!=None: rule["iter"] = iter lines = L_System(rule).get_lines() linecollections = collections.LineCollection(lines) ax.add_collection(linecollections, autolim=True) ax.axis("equal") ax.set_axis_off() ax.set_xlim(ax.dataLim.xmin, ax.dataLim.xmax) ax.invert_yaxis() fig = pl.figure(figsize=(7,4.5))fig.patch.set_facecolor("papayawhip")for i in xrange(6): ax = fig.add_subplot(231+i) draw(ax, rules[i])fig.subplots_adjust(left=0,right=1,bottom=0,top=1,wspace=0,hspace=0)pl.show()