在Matlab中，有没有更简单的构造Mandelbrot集合的方法？

Question

下面是绘制曼德博集合的代码，我认为我的代码在构建Matrix M时有些冗余。在Python中，我知道有一种简洁的方法来做到这一点， M = [[mandel(complex(r, i)) for r in np.arange(-2, 0.5,0.005) ] for i in np.range(-1,1,0.005)]] 在Matlab中是否有类似的方法可以实现？

function M=mandelPerf()
rr=-2:0.005:0.5;
ii=-1:0.005:1;
M = zeros(length(ii), length(rr));
id1 = 1;
for i =ii
    id2 = 1;
    for r = rr
        M(id1, id2) = mandel(complex(r,i));
        id2 = id2 + 1;
    end
    id1 = id1 + 1;
end
end

function n = mandel(z)
n = 0;
c = z;
for n=0:100
    if abs(z)>2
        break
    end
    z = z^2+c;
end
end

- matrix42

---

请注意，由于Madelbrot集合在y=0处具有线性对称性，因此您至少需要将计算量加倍。 - alko

Answer 1

您可以完全避免循环。您可以以矢量化的方式进行迭代 z = z.^2 + c 。为避免不必要的操作，在每次迭代中，跟踪已经超过阈值的点c，并仅使用剩余点进行迭代（这是下面代码中索引ind和ind2的目的）：

rr =-2:0.005:0.5;
ii =-1:0.005:1;
max_n = 100;
threshold = 2;
c = bsxfun(@plus, rr(:).', 1i*ii(:)); %'// generate complex grid
M = max_n*ones(size(c)); %// preallocation.
ind = 1:numel(c); %// keeps track of which points need to be iterated on
z = zeros(size(c)); %// initialization
for n = 0:max_n;
    z(ind) = z(ind).^2 + c(ind);
    ind2 = abs(z(ind)) > threshold;
    M(ind(ind2)) = n; %// store result for these points...
    ind = ind(~ind2); %// ...and remove them from further consideration
end

imagesc(rr,ii,M)
axis equal

- Luis Mendo

---

更改颜色映射:P - Ander Biguri

@AnderBiguri 我知道，我知道... 那时候还没有hg2 :-) - Luis Mendo

@AnderBiguri，实际上我刚刚尝试了parula和viridis，但我发现旧的彩虹色图更适合这个问题，因为它可以突出边缘。 - Luis Mendo

你可能是对的。最终，这个图像只有两种颜色和边缘。 - Ander Biguri

Answer 2

你至少可以避免使用for循环：

function M=mandelPerf()

rr = -2:0.005:0.5;
ii = -1:0.005:1;

[R,I] = meshgrid(rr,ii);
M = arrayfun(@(x) mandel(x), R+1i*I);

end

function n = mandel(z)
n = 0;
c = z;
for n=0:100
    if abs(z)>2
        break
    end
    z = z^2+c;
end
end

- Robert Seifert

---

请注意，arrayfun通常只是一个伪装成for循环的东西。它不一定比for循环节省时间。例如，请参见这里。 - Luis Mendo

@LuisMendo 我知道这一点，但对于曼德博集背后的理论不够熟悉，因此不能跳过不必要的步骤。OP给出的Python示例也只是在循环执行，对吧？ - Robert Seifert

关于Python，我无话可说... 这里的for循环可以轻松向量化；请参见我的答案。 - Luis Mendo