Mandelbrot set escape value is complex?

My codes are as follows,

clear all, clc
[re,im] = meshgrid(-3:0.5:3,-3:0.5:2);                                       % create a 2D grid ranging from -5 to 5 in both x and y
c = re+ j*im;                                                                % turn the c value into complex values.

row = size(c,1);                                                             % gives the row dimension of the matrix
column = size(c,2);                                                          % gives the column dimension of the matrix
k=0;
Lmat = zeros(row,column);                                                    % a zero matrix that stores the lastval
% The nested for loop here prints the value of each
for i = 1:1:row
    for j = 1:1:column
        cval = c(i,j);
        k = k+1;
        % Run the iteration loop here
        zn = 0;
        iteration = 0;                                                       % i stands for iteration
        n = 10;

        while iteration <(n+1)                                               % (n + 1) to denote the number of iterations
            zn = (zn)*(zn')+cval;
            lastval = zn;
            iteration = iteration+1;

            %fprintf('Iteration #%5.0f.lastval: %4.5f\n',iteration,lastval)
            if iteration == n+1
                %fprintf('loop ended.\n')
                break
            end
            Lmat(i,j) = lastval;
        end
        % Lmat(i,j) = lastval;
        % fprintf('j=%3.0f.\n',j)
    end
    % fprintf('i=%3.0f.\n',i)
end
surf(re,im,real(Lmat))

And when I generated the plot it just doesn't look like what it's supposed to.

Is the escape value supposed to be complex? I think it should? because escapevalue = zn*zn' + cval and cval by itself is complex. and therefore the cval alone has a contribution to the escape value.

I am also thinking that my mathematical logic in this code is simply incorrect. Since I attempted to plot every cval, but I thought the purpose of mandelbrot set is to test out which cval can escape and which cval cannot escape.

Please if you could let me know where I got it wrong I would greatly appreciate it!

Thanks,

Michael.

[Jan]

Plotting the positions of the poiunt over the iterations is the Julia set, not the Mandelbrot set. For the Mandelbrot set count the number of iterations you need to leave a radius of 2.

Why do you use 2 stopping methods to limit the loop?

while iteration <(n+1)
    ...
        if iteration == n+1
        

[Michael Lam]

Did not see your comment here until just now.

It is true that my while loop is set up in a not very smart way, but when I took out

iteration == n + 1

then it stops working. I haven't bothered to fix the code there yet. But that's probably one of the errors that I have to fix.

Just earlier today I was able to plot the mandelbrot set already.

The error that I ran into yesterday was that I plotted the not , of which I would need to take the norm of it before storing it to a zero matrix.

Check the code here.

clear all, close all, clc
dr = 500;
xmin = -2; ymin = -1;
xmax = 1; ymax = xmax;
[re,im] = meshgrid(linspace(xmin,xmax,dr),linspace(ymin,ymax,dr));          % create a 2D grid ranging from -5 to 5 in both x and y
c = re + im*1i;                                                             % turn the c value into complex values.

% declaring empty matrices to store values
row = size(c,1); column = size(c,2);                                        % gives the row and columns of the matrix 
L = zeros(row,column);                                                      % a zero matrix that stores the Last value of zn
E = zeros(row,column);                                                      % a zero matrix that stores the Escape values, which is mathematically norm(z).

% This here actually plots it 
for a = 1:1:row
    for b = 1:1:column
        cval = c(a,b);                                                      % Run the iteration loop here
        zn = 0;
        iteration = 0;                                                      
        n = 10;
        % Now the error is less of an error but something's still not right
        while iteration < (n+1)                                             % (n + 1) to denote the number of iterations
            zn = (zn).^2+ cval;                                             % the mandelbrot equation which tells you if it's real or complex. 
            E(a,b) = norm(zn);
            iteration = iteration+1;                                        %This line is unnecesary.
            if norm(zn) > 2
                break
            end
            L(a,b) = zn;
        end
    end
    fprintf('Working on it.%3.2f%%.\n',a.*100/dr)
end
%contour(re,im,E)
% plot(re,im,E)
imagesc(E)
colormap copper
axis normal
fprintf('done\n')
% savefig('Mandelbrot_set_contour.fig')

The graph is generated decently enough, but the details (not the resolution) of the graph is not as well as I wanted it, and at this point I'm not sure what I can modify in my code to improve it.

As with Julia set, plotting over the "c space" makes for the Julia set that is correct. Mandelbrot set plots the varying c values, which I simply have c be a matrix by itself, and called the respective values in the c matrix c(a,b), following the format that .

I guess I'm not really plotting over the c space, I'm just using the c matrix to extract the varying c values and then plot it. I do plan on making a Julia set and finish it in the next few days. So I'll probably keep posting here or create another post.

Thanks for the help!

Answer 1

[Jan]

Here a cleaned version of your code:

• Calculating norm() twice is a waste of time.
• abs() is faster than norm().
• No need to create the c matrix explicitly.
• Run 3 times faster
• FOR instead of WHILE for more compact code.

dr   = 500;
xmin = -2; ymin = -1;
xmax =  1; ymax = xmax;
re = linspace(xmin, xmax, dr);
im = linspace(ymin, ymax, dr);

row    = numel(re);
column = numel(im);

E = zeros(row, column);
n = 100;
for a = 1:row
   for b = 1:column
      c = re(b) + 1i* im(a);                                                      
      z = 0;
      for iter = 0:n                                          
         z = z ^ 2 + c;
         if abs(z) > 2
            break
         end
      end
      E(a, b) = iter;
   end
end

imagesc(re, im, E)
colormap copper
axis normal

[Michael Lam]

Thank you for the correction! This helps me to improve the efficiency and details in the way I code. I'll make sure to compare between the two and figure out the differences.

[Jan]

I've played with the code a little bit and found some interesting timings:

dr   = 500;
xmin = -2;
xmax =  1;
ymin = -1;
ymax =  1;
re   = linspace(xmin, xmax, dr);
im   = linspace(ymin, ymax, dr);

tic; E1 = mandel_1(re, im); toc

Elapsed time is 0.203943 seconds.

tic; E2 = mandel_2(re, im); toc

Elapsed time is 0.047390 seconds.

% ==================================
function E = mandel_1(re, im)
row = numel(re);
col = numel(im);
E   = zeros(row, col);
n   = 100;
for b = 1:col
   reb = re(b);
   for a = 1:row
      c = reb + 1i * im(a);                                                      
      z = 0;
      for iter = 0:n                                          
         z = z ^ 2 + c;
         if abs(z) > 2
            break
         end
      end
      E(a, b) = iter;
   end
end
end

% ==================================
function E = mandel_2(re, im)
row = numel(re);
col = numel(im);
E   = zeros(row, col);
n   = 100;
% Swap loops to write to E columnwise - some percent faster
for b = 1:col  % Use PARFOR for huge input
   cr = re(b);
   for a = 1:row
      ci = im(a);                                                      
      zr = cr;
      zi = ci;
      
      % FOR ==> WHILE: just some percent faster
      % Avoid SQRT or ABS and measure squared distance: 50% faster
      iter = 0;
      while iter < n && (zr * zr + zi * zi) <= 4
         % Working with real numbers only: 50% faster:
         zr_ = zr * zr - zi * zi + cr;
         zi  = 2 * zr * zi + ci;
         zr  = zr_;
         
         iter = iter + 1;
      end
      E(a, b) = iter;
   end
end
end