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

大数据文摘作品,转载具体要求见文末

编译团队:Aileen,徐凌霄

用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)

后台回复“数学”查看完整代码集哦

一人一笔 | 数据团队建设“全景报告”

清华数据科学研究院联合大数据文摘,发起一次数据团队全行业调研。本次调研将对国内外数据团队发展现状进行盘点和趋势预测,同时探索数据团队应如何建设。我们将结合一系列专访与调查问卷内容,在7月初发布《数据团队建设全景报告》。

如果你是数据团队的一员、和数据团队一起工作,或者希望了解其他数据团队的发展现状和未来。

那么恳请你花费5分钟时间,点击“阅读原文”填写问卷,帮助我们完成这次调研。

原文链接:https://github.com/neozhaoliang/pywonderland/blob/master/README.md

原文发布于微信公众号 - 大数据文摘(BigDataDigest)

原文发表时间:2017-04-15

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏机器之心

入门 | 自然语言处理是如何工作的?一步步教你构建 NLP 流水线

计算机非常擅长使用结构化数据,例如电子表格和数据库表。但是我们人类通常用文字交流,而不是使用电子表格来交流。这对计算机来说不是一件好事。

1093
来自专栏机器学习原理

NLP(6)——命名实体识别

普通的工具如hanlp,htp,不能识别特定领域的专有名词,所以需要实体识别的算法。下面就以医疗专业为例子来谈一下医疗专业的命名实体识别。

2963
来自专栏深度学习自然语言处理

这么好的视频不看吗?深度学习和线代,微积分

大家盼望的中秋节和十一已经基本都要过去了,大家是不是都玩的挺开心呀?(哎,我可没0.0,基本没离开过实验室,别认为我在学习

2312
来自专栏java一日一条

为什么这段代码输出的是”Hello World”

明明是在程序里使用了java.util.Ramdom()函数产生随机数,为什么每次打出的结果都是Hello world? 各位程序员,你们怎么看?请务必要独立思...

642
来自专栏机器学习之旅

动态最优化经典面试题

最近看到了一条史前的算法面试题,觉得挺有意思的,虽然网上已经有了很多完善的答案,但是我还是想自己整理一遍,强化印象,同时也和大家分享一下这道12年的Google...

952
来自专栏机器之心

资源 | 2017年最流行的15个数据科学Python库

选自Medium 作者:Igor Bobriakov 机器之心编译 参与:朱朝阳、吴攀 Python 近几年在数据科学行业获得了人们的极大青睐,各种资源也层出不...

2914
来自专栏牛客网

2018年5月份找实习经历(计算机视觉与深度学习岗)

2043
来自专栏大数据文摘

人类对随机数的探索:如何才能生成一个均匀的随机数列

2597
来自专栏java一日一条

为什么这段代码输出的是”Hello World”

明明是在程序里使用了java.util.Ramdom()函数产生随机数,为什么每次打出的结果都是Hello world? 各位程序员,你们怎么看?请务必要独立思...

982
来自专栏大数据钻研

大数据入门之路 献给迷茫的你

假如你想成为一个数据科学家,或者已经是数据科学家的你想扩展你的技能,那么你已经来对地方了。本文的目的就是给数据分析方面的Python新手提供一个完整的学习路径。...

3304

扫码关注云+社区

领取腾讯云代金券