如何用Python画各种著名数学图案 | 附图+代码

xiaoxiao2023-12-19 163

代码：46 lines (34 sloc) 1.01 KB

'''A fast Mandelbrot set wallpaper rendererreddit discussion: https://www.reddit.com/r/math/comments/2abwyt/smooth_colour_mandelbrot/''' import numpy as np from PIL import Image from numba import jit MAXITERS = 200 RADIUS = 100 @jit def color(z, i):v = np.log2(i + 1 - np.log2(np.log2(abs(z)))) / 5 if v < 1.0: return v**4, v**2.5, v else:v = max(0, 2-v) return v, v**1.5, v**3 @jit def iterate(c):z = 0j for i in range(MAXITERS): if z.real*z.real + z.imag*z.imag > RADIUS: return color(z, i)z = z*z + c return 0, 0 ,0 def main(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.''' import numpy as np from PIL import Image from matplotlib.colors import hsv_to_rgb def Klein(z):'''Klein's j-function''' return 1728 * (z * (z**10 + 11 * z**5 - 1))**5 / \(-(z**20 + 1) + 228 * (z**15 - z**5) - 494 * z**10)**3 def RiemannSphere(z):''' map the complex plane to Riemann's sphere via stereographic projection '''t = 1 + z.real*z.real + z.imag*z.imag return 2*z.real/t, 2*z.imag/t, 2/t-1 def Mobius(z):''' distort the result image by a mobius transformation ''' return (z - 20)/(3*z + 1j) def main(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 spaceH = 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 spaceimg = hsv_to_rgb(HSV)Image.fromarray(np.uint8(img*255)).save('kaleidoscope.png') if __name__ == '__main__': import timestart = time.time()main(imgsize=800)end = time.time() print('runtime: {:3f} seconds'.format(end - start))

Newton 迭代分形

代码：46 lines (35 sloc) 1.05 KB

import numpy as np import matplotlib.pyplot as plt from numba import jit# define functions manually, do not use numpy's poly1d funciton! @jit('complex64(complex64)', nopython=True) def f(z):# z*z*z is faster than z**3 return z*z*z - 1 @jit('complex64(complex64)', nopython=True) def df(z): return 3*z*z @jit('float64(complex64)', nopython=True) def iterate(z):num = 0 while abs(f(z)) > 1e-4:w = z - f(z)/df(z)num += np.exp(-1/abs(w-z))z = w return num def render(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__': import timestart = 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 from itertools import combinations, product import numpy as np from vapory import * class Penrose(object): GRIDS = [np.exp(2j * np.pi * i / 5) for i in range(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) def compute_pov_objs(self, **config):objects_pool = [] for rhombi, color in self.tile():p1, p2, p3, p4 = rhombipolygon = Polygon(5, p1, p2, p3, p4, p1,Texture(Pigment('color', color), config['default']))objects_pool.append(polygon) for p, q in zip(rhombi, [p2, p3, p4, p1]):cylinder = Cylinder(p, q, config['edge_thickness'], config['edge_texture'])objects_pool.append(cylinder) for point in rhombi:x, y = pointsphere = Sphere((x, y, 0), config['vertex_size'], config['vertex_texture'])objects_pool.append(sphere) return Object(Union(*objects_pool)) def rhombus(self, r, s, kr, ks): if (s - r)**2 % 5 == 1:color = self.thin_color else:color = self.fat_colorpoint = (Penrose.GRIDS[r] * (ks - self.shift[s]) - Penrose.GRIDS[s] * (kr - self.shift[r])) *1j / Penrose.GRIDS[s-r].imagindex = [np.ceil((point/grid).real + shift) for grid, shift in zip(Penrose.GRIDS, self.shift)]vertices = [] for index[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) for z in vertices] return vertices_real, color def tile(self): for r, s in combinations(range(5), 2): for kr, ks in product(range(-self.num_lines, self.num_lines+1), repeat=2): yield self.rhombus(r, s, kr, ks) def put_objs(self, *args): return Object(self.objs, *args)

原文发布时间为：2017-04-15

本文来自云栖社区合作伙伴“大数据文摘”，了解相关信息可以关注“BigDataDigest”微信公众号