% ### EXfractal1.m ### 10.08.14
% Creates a (color) fractal based upon the Mandelbrot set
% [modified from source code of Steve Lord/Birkhauser]
% Uses the difference map (i.e., discrete ODE)
% z_{n+1} = (z_n)^2 + c (z,c are both complex)
% where c = x+i*y (i.e., some point in the complex plane)
% Basic idea is to choose some initial z (e.g., z_0= 0) and 'iterate' the
% map (i.e., increment n some # of times). For that c, see if the value
% 'diverges' (i.e., blows up towards infinity). If not, we say that that
% point is part of the Mandelbrot set. If so, we can color code the
% divergence rate (so to get a prettier color image)
% NOTE: At best, we only have an approximation of the Mandelbrot set, since
% we can only iterate the map a finite # of times (and a limited # at
% that!). Basically, if |z|>2, the point will likely diverge. In general,
% this can be seen after ~30 iterations.
% Some values to try (successively zooming in):
% - xB=[-2.4 1.2]; yB=[-1.5 1.5]; scale=0.005; % {BROAD VIEW}
% - xB=[-1.4 -1.05]; yB=[0.15 0.5]; scale=0.0005;
% - xB=[-1.31 -1.22]; yB=[0.3 0.45]; scale=0.0001;
clear; figure(1); clf;
% -------------------------------------
Nmax = 30; % # of iterations to run to check for divergence {30}
maxZ= 2; % threshold |z| to determine divergence {2}
% spatial scale of complex plane (x+i*y) to consider
xB= [-2.4 1.2]; % x-range
yB= [-1.5 1.5]; % y-range
scale = 0.005; % spacing between adjacent points
fileN= './mandelbrot1.jpg'; % filename to save image to
% -------------------------------------
% generate a grid of points in the specified complex plane, and determine
% the associated c value (i.e., simply x+i*y)
[x,y]=meshgrid(xB(1):scale:xB(2), yB(1):scale:yB(2));
c = x+ i*y;
% Initialize matricies for iterated value (z) and number of iterations
% before blowing up (k); thus everything will be done at once in a single
% while loop (rather than nested for loops for each x,y pair)
z = zeros(size(c));
k = zeros(size(c));
% +++++
% loop to iterate the grid and check for divergence (flagging the array k);
% loop stops once N reaches Nmax or all points have been determined to have diverged
N = 1; % initialize index
while NmaxZ at this iteration and no
% previous iteration get assigned the value of N (indicates rate of
% divergence and therefore means to color code)
k(~k & abs(z)>maxZ) = N; % clever way to update values conditionally
N = N+1; % Increment iteration count
end
% If any k's are equal to 0 (i.e. the corresponding w's never blew up) set
% them to the final iteration number
k(k==0) = Nmax;
% Display the matrix as a surface (color by default)
s=pcolor(x,y,k);
tb = colorbar('peer',gca); set(get(tb,'ylabel'),'String', 'k'); % create colorbar and label
% plot in B&W??
if (1==0), colormap(bone); end
% Turn off the edges of the surface (because the cells are so small, the
% edges would drown out any useful information if we left them black)
set(s,'edgecolor','none')
% Adjust axis limits, ticks, and tick labels
axis([xB(1) xB(2) yB(1) yB(2)])
fs=15; xlabel('Re z','FontSize',fs); ylabel('Im z','FontSize',fs);
% +++++
% save picture to file as a jpg w/ a user=specified resolution
REZ= '-r180'; % resolution for exporting colormaps to jpg
print('-djpeg',REZ,[fileN]);