基于Python画曼德勃罗集

最近有点迷分型几何，看到“上帝指纹”曼德勃罗集，想用Python实现一下。
源码很简单20行不到。

import matplotlib.pyplot as plt
import numpy

def mb(x,y):
    C = complex(x,y)
    Z = 0
    for i in range(100):
        Z = Z*Z +C
        if abs(Z)>2:
            return False
    return True

for x in numpy.arange(-2,2,0.01):
    for y in numpy.arange(-1,1,0.01):
        if mb(x,y):
            plt.scatter(x,y,s=0.5)

plt.show()

运行结果↓

由于电脑性能受限，步长实在是不能取的太精细。。。实在有点丑。附上网图一张。

从运行结果来看（也是一步一步探索出来的），曼德勃罗集的C的范围大概是：实部∈(-2,0.5)，虚部∈(-1,1)。（这里不禁吐槽下，关于曼德罗集的图片很多，但是几乎没有人标上坐标轴）。
既然有了代码，我们就可以再看下细节。

继续放大（当然毫无疑问，肯定是无限细节）

最后再放大一下

------------------更新-----------------
发现用imshow的方式画出的图会好看很多，直接再附源码

from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
import math


def mb(x, y):
    C = complex((x - 1500) /600, (500 - y) / 500)
    Z = 0
    for i in range(200):
        Z = Z**2 + C
        if abs(Z)>2:
            return False
    return True


image = Image.new('RGB', (2000, 1000))
image_array = np.array(image)

for k in range(1000):
    for v in range(2000):
        if mb(v, k):
            image_array[k, v, :] = 255
        else:
            image_array[k, v, 0] = abs((v - 1500) /600 * 255/(1500/600))
            image_array[k, v, 1] = abs((500 - k) / 500 * 255)
            image_array[k, v, 2] = abs((v - 1500) /600 * 255/(1500/600)/2 + (500 - k) / 500 * 255/2)

plt.imshow(image_array)
plt.show()

运行如下↓

既然有了美美的代码，就可以玩点花的了：稍微改下曼德勃罗集的公式。
Zn+1 = Zn^3 + C 结果↓

Zn+1 = Zn^2+Zn + C 结果↓（这个挺好看的）

-----------------再更新-------------------
（没想到这个还能再更一次。。。）
根据迭代公式，实现了曼德勃罗集的绘制。
那么如果不仅限于二维，如果扩展到三维会是怎么样的？
很遗憾，曼德勃罗集在三维上是不存在的。但是可以在四维上用双复数，按照曼德勃罗集的迭代方式构建四维的曼德勃罗集：曼德尔球

思维再扩展一下，那么能否也用简单的迭代规则，构建一个三维的分型呢？

我自己构建了一个规则，但是受制于电脑性能，并未看出它的分型结构（而且从初步的结果来看，它很有可能还不是分型的，看起来似乎是有解析解）

构建这个三维图形的规则如下：

这里我们要找的集合就是这样的 Z1 的集合，使得在有限次迭代后，Zn能够收敛

说明：

• 这里向量的每次迭代都是用外积，理由也很简单，两个向量外积之后仍然是向量，能够迭代下去；
• 这里指定了Z0 是单位向量（1，1，1），实际上Z0也可以取任意值，然后找出Z0的集合（这个集合也就和Z1一样），这里为了减少计算机的工作量，指定了Z0为单位向量。

运行结果

受制于电脑的性能，点云没办法取得太密集。
从下面Z1的集合来看，并不能看出分型的结构。

既然三维结构看不出来的话，只能取截面看看它的形状。
奇怪且有趣的是，这个集合似乎在xoy, yoz, zox三个面上的截面图形都是一样的。

最后，如果有代码大神能绘制出这个集合的三维图像，还请指教一下。

最后的最后，附上源码

import matplotlib.pyplot as plt
from PIL import Image
import numpy
import warnings
warnings.filterwarnings("ignore")

def iter(x,y,z):
    x1, y1, z1=1, 1, 1
    x2 = z-y
    y2 = x-z
    z2 = y-x
    for i in range(10):
        x_mid = y1*z2 - y2*z1
        y_mid = z1*x2 - z2*x1
        z_mid = x1*y2 - x2*y1
        x1 = x2
        y1 = y2
        z1 = z2
        x2 = x_mid
        y2 = y_mid
        z2 = z_mid
        # print(x2,y2,z2)
    if abs(x2) and abs(y2) and abs(z2) >1:
        return False
    else:
        return True

#-------------------------三维点云---------------------------------
# fig = plt.figure()
# ax = fig.add_subplot(111, projection='3d')
#
# k =10
# for x in range(-k,k):
#     for y in range(-k,k):
#         for z in range(-k, k):
#             if iter(x/k, y/k, z/k):
#                 ax.scatter(x/k, y/k, z/k,s=2, c='r')
#
# plt.show()
#----------------------------------------------------------------

#----------------xoy二维截面------------------------------------------
# image = Image.new('RGB', (1000, 1000))
# image_array = numpy.array(image)
#
# k=500
# for x in range(-k,k):
#     for y in range(-k,k):
#         if iter(x/k, y/k, 0):
#             image_array[x+k,y+k]= 255
#         else:
#             image_array[x+k, y+k] = 0
#
# plt.imshow(image_array)
# plt.show()
#----------------------------------------------------------------

#----------------xoz二维截面------------------------------------------
# image = Image.new('RGB', (1000, 1000))
# image_array = numpy.array(image)
#
# k=500
# for x in range(-k,k):
#     for z in range(-k,k):
#         if iter(x/k, 0, z/k):
#             image_array[x+k,z+k]= 255
#         else:
#             image_array[x+k, z+k] = 0
#
# plt.imshow(image_array)
# plt.show()
#----------------------------------------------------------------

#----------------yoz二维截面------------------------------------------
image = Image.new('RGB', (1000, 1000))
image_array = numpy.array(image)

k=500
for y in range(-k,k):
    for z in range(-k,k):
        if iter(0, y/k, z/k):
            image_array[y+k,z+k]= 255
        else:
            image_array[y+k, z+k] = 0

plt.imshow(image_array)
plt.show()
#----------------------------------------------------------------