在Python中设置简单的Mandelbrot

Question

我已经看过其他关于这个问题，但我似乎无法弄清楚我要去哪里错了。我们的目标是“编写一个程序来生成Mandelbrot集的图像，方法是对跨越-2≤x≤2和-2≤y区域的N×N网格执行c = x + iy的所有值的迭代≤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() 

对于任何未来的读者，这是我新的代码是如何结束了寻找：

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做大部分循环的版本。它使用更多的内存，但比前一版本快大约2.5倍，并且略短一些。

此代码使用几个二维数组。 c包含所有复杂的种子编号，我们在z执行核心Mandelbrot计算，它被初始化为零。 mask是一个布尔数组，用于控制需要执行Mandelbrot计算的位置。其所有项目初始设置为True，并且在每个迭代True项目mask中对应于z中已从Mandelbrot集合中跳出的项目设置为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)