python - Python 中的简单 Mandelbrot 集

Question

我已经查看了与此相关的其他问题，但我似乎无法弄清楚我哪里出了问题。目标是“编写一个程序，通过在跨越 −2 ≤ x ≤ 2 和 −2 ≤ y 区域的 N × N 网格上对 c = x + iy 的所有值执行迭代来制作 Mandelbrot 集的图像≤ 2. 绘制密度图，其中 Mandelbrot 集内的网格点涂为黑色，外部的网格点涂为白色。在 Mandelbrot 集合中，当迭代为 z' = z^2 + c 时，z 的大小永远不会大于 2。

我的代码生成了一个几乎全黑的网格，带有一小部分白色。

from pylab import imshow,show,gray
from numpy import zeros,linspace

z = 0 + 0j
n=100

M = zeros([n,n],int)
xvalues = linspace(-2,2,n)
yvalues = linspace(-2,2,n)

for x in xvalues:
    for y in yvalues:
        c = complex(x,y)
        for i in range(100):
            z = z*z + c
            if abs(z) > 2.0:
                M[y,x] = 1
                break

imshow(M,origin="lower")
gray()
show()

对于任何 future 的读者，以下是我的新代码的最终外观:

from pylab import imshow,show,gray
from numpy import zeros,linspace

n=1000

M = zeros([n,n],int)
xvalues = linspace(-2,2,n)
yvalues = linspace(-2,2,n)

for u,x in enumerate(xvalues):
    for v,y in enumerate(yvalues):
        z = 0 + 0j
        c = complex(x,y)
        for i in range(100):
            z = z*z + c
            if abs(z) > 2.0:
                M[v,u] = 1
                break

imshow(M,origin="lower")
gray()
show()

Answer 1

您的代码存在一些问题。

首先，您使用xvalues和yvalues来索引M，但这些索引应该是0..范围内的像素索引整数。 (n-1)。 Numpy 会将 xvalues 和 yvalues 中的 float 转换为整数，但结果数字将为 -2..2，因此绘制的像素不会很多，图像会很小，而且由于 Python 中负索引的工作方式，你会得到换行。

获取所需像素索引的一个简单方法是使用内置的 Python 函数 enumerate，但可能有一种方法可以重新组织代码以使用 Numpy 函数来执行此操作。

另一个问题是您需要将每个像素的 z 重置为零。目前，您的代码会重用前一个像素的最后一个 z 内容，如果该像素位于 Mandelbrot 集中，则 z 将太大。

这是代码的修复版本。我没有pylab，所以我使用PIL编写了一个简单的位图查看器。您可以通过在我的 show 函数中调用 img.save(filename) 将图像保存到文件中； PIL 将从文件扩展名中找出正确的文件格式。

import numpy as np
from PIL import Image

def show(data):
    img = Image.frombytes('1', data.shape[::-1], np.packbits(data, 1))
    img.show()

n = 100
maxiter = 100

M = np.zeros([n, n], np.uint8)
xvalues = np.linspace(-2, 2, n)
yvalues = np.linspace(-2, 2, n)

for u, x in enumerate(xvalues):
    for v, y in enumerate(yvalues):
        z = 0
        c = complex(x, y)
        for i in range(maxiter):
            z = z*z + c
            if abs(z) > 2.0:
                M[v, u] = 1
                break

show(M)

这是输出图像:

<小时/>

当然，每当您发现自己迭代 Numpy 数组索引时，就表明您做错了。使用 Numpy 的要点是它可以通过以 C 语言的速度在内部迭代整个数组来一次性对整个数组执行操作；上面的代码也可能使用普通的 Python 列表而不是 Numpy 数组。

这是一个让 Numpy 完成大部分循环的版本。它使用更多的 RAM，但速度比以前的版本大约快 2.5 倍，而且长度也更短。

此代码使用多个二维数组。 c 包含所有复数种子数，我们在 z 中执行核心 Mandelbrot 计算，其初始化为零。 mask 是一个 bool 数组，用于控制需要执行 Mandelbrot 计算的位置。其所有项目最初都设置为 True，并且在每次迭代时，mask 中的 True 项目对应于 z 中的项目> 已从 Mandelbrot 集中转义的设置为 False。

为了测试某个点是否已逃逸，我们使用 z.real**2 + z.imag**2 > 4.0 而不是 abs(z) > 2.0，这节省了一点时间，因为它避免了昂贵的平方根计算和 abs 函数调用。

我们可以使用 mask 的最终值来绘制 Mandelbrot 集合，但是为了使 Mandelbrot 集合中的点变为黑色，我们需要反转其值，我们可以使用 1 来实现- 掩码。

import numpy as np
from PIL import Image

def show(data):
    img = Image.frombytes('1', data.shape[::-1], np.packbits(data, 1))
    img.show()
    img.save('mset.png')

n = 100
maxiter = 100

a = np.linspace(-2, 2, n)
c = a + 1.j * a[:, None]
z = np.zeros((n, n), np.complex128)
mask = np.ones((n, n), np.bool)

for i in range(maxiter):
    mask[mask] = z[mask].real**2 + z[mask].imag**2 < 4.0
    z[mask] = z[mask]**2 + c[mask]

show(1 - mask)