如何用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/ | |
''' | |
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 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__': | |
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 | |
# 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 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
本文转载自大数据文摘
1. 时间的方向 || Stephen W. Hawking
6. 数学是什么?
10. 漫谈经济物理学
END
更多精彩文章,请关注微信号:cpsjournals