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 C 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: (mandel.c) #include #include #include "complex.c" #define OUTFILE "fractal.bin" #define NMAX 90000 #define PERMS 0666 main() { int fp1; /* not a pointer for binary files */ int nx, ny, ix, iy, it, idat, ndat, nbytes, status; float x1, x2, y1, y2, dx, dy; fcomplex c, z, z2; int dat[NMAX],nxny[2]; float wind[4]; printf("Enter x1,x2,y1,y2,nx,ny \n"); scanf("%f %f %f %f %d %d",&x1,&x2,&y1,&y2,&nx,&ny); dx = (x2 - x1) / (float)nx; dy = (y2 - y1) / (float)ny; if (nx*ny > NMAX) { printf("Too big for current array size \n"); return 0; } idat=0; for (ix=0; ix 4.0) break; } dat[idat] =it; idat++; } } printf("Finished calculation \n"); ndat = idat; printf("ndat = %d \n",ndat); fp1 = creat(OUTFILE,PERMS); if (fp1 == -1) { printf("***Can't open file: %s\n", OUTFILE); return 1; } nxny[0] = nx; nxny[1] = ny; nbytes = 8; status = write (fp1,nxny,nbytes); if (status != nbytes) { printf(" *** ERROR1 *** disk file write error: %d \n",status); } wind[0]=x1; wind[1]=x2; wind[2]=y1; wind[3]=y2; nbytes = 16; status = write (fp1,wind,nbytes); if (status != nbytes) { printf(" *** ERROR2 *** disk file write error: %d \n",status); } nbytes = 4*ndat; printf("writing nbytes = %d \n",nbytes); status = write (fp1,dat,nbytes); if (status != nbytes) { printf(" *** ERROR3 *** disk file write error: %d \n",status); } return 0; } 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. The maximum allowed resolution is set by NMAX where nx*ny <= NMAX. Although the image is really two-dimensional, I store the results in the 1-D array, dat[NMAX}, because this simplies the output of the results. The spacing in x and y are defined by: dx = (x2 - x1) / (float)nx; dy = (y2 - y1) / (float)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; for (ix=0; ix 4.0) break; } dat[idat] =it; idat++; } } The real and imaginary parts of c are set to the center of the "pixel" with sides dx and dy. The calculation is done with the Numerical Recipes routines Cadd and Cmul for a maximum of 1000 iterations. The calculation is halted if |z| > 2 (computed as z^2 > 4 because it is faster to avoid the square root) by exiting the for loop. The iteration number is saving 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. Note that the file number, fp1, in this case is NOT defined as a pointer as we have seen in earlier IO examples. IN C BINARY FILE NUMBERS ARE INTEGERS, NOT POINTERS!!!!!! We create and open the file with fp1 = creat(OUTFILE,PERMS); where we have defined OUTFILE as "fractal.bin" and PERMS as 0666. PERMS is supposed to set the permissions for the file but I have not been able to get this to work right. If there is an error in opening the file, then fp1 = -1 so the program checks for this and prints out a suitable message if there is a problem. 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 C is done using the "write" statement which has the form: write (fp1, xarray, nbytes); where fp1 is the file number, xarray is the name of an array, and nbytes is the number of bytes of data to write. For int or float arrays, each number is stored with 4 bytes, so normally nbytes should be set to 4 times the number of data points. (In contrast, char arrays have 1 byte per element, double arrays have 8 bytes per element.) So we set up a 2-element int array, nxny, to handle the output of nx and ny: nxny[0] = nx; nxny[1] = ny; nbytes = 8; status = write (fp1,nxny,nbytes); if (status != nbytes) { printf(" *** ERROR1 *** disk file write error: %d \n",status); } and we set up a 4-element float array, wind, to handle the output of x1, etc. wind[0]=x1; wind[1]=x2; wind[2]=y1; wind[3]=y2; nbytes = 16; status = write (fp1,wind,nbytes); if (status != nbytes) { printf(" *** ERROR2 *** disk file write error: %d \n",status); } The write function returns an integer value set to the number of bytes that were successfully written. The program checks to make sure that this equals what was supposed to be written and prints an error message if there is a problem. Finally, we write the output array to the file" nbytes = 4*ndat; printf("writing nbytes = %d \n",nbytes); status = write (fp1,dat,nbytes); if (status != nbytes) { printf(" *** ERROR3 *** disk file write error: %d \n",status); } For this program, it is a good idea to compile with the -fast compiler option to make it run faster: cc -fast -lm mandel.c -o mandel Of course, we also need to explicitly load the math library in the case of the Sun compiler. 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: (plotfractal.m) clf reset clear [fid,message] = fopen('fractal.bin') [size,count] = fread(fid,[2],'int'); nx=size(1) ny=size(2) [wind,count] = fread(fid,[4],'float') x1=wind(1) x2=wind(2) y1=wind(3) y2=wind(4) [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; In this case, 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. SPECIAL PROJECT 1 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 C and Matlab source code. Have fun!