用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/'''
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))))
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))
img = hsv_to_rgb(HSV)
Image.fromarray(np.uint8(img*255)).save('kaleidoscope.png')
if __name__ == '__main__':
import time
start = 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
@jit('complex64(complex64)', nopython=True)
def f(z):
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 time
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
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 = rhombi
polygon = 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 = point
sphere = 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_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) 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)
原文链接：
https://github.com/neozhaoliang/pywonderland/blob/master/README.md