mandelbrot集合代码opencv c语言_来画分形吧！用30行代码画个Mandelbrot集合

这是什么炫技吗？完全没有，因为这是只用最基本的砖块就能达成的事情，我只是以自己写的代码为例复述一遍。

先来重复介绍一遍什么是Mandelbrot集合。简单来说，这是把实数域上的二次迭代搬到复数上形成的。考虑如下形式的映射：

这个复数“数列”的性质如何？

先限定到实数域，很容易观察到，c = 0时数列为不动点，c = 1时数列则快速发散。先不谈收敛到什么值，什么情况下这个数列是有界的？容易证明，这需要初值满足 -2 ≤ c ≤ 1/4。那么对于整个复平面如何呢？是不是会是一个较为均匀分布的图形，例如一个椭圆？

完全不是这样——正如容易猜到的，上面的Mandelbrot集合就是一个在复平面上标记出使数列不向无穷远发散的初值 c 点组成的集合，它看起来一点也不规则。(好消息是，虽然图上由于分辨率的原因看起来像有什么孤岛，但它的确是个连通集。)图中心脏形状的大区域对应于实部-3/4 ≤ Re(c) ≤ 1/4，而略小一号的圆盘则是以-1为圆心，1/4为半径的圆。

比起略显单调的雪花曲线，Mandelbrot集合可能是一类被称为“分形”的物体中最为熟知的代表。虽然它没有那么严格的自相似性，却增加了多得多的美感。下面是一段放大该集合局部的展示，可以看到它是如何改变了自己，又是如何与原来一样。

(为什么这图片有颜色？稍后解释。)

介绍完毕，来贴上代码：

import turtleimport numpy as npdef mandelbrot(x, y, max_repeat):    z = c = complex(x, y)    for i in range(max_repeat):        z = z ** 2 + c        if abs(z) > 2:            return False    return Trueturtle.setworldcoordinates(-2, -1.5, 1, 1.5)t = turtle.Turtle()t.hideturtle()turtle.tracer(0)nx, ny = 750, 750max_repeat = 100x = np.linspace(-2, 1, nx, endpoint=False)y = np.linspace(-1.5, 1.5, ny, endpoint=False)xv, yv = np.meshgrid(x, y, indexing='ij')for i in range(nx):    for j in range(ny):        t.penup()        t.goto(xv[i, j], yv[i, j])        if mandelbrot(xv[i, j], yv[i, j], max_repeat):            t.pendown()            t.dot(size=2)    turtle.update()turtle.done()

比起具体的代码，我认为，更重要的是具体思路的讲解。

首先，如何“画”一个函数图像？我们当然有matplotlib等高级的作图库，但其实只用Python自带的库已经能解决不少问题了，比如这个问题通常就是用turtle库解决。

(这里我想起高中时老师的一个说法：一个好律师未必能熟背所有的常用法规，但遇到一个问题时，ta一定知道在什么样的地方能找到什么样的条款，即使会找到什么样的条款还不清楚。当然，如果读者中有学法的，那见笑了。)

然后，如何画一个“函数图像”？注意到这里画的是一个连通的集合而不是一条曲线，因此描点法严格来说是不可行的。然而电脑的分辨率有自己的限制，我们可以利用这点，用turtle库的dot()函数“糊弄”过去——在这个集合可能的范围里采样选点，然后如果判定为在集合内就点一个黑色的点，否则就移动到下一个点。

总结一下目前的思路：在必要的模式化代码之外，需要写一个循环体，其中的顺序流程是选点–移动画笔–判定–画图；由于移动和画图都有现成的代码，需要解决的是选点和判定。

(我曾经以为这是很自然而然的事情，但后来才意识到，不是所有人都能跨过编程思维的门槛的。会有人问：可是我还不知道怎么判定，怎么就把循环体的代码写好了？于是因为不会写判定的代码而纠结半天。实际上，编程需要一种模块化的思维，或者用另外一本书上的说法叫做“假设事情已经做成了”的“算盘思维”。先穿鞋戴帽，写好必要的代码，再在此基础上实现各种功能，完全没问题。)

选点的话，当然可以用系统自带的循环，但为了数据处理的效率，可以用numpy这个库生成网格(为什么知道这里可以用numpy？这是熟练使用Python与否的问题了)，正如代码里展示的那样。

判定怎么办？我们不可能真的计算无穷多次来确定数列会不会发散，而且实际上有定理说明，Mandelbrot集合本身是“不可计算”的，其补集则是“递归可枚举”(可以粗略理解成存在一个算法，使得只要点属于该集合，算法一定能在有限时间判定出来)的。

我们的做法是，只迭代计算若干次(用一个变量max_repeat控制)，如果届时还不发散，便判定为在集合内。更妙的是，因为有以下定理，判断发散的上限也变得方便了：如果某时刻 n 有

则该数列最终将发散到无穷远。

所以我们的判断方法是：不断用那个映射迭代，发现得到的绝对值“不对劲”就立刻抛弃；如果在max_repeat之后依然没有逃出这个圆盘(因为上面的判据的几何表示恰好是一个以原点为圆心的圆盘)，我们就谨慎地判定这个点合格(，并把它交给我们的画笔去描点)。

总结一下，上面这段代码顺序做了以下事情：

基本的调用库；

定义判定用的函数，后面再调用；(有趣的是这段代码在开发过程中是最晚写上的)

基本的画布设定；

生产接下来描点步骤要用的点；

循环体描点，其中调用前面的判定函数；

done()。

再强调一遍，代码的“物理顺序”与逻辑顺序，甚至是开发时写作的顺序可以有很大的差别，需要习惯这一点。

最后说一下这段代码还可以如何改良吧。

就以前面的有色图片为例，如果不使用简单的判定函数，而是使判定函数输出判定为发散时的迭代次数(判定为属于集合的输出max_repeat)，则就可以根据迭代次数选择不同的画笔颜色，表现出介于有界与发散边界的值发散的速度。(实际上还有选择适当的调色板等细节，这里从略。)

如果想把这段代码继续改短，其实可以考虑把所有的常数变量名都用对应的值代替(例如那个max_repeat实际上可以不定义)。但我强烈反对这种做法，因为这样一来会让人不理解函数传入参数的含义，在试图修改代码以适合新环境时也会遇到阅读的困难。

以上。