分形几何python代码_Python, Cython绘制美妙绝伦的Mandelbrot集, 曼德博集分形图案

上世纪60-70年代，美籍数学家曼德博 - Benoit B. Mandelbrot几乎单枪匹马的创立了一个新的数学分支，即分形几何学 - fractal geometry。这个新的数学分支有助于人类探索物理现象背后的数学规律，分形混沌之旋风，横扫数学、理化、生物、大气、海洋以至社会学科，在音乐、美术领域也产生了一定的影响。

分形艺术 - fractal art不同于普通的电脑绘画，它利用分形几何学和计算机强大的运算能力，将数学公式反复迭代运算，再结合作者的审美及艺术性的塑造，从而将抽象神秘的数学公式变成一幅幅精美绝伦的艺术画作。

本文节选自作者的《Python编程基础及应用》视频教程。想完整零基础学习Python程序设计，欢迎使用此免费视频教程。Python编程基础及应用_哔哩哔哩 (゜-゜)つロ 干杯~-bilibili​www.bilibili.com

1. Mandelbrot集

1.1 定义

​ 说明：阅读Mandelbrot集相关内容需要一点复数及复平面的知识，缺乏相关背景知识的读者可以略过，不影响后续内容的理解。

​ Mandelbrot(曼德布洛特)集是在复平面上构成分形图案的点的集合。它可以用下述复二次多项式定义：

这里的c = x + yi是一个复数，如果把实部x和虚部y视为复平面上的横纵坐标的话，那么c对应该复平面上的一个点。从z=0开始，反复应用上述二次多项式进行迭代计算，z值将不断变化，其或者延伸至无限大，或者停留在有限半径的圆内部。如果是后者，我们称参数c导致了不发散的迭代z序列，此时，c点属于Mandelbrot集合，否则不属于。

1.2 计算

​ 计算机显然没有办法进行无穷尽的计算才确认c所导致的迭代z值是否发散，所以我们简单地计算c点的逃逸时间 - 也就是经过多少轮迭代计算，z值的模超过2。在绘图过程中，该逃逸时间将作为c点在平面上的颜色。

def getEscapeTime(c):

"计算参数c的逃逸时间，该逃逸时间将用作点的颜色"

z = 0

for i in range(100):

if abs(z) > 2:

return i

z = z*z + c

return i

​ 如代码所示，abs(z)用于计算复数的模，在该代码中，逃逸时间最大为100-1 = 99。同时，应注意到，Python原生支持复数运算，z*z + c就是在复数之间的运算，其代码与普通实数间的运算没有什么不同。

​ 接下来，就可以计算Mandelbrot集合了。广义的Mandelbrot集合是个无穷大的复平面，显然，我们只可能画其中一小块。下述代码计算以(xCenter,yCenter)为中心点，semiWidth为半宽的正方形区域内的Mandelbrot集。这个正方形区域x坐标范围为xCenter±semiWidth，y坐标范围为yCenter±semiWidth。考虑到计算机只能计算离散点，所以该正方形区域将被分成N*N个点进行计算。当N=600时，下述函数的结果矩阵就是一个600x600的图像，矩形元素值代表该点的颜色。

import numpy as np

def computeMandelbrot(xCenter,yCenter,semiWidth,N):

xFrom,xTo,yFrom,yTo = xCenter-semiWidth,xCenter+semiWidth,\

yCenter-semiWidth,yCenter+semiWidth

y,x = np.ogrid[yFrom:yTo:N*1j,xFrom:xTo:N*1j]

c = x + y*1j

print("c.shape:",c.shape,"x.shape:",x.shape,"y.shape:",y.shape)

return np.frompyfunc(getEscapeTime,1,1)(c).astype(np.float)

​ 上述代码的理解并不容易，需要先从多维数据np.ndarray的广播运算开始。对于复杂的问题，做一点小规模的模拟运算可以帮助理解。

import numpy as np

y,x = np.ogrid[1:4:4j,1:3:3j]

print("y=")

print(y)

print("x=",x)

print("shape of y:",y.shape,"shape of x:",x.shape)

执行如果如下：

y=

[[1.]

[2.]

[3.]

[4.]]

x= [[1. 2. 3.]]

shape of y: (4, 1) shape of x: (1, 3)

​ np.ogrid是一个特殊对象，它使用切片元组作为下标，专门用于生成用来广播的数组。1:4:4j有点类似于np.linspace(1,4,4)，用于生成从1到4(含4)的有4个元素的等差数列，这里的j并不表示复数，只是一种表达形式。可以看到，ogrid生成的y是一个4行1列的二维数组，x则是一个1行3列的二维数组。其值见上。

z = y*1j

print("z = y*1j =",z)

执行结果：

z = y*1j = [[0.+1.j]

[0.+2.j]

[0.+3.j]

[0.+4.j]]

​ y*1j是把复数1j同y里的每一个元素相乘，并返回相同形状(维度)的数组。由于numpy模块并不是Python的原生内容，所以这里的*号并不是普通意义的Python乘法，它事实上是一个被重载的操作符函数，这个函数位于numpy模块，是用C语言书写的。可见，y*1j的结果元素的实部全是0，虚部则等于原元素值乘以1j。<