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里的每一个元素相乘，并返回相同形状