Python 繪製分形圖（曼德勃羅集、分形樹葉、科赫曲線、分形龍、謝爾賓斯基三角等）附代碼

1. 曼德勃羅集

import numpy as npnimport pylab as plnimport timenfrom matplotlib import cmnndef iter_point(c):n z = cn for i in xrange(1, 100): # 最多迭代100次n if abs(z)>2: break # 半徑大於2則認為逃逸n z = z*z+cn return i # 返回迭代次數nndef draw_mandelbrot(cx, cy, d):n """n 繪製點(cx, cy)附近正負d的範圍的Mandelbrotn """n x0, x1, y0, y1 = cx-d, cx+d, cy-d, cy+d n y, x = np.ogrid[y0:y1:200j, x0:x1:200j]n c = x + y*1jn start = time.clock()n mandelbrot = np.frompyfunc(iter_point,1,1)(c).astype(np.float)n print "time=",time.clock() - startn pl.imshow(mandelbrot, cmap=cm.jet, extent=[x0,x1,y0,y1])n pl.gca().set_axis_off()nnx,y = 0.27322626, 0.595153338nnpl.subplot(231)ndraw_mandelbrot(-0.5,0,1.5)nfor i in range(2,7): n pl.subplot(230+i)n draw_mandelbrot(x, y, 0.2**(i-1))npl.subplots_adjust(0.02, 0, 0.98, 1, 0.02, 0)npl.show()n

2. 分形樹葉

import numpy as npnimport matplotlib.pyplot as plnimport timenn# 蕨類植物葉子的迭代函數和其概率值neq1 = np.array([[0,0,0],[0,0.16,0]])np1 = 0.01nneq2 = np.array([[0.2,-0.26,0],[0.23,0.22,1.6]])np2 = 0.07nneq3 = np.array([[-0.15, 0.28, 0],[0.26,0.24,0.44]])np3 = 0.07nneq4 = np.array([[0.85, 0.04, 0],[-0.04, 0.85, 1.6]])np4 = 0.85nndef ifs(p, eq, init, n):n"""n 進行函數迭代n p: 每個函數的選擇概率列表n eq: 迭代函數列表n init: 迭代初始點n n: 迭代次數nn 返回值： 每次迭代所得的X坐標數組， Y坐標數組， 計算所用的函數下標 n """nn# 迭代向量的初始化npos = np.ones(3, dtype=np.float)npos[:2] = initnn# 通過函數概率，計算函數的選擇序列np = np.add.accumulate(p) nrands = np.random.rand(n)nselect = np.ones(n, dtype=np.int)*(n-1)nfor i, x in enumerate(p[::-1]):nselect[rands<x] = len(p)-i-1nn# 結果的初始化nresult = np.zeros((n,2), dtype=np.float)nc = np.zeros(n, dtype=np.float)nnfor i in range(n):neqidx = select[i] # 所選的函數下標ntmp = np.dot(eq[eqidx], pos) # 進行迭代npos[:2] = tmp # 更新迭代向量nn# 保存結果nresult[i] = tmpnc[i] = eqidxnnreturn result[:,0], result[:, 1], cnnstart = time.clock()nx, y, c = ifs([p1,p2,p3,p4],[eq1,eq2,eq3,eq4], [0,0], 100000)ntime.clock() - startnpl.figure(figsize=(6,6))npl.subplot(121)npl.scatter(x, y, s=1, c="g", marker="s", linewidths=0)npl.axis("equal")npl.axis("off")npl.subplot(122)npl.scatter(x, y, s=1,c = c, marker="s", linewidths=0)npl.axis("equal")npl.axis("off")npl.subplots_adjust(left=0,right=1,bottom=0,top=1,wspace=0,hspace=0)npl.gcf().patch.set_facecolor("#D3D3D3")npl.show()n

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

from math import sin, cos, pinimport matplotlib.pyplot as plnfrom matplotlib import collectionsnnclass L_System(object):ndef __init__(self, rule):ninfo = rule[S]nfor i in range(rule[iter]):nninfo = []nfor c in info:nif c in rule:nninfo.append(rule[c])nelse:nninfo.append(c)ninfo = "".join(ninfo)nself.rule = rulenself.info = infonndef get_lines(self):nd = self.rule[direct]na = self.rule[angle]np = (0.0, 0.0)nl = 1.0nlines = []nstack = []nfor c in self.info:nif c in "Ff":nr = d * pi / 180nt = p[0] + l*cos(r), p[1] + l*sin(r)nlines.append(((p[0], p[1]), (t[0], t[1])))np = tnelif c == "+":nd += anelif c == "-":nd -= anelif c == "[":nstack.append((p,d))nelif c == "]":np, d = stack[-1]ndel stack[-1]nreturn linesnnrules = [n{n"F":"F+F--F+F", "S":"F",n"direct":180,n"angle":60,n"iter":5,n"title":"Koch"n},n{n"X":"X+YF+", "Y":"-FX-Y", "S":"FX",n"direct":0,n"angle":90,n"iter":13,n"title":"Dragon"n},n{n"f":"F-f-F", "F":"f+F+f", "S":"f",n"direct":0,n"angle":60,n"iter":7,n"title":"Triangle"n},n{n"X":"F-[[X]+X]+F[+FX]-X", "F":"FF", "S":"X",n"direct":-45,n"angle":25,n"iter":6,n"title":"Plant"n},n{n"S":"X", "X":"-YF+XFX+", "Y":"+XF-YFY-FX+",n"direct":0,n"angle":90,n"iter":6,n"title":"Hilbert"n},n{n"S":"L--F--L--F", "L":"+R-F-R+", "R":"-L+F+",n"direct":0,n"angle":45,n"iter":10,n"title":"Sierpinski"n},nn]nndef draw(ax, rule, iter=None):nif iter!=None:nrule["iter"] = iternlines = L_System(rule).get_lines()nlinecollections = collections.LineCollection(lines)nax.add_collection(linecollections, autolim=True)nax.axis("equal")nax.set_axis_off()nax.set_xlim(ax.dataLim.xmin, ax.dataLim.xmax)nax.invert_yaxis()nnfig = pl.figure(figsize=(7,4.5))nfig.patch.set_facecolor("papayawhip")nnfor i in xrange(6):nax = fig.add_subplot(231+i)ndraw(ax, rules[i])nnfig.subplots_adjust(left=0,right=1,bottom=0,top=1,wspace=0,hspace=0)npl.show()n

------------------

作者：霍開拓

博客專欄：霍開拓

知乎專欄：繁星悟語(Physics, Math & Philosophy)

大家也可以加小編微信：tszhihu （備註：Python），拉大家到 Python愛好者社區 微信群，可以跟各位老師互相交流。謝謝。第一時間獲取視頻更新動態。

也可以關注官網微信公眾號：Python愛好者社區