FUN WITH FRACTALS Generating beautiful images of fractals, such as the Mandelbrot set, is surprisingly easy on the computer and provides good practice in using complex numbers. The basic idea behind most fractal calculations is that there are certain complex functions that, when computed over and over again, will either diverge or remain bounded. Whether or not they diverge turns out to be extremely sensitive to small changes in the initial complex value that starts the calculation, at least in certain regions of the complex plane. The most famous of all fractal images is the Mandelbrot set. To generate the Mandelbrot set, we begin by considering a complex number, c. We then apply the following algorithm: set z = 0 to start then repeatedly compute z = z*z + c until |z| > 2 OR the number of iterations exceeds some threshold then output the number of iterations For example, if c = 0.3 + 0.3i then 1st iteration: z = 0.30 + 0.30i |z| = 0.42 2nd iteration: z = 0.30 + 0.48i |z| = 0.57 3rd iteration: z = 0.16 + 0.59i |z| = 0.61 4th iteration: z = -0.02 + 0.49i |z| = 0.49 In this case, z will remain bounded even after an infinite number of iterations. However, if c = 0.5 + 1.0i then 1st iteration: z = 0.50 + 1.00i |z| = 1.1 2nd iteration: z = -0.25 + 2.00i |z| = 2.0 3rd iteration: z = -3.44 + 0.00i |z| = 3.4 4th iteration: z = 12.32 + 1.00i |z| = 12.4 5th iteration: z = 151.19 + 25.63i |z| = 153.4 and the size of z explodes toward infinite values. By the 10th iteration it will exceed the floating point range of most computers. However, we can stop computing z once its absolute value (the complex absolute value or modulus is defined as the distance of a point in the complex plane from the origin, i.e. |z| = abs(a+bi) = sqrt(a*a+b*b) ) exceeds 2 because it can be shown that divergence is guaranteed at this point. We repeat this calulation for a series of different values of c and plot the resulting output as a function of the location of c on the complex plane. The Mandelbrot set is the set of all complex numbers c for which the size of z^2 + c is finite even after an infinite number of iterations. We may obtain a good approximation, however, by checking after some large number of iterations (1000 is a common choice). We also may take advantage of the fact that it is known that z^2+c will diverge whenever |z| > 2. Here is an example F90 program that performs this calculation for user-specified parts of the complex plane and outputs a binary file that can be read and plotted using MatLab: program mandel implicit none real, dimension(:), allocatable :: dat real :: x1, x2, y1, y2, dx, dy, cr, ci integer :: nx, ny, ix, iy, it, idat, ndat, nbytes, status integer, dimension(2) :: nxny complex :: c, z, z2 real, dimension(4) :: wind print *, 'Enter x1, x2, y1, y2, nx, ny' read *, x1, x2, y1, y2, nx, ny dx = (x2 - x1) / real(nx) dy = (y2 - y1) / real(ny) allocate(dat(nx*ny)) idat=0 do ix = 1, nx do iy = 1, ny cr = x1 + dx/2. + dx*real(ix-1) ci = y1 + dy/2. + dy*real(iy-1) c = cmplx(cr, ci) z = cmplx(0.0, 0.0) do it = 1, 1000 z = c + z*z if (cabs(z) > 2) exit enddo idat = idat + 1 dat(idat) = it enddo enddo print *, 'Finished calculation' ndat = idat print *, 'ndat = ',ndat open (12, file='fractal.bin', form='unformatted') nxny(1) = nx nxny(2) = ny write (12) nxny wind(1) = x1 wind(2) = x2 wind(3) = y1 wind(4) = y2 write (12) wind write (12) dat end program mandel The program computes the result for a grid of points in the complex plane, defined by the user-defined limits x1 and x2 for the real part and y1 and y2 for the complex part. The number of points computed in x and y are input as nx and ny. This will determine the resolution of the plot. For example, nx=50 and ny=50 will result in an image 50 pixels wide by 50 pixels tall. Although the image is really two-dimensional, I store the results in the 1-D array, dat, because this simplies the output of the results. Notice how allocate is used to set dat to dimension(nx*ny) on the fly. The spacing in x and y are defined by: dx = (x2 - x1) / real(nx) dy = (y2 - y1) / real(ny) Note the conversion of the integer variables nx and ny to floats prior to performing the division. Here are the main loops that actually do the calculation: idat=0 do ix = 1, nx do iy = 1, ny cr = x1 + dx/2. + dx*real(ix-1) ci = y1 + dy/2. + dy*real(iy-1) c = cmplx(cr, ci) z = cmplx(0.0, 0.0) do it = 1, 1000 z = c + z*z if (cabs(z) > 2) exit enddo idat = idat + 1 dat(idat) = it enddo enddo The real and imaginary parts of c are set to the center of the "pixel" with sides dx and dy. The calculation is done for a maximum of 1000 iterations. The calculation is halted if |z| > 2 by exiting the for loop. The iteration number is saved in the dat array, using idat as a counter. The next step is to save this information to a file. For speed in the case of a large dat array, we use a binary file. We create and open the file with open (12, file='fractal.bin', form='unformatted') We use form='unformatted' to specify that this is a binary file, not an ascii file. Before writing the dat array to the file, we write some header information that will help MatLab to know the size of the array and the values of x1, x2, y1, and y2. Writing to a binary file in Fortran is done using the "write" statement without a format. We set up the 2-element int array, nxny, to handle the output of nx and ny: nxny(1) = nx nxny(2) = ny write (12) nxny and we set up a 4-element float array, wind, to handle the output of x1, etc. wind(1) = x1 wind(2) = x2 wind(3) = y1 wind(4) = y2 write (12) wind Finally, we write the output array to the file write (2) dat For this program, it is a good idea to compile with the -O compiler option to make it run faster: The Mandelbrot set is approximated in our program by places on the complex plane where the output value (contained in the dat array) is 1000, i.e., the z function has not diverged even after 1000 iterations. Output values less than 1000 provide a message of how quickly the function diverges (small values indicate rapid divergence, larger values indicate slower divergence). Plots of the output usually use colors to indicate the different output values. Here is how to plot our output using Matlab: clf reset clear [fid, message] = fopen('fractal.bin') % Fortran I/O has extra four bytes at beginning and end of writes [dum, count] = fread(fid, [1], 'int'); [size, count] = fread(fid, [2], 'int'); [dum, count] = fread(fid, [1], 'int'); nx = size(1) ny = size(2) [dum, count] = fread(fid, [1], 'int'); [wind, count] = fread(fid, [4], 'float') [dum, count] = fread(fid, [1], 'int'); x1 = wind(1) x2 = wind(2) y1 = wind(3) y2 = wind(4) [dum, count] = fread(fid, [1], 'int'); [z, count] = fread(fid, [nx,ny], 'int'); fclose(fid); z2 = log10(z); axis( [x1, x2, y1, y2] ); axis('square'); hold on; imagesc( [x1, x2], [y1, y2], z2 ); colorbar; One detail is that Fortran (unlike C) writes an extra 4 bytes before and after each of the binary writes. So we need to read these extra btyes with our Matlab script into a dummy variable or we will not be properly aligned with the data fields. We set the output array to z and use imagesc to plot the result using the default color scheme. We plot log(z) rather than z because it shows some of the details a little better. The entire Mandelbrot set may be imaged in the window real -2 to 0.5, imag -1.25 to 1.25 (i.e., -2, 0.5, -1.25, 1.25) It is also fun to zoom in on details of the image. Here are some suggestions for places to look: -0.95 -0.90 0.25 0.30 -0.75 -0.65 0.25 0.35 -0.73 -0.72 0.282 0.292 -0.703 -0.693 0.26 0.27 Further details about generating fractals may be found in a series of articles by A.K. Dewdney in the Computer Recreations section that used to appear in Scientific American. The August 1985 article describes the Mandelbrot set, the November 1987 article talks about the related Julia sets, and the July 1989 article is about fractals called "biomorphs" that look like little creatures. There are also lots of books on fractals that have many beautiful images. The study of fractal is related to chaos theory and the fact that inphysical systems (weather is one example) are inherently unpredictable on large time scales because small perturbations to the starting conditions will cause large changes over time. FRACTALS HW1 Rewrite the mandel.c and plotfractal.m programs to plot examples of the Julia sets and the biomorphs. You will want to first read the appropriate Scientific American articles. (a) Produce hard copies of 3 especially nice images of Julia sets. (b) Reproduce and print hard copies of as many of the biomorphs as you can that are shown on p. 110 of the biomorph article. Part (b) is harder because you will have to generate more complicated complex functions than those that are used for the Mandelbrot and Julia sets. You may have to consult a book on complex numbers to find suitable formulae to compute sin(z), z^z, etc. Please turn in the hard copies and e-mail me copies of your F90 and Matlab source code. Have fun!