% ### 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]);