Python 中的简单 Mandelbrot 集

Question

我已经查看了与此相关的其他问题，但我似乎无法弄清楚我哪里出错了。目标是 "Write a program to make an image of the Mandelbrot set by performing the iteration for all values of c = x + iy on an N × N grid spanning the region where −2 ≤ x ≤ 2 and −2 ≤ y ≤ 2. Make a density plot in which grid points inside the Mandelbrot set are colored black and those outside are colored white." 当迭代为 z' = z^2 + c 时，在 Mandelbrot 集中被认为是 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()

对于任何未来的读者，这是我的新代码最终的样子：

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 中的浮点数转换为整数，但结果数字将在 -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 是一个布尔数组，控制需要在哪里进行Mandelbrot计算。它的所有项目最初都设置为 True，并且在每次迭代中 mask 中对应于 z 中已从 Mandelbrot 集合中逃脱的项目的 True 项目被设置为 False。

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

我们可以使用 mask 的最终值来绘制 Mandelbrot 集，但是要使 Mandelbrot 集中的点变黑，我们需要反转它的值，我们可以使用 1 - mask .

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)

Answer 2

from pylab import imshow,show,gray
from numpy import zeros,linspace
N=1000
A=zeros([N+1,N+1],float)
def mandler(x,y):
    c=complex(x,y)
    z=c
    n=0
    while n<N:
        z=z*z+c
        if abs(z)>2:
            return 0
        n+=1
    else:
        return 1
for x in linspace(-2,2,N):
    for y in linspace(2,-2,N):
        a=mandler(x,y)
        s=int(((x+2)*N)/4)
        t=int(((y+2)*N)/4)
        A[s,t]=a
imshow(A)
gray()
show()

Python 中的简单 Mandelbrot 集

Simple Mandelbrot set in Python

python

numpy

mandelbrot