原标题：如何用Python画各种著名数学图案 | 附图+代码
用Python绘制著名的数学图片或动画，展示数学中的算法魅力。
Mandelbrot 集
代码：46 lines (34 sloc) 1.01 KB
”’
A fast Mandelbrot set wallpaper renderer
reddit discussion: https://www.reddit.com/r/math/comments/2abwyt/smooth_colour_mandelbrot/
”’
importnumpy asnp
fromPILimportImage
fromnumba importjit
MAXITERS=200
RADIUS=100
@jit
defcolor(z, i):
v =np.log2(i +1-np.log2(np.log2(abs(z)))) /5
ifv <1.0:
returnv**4, v**2.5, v
else:
v =max(0, 2-v)
returnv, v**1.5, v**3
@jit
defiterate(c):
z =0j
fori inrange(MAXITERS):
ifz.real*z.real +z.imag*z.imag >RADIUS:
returncolor(z, i)
z =z*z +c
return0, 0,0
defmain(xmin, xmax, ymin, ymax, width, height):
x =np.linspace(xmin, xmax, width)
y =np.linspace(ymax, ymin, height)
z =x[None, :] +y[:, None]*1j
red, green, blue =np.asarray(np.frompyfunc(iterate, 1, 3)(z)).astype(np.float)
img =np.dstack((red, green, blue))
Image.fromarray(np.uint8(img*255)).save(‘mandelbrot.png’)
if__name__==’__main__’:
main(-2.1, 0.8, -1.16, 1.16, 1200, 960)
多米诺洗牌算法
代码链接：https://github.com/neozhaoliang/pywonderland/tree/master/src/domino
正二十面体万花筒
代码：53 lines (40 sloc) 1.24 KB
”’
A kaleidoscope pattern with icosahedral symmetry.
”’
importnumpy asnp
fromPILimportImage
frommatplotlib.colors importhsv_to_rgb
defKlein(z):
”’Klein’s j-function”’
return1728*(z *(z**10+11*z**5-1))**5/
(-(z**20+1) +228*(z**15-z**5) -494*z**10)**3
defRiemannSphere(z):
”’
map the complex plane to Riemann’s sphere via stereographic projection
”’
t =1+z.real*z.real +z.imag*z.imag
return2*z.real/t, 2*z.imag/t, 2/t-1
defMobius(z):
”’
distort the result image by a mobius transformation
”’
return(z -20)/(3*z +1j)
defmain(imgsize):
x =np.linspace(-6, 6, imgsize)
y =np.linspace(6, -6, imgsize)
z =x[None, :] +y[:, None]*1j
z =RiemannSphere(Klein(Mobius(Klein(z))))
#define colors in hsv space
H =np.sin(z[0]*np.pi)**2
S =np.cos(z[1]*np.pi)**2
V =abs(np.sin(z[2]*np.pi) *np.cos(z[2]*np.pi))**0.2
HSV=np.dstack((H, S, V))
#transform to rgb space
img =hsv_to_rgb(HSV)
Image.fromarray(np.uint8(img*255)).save(‘kaleidoscope.png’)
if__name__==’__main__’:
importtime
start =time.time()
main(imgsize=800)
end =time.time()
print(‘runtime: {:3f}seconds’.format(end -start))
Newton 迭代分形
代码：46 lines (35 sloc) 1.05 KB
importnumpy asnp
importmatplotlib.pyplot asplt
fromnumba importjit
#define functions manually, do not use numpy’s poly1d funciton!
@jit(‘complex64(complex64)’, nopython=True)
deff(z):
#z*z*z is faster than z**3
returnz*z*z -1
@jit(‘complex64(complex64)’, nopython=True)
defdf(z):
return3*z*z
@jit(‘float64(complex64)’, nopython=True)
defiterate(z):
num =0
whileabs(f(z)) >1e-4:
w =z -f(z)/df(z)
num +=np.exp(-1/abs(w-z))
z =w
returnnum
defrender(imgsize):
x =np.linspace(-1, 1, imgsize)
y =np.linspace(1, -1, imgsize)
z =x[None, :] +y[:, None] *1j
img =np.frompyfunc(iterate, 1, 1)(z).astype(np.float)
fig =plt.figure(figsize=(imgsize/100.0, imgsize/100.0), dpi=100)
ax =fig.add_axes([0, 0, 1, 1], aspect=1)
ax.axis(‘off’)
ax.imshow(img, cmap=’hot’)
fig.savefig(‘newton.png’)
if__name__==’__main__’:
importtime
start =time.time()
render(imgsize=400)
end =time.time()
print(‘runtime: {:03f}seconds’.format(end -start))
李代数E8 的根系
代码链接：https://github.com/neozhaoliang/pywonderland/blob/master/src/misc/e8.py
模群的基本域
代码链接：
https://github.com/neozhaoliang/pywonderland/blob/master/src/misc/modulargroup.py
彭罗斯铺砌
代码链接：
https://github.com/neozhaoliang/pywonderland/blob/master/src/misc/penrose.py
Wilson 算法
代码链接：https://github.com/neozhaoliang/pywonderland/tree/master/src/wilson
反应扩散方程模拟
代码链接：https://github.com/neozhaoliang/pywonderland/tree/master/src/grayscott
120 胞腔
代码：69 lines (48 sloc) 2.18 KB
#pylint: disable=unused-import
#pylint: disable=undefined-variable
fromitertools importcombinations, product
importnumpy asnp
fromvapory import*
classPenrose(object):
GRIDS=[np.exp(2j*np.pi *i /5) fori inrange(5)]
def__init__(self, num_lines, shift, thin_color, fat_color, **config):
self.num_lines =num_lines
self.shift =shift
self.thin_color =thin_color
self.fat_color =fat_color
self.objs =self.compute_pov_objs(**config)
defcompute_pov_objs(self, **config):
objects_pool =[]
forrhombi, color inself.tile():
p1, p2, p3, p4 =rhombi
polygon =Polygon(5, p1, p2, p3, p4, p1,
Texture(Pigment(‘color’, color), config[‘default’]))
objects_pool.append(polygon)
forp, q inzip(rhombi, [p2, p3, p4, p1]):
cylinder =Cylinder(p, q, config[‘edge_thickness’], config[‘edge_texture’])
objects_pool.append(cylinder)
forpoint inrhombi:
x, y =point
sphere =Sphere((x, y, 0), config[‘vertex_size’], config[‘vertex_texture’])
objects_pool.append(sphere)
returnObject(Union(*objects_pool))
defrhombus(self, r, s, kr, ks):
if(s -r)**2%5==1:
color =self.thin_color
else:
color =self.fat_color
point =(Penrose.GRIDS[r] *(ks -self.shift[s])
-Penrose.GRIDS[s] *(kr -self.shift[r])) *1j/Penrose.GRIDS[s-r].imag
index =[np.ceil((point/grid).real +shift)
forgrid, shift inzip(Penrose.GRIDS, self.shift)]
vertices =[]
forindex[r], index[s] in[(kr, ks), (kr+1, ks), (kr+1, ks+1), (kr, ks+1)]:
vertices.append(np.dot(index, Penrose.GRIDS))
vertices_real =[(z.real, z.imag) forz invertices]
returnvertices_real, color
deftile(self):
forr, s incombinations(range(5), 2):
forkr, ks inproduct(range(-self.num_lines, self.num_lines+1), repeat=2):
yieldself.rhombus(r, s, kr, ks)
defput_objs(self, *args):
returnObject(self.objs, *args)
1.
免责声明：本文系网络转载，版权归原作者所有。如涉及作品版权问题，请与我们联系，我们将根据您提供的版权证明材料确认版权并支付稿酬或者删除内容。返回搜狐，查看更多
责任编辑：