【问题标题】:Is there a simpler way to construct Mandelbrot set in Matlab?有没有更简单的方法在 Matlab 中构造 Mandelbrot 集?
【发布时间】:2013-12-24 20:27:49
【问题描述】:

下面显示的代码用于绘制Mandelbrot set,我认为我的代码对于构造矩阵 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

【问题讨论】:

  • 请注意,由于madelbrot集相对于y=0是线性对称的,因此您至少要加倍计算

标签: python matlab matrix


【解决方案1】:

您可以完全避免循环。您可以以矢量化方式进行迭代z = z.^2 + c。为避免不必要的操作,在每次迭代时跟踪哪些点 c 已经超过了您的阈值,并仅使用剩余点继续迭代(这就是下面代码中索引 indind2 的目的):

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

【讨论】:

  • 更改颜色图:P
  • @AnderBiguri 我知道,我知道...hg2 那时还不存在 :-)
  • @AnderBiguri 实际上我刚刚尝试过parulaviridis,我发现旧的彩虹色图更适合这个问题,其中所有的重点都是突出边缘
  • 你可能是对的。最后,这张图片是2种颜色,然后是边缘
【解决方案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

【讨论】:

  • 请注意,arrayfun 通常只是伪装的for 循环。与for 循环相比,它不一定会节省时间。以here 为例
  • @LuisMendo 我知道这一点,但对 Mandelbrot 背后的理论不够熟悉,无法跳过不必要的步骤。 OP给出的python示例也只是循环它,对吧?
  • 关于Python,我无话可说……这里的for循环可以很容易地向量化;看我的回答
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 2023-02-21
  • 1970-01-01
  • 1970-01-01
  • 2022-12-28
  • 1970-01-01
  • 1970-01-01
  • 2010-12-15
相关资源
最近更新 更多