Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Is there a simpler way to construct Mandelbrot set in Matlab?

The code shown below for drawing a Mandelbrot set, I think my code a bit redundancy to construction the Matrix M. In Python I know there is a clean way do this,

M = [[mandel(complex(r, i)) for r in np.arange(-2, 0.5,0.005) ] for i in np.range(-1,1,0.005)]

Is there a similar way do this in 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
like image 882
matrix42 Avatar asked Jan 11 '23 18:01

matrix42


2 Answers

You can avoid the loop altogether. You can do the iteration z = z.^2 + c in a vectorized manner. To avoid unnecessary operations, at each iteration keep track of which points c have already surpassed your threshold, and continue iterating only with the remaining points (that's the purpose of indices ind and ind2 in the code below):

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

enter image description here

like image 162
Luis Mendo Avatar answered Jan 14 '23 09:01

Luis Mendo


You could at least avoid the for-loop:

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
like image 43
Robert Seifert Avatar answered Jan 14 '23 09:01

Robert Seifert