Numerical Computing with MATLAB - Cleve Moler, Notas de estudo de Engenharia Elétrica

Baixe Numerical Computing with MATLAB - Cleve Moler e outras Notas de estudo em PDF para Engenharia Elétrica, somente na Docsity! Preface Numerical Computing with MATLAB is a textbook for an introductory course in numerical methods, Matlab, and technical computing. The emphasis is on informed use of mathematical software. We want you to learn enough about the mathematical functions in MATLAB that you will be able to use them correctly, appreciate their limitations, and modify them when necessary to suit your own needs. The topics include: • introduction to MATLAB • linear equations • interpolation • zero finding • least squares • quadrature • ordinary differential equations • random numbers • Fourier analysis • eigenvalues and singular values • partial differential equations George Forsythe initiated a software-based numerical methods course at Stan- ford University in the late 1960s. The textbooks by Forsythe, Malcolm, and Moler [1] and Kahaner, Moler, and Nash [2] that evolved from the Stanford course were based upon libraries of Fortran subroutines. This textbook is based upon MATLAB. NCM, a collection of over 70 M- files, forms an essential part of the book. Many of the over 200 exercises involve modifying and extending the programs in NCM. The book also makes extensive use of computer graphics, including interactive graphical expositions of numerical algorithms. The prerequisites for the course, and the book, include: 1 2 Preface • calculus • some familiarity with ordinary differential equations • some familiarity with matrices • some computer programming experience If you’ve never used Matlab before, the first chapter will help you get started. If you’re already familiar with Matlab you can glance over most of the first chapter quickly. Everyone should read the section in the first chapter about floating point arithmetic. There is probably too much material here for a one-quarter or one-semester course. Plan to cover the first several chapters and then choose the portions of the last four chapters that interest you. Make sure that the NCM collection is installed on your network or your per- sonal computer as you read the book. The software is available from a Web site devoted to the book, http://www.mathworks.com/moler There are three types of NCM files: • gui files. Interactive graphical demonstrations. • tx files. Textbook implementations of built-in Matlab functions. • Others. Miscellaneous files, primarily associated with exercises. When you have NCM available, ncmgui produces the figure shown on the next page. Each thumbnail plot is actually a push button that launches the corresponding gui. This book would not have been possible without the staff of The MathWorks. They are a terrific group of people and have been especially supportive of this book project. Out of the many friends and colleagues who have made specific contributions, I want to mention five in particular. Kathryn Ann Moler has used early drafts of the book several times in courses at Stanford and has been my best critic. Tim Davis and Charlie Van Loan wrote especially helpful reviews. Lisl Urban did an immaculate editing job. My wife Patsy has lived with my work habits and my laptop and loves me anyway. Thanks, everyone. – Cleve Moler, January 5, 2004 Bibliography [1] G. Forsythe, M. Malcolm, and C. Moler, Computer Methods for Math- ematical Computations, Prentice Hall, Englewood Cliffs, 1977. [2] D. Kahaner, C. Moler, and S. Nash, Numerical Methods and Software, Prentice Hall, Englewood Cliffs, 1989. [3] The MathWorks, Inc., Numerical Computing with MATLAB, http://www.mathworks.com/moler 5 Chapter 1 Introduction to MATLAB This book is an introduction to two subjects, Matlab and numerical computing. This first chapter introduces Matlab by presenting several programs that inves- tigate elementary, but interesting, mathematical problems. If you already have some experience programming in another language, we hope that you can see how Matlab works by simply studying these programs. If you want a more comprehensive introduction, an on-line manual from The MathWorks is available. Select Help in the toolbar atop the Matlab command window, then select MATLAB Help and Getting Started. A PDF version is available under Printable versions. The document is also available from The MathWorks Web site [10]. Many other manuals produced by The MathWorks are available on line and from the Web site. A list of over 600 Matlab based books by other authors and publishers, in sev- eral languages, is available at [11]. Three introductions to Matlab are of particular interest here, a relatively short primer by Sigmon and Davis [8], a medium-sized, mathematically oriented text by Higham and Higham [3], and a large, comprehen- sive manual by Hanselman and Littlefield [2]. You should have a copy of Matlab close at hand so you can run our sample programs as you read about them. All of the programs used in this book have been collected in a directory (or folder) named NCM (The directory name is the initials of the book title.) You can either start Matlab in this directory, or use pathtool to add the directory to the Matlab path. 1.1 The Golden Ratio What is the world’s most interesting number? Perhaps you like π, or e, or 17. Some people might vote for φ, the golden ratio, computed here by our first Matlab 1 2 Chapter 1. Introduction to MATLAB statement phi = (1 + sqrt(5))/2 This produces phi = 1.6180 Let’s see more digits. format long phi phi = 1.61803398874989 This didn’t recompute φ, it just displayed 15 significant digits instead of five. The golden ratio shows up in many places in mathematics; we’ll see several in this book. The golden ratio gets its name from the golden rectangle, shown in figure 1.1. The golden rectangle has the property that removing a square leaves a smaller rectangle with the same shape. φ φ − 1 1 1 Figure 1.1. The golden rectangle Equating the aspect ratios of the rectangles gives a defining equation for φ. 1 φ = φ− 1 1 This equation says that you can compute the reciprocal of φ by simply subtracting one. How many numbers have that property? Multiplying the aspect ratio equation by φ produces a polynomial equation φ2 − φ− 1 = 0 The roots of this equation are given by the quadratic formula. φ = 1±√5 2 1.1. The Golden Ratio 5 0 0.5 1 1.5 2 2.5 3 3.5 4 −3 −2 −1 0 1 2 3 4 5 6 7 x 1/x − (x−1) Figure 1.2. f(φ) = 0 hold on plot(phi,0,’o’) The following Matlab program produces the picture of the golden rectangle shown in figure 1.1. The program is contained in an M-file named goldrect.m, so issuing the command goldrect runs the script and creates the picture. % GOLDRECT Plot the golden rectangle phi = (1+sqrt(5))/2; x = [0 phi phi 0 0]; y = [0 0 1 1 0]; u = [1 1]; v = [0 1]; plot(x,y,’b’,u,v,’b--’) text(phi/2,1.05,’\phi’) text((1+phi)/2,-.05,’\phi - 1’) text(-.05,.5,’1’) text(.5,-.05,’1’) axis equal axis off set(gcf,’color’,’white’) 6 Chapter 1. Introduction to MATLAB The vectors x and y each contain five elements. Connecting consecutive (xk, yk) pairs with straight lines produces the outside rectangle. The vectors u and v each contain two elements. The line connecting (u1, v1) with (u2, v2) sepa- rates the rectangle into the square and the smaller rectangle. The plot command draws these lines, the x−y lines in solid blue and the u−v line in dashed blue. The next four statements place text at various points; the string ’\phi’ denotes the Greek letter. The two axis statements cause the scaling in the x and y directions to be equal and then turn off the display of the axes. The last statement sets the background color of gcf, which stands for get current figure, to white. A continued fraction is an infinite expression of the form a0 + 1 a1 + 1a2+ 1a3+... If all the ak’s are equal to 1, the continued fraction is another representation of the golden ratio. φ = 1 + 1 1 + 1 1+ 11+... The following Matlab function generates and evaluates truncated continued frac- tion approximations to φ. The code is stored in an M-file named goldfract.m. function goldfract(n) %GOLDFRACT Golden ratio continued fraction. % GOLDFRACT(n) displays n terms. p = ’1’; for k = 1:n p = [’1+1/(’ p ’)’]; end p p = 1; q = 1; for k = 1:n s = p; p = p + q; q = s; end p = sprintf(’%d/%d’,p,q) format long p = eval(p) format short err = (1+sqrt(5))/2 - p 1.1. The Golden Ratio 7 The statement goldfract(6) produces p = 1+1/(1+1/(1+1/(1+1/(1+1/(1+1/(1)))))) p = 21/13 p = 1.61538461538462 err = 0.0026 The three p’s are all different representations of the same approximation to φ. The first p is the continued fraction truncated to six terms. There are six right parentheses. This p is a string generated by starting with a single ’1’ (that’s goldfract(0)) and repeatedly inserting the string ’1+1/(’ in front and the string ’)’ in back. No matter how long this string becomes, it is a valid Matlab expression. The second p is an “ordinary” fraction with a single integer numerator and denominator obtained by collapsing the first p. The basis for the reformulation is 1 + 1 p q = p + q p So the iteration starts with 1 1 and repeatedly replaces the fraction p q by p + q p The statement p = sprintf(’%d/%d’,p,q) prints the final fraction by formatting p and q as decimal integers and placing a ’/’ between them. The third p is the same number as the first two p’s, but is represented as a conventional decimal expansion, obtained by having the Matlab eval function actually do the division expressed in the second p. 10 Chapter 1. Introduction to MATLAB The name of the function is in uppercase because historically Matlab was case insensitive and ran on terminals with only a single font. The use of capital letters may be confusing to some first-time Matlab users, but the convention persists. It is important to repeat the input and output arguments in these comments because the first line is not displayed when you ask for help on the function. The next line f = zeros(n,1); creates an n-by-1 matrix containing all zeros and assigns it to f. In Matlab, a matrix with only one column is a column vector and a matrix with only one row is a row vector. The next two lines, f(1) = 1; f(2) = 2; provide the initial conditions. The last three lines are the for statement that does all the work. for k = 3:n f(k) = f(k-1) + f(k-2); end We like to use three spaces to indent the body of for and if statements, but other people prefer two or four spaces, or a tab. You can also put the entire construction on one line if you provide a comma after the first clause. This particular function looks a lot like functions in other programming lan- guages. It produces a vector, but it does not use any of the Matlab vector or matrix operations. We will see some of these operations soon. Here is another Fibonacci function, fibnum.m. Its output is simply the nth Fibonacci number. function f = fibnum(n) % FIBNUM Fibonacci number. % FIBNUM(n) generates the nth Fibonacci number. if n <= 1 f = 1; else f = fibnum(n-1) + fibnum(n-2); end The statement fibnum(12) produces ans = 233 1.2. Fibonacci Numbers 11 The fibnum function is recursive. In fact, the term recursive is used in both a mathematical and a computer science sense. The relationship fn = fn−1 + fn−2 is known as a recursion relation and a function that calls itself is a recursive function. A recursive program is elegant, but expensive. You can measure execution time with tic and toc. Try tic, fibnum(24), toc Do not try tic, fibnum(50), toc Now compare the results produced by goldfract(6) and fibonacci(7). The first contains the fraction 21/13 while the second ends with 13 and 21. This is not just a coincidence. The continued fraction is collapsed by repeating the statement p = p + q; while the Fibonacci numbers are generated by f(k) = f(k-1) + f(k-2); In fact, if we let φn denote the golden ratio continued fraction truncated at n terms, then fn+1 fn = φn In the infinite limit, the ratio of successive Fibonacci numbers approaches the golden ratio. lim n→∞ fn+1 fn = φ To see this, compute 40 Fibonacci numbers: n = 40; f = fibonacci(n); Then compute their ratios: f(2:n)./f(1:n-1) This takes the vector containing f(2) through f(n) and divides it, element by element, by the vector containing f(1) through f(n-1). The output begins with 2.00000000000000 1.50000000000000 1.66666666666667 1.60000000000000 1.62500000000000 1.61538461538462 1.61904761904762 1.61764705882353 1.61818181818182 12 Chapter 1. Introduction to MATLAB and ends with 1.61803398874990 1.61803398874989 1.61803398874990 1.61803398874989 1.61803398874989 Do you see why we chose n = 40? Use the up arrow key on your keyboard to bring back the previous expression. Change it to f(2:n)./f(1:n-1) - phi and then press the Enter key. What is the value of the last element? The population of Fibonacci’s rabbit pen doesn’t double every month; it is multiplied by the golden ratio every month. It is possible to find a closed-form solution to the Fibonacci number recurrence relation. The key is to look for solutions of the form fn = cρn for some constants c and ρ. The recurrence relation fn = fn−1 + fn−2 becomes ρ2 = ρ + 1 We’ve seen this equation before. There are two possible values of ρ, namely φ and 1− φ. The general solution to the recurrence is fn = c1φn + c2(1− φ)n The constants c1 and c2 are determined by initial conditions, which are now conveniently written f0 = c1 + c2 = 1 f1 = c1φ + c2(1− φ) = 1 An exercise will ask you to use the Matlab backslash operator to solve this 2-by-2 system of simultaneous linear equations, but it is actually easier to solve the system by hand. c1 = φ 2φ− 1 c2 = − (1− φ)2φ− 1 Inserting these in the general solution gives fn = 1 2φ− 1(φ n+1 − (1− φ)n+1) 1.3. Fractal Fern 15 If you like the image, you might even choose to make it your computer desktop background. However, you should really run fern on your own computer to see the dynamics of the emerging fern in high resolution. The fern is generated by repeated transformations of a point in the plane. Let x be a vector with two components, x1 and x2, representing the point. There are four different transformations, all of them of the form x → Ax + b with different matrices A and vectors b. These are known as affine transformations. The most frequently used transformation has A = ( .85 .04 −.04 .85 ) , b = ( 0 1.6 ) This transformation shortens and rotates x a little bit, then adds 1.6 to its second component. Repeated application of this transformation moves the point up and to the right, heading towards the upper tip of the fern. Every once in a while, one of the other three transformations is picked at random. These transformations move the point into the lower subfern on the right, the lower subfern on the left, or the stem. Here is the complete fractal fern program. function fern %FERN MATLAB implementation of the Fractal Fern % Michael Barnsley, Fractals Everywhere, Academic Press,1993 % This version runs forever, or until stop is toggled. % See also: FINITEFERN. shg clf reset set(gcf,’color’,’white’,’menubar’,’none’, ... ’numbertitle’,’off’,’name’,’Fractal Fern’) x = [.5; .5]; h = plot(x(1),x(2),’.’); darkgreen = [0 2/3 0]; set(h,’markersize’,1,’color’,darkgreen,’erasemode’,’none’); axis([-3 3 0 10]) axis off stop = uicontrol(’style’,’toggle’,’string’,’stop’, ... ’background’,’white’); drawnow p = [ .85 .92 .99 1.00]; A1 = [ .85 .04; -.04 .85]; b1 = [0; 1.6]; A2 = [ .20 -.26; .23 .22]; b2 = [0; 1.6]; A3 = [-.15 .28; .26 .24]; b3 = [0; .44]; A4 = [ 0 0 ; 0 .16]; 16 Chapter 1. Introduction to MATLAB cnt = 1; tic while ~get(stop,’value’) r = rand; if r < p(1) x = A1*x + b1; elseif r < p(2) x = A2*x + b2; elseif r < p(3) x = A3*x + b3; else x = A4*x; end set(h,’xdata’,x(1),’ydata’,x(2)); cnt = cnt + 1; drawnow end t = toc; s = sprintf(’%8.0f points in %6.3f seconds’,cnt,t); text(-1.5,-0.5,s,’fontweight’,’bold’); set(stop,’style’,’pushbutton’,’string’,’close’, ... ’callback’,’close(gcf)’) Let’s examine this program a few statements at a time. shg stands for “show graph window.” It brings an existing graphics window forward, or creates a new one if necessary. clf reset resets most of the figure properties to their default values. set(gcf,’color’,’white’,’menubar’,’none’, ... ’numbertitle’,’off’,’name’,’Fractal Fern’) changes the background color of the figure window from the default gray to white and provides a customized title for the window. x = [.5; .5]; provides the initial coordinates of the point. h = plot(x(1),x(2),’.’); plots a single dot in the plane and saves a handle, h, so that we can later modify the properties of the plot. darkgreen = [0 2/3 0]; 1.3. Fractal Fern 17 defines a color where the red and blue components are zero and the green component is two-thirds of its full intensity. set(h,’markersize’,1,’color’,darkgreen,’erasemode’,’none’); makes the dot referenced by h smaller, changes its color, and specifies that the image of the dot on the screen should not be erased when its coordinates are changed. A record of these old points is kept by the computer’s graphics hardware (until the figure is reset), but Matlab itself does not remember them. axis([-3 3 0 10]) axis off specifies that the plot should cover the region −3 ≤ x1 ≤ 3, 0 ≤ x2 ≤ 10 but that the axes should not be drawn. stop = uicontrol(’style’,’toggle’,’string’,’stop’, ... ’background’,’white’); creates a toggle user interface control, labeled ’stop’ and colored white, in the default position near the lower left corner of the figure. The handle for the control is saved in the variable stop. drawnow causes the initial figure, including the initial point, to actually be plotted on the computer screen. The statement p = [ .85 .92 .99 1.00]; sets up a vector of probabilities. The statements A1 = [ .85 .04; -.04 .85]; b1 = [0; 1.6]; A2 = [ .20 -.26; .23 .22]; b2 = [0; 1.6]; A3 = [-.15 .28; .26 .24]; b3 = [0; .44]; A4 = [ 0 0 ; 0 .16]; define the four affine transformations. The statement cnt = 1; initializes a counter that keeps track of the number of points plotted. The statement tic initializes a stopwatch timer. The statement while ~get(stop,’value’) 20 Chapter 1. Introduction to MATLAB 15 The opposite diagonal, which runs from upper right to lower left, is less important in linear algebra, so finding its sum is a little trickier. One way to do it makes use of the function that “flips” a matrix “up-down”: sum(diag(flipud(A))) produces 15 This verifies that A has equal row, column, and diagonal sums. Why is the magic sum equal to 15? The command sum(1:9) tells us that the sum of the integers from 1 to 9 is 45. If these integers are allocated to three columns with equal sums, that sum must be sum(1:9)/3 which is 15. There are eight possible ways to place a transparency on an overhead projec- tor. Similarly, there are eight magic squares of order three that are rotations and reflections of A. The statements for k = 0:3 rot90(A,k) rot90(A’,k) end display all eight of them. 8 1 6 8 3 4 3 5 7 1 5 9 4 9 2 6 7 2 6 7 2 4 9 2 1 5 9 3 5 7 8 3 4 8 1 6 2 9 4 2 7 6 7 5 3 9 5 1 6 1 8 4 3 8 4 3 8 6 1 8 9 5 1 7 5 3 2 7 6 2 9 4 This is all the magic squares of order three. Now for some linear algebra. The determinant of our magic square, 1.4. Magic Squares 21 det(A) is -360 The inverse, X = inv(A) is X = 0.1472 -0.1444 0.0639 -0.0611 0.0222 0.1056 -0.0194 0.1889 -0.1028 The inverse looks better if it is displayed with a rational format. format rat X shows that the elements of X are fractions with det(A) in the denominator. X = 53/360 -13/90 23/360 -11/180 1/45 19/180 -7/360 17/90 -37/360 The statement format short restores the output format to its default. Three other important quantities in computational linear algebra are matrix norms, eigenvalues, and singular values. The statements r = norm(A) e = eig(A) s = svd(A) produce r = 15 e = 15.0000 4.8990 -4.8990 s = 15.0000 6.9282 3.4641 22 Chapter 1. Introduction to MATLAB The magic sum occurs in all three because the vector of all ones is an eigenvector, and is also a left and right singular vector. So far, all the computations in this section have been done using floating-point arithmetic. This is the arithmetic used for almost all scientific and engineering computation, especially for large matrices. But for a 3-by-3 matrix, it is easy to repeat the computations using symbolic arithmetic and the Symbolic Toolbox connection to Maple. The statement A = sym(A) changes the internal representation of A to a symbolic form that is displayed as A = [ 8, 1, 6] [ 3, 5, 7] [ 4, 9, 2] Now commands like sum(A), sum(A’)’, det(A), inv(A), eig(A), svd(A) produce symbolic results. In particular, the eigenvalue problem for this matrix can be solved exactly, and e = [ 15] [ 2*6^(1/2)] [ -2*6^(1/2)] A 4-by-4 magic square is one of several mathematical objects on display in Melancolia, a Renaissance etching by Albrect Durer. An electronic copy of the etching is available in a Matlab data file. load durer whos produces X 648x509 2638656 double array caption 2x28 112 char array map 128x3 3072 double array The elements of the matrix X are indices into the gray-scale color map named map. The image is displayed with image(X) colormap(map) axis image Click the magnifying glass with a “+” in the toolbar and use the mouse to zoom in on the magic square in the upper right-hand corner. The scanning resolution becomes evident as you zoom in. The commands 1.4. Magic Squares 25 0 5 10 15 20 25 0 5 10 15 20 25 Rank of magic squares Figure 1.4. Rank of magic squares • Singly even order; n is a multiple of 2, but not 4. • Doubly even order; n is a multiple of 4. Odd-ordered magic squares, n = 3, 5, 7, ..., have full rank n. They are nonsingular and have inverses. Doubly even magic squares, n = 4, 8, 12, ..., have rank three no matter how large n is. They might be called very singular. Singly even magic squares, n = 6, 10, 14, ..., have rank n/2 + 2. They are also singular, but have fewer row and column dependencies than the doubly even squares. If you have Matlab Version 6 or later, you can look at the M-file that gener- ates magic squares with edit magic.m or type magic.m You will see the three different cases in the code. The different kinds of magic squares also produce different three-dimensional surface plots. Try the following for various values of n. surf(magic(n)) axis off set(gcf,’doublebuffer’,’on’) cameratoolbar Double buffering prevents flicker when you use the various camera tools to move the viewpoint. The following code produces figure 1.5. 26 Chapter 1. Introduction to MATLAB 8 9 10 11 Figure 1.5. Surface plots of magic squares for n = 8:11 subplot(2,2,n-7) surf(magic(n)) title(num2str(n)) axis off view(30,45) axis tight end 1.5 Cryptography This section uses a cryptography example to show how Matlab deals with text and character strings. The cryptographic technique, which is known as a Hill cipher, involves arithmetic in a finite field. Almost all modern computers use the ASCII character set to store basic text. ASCII stands for American Standard Code for Information Interchange. The char- acter set uses seven of the eight bits in a byte to encode 128 characters. The first 32 characters are nonprinting control characters, such as tab, backspace, and end-of- line. The 128th character is another nonprinting character that corresponds to the Delete key on your keyboard. In between these control characters are 95 printable characters, including a space, 10 digits, 26 lowercase letters, 26 uppercase letters, 1.5. Cryptography 27 and 32 punctuation marks. Matlab can easily display all the printable characters, in the order deter- mined by their ASCII encoding. Start with x = reshape(32:127,32,3)’ This produces a 3-by-32 matrix. x = 32 33 34 ... 61 62 63 64 65 66 ... 93 94 95 96 97 98 ... 125 126 127 The char function converts numbers to characters. The statement c = char(x) produces c = !"#$%&’()*+,-./0123456789:;<=>? @ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_ ‘abcdefghijklmnopqrstuvwxyz{|}~ We have cheated a little bit because the last element of x is 127, which corresponds to the nonprinting delete character, and we have not shown the last character in c. You can try this on your computer and see what is actually displayed. The first character in c is blank, indicating that char(32) is the same as ’ ’ The last printable character in c is the tilde, indicating that char(126) is the same as ’~’ The characters representing digits are in the first line of c. In fact d = char(48:57) displays a ten-character string d = 0123456789 This string can be converted to the corresponding numerical values with double or real. The statement 30 Chapter 1. Introduction to MATLAB A comment precedes the statement that assigns the prime p. % Use a two-character Hill cipher with arithmetic % modulo 97, a prime. p = 97; Choose two characters above ASCII 128 to expand the size of the character set from 95 to 97. c1 = char(169); c2 = char(174); x(x==c1) = 127; x(x==c2) = 128; The conversion from characters to numerical values is done by x = mod(real(x-32),p); Prepare for the matrix-vector product by forming a matrix with two rows and lots of columns. n = 2*floor(length(x)/2); X = reshape(x(1:n),2,n/2); All this preparation has been so that we can do the actual finite field arithmetic quickly and easily. % Encode with matrix multiplication modulo p. A = [71 2; 2 26]; Y = mod(A*X,p); Reshape into a single row. y = reshape(Y,1,n); If length(x) is odd, encode the last character. if length(x) > n y(n+1) = mod((p-1)*x(n+1),p); end Finally, convert the numbers back to characters. y = char(y+32); y(y==127) = c1; y(y==128) = c2; Let’s follow the computation of y = crypto(’Hello world’). We begin with a character string. x = ’Hello world’ This is converted to an integer vector. 1.6. The 3n+1 Sequence 31 x = 40 69 76 76 79 0 87 79 82 76 68 The length(x) is odd, so the reshaping temporarily ignores the last element. X = 40 76 79 87 82 69 76 0 79 76 A conventional matrix-vector multiplication A*X produces an intermediate matrix. 2978 5548 5609 6335 5974 1874 2128 158 2228 2140 Then the mod(.,p) operation produces Y = 68 19 80 30 57 31 91 61 94 6 This is rearranged to a row vector. y = 68 31 19 91 80 61 30 94 57 6 Now the last element of x is encoded by itself and attached to the end of y. y = 68 31 19 91 80 61 30 94 57 6 29 Finally, y is converted back to a character string to produce the encrypted result. y = ’d?3{p]>~Y&=’ If we now compute crypto(y), we get back our original ’Hello world’. 1.6 The 3n+1 Sequence This section describes a famous unsolved problem in number theory. Start with any positive integer n. Repeat the following steps: • If n = 1, stop. • If n is even, replace it with n/2. • If n is odd, replace it with 3n + 1. For example, starting with n = 7 produces 7, 22, 11, 34, 17, 52, 26, 13, 40, 20, 10, 5, 16, 8, 4, 2, 1 The sequence terminates after 17 steps. Note that whenever n reaches a power of 2, the sequence terminates in log2 n more steps. 32 Chapter 1. Introduction to MATLAB The unanswered question is, does the process always terminate? Or is there some starting value that causes the process to go on forever, either because the numbers get larger and larger, or because some periodic cycle is generated? This problem is known as the 3n + 1 problem. It has been studied by many eminent mathematicians, including Collatz, Ulam, and Kakatani, and is discussed in a survey paper by Jeffrey Lagarias [5]. The following Matlab code fragment generates the sequence starting with any specified n. y = n; while n > 1 if rem(n,2)==0 n = n/2; else n = 3*n+1; end y = [y n]; end We don’t know ahead of time how long the resulting vector y is going to be. But the statement y = [y n]; automatically increases length(y) each time it is executed. In principle, the unsolved mathematical problem is: can this code fragment run forever? In actual fact, floating-point roundoff error causes the calculation to misbehave whenever 3n + 1 becomes greater than 253, but it is still interesting to investigate modest values of n. Let’s embed our code fragment in a GUI. The complete function is in M-file threenplus1.m. For example, the statement threenplus1(7) produces figure 1.6. The M-file begins with a preamble containing the function header and the help information. function threenplus1(n) % ‘‘Three n plus 1’’. % Study the 3n+1 sequence. % threenplus1(n) plots the sequence starting with n. % threenplus1 with no arguments starts with n = 1. % uicontrols decrement or increment the starting n. % Is it possible for this to run forever? The next section of code brings the current graphics window forward and resets it. Two pushbuttons, which are the default uicontrols, are positioned near the bottom center of the figure at pixel coordinates [260,5] and [300,5]. Their size is 25-by- 22 pixels and they are labeled with ’<’ and ’>’. If either button is subsequently 1.7. Floating-Point Arithmetic 35 underflow, and overflow. Most of the time, it is possible to use Matlab effectively without worrying about these details, but every once in a while, it pays to know something about the properties and limitations of floating-point numbers. Twenty years ago, the situation was far more complicated than it is today. Each computer had its own floating-point number system. Some were binary; some were decimal. There was even a Russian computer that used trinary arithmetic. Among the binary computers, some used 2 as the base; others used 8 or 16. And everybody had a different precision. In 1985, the IEEE Standards Board and the American National Standards Institute adopted the ANSI/IEEE Standard 754-1985 for Binary Floating-Point Arithmetic. This was the culmination of almost a decade of work by a 92-person working group of mathematicians, computer scientists, and engineers from universities, computer manufacturers, and microprocessor compa- nies. All computers designed since 1985 use IEEE floating-point arithmetic. This doesn’t mean that they all get exactly the same results, because there is some flexibility within the standard. But it does mean that we now have a machine- independent model of how floating-point arithmetic behaves. Matlab has traditionally used the IEEE double-precision format. There is a single-precision format that saves space, but that isn’t much faster on modern machines. Matlab 7 will have support for single-precision arithmetic, but we will deal exclusively with double-precision in this book. There is also an extended precision format, which is optional and therefore is one of the reasons for lack of uniformity among different machines. Most nonzero floating-point numbers are normalized. This means they can be expressed as x = ±(1 + f) · 2e The quantity f is the fraction or mantissa and e is the exponent. The fraction satisfies 0 ≤ f < 1 and must be representable in binary using at most 52 bits. In other words, 252f is an integer in the interval 0 ≤ 252f < 252 The exponent e is an integer in the interval −1022 ≤ e ≤ 1023 The finiteness of f is a limitation on precision. The finiteness of e is a limitation on range. Any numbers that don’t meet these limitations must be approximated by ones that do. Double-precision floating-point numbers are stored in a 64 bit word, with 52 bits for f , 11 bits for e, and one bit for the sign of the number. The sign of e is accommodated by storing e + 1023, which is between 1 and 211 − 2. The two extreme values for the exponent field, 0 and 211 − 1, are reserved for exceptional floating-point numbers that we will describe later. 36 Chapter 1. Introduction to MATLAB The entire fractional part of a floating-point number is not f , but 1+f , which has 53 bits. However the leading 1 doesn’t need to be stored. In effect, the IEEE format packs 65 bits of information into a 64 bit word. The program floatgui shows the distribution of the positive numbers in a model floating-point system with variable parameters. The parameter t specifies the number of bits used to store f . In other words, 2tf is an integer. The parameters emin and emax specify the range of the exponent, so emin ≤ e ≤ emax. Initially, floatgui sets t = 3, emin = −4, and emax = 3 and produces the distribution shown in figure 1.7. 1/16 1/2 1 2 4 8−1/2 ||||||||||||||||||||| | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | | Figure 1.7. floatgui Within each binary interval, 2e ≤ x ≤ 2e+1, the numbers are equally spaced, with an increment of 2e−t. If e = 0 and t = 3, for example, the spacing of the numbers between 1 and 2 is 1/8. As e increases, the spacing increases. It is also instructive to display the floating-point numbers with a logarithmic scale. Figure 1.8 shows floatgui with logscale checked and t = 5, emin = −4 and emax = 3. With this logarithmic scale, it is more apparent that the distribution in each binary interval is the same. 1/16 1/8 1/4 1/2 1 2 4 8 16−1/4 |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| Figure 1.8. floatgui(log scale) A very important quantity associated with floating-point arithmetic is high- lighted in red by floatgui. Matlab calls this quantity eps, which is short for machine epsilon. eps is the distance from 1 to the next larger floating-point number. For the floatgui model floating-point system, eps = 2^(-t). Before the IEEE standard, different machines had different values of eps. Now, for IEEE double-precision, eps = 2^(-52) 1.7. Floating-Point Arithmetic 37 The approximate decimal value of eps is 2.2204 ·10−16. Either eps/2 or eps can be called the roundoff level. The maximum relative error incurred when the result of an arithmetic operation is rounded to the nearest floating-point number is eps/2. The maximum relative spacing between numbers is eps. In either case, you can say that the roundoff level is about 16 decimal digits. A frequent instance of roundoff occurs with the simple Matlab statement t = 0.1 The mathematical value t stored in t is not exactly 0.1 because expressing the decimal fraction 1/10 in binary requires an infinite series. In fact, 1 10 = 1 24 + 1 25 + 0 26 + 0 27 + 1 28 + 1 29 + 0 210 + 0 211 + 1 212 + ... After the first term, the sequence of coefficients 1, 0, 0, 1 is repeated infinitely often. Grouping the resulting terms together four at a time expresses 1/10 in a base 16, or hexadecimal, series. 1 10 = 2−4 · (1 + 9 16 + 9 162 + 9 163 + 9 164 + . . .) Floating-point numbers on either side of 1/10 are obtained by terminating the fractional part of this series after 52 binary terms, or 13 hexadecimal terms, and rounding the last term up or down. Thus t1 < 1/10 < t2 Where t1 = 2−4 · (1 + 916 + 9 162 + 9 163 + ... + 9 1612 + 9 1613 ) t2 = 2−4 · (1 + 916 + 9 162 + 9 163 + ... + 9 1612 + 10 1613 ) It turns out that 1/10 is closer to t2 than to t1, so t is equal to t2. In other words, t = (1 + f) · 2e where f = 9 16 + 9 162 + 9 163 + ... + 9 1612 + 10 1613 e = −4 The Matlab command format hex causes t to be displayed as 3fb999999999999a 40 Chapter 1. Introduction to MATLAB elegant way to handle underflow, but their practical importance for Matlab style computation is very rare. Denormal numbers are represented by taking e = −1023, so the biased exponent e + 1023 is 0. Matlab uses the floating-point system to handle integers. Mathematically, the numbers 3 and 3.0 are the same, but many programming languages would use different representations for the two. Matlab does not distinguish between them. We sometimes use the term flint to describe a floating-point number whose value is an integer. Floating-point operations on flints do not introduce any roundoff error, as long as the results are not too large. Addition, subtraction, and multiplication of flints produce the exact flint result, if it is not larger than 253. Division and square root involving flints also produce a flint if the result is an integer. For example, sqrt(363/3) produces 11, with no roundoff. Two Matlab functions that take apart and put together floating-point num- bers are log2 and pow2. help log2 help pow2 produces [F,E] = LOG2(X) for a real array X, returns an array F of real numbers, usually in the range 0.5 <= abs(F) < 1, and an array E of integers, so that X = F .* 2.^E. Any zeros in X produce F = 0 and E = 0. X = POW2(F,E) for a real array F and an integer array E computes X = F .* (2 .^ E). The result is computed quickly by simply adding E to the floating-point exponent of F. The quantities F and E used by log2 and pow2 predate the IEEE floating-point standard and so are slightly different from the f and e we are using in this section. In fact, f = 2*F-1 and e = E-1. [F,E] = log2(phi) produces F = 0.80901699437495 E = 1 Then phi = pow2(F,E) gives back phi = 1.61803398874989 1.7. Floating-Point Arithmetic 41 As an example of how roundoff error affects matrix computations, consider the two-by-two set of linear equations 17x1 + 5x2 = 22 1.7x1 + 0.5x2 = 2.2 The obvious solution is x1 = 1, x2 = 1. But the Matlab statements A = [17 5; 1.7 0.5] b = [22; 2.2] x = A\b produce x = -1.0588 8.0000 Where did this come from? Well, the equations are singular, but consistent. The second equation is just 0.1 times the first. The computed x is one of infinitely many possible solutions. But the floating-point representation of the matrix A is not exactly singular because A(2,1) is not exactly 17/10. The solution process subtracts a multiple of the first equation from the second. The multiplier is mu = 1.7/17, which turns out to be the floating-point number obtained by truncating, rather than rounding, the binary expansion of 1/10. The matrix A and the right-hand side b are modified by A(2,:) = A(2,:) - mu*A(1,:) b(2) = b(2) - mu*b(1) With exact computation, both A(2,2) and b(2) would become zero, but with floating-point arithmetic, they both become nonzero multiples of eps. A(2,2) = (1/4)*eps = 5.5511e-17 b(2) = 2*eps = 4.4408e-16 Matlab notices the tiny value of the new A(2,2) and displays a message warning that the matrix is close to singular. It then computes the solution of the modified second equation by dividing one roundoff error by another. x(2) = b(2)/A(2,2) = 8 This value is substituted back into the first equation to give x(1) = (22 - 5*x(2))/17 = -1.0588 42 Chapter 1. Introduction to MATLAB 0.985 0.99 0.995 1 1.005 1.01 1.015 −5 −4 −3 −2 −1 0 1 2 3 4 5 x 10 −14 Figure 1.9. Is this a polynomial? The details of the roundoff error lead Matlab to pick out one particular solution from among the infinitely many possible solutions to the singular system. Our final example plots a seventh-degree polynomial. x = 0.988:.0001:1.012; y = x.^7-7*x.^6+21*x.^5-35*x.^4+35*x.^3-21*x.^2+7*x-1; plot(x,y) The resulting plot doesn’t look anything like a polynomial. It isn’t smooth. You are seeing roundoff error in action. The y-axis scale factor is tiny, 10−14. The tiny values of y are being computed by taking sums and differences of numbers as large as 35 · 1.0124 . There is severe subtractive cancellation. The example was contrived by using the Symbolic Toolbox to expand (x− 1)7 and carefully choosing the range for the x-axis to be near x = 1. If the values of y are computed instead by y = (x-1).^7; then a smooth (but very flat) plot results. 1.8 Further Reading Additional information about floating-point arithmetic and roundoff error can be found in Higham [4] and Overton [6]. Exercises 45 plot(...) hold on plot(...) plot(...) hold off 1.16. Use the following code to make your own Portable Network Graphics file from the fern. Then compare your image with one obtained from ncm/fern.png. bg = [0 0 85]; % Dark blue background fg = [255 255 255]; % White dots sz = get(0,’screensize’); rand(’state’,0) X = finitefern(500000,sz(4),sz(3)); d = fg - bg; R = uint8(bg(1) + d(1)*X); G = uint8(bg(2) + d(2)*X); B = uint8(bg(3) + d(3)*X); F = cat(3,R,G,B); imwrite(F,’myfern.png’,’png’,’bitdepth’,8) 1.17. Modify fern.m or finitefern.m so that it produces Sierpinski’s triangle. Start at x = ( 0 0 ) At each iterative step the current point x is replaced by Ax + b where the matrix A is always A = ( 1/2 0 0 1/2 ) and the vector b is chosen at random with equal probability from among the three vectors b = ( 0 0 ) , b = ( 1/2 0 ) , or b = ( 1/4√ 3/4 ) 1.18. greetings(phi) generates a seasonal holiday fractal that depends upon the parameter φ. The default value of φ is the golden ratio. What happens for other values of φ? Try both simple fractions and floating point approxima- tions to irrational values. 1.19. A = magic(4) is singular. Its columns are linearly dependent. What do null(A), null(A,’r’), null(sym(A)), and rref(A) tell you about that de- pendence? 46 Chapter 1. Introduction to MATLAB 1.20. Let A = magic(n) for n = 3, 4, or 5. What does p = randperm(n); q = randperm(n); A = A(p,q); do to sum(A) sum(A’)’ sum(diag(A)) sum(diag(flipud(A))) rank(A) 1.21. The character char(7) is a control character. What does it do? 1.22. What does char([169 174]) display on your computer? 1.23. What fundamental physical law is hidden in this string? s = ’/b_t3{$H~MO6JTQI>v~#3GieW*l(p,nF’ 1.24. Find the two files encrypt.m and gettysburg.txt. Use encrypt to encrypt gettysburg.txt. Then decrypt the result. Use encrypt to encrypt itself. 1.25. With the NCM directory on you path, you can read the text of Lincoln’s Gettysburg address with fp = fopen(’gettysburg.txt’); G = char(fread(fp))’ fclose(fp); (a) How many characters are in the text? (b) Use the unique function to find the unique characters in the text. (c) How many blanks are in the text? What punctuation characters, and how many of each, are there? (d) Remove the blanks and the punctuation and convert the text to all upper or lower case. Use the histc function to count the number of letters. What is the most frequent letter? What letters are missing? (e) Use the bar function as described in help histc to plot a histogram of the letter frequencies. (f) Use get(gca,’xtick’) and get(gca,’xticklabel’) to see how the x- axis of the histogram is labeled. Then use set(gca,’xtick’,...,’xticklabel’,...) to relabel the x-axis with the letters in the text. 1.26. If x is the character string consisting of just two blanks, x = ’ ’ Exercises 47 then crypto(x) is actually equal to x. Why does this happen? Are there any other two-character strings that crypto does not change? 1.27. Find another 2-by-2 integer matrix A for which mod(A*A,97) is the identity matrix. Replace the matrix in crypto.m with your matrix and verify that the function still works correctly. 1.28. The function crypto works with 97 characters instead of 95. It can produce output, and correctly handle input, that contains two characters with ASCII values greater than 127. What are these characters? Why are they necessary? What happens to other characters with ASCII values greater than 127? 1.29. Create a new crypto function that works with just 29 characters, the 26 lowercase letters, plus blank, period, and comma. You will need to find a 2-by-2 integer matrix A for which mod(A*A,29) is the identity matrix. 1.30. The graph of the 3n + 1 sequence has a particular characteristic shape if the starting n is 5, 10, 20, 40, . . ., that is, n is five times a power of 2. What is this shape and why does it happen? 1.31. The graphs of the 3n+1 sequences starting at n = 108, 109, and 110 are very similar to each other. Why? 1.32. Let L(n) be the number of terms in the 3n + 1 sequence that starts with n. Write a Matlab function that computes L(n) without using any vectors or unpredictable amounts of storage. Plot L(n) for 1 ≤ n ≤ 1000. What is the maximum value of L(n) for n in this range, and for what value of n does it occur? Use threenplus1 to plot the sequence that starts with this particular value of n. 1.33. Modify floatgui.m by changing its last line from a comment to an executable statement and changing the question mark to a simple expression that counts the number of floating-point numbers in the model system. 1.34. Explain the output produced by t = 0.1 n = 1:10 e = n/10 - n*t 1.35. What does each of these programs do? How many lines of output does each program produce? What are the last two values of x printed? x = 1; while 1+x > 1, x = x/2, pause(.02), end x = 1; while x+x > x, x = 2*x, pause(.02), end x = 1; while x+x > x, x = x/2, pause(.02), end 50 Chapter 1. Introduction to MATLAB 43 44 45 46 47 48 49 42 21 22 23 24 25 26 41 20 7 8 9 10 27 40 19 6 1 2 11 28 39 18 5 4 3 12 29 38 17 16 15 14 13 30 37 36 35 34 33 32 31 Figure 1.10. primespiral(7) piecewise quadratic function of i and j, so each diagonal segment represents a little mini-theorem about the distribution of primes. The phenomenon was discovered by Stanislaw Ulam in 1963 and appeared on the cover of Scientific American in 1964. There are a number of interesting Web pages devoted to prime spirals. Start with [7] and [9]. (a) The Matlab demos directory contains an M-file spiral.m. The integers from 1 to n2 are arranged in a spiral pattern, starting in the center of the matrix. The code in demos/spiral.m is not very elegant. Here is a better version. function S = spiral(n) %SPIRAL SPIRAL(n) is an n-by-n matrix with elements % 1:n^2 arranged in a rectangular spiral pattern. S = []; for m = 1:n S = rot90(S,2); S(m,m) = 0; p = ??? v = (m-1:-1:0); S(:,m) = p-v’; S(m,:) = p+v; end if mod(n,2)==1 S = rot90(S,2); end What value should be assigned to p each time through the loop so that this function generates the same matrices as the spiral.m in the demos directory? (b) Why do half of the diagonals of spiral(n) contain no primes? (c) Let S = spiral(2*n) and let r1 and r2 be rows that go nearly halfway Exercises 51 0 50 100 150 200 250 0 50 100 150 200 250 nz = 6275 Figure 1.11. primespiral(250) across the middle of the matrix r1 = S(n+1,1:n-2) r2 = S(n-1,n+2:end). Why do these rows contain no primes? (d) What is particularly remarkable about primespiral(17,17) primespiral(41,41) (e) Find values of n and c, both less than 50, and not equal to 17 or 41, so that [S,P] = primespiral(n,c) contains a diagonal segment with 8 or more primes. 1.42. Triangular numbers are integers of the form n(n + 1)/2. The term comes from the fact that a triangular grid with n points on a side has a total of n(n + 1)/2 points. Write a function trinums(m) that generates all the triangular numbers less than or equal to m. Modify primespiral to use your trinums and become trinumspiral. 52 Chapter 1. Introduction to MATLAB 1.43. Here is a puzzle that does not have much to do with this chapter, but you might find it interesting nevertheless. What familiar property of the integers is represented by the following plot? 0 10 20 30 40 50 60 70 80 90 100 0 2 4 6 8 1.44. In the Gregorian calendar, a year y is a leap year if and only if (mod(y,4) == 0) & (mod(y,100) ~= 0) | (mod(y,400) == 0) Thus, 2000 was a leap year, but 2100 will not be a leap year. This rule implies that the Gregorian calendar repeats itself every 400 years. In that 400 year period, there are 97 leap years, 4800 months, 20871 weeks and 146097 days. The Matlab functions datenum, datevec, datestr, and weekday use these facts to facilitate computations involving calendar dates. For example, either of the statements [d,w] = weekday(’Aug. 17, 2003’) or [d,w] = weekday(datenum([2003 8 17])) tells me that my birthday was on a Sunday in 2003. Use Matlab to answer the following questions. (a) On which day of the week were you born? (b) In a 400 year Gregorian calendar cycle, which week day is the most likely for your birthday? (c) What is the probability that the 13th of any month falls on a Friday? The answer is close to, but not exactly equal to, 1/7. 1.45. Biorhythms were very popular in the ’60’s. You can still find many Web sites today that offer to prepare personalized biorhythms, or that sell software to compute them. Biorhythms are based on the notion that three sinusoidal cycles influence our lives. The physical cycle has a period of 23 days, the emotional cycle has a period of 28 days, and the intellectual cycle has a period of 33 days. For any individual, the cycles are initialized at birth. Figure 1.12 is my biorhythm, which begins on Aug. 17, 1939, plotted for a eight-week period centered around the date this is being written, Oct. 19, 2003. It shows that my intellectual power reached a peak yesterday, that my physical strength and Bibliography [1] M. Barnsley, Fractals Everywhere, Academic Press, 1993. [2] D. C. Hanselman and B. Littlefield, Mastering MATLAB 6, A Compre- hensive Tutorial and Reference, Prentice-Hall, 2000, 832 pages. [3] D. J. Higham and N. J. Higham, MATLAB Guide, SIAM, 2000, 283 pages. [4] N. J. Higham, Accuracy and Stability of Numerical Algorithms, SIAM, 2002, 680 pages. [5] J. Lagarias, The 3x+1 problem and its generalizations, American Mathemat- ical Monthly, 92 (1985), pp. 3–23. http://www.cecm.sfu.ca/organics/papers/lagarias [6] M. Overton, Numerical Computing with IEEE Floating Point Arithmetic, SIAM, 2001, 104 pages. [7] Ivars Peterson, Prime Spirals, Science News Online, 161 (2002). http://www.sciencenews.org/20020504/mathtrek.asp [8] K. Sigmon and T. A. Davis, MATLAB Primer, Sixth Edition, Chapman and Hall/CRC, 2002, 176 pages. [9] Eric Weisstein, World of Mathematics, Prime Spiral, http://mathworld.wolfram.com/PrimeSpiral.html [10] MathWorks, The., Getting Started with MATLAB. http://www.mathworks.com/access/helpdesk/help/techdoc ... /learn_matlab/learn_matlab.shtml [11] MathWorks, The., List of Matlab based books. http://www.mathworks.com/support/books/index.jsp 55 Chapter 2 Linear Equations One of the problems encountered most frequently in scientific computation is the solution of systems of simultaneous linear equations. This chapter covers the solu- tion of linear systems by Gaussian elimination and the sensitivity of the solution to errors in the data and roundoff errors in the computation. 2.1 Solving Linear Systems With matrix notation, a system of simultaneous linear equations is written Ax = b In the most frequent case when there are as many equations as unknowns, A is a given square matrix of order n, b is a given column vector of n components, and x is an unknown column vector of n components. Students of linear algebra learn that the solution to Ax = b can be written x = A−1b where A−1 is the inverse of A. However, in the vast majority of practical computational problems, it is unnecessary and inadvisable to actually compute A−1. As an extreme but illustrative example, consider a system consisting of just one equation, such as 7x = 21 The best way to solve such a system is by division, x = 21 7 = 3 Use of the matrix inverse would lead to x = 7−1 × 21 = .142857× 21 = 2.99997 The inverse requires more arithmetic — a division and a multiplication instead of just a division — and produces a less accurate answer. Similar considerations apply to systems of more than one equation. This is even true in the common situation 1 2 Chapter 2. Linear Equations where there are several systems of equations with the same matrix A but different right-hand sides b. Consequently, we shall concentrate on the direct solution of systems of equations rather than the computation of the inverse. 2.2 The MATLAB Backslash Operator To emphasize the distinction between solving linear equations and computing in- verses, Matlab has introduced nonstandard notation using backward slash and forward slash operators, “\” and “/”. If A is a matrix of any size and shape and B is a matrix with as many rows as A, then the solution to the system of simultaneous equations AX = B is denoted by X = A\B Think of this as dividing both sides of the equation by the coefficient matrix A. Because matrix multiplication is not commutative and A occurs on the left in the original equation, this is left division. Similarly, the solution to a system with A on the right and B with as many columns as A, XA = B is obtained by right division, X = B/A This notation applies even if A is not square, so that the number of equations is not the same as the number of unknowns. However, in this chapter, we limit ourselves to systems with square coefficient matrices. 2.3 A 3-by-3 Example To illustrate the general linear equation solution algorithm, consider an example of order three:  10 −7 0 −3 2 6 5 −1 5     x1 x2 x3   =   7 4 6   This of course, represents the three simultaneous equations 10x1 − 7x2 = 7 −3x1 + 2x2 + 6x3 = 4 5x1 − x2 + 5x3 = 6 The first step of the solution algorithm uses the first equation to eliminate x1 from the other equations. This is accomplished by adding 0.3 times the first equation 2.5. LU Factorization 5 An upper triangular matrix has all its nonzero elements above or on the main diagonal. A unit lower triangular matrix has ones on the main diagonal and all the rest of its nonzero elements below the main diagonal. For example, U =   1 2 3 4 0 5 6 7 0 0 8 9 0 0 0 10   is upper triangular, and L =   1 0 0 0 2 1 0 0 3 5 1 0 4 6 7 1   is unit lower triangular. Linear equations involving triangular matrices are also easily solved. There are two variants of the algorithm for solving an n-by-n upper triangular system, Ux = b. Both begin by solving the last equation for the last variable, then the next to last equation for the next to last variable, and so on. One subtracts multiples of the columns of U from b: x = zeros(n,1); for k = n:-1:1 x(k) = b(k)/U(k,k); i = (1:k-1)’; b(i) = b(i) - x(k)*U(i,k); end The other uses inner products between the rows of U and portions of the emerging solution x: x = zeros(n,1); for k = n:-1:1 j = k+1:n; x(k) = (b(k) - U(k,j)*x(j))/U(k,k); end 2.5 LU Factorization The algorithm that is almost universally used to solve square systems of simultane- ous linear equations is one of the oldest numerical methods, the systematic elimi- nation method, generally named after C. F. Gauss. Research in the period 1955 to 1965 revealed the importance of two aspects of Gaussian elimination that were not emphasized in earlier work: the search for pivots and the proper interpretation of the effect of rounding errors. In general, Gaussian elimination has two stages, the forward elimination and the back substitution. The forward elimination consists of n − 1 steps. At the kth 6 Chapter 2. Linear Equations step, multiples of the kth equation are subtracted from the remaining equations to eliminate the kth variable. If the coefficient of xk is “small,” it is advisable to interchange equations before this is done. The elimination steps can be simultane- ously applied to the right-hand side, or the interchanges and multipliers saved and applied to the right-hand side later. The back substitution consists of solving the last equation for xn, then the next-to-last equation for xn−1, and so on, until x1 is computed from the first equation. Let Pk, k = 1, · · · , n − 1, denote the permutation matrix obtained by in- terchanging the rows of the identity matrix in the same way the rows of A are interchanged at the kth step of the elimination. Let Mk denote the unit lower tri- angular matrix obtained by inserting the negatives of the multipliers used at the kth step below the diagonal in the kth column of the identity matrix. Let U be the final upper triangular matrix obtained after the n− 1 steps. The entire process can be described by one matrix equation, U = Mn−1Pn−1 · · ·M2P2M1P1A It turns out that this equation can be rewritten L1L2 · · ·Ln−1U = Pn−1 · · ·P2P1A where Lk is obtained from Mk by permuting and changing the signs of the multi- pliers below the diagonal. So, if we let L = L1L2 · · ·Ln−1 P = Pn−1 · · ·P2P1 then we have LU = PA The unit lower triangular matrix L contains all the multipliers used during the elimination and the permutation matrix P accounts for all the interchanges. For our example A =   10 −7 0 −3 2 6 5 −1 5   the matrices defined during the elimination are P1 =   1 0 0 0 1 0 0 0 1   , M1 =   1 0 0 0.3 1 0 −0.5 0 1   , P2 =   1 0 0 0 0 1 0 1 0   , M2 =   1 0 0 0 1 0 0 0.04 1   , 2.6. Why Is Pivoting Necessary? 7 The corresponding L’s are L1 =   1 0 0 0.5 1 0 −0.3 0 1   , L2 =   1 0 0 0 1 0 0 −0.04 1   , The relation LU = PA is called the LU factorization or the triangular de- composition of A. It should be emphasized that nothing new has been introduced. Computationally, elimination is done by row operations on the coefficient matrix, not by actual matrix multiplication. LU factorization is simply Gaussian elimina- tion expressed in matrix notation. With this factorization, a general system of equations Ax = b becomes a pair of triangular systems Ly = Pb Ux = y 2.6 Why Is Pivoting Necessary? The diagonal elements of U are called pivots. The kth pivot is the coefficient of the kth variable in the kth equation at the kth step of the elimination. In our 3-by-3 example, the pivots are 10, 2.5, and 6.2. Both the computation of the multipliers and the back substitution require divisions by the pivots. Consequently, the algorithm cannot be carried out if any of the pivots are zero. Intuition should tell us that it is a bad idea to complete the computation if any of the pivots are nearly zero. To see this, let us change our example slightly to   10 −7 0 −3 2.099 6 5 −1 5     x1 x2 x3   =   7 3.901 6   The (2,2) element of the matrix has been changed from 2.000 to 2.099, and the right-hand side has also been changed so that the exact answer is still (0,−1, 1)T . Let us assume that the solution is to be computed on a hypothetical machine that does decimal floating-point arithmetic with five significant digits. The first step of the elimination produces   10 −7 0 0 −0.001 6 0 2.5 5     x1 x2 x3   =   7 6.001 2.5   The (2,2) element is now quite small compared with the other elements in the ma- trix. Nevertheless, let us complete the elimination without using any interchanges. The next step requires adding 2.5 · 103 times the second equation to the third. (5 + (2.5 · 103)(6))x3 = (2.5 + (2.5 · 103)(6.001)) 10 Chapter 2. Linear Equations Study this function carefully. Almost all the execution time is spent in the statement A(i,j) = A(i,j) - A(i,k)*A(k,j); At the kth step of the elimination, i and j are index vectors of length n-k. The operation A(i,k)*A(k,j) multiplies a column vector by a row vector to produce a square, rank one matrix of order n-k. This matrix is then subtracted from the submatrix of the same size in the bottom right corner of A. In a programming language without vector and matrix operations, this update of a portion of A would be done with doubly nested loops on i and j. The second function, bslashtx, is a simplified version of the built-in Matlab backslash operator. It begins by checking for three important special cases, lower triangular, upper triangular, and symmetric positive definite. Linear systems with these properties can be solved is less time than a general system. function x = bslashtx(A,b) % BSLASHTX Solve linear system (backslash) % x = bslashtx(A,b) solves A*x = b [n,n] = size(A); if isequal(triu(A,1),zeros(n,n)) % Lower triangular x = forward(A,b); return elseif isequal(tril(A,-1),zeros(n,n)) % Upper triangular x = backsubs(A,b); return elseif isequal(A,A’) [R,fail] = chol(A); if ~fail % Positive definite y = forward(R’,b); x = backsubs(R,y); return end end If none of the special cases is detected, bslashtx calls lutx to permute and fac- tor the coefficient matrix, then uses the permutation and factors to complete the solution of a linear system. % Triangular factorization [L,U,p] = lutx(A); % Permutation and forward elimination y = forward(L,b(p)); 2.8. Effect of Roundoff Errors 11 % Back substitution x = backsubs(U,y); The bslashtx function employs subfunctions to carry out the solution of lower and upper triangular systems. function x = forward(L,x) % FORWARD. Forward elimination. % For lower triangular L, x = forward(L,b) solves L*x = b. [n,n] = size(L); for k = 1:n j = 1:k-1; x(k) = (x(k) - L(k,j)*x(j))/L(k,k); end function x = backsubs(U,x) % BACKSUBS. Back substitution. % For upper triangular U, x = backsubs(U,b) solves U*x = b. [n,n] = size(U); for k = n:-1:1 j = k+1:n; x(k) = (x(k) - U(k,j)*x(j))/U(k,k); end A third function, lugui, shows the steps in LU decomposition by Gaussian elimination. It is a version of lutx that allows you to experiment with various pivot selection strategies. At the kth step of the elimination, the largest element in the unreduced portion of the kth column is shown in magenta. This is the element that partial pivoting would ordinarily select as the pivot. You can then choose among four different pivoting strategies: • Pick a pivot. Use the mouse to pick the magenta element, or any other element, as pivot. • Diagonal pivoting. Use the diagonal element as the pivot. • Partial pivoting. Same strategy as lu and lutx. • Complete pivoting. Use the largest element in the unfactored submatrix as the pivot. The chosen pivot is shown in red and the resulting elimination step is taken. As the process proceeds, the emerging columns of L are shown in green, and the emerging rows of U in blue. 2.8 Effect of Roundoff Errors The rounding errors introduced during the solution of a linear system of equations almost always cause the computed solution — which we now denote by x∗ — to 12 Chapter 2. Linear Equations differ somewhat from the theoretical solution, x = A−1b. In fact, if the elements of x are not floating-point numbers, then x∗ cannot equal x. There are two common measures of the discrepancy in x∗, the error, e = x− x∗ and the residual, r = b−Ax∗ Matrix theory tells us that, because A is nonsingular, if one of these is zero, the other must also be zero. But they are not necessarily both “small” at the same time. Consider the following example: ( 0.780 0.563 0.913 0.659 )( x1 x2 ) = ( 0.217 0.254 ) What happens if we carry out Gaussian elimination with partial pivoting on a hypothetical three-digit decimal computer? First, the two rows (equations) are interchanged so that 0.913 becomes the pivot. Then the multiplier 0.780 0.913 = 0.854 (to three places) is computed. Next, 0.854 times the new first row is subtracted from the new second row to produce the system ( 0.913 0.659 0 0.001 )( x1 x2 ) = ( 0.254 0.001 ) Finally, the back substitution is carried out: x2 = 0.001 0.001 = 1.00 (exactly), x1 = 0.254− 0.659x2 0.913 = −0.443 (to three places). Thus the computed solution is x∗ = (−0.443 1.000 ) To assess the accuracy without knowing the exact answer, we compute the residuals (exactly): r = b−Ax∗ = ( 0.217− ((0.780)(−0.443) + (0.563)(1.00)) 0.254− ((0.913)(−0.443) + (0.659)(1.00)) ) = (−0.000460 −0.000541 ) 2.9. Norms and Condition Numbers 15 1 ≤ p ≤ ∞. ‖x‖p = ( n∑ i=1 |xi|p)1/p We almost always use p = 1, p = 2 or lim p →∞; ‖x‖1 = n∑ i=1 |xi| ‖x‖2 = ( n∑ i=1 |xi|2)1/2 ‖x‖∞ = max i |xi| The l1 norm is also known as the Manhattan norm because it corresponds to the distance traveled on a grid of city streets. The l2 norm is the familiar Euclidean distance. The l∞ norm is also known as the Chebyshev norm. The particular value of p is often unimportant and we simply use ‖x‖. All vector norms have the following basic properties associated with the notion of dis- tance. ‖x‖ > 0 if x 6= 0 ‖0‖ = 0 ‖cx‖ = |c|‖x‖ for all scalars c ‖x + y‖ ≤ ‖x‖+ ‖y‖, (the triangle inequality) In Matlab, ‖x‖p is computed by norm(x,p) and norm(x) is the same as norm(x,2). For example: x = (1:4)/5 norm1 = norm(x,1) norm2 = norm(x) norminf = norm(x,inf) produces x = 0.2000 0.4000 0.6000 0.8000 norm1 = 2.0000 norm2 = 1.0954 norminf = 0.8000 16 Chapter 2. Linear Equations Multiplication of a vector x by a matrix A results in a new vector Ax that can have a very different norm from x. This change in norm is directly related to the sensitivity we want to measure. The range of the possible change can be expressed by two numbers, M = max ‖Ax‖ ‖x‖ m = min ‖Ax‖ ‖x‖ The max and min are taken over all nonzero vectors, x. Note that if A is singular, then m = 0. The ratio M/m is called the condition number of A, κ(A) = max ‖Ax‖‖x‖ min ‖Ax‖‖x‖ The actual numerical value of κ(A) depends on the vector norm being used, but we are usually only interested in order of magnitude estimates of the condition number, so the particular norm is usually not very important. Consider a system of equations Ax = b and a second system obtained by altering the right-hand side: A(x + δx) = b + δb We think of δb as being the error in b and δx as being the resulting error in x, although we need not make any assumptions that the errors are small. Because A(δx) = δb, the definitions of M and m immediately lead to ‖b‖ ≤ M‖x‖ and ‖δb‖ ≥ m‖δx‖ Consequently, if m 6= 0, ‖δx‖ ‖x‖ ≤ κ(A) ‖δb‖ ‖b‖ The quantity ‖δb‖/‖b‖ is the relative change in the right-hand side, and the quantity ‖δx‖/‖x‖ is the relative error caused by this change. The advantage of using relative changes is that they are dimensionless, that is, they are not affected by overall scale factors. This shows that the condition number is a relative error magnification factor. Changes in the right-hand side can cause changes κ(A) times as large in the solution. It turns out that the same is true of changes in the coefficient matrix itself. 2.9. Norms and Condition Numbers 17 The condition number is also a measure of nearness to singularity. Although we have not yet developed the mathematical tools necessary to make the idea pre- cise, the condition number can be thought of as the reciprocal of the relative distance from the matrix to the set of singular matrices. So, if κ(A) is large, A is close to singular. Some of the basic properties of the condition number are easily derived. Clearly, M ≥ m, and so κ(A) ≥ 1 If P is a permutation matrix, then the components of Px are simply a rearrangement of the components of x. It follows that ‖Px‖ = ‖x‖ for all x, and so κ(P ) = 1 In particular, κ(I) = 1. If A is multiplied by a scalar c, then M and m are both multiplied by the same scalar, and so κ(cA) = κ(A) If D is a diagonal matrix, then κ(D) = max |dii| min |dii| These last two properties are two of the reasons that κ(A) is a better measure of nearness to singularity than the determinant of A. As an extreme example, consider a 100-by-100 diagonal matrix with 0.1 on the diagonal. Then det(A) = 10−100, which is usually regarded as a small number. But κ(A) = 1, and the components of Ax are simply 0.1 times the corresponding components of x. For linear systems of equations, such a matrix behaves more like the identity than like a singular matrix. The following example uses the l1 norm. A = ( 4.1 2.8 9.7 6.6 ) b = ( 4.1 9.7 ) x = ( 1 0 ) Clearly, Ax = b, and ‖b‖ = 13.8, ‖x‖ = 1 If the right-hand side is changed to b̃ = ( 4.11 9.70 ) the solution becomes x̃ = ( 0.34 0.97 ) 20 Chapter 2. Linear Equations The residual involves the product Ax∗ so it is appropriate to consider the relative residual, which compares the norm of b− Ax to the norms of A and x∗. It follows directly from the above inequalities that ‖b−Ax∗‖ ‖A‖‖x∗‖ ≤ ρ² If A is nonsingular, the error can be expressed using the inverse of A by x− x∗ = A−1(b−Ax∗) and so ‖x− x∗‖ ≤ ‖A−1‖‖E‖‖x∗‖ It is simplest to compare the norm of the error with the norm of the computed solution. Thus the relative error satisfies ‖x− x∗‖ ‖x∗‖ ≤ ρ‖A‖‖A −1‖² Hence ‖x− x∗‖ ‖x∗‖ ≤ ρκ(A)² The actual computation of κ(A) requires knowing ‖A−1‖. But computing A−1 requires roughly three times as much work as solving a single linear system. Computing the l2 condition number requires the singular value decomposition and even more work. Fortunately, the exact value of κ(A) is rarely required. Any reasonably good estimate of it is satisfactory. Matlab has several functions for computing or estimating condition numbers. • cond(A) or cond(A,2) computes κ2(A). Uses svd(A). Suitable for smaller matrices where the geometric properties of the l2 norm are important. • cond(A,1) computes κ1(A). Uses inv(A). Less work than cond(A,2). • cond(A,inf) computes κ∞(A). Uses inv(A). Same as cond(A’,1). • condest(A) estimates κ1(A). Uses lu(A) and a recent algorithm of Higham and Tisseur [9]. Especially suitable for large, sparse matrices. • rcond(A) estimates 1/κ1(A). Uses lu(A) and an older algorithm developed by the LINPACK and LAPACK projects. Primarily of historical interest. 2.10 Sparse Matrices and Band Matrices Sparse matrices and band matrices occur frequently in technical computing. The sparsity of a matrix is the fraction of its elements that are zero. The Matlab function nnz counts the number of nonzeros in a matrix, so the sparsity of A is given by 2.10. Sparse Matrices and Band Matrices 21 density = nnz(A)/prod(size(A)) sparsity = 1 - density A sparse matrix is a matrix whose sparsity is nearly equal to 1. The bandwidth of a matrix is the maximum distance of the nonzero elements from the main diagonal. [i,j] = find(A) bandwidth = max(abs(i-j)) A band matrix is a matrix whose bandwidth is small. As you can see, both sparsity and bandwidth are matters of degree. An n-by-n diagonal matrix with no zeros on the diagonal has sparsity 1− 1/n and bandwidth 0, so it is an extreme example of both a sparse matrix and a band matrix. On the other hand, an n-by-n matrix with no zero elements, such as the one created by rand(n,n), has sparsity equal to zero, bandwidth equal to n− 1, and so is far from qualifying for either category. The Matlab sparse data structure stores the nonzero elements together with information about their indices. The sparse data structure also provides efficient handling of band matrices, so Matlab does not have a separate band matrix storage class. The statement S = sparse(A) converts a matrix to its sparse representation. The statement A = full(S) reverses the process. However, most sparse matrices have orders so large that it is impractical to store the full representation. More frequently, sparse matrices are created by S = sparse(i,j,x,m,n) This produces a matrix S with [i,j,x] = find(S) [m,n] = size(S) Most Matlab matrix operations and functions can be applied to both full and sparse matrices. The dominant factor in determining the execution time and mem- ory requirements for sparse matrix operations is the number of nonzeros, nnz(S), in the various matrices involved. A matrix with bandwidth equal to 1 is known as a tridiagonal matrix. It is worthwhile to have a specialized function for one particular band matrix operation, the solution of a tridiagonal system of simultaneous linear equations.   b1 c1 a1 b2 c2 a2 b3 c3 . . . . . . . . . an−2 bn−1 cn−1 an−1 bn     x1 x2 x3 ... xn−1 xn   =   d1 d2 d3 ... dn−1 dn   22 Chapter 2. Linear Equations The function tridisolve is included in the NCM directory. The statement x = tridisolve(a,b,c,d) solves the tridiagonal system with subdiagonal a, diagonal b, superdiagonal c, and right-hand side d. We have already seen the algorithm that tridisolve uses; it is Gaussian elimination. In many situations involving tridiagonal matrices, the diagonal elements dominate the off-diagonal elements, so pivoting is unnecessary. Furthermore, the right-hand side is processed at the same time as the matrix itself. In this context, Gaussian elimination without pivoting is also known as the Thomas algorithm. The body of tridisolve begins by copying the right-hand side to a vector that will become the solution. x = d; n = length(x); The forward elimination step is a simple for loop. for j = 1:n-1 mu = a(j)/b(j); b(j+1) = b(j+1) - mu*c(j); x(j+1) = x(j+1) - mu*x(j); end The mu’s would be the multipliers on the subdiagonal of L if we were saving the LU factorization. Instead, the right-hand side is processed in the same loop. The back substitution step is another simple loop. x(n) = x(n)/b(n); for j = n-1:-1:1 x(j) = (x(j)-c(j)*x(j+1))/b(j); end Because tridisolve does not use pivoting, the results might be inaccurate if abs(b) is much smaller than abs(a)+abs(c). More robust, but slower, alternatives that do use pivoting include generating a full matrix with diag, T = diag(a,-1) + diag(b,0) + diag(c,1); x = T\d or generating a sparse matrix with spdiags, S = spdiags([a b c],[-1 0 1],n,n); x = S\d 2.11 PageRank and Markov Chains One of the reasons why GoogleTM is such an effective search engine is the PageRankTM algorithm developed by Google’s founders, Larry Page and Sergey Brin, when they 2.11. PageRank and Markov Chains 25 to divide each column of G by the corresponding element of c. Until this is available, it is best to use the spdiags function to create a sparse diagonal matrix, D = spdiags(1./c’,0,n,n) The sparse matrix product G*D will then be computed efficiently. The statements p = .85 delta = (1-p)/n e = ones(n,1) I = speye(n,n) x = (I - p*G*D)\(delta*e) compute PageRank by solving the sparse linear system with Gaussian elimination. indexbackslash indexinverse iteration It is also possible to use an algorithm known as inverse iteration. A = p*G*D + delta x = (I - A)\e x = x/sum(x) At first glance, this appears to be a very dangerous idea. Because I − A is the- oretically singular, with exact computation some diagonal element of the upper triangular factor of I − A should be zero and this computation should fail. But with roundoff error, the computed matrix I - A is probably not exactly singular. Even if it is singular, roundoff during Gaussian elimination will most likely pre- vent any exact zero diagonal elements. We know that Gaussian elimination with partial pivoting always produces a solution with a small residual, relative to the computed solution, even if the matrix is badly conditioned. The vector obtained with the backslash operation, (I - A)\e, usually has very large components. If it is rescaled by its sum, the residual is scaled by the same factor and becomes very small. Consequently, the two vectors x and A*x almost equal each other to within roundoff error. In this setting, solving the singular system with Gaussian elimination blows up, but it blows up in exactly the right direction. Figure 2.1 is the graph for a tiny example, with n = 6 instead of n = 3 · 109. Pages on the Web are identified by strings known as uniform resource locators, or URLs. Most URLs begin with http because they use the hypertext transfer protocol. In Matlab we can store the URLs as an array of strings in a cell array. This example involves a 6-by-1 cell array. U = {’http://www.alpha.com’ ’http://www.beta.com’ ’http://www.gamma.com’ ’http://www.delta.com’ ’http://www.rho.com’ ’http://www.sigma.com’} Two different kinds of indexing into cell arrays are possible. Parentheses denote subarrays, including individual cells, and curly braces denote the contents of the 26 Chapter 2. Linear Equations alpha beta gamma delta sigma rho Figure 2.1. A tiny Web cells. If k is a scalar, then U(k) is a 1-by-1 cell array consisting of the kth cell in U, while U{k} is the string in that cell. Thus U(1) is a single cell and U{1} is the string ’http://www.alpha.com’. Think of mail boxes with addresses on a city street. B(502) is the box at number 502, while B{502} is the mail in that box. We can generate the connectivity matrix by specifying the pairs of indices (i,j) of the nonzero elements. Because there is a link to beta.com from alpha.com, the (2,1) element of G is nonzero. The nine connections are described by i = [ 2 3 4 4 5 6 1 6 1] j = [ 1 2 2 3 3 3 4 5 6] A sparse matrix is stored in a data structure that requires memory only for the nonzero elements and their indices. This is hardly necessary for a 6-by-6 matrix with only 27 zero entries, but it becomes crucially important for larger problems. The statements n = 6 G = sparse(i,j,1,n,n); full(G) generate the sparse representation of an n-by-n matrix with ones in the positions specified by the vectors i and j and display its full representation. 0 0 0 1 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 1 1 0 0 0 0 0 1 0 0 0 0 0 1 0 1 0 2.11. PageRank and Markov Chains 27 The statement c = full(sum(G)) computes the column sums c = 1 2 3 1 1 1 The statement x = (I - p*G*D)\(delta*e) solves the sparse linear system to produce x = 0.2675 0.2524 0.1323 0.1697 0.0625 0.1156 For this tiny example, the smallest element of the Markov transition matrix is δ = .15/6 = .0250. A = p*G*D + delta A = 0.0250 0.0250 0.0250 0.8750 0.0250 0.8750 0.8750 0.0250 0.0250 0.0250 0.0250 0.0250 0.0250 0.4500 0.0250 0.0250 0.0250 0.0250 0.0250 0.4500 0.3083 0.0250 0.0250 0.0250 0.0250 0.0250 0.3083 0.0250 0.0250 0.0250 0.0250 0.0250 0.3083 0.0250 0.8750 0.0250 Notice that the column sums of A are all equal to one. Computing PageRank with inverse iteration, x = (I - A)\e produces a warning message about ill conditioning and a vector with elements on the order of 1016. On one particular computer the elements of x happen to be negative and their sum is s = sum(x) = -6.6797e+016 Other computers with different roundoff error might give other results. But in all cases, the rescaled solution x = x/sum(x) 30 Chapter 2. Linear Equations 0 50 100 150 200 250 300 350 400 450 500 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 Page Rank Figure 2.4. Page Rank of the harvard500 graph page-rank in out url 1 0.0823 195 26 http://www.harvard.edu 10 0.0161 21 18 http://www.hbs.edu 42 0.0161 42 0 http://search.harvard.edu:8765/custom/query.html 130 0.0160 24 12 http://www.med.harvard.edu 18 0.0135 45 46 http://www.gse.harvard.edu 15 0.0129 16 49 http://www.hms.harvard.edu 9 0.0112 21 27 http://www.ksg.harvard.edu 17 0.0109 13 6 http://www.hsph.harvard.edu 46 0.0097 18 21 http://www.gocrimson.com 13 0.0084 9 1 http://www.hsdm.med.harvard.edu 260 0.0083 26 1 http://search.harvard.edu:8765/query.html 19 0.0081 23 21 http://www.radcliffe.edu The URL where the search began, www.harvard.edu, dominates. Like most uni- versities, Harvard is organized into various colleges and institutes, including the Kennedy School of Government, the Harvard Medical School, the Harvard Busi- ness School, and the Radcliffe Institute. You can see that the home pages of these schools have high PageRank. With a different sample, such as the one generated by Google itself, the ranks would be different. 2.12 Further Reading Further reading on matrix computation includes books by Demmel [2], Golub and Van Loan [3], Stewart [4, 5], and Trefethen and Bau [6]. The definitive references on Fortran matrix computation software are the LAPACK Users’ Guide and Web site [1]. The Matlab sparse matrix data structure and operations are described in [8]. Exercises 31 Information available on Web sites about PageRank includes a brief explanation at Google [7], a technical report by Page, Brin, and colleagues [10], and a paper by John Tomlin and colleagues [11]. Exercises 2.1. Alice buys three apples, a dozen bananas and one cantaloupe for $2.36. Bob buys a dozen apples and two cantaloupes for $5.26. Carol buys two bananas and three cantaloupes for $2.77. How much do single pieces of each fruit cost? (You might want to set format bank.) 2.2. What Matlab function computes the reduced row echelon form of a ma- trix? What Matlab function generates magic square matrices? What is the reduced row echelon form of the magic square of order six? 2.3. Figure 2.5 depicts a plane truss having 13 members (the numbered lines) connecting 8 joints (the numbered circles). The indicated loads, in tons, are applied at joints 2, 5, and 6, and we want to determine the resulting force on each member of the truss. 1 2 5 6 8 3 4 7 1 3 5 7 9 11 12 2 6 10 13 4 8 10 15 20 Figure 2.5. A plane truss For the truss to be in static equilibrium, there must be no net force, hor- izontally or vertically, at any joint. Thus, we can determine the member forces by equating the horizontal forces to the left and right at each joint, and similarly equating the vertical forces upward and downward at each joint. For the eight joints, this would give 16 equations, which is more than the 13 unknown factors to be determined. For the truss to be statically determi- nate, that is, for there to be a unique solution, we assume that joint 1 is rigidly fixed both horizontally and vertically, and that joint 8 is fixed verti- cally. Resolving the member forces into horizontal and vertical components and defining α = 1/ √ 2, we obtain the following system of equations for the 32 Chapter 2. Linear Equations member forces fi: Joint 2: f2 = f6 f3 = 10 Joint 3: αf1 = f4 + αf5 αf1 + f3 + αf5 = 0 Joint 4: f4 = f8 f7 = 0 Joint 5: αf5 + f6 = αf9 + f10 αf5 + f7 + αf9 = 15 Joint 6: f10 = f13 f11 = 20 Joint 7: f8 + αf9 = αf12 αf9 + f11 + αf12 = 0 Joint 8: f13 + αf12 = 0 Solve this system of equations to find the vector f of member forces. 2.4. Figure 2.6 is the circuit diagram for a small network of resistors. 1 2 3 4 5 r 23 r 34 r 45 r 25 r 13 r 12 r 14 r 35 v s i 1 i 2 i 3 i 4 Figure 2.6. A resistor network There are five nodes, eight resistors, and one constant voltage source. We want to compute the voltage drops between the nodes and the currents around each of the loops. Several different linear systems of equations can be formed to describe this circuit. Let vk, k = 1, . . . , 4 denote the voltage difference between each of the first four nodes and node number 5 and let ik, k = 1, . . . , 4 denote the clockwise current around each of the loops in the diagram. Ohm’s law says that the voltage drop across a resistor is the resistance times the current. For example, the branch between nodes 1 and 2 gives v1 − v2 = r12(i2 − i1) Exercises 35 M = magic(n) H = hilb(n) P = pascal(n) I = eye(n,n) R = randn(n,n) R = randn(n,n); A = R’ * R R = randn(n,n); A = R’ + R R = randn(n,n); I = eye(n,n); A = R’ + R + n*I (b) If the matrix R is upper triangular, then equating individual elements in the equation A = RT R gives akj = k∑ i=1 rikrij , k ≤ j Using these equations in different orders yields different variants of the Cholesky algorithm for computing the elements of R. What is one such algorithm? 2.6. This example shows that a badly conditioned matrix does not necessarily lead to small pivots in Gaussian elimination. The matrix is the n-by-n upper triangular matrix A with elements aij =    −1, i < j 1, i = j 0, i > j Show how to generate this matrix in Matlab with eye, ones, and triu. Show that κ1(A) = n2n−1 For what n does κ1(A) exceed 1/eps? This matrix is not singular, so Ax cannot be zero unless x is zero. However, there are vectors x for which ‖Ax‖ is much smaller than ‖x‖. Find one such x. Because this matrix is already upper triangular, Gaussian elimination with partial pivoting has no work to do. What are the pivots? Use lugui to design a pivot strategy that will produce smaller pivots than partial pivoting. (Even these pivots do not completely reveal the large con- dition number.) 2.7. The matrix factorization LU = PA can be used to compute the determinant of A. We have det(L)det(U) = det(P )det(A) Because L is triangular with ones on the diagonal, det(L) = 1. Because U is triangular, det(U) = u11u22 · · ·unn. Because P is a permutation, det(P ) = +1 if the number of interchanges is even and −1 if it is odd. So det(A) = ±u11u22 · · ·unn 36 Chapter 2. Linear Equations Modify the lutx function so that it returns four outputs: function [L,U,p,sig] = lutx(A) %LU Triangular factorization % [L,U,p,sig] = lutx(A) computes a unit lower triangular % matrix L, an upper triangular matrix U, a permutation % vector p and a scalar sig, so that L*U = A(p,:) and % sig = +1 or -1 if p is an even or odd permutation. Write a function mydet(A) that uses your modified lutx to compute the determinant of A. In Matlab, the product u11u22 · · ·unn can be computed with prod(diag(U)). 2.8. Modify the lutx function so that it uses explicit for loops instead of Matlab vector notation. For example, one section of your modified program will read % Compute the multipliers for i = k+1:n A(i,k) = A(i,k)/A(k,k); end Compare the execution time of your modified lutx program with the original lutx program and with the built-in lu function by finding the order of the matrix for which each of the three programs takes about 10 seconds on your computer. 2.9. Let A =   1 2 3 4 5 6 7 8 9   , b =   1 3 5   (a) Show that the set of linear equations Ax = b has infinitely many solutions. Describe the set of possible solutions. (b) Suppose Gaussian elimination is used to solve Ax = b using exact arith- metic. Because there are infinitely many solutions, it is unreasonable to expect one particular solution to be computed. What does happen? (c) Use bslashtx to solve Ax = b on an actual computer with floating-point arithmetic. What solution is obtained? Why? In what sense is it a “good” solution? In what sense is it a “bad” solution? indexbackslash (d) Explain why the built-in backslash operator, x = A\b, gives a different solution from x = bslashtx(A,b). 2.10. Section 2.4 gives two algorithms for solving triangular systems. One subtracts columns of the triangular matrix from the right hand side; the other uses inner products between the rows of the triangular matrix and the emerging solution. (a) Which of these two algorithms does bslashtx use? (b) Write another function, bslashtx2, that uses the other algorithm. Exercises 37 2.11. The inverse of a matrix A can be defined as the matrix X whose columns xj solve the equations Axj = ej where ej is the jth column of the identity matrix. (a) Starting with the function bslashtx, write a Matlab function X = myinv(A) that computes the inverse of A. Your function should call lutx only once and should not use the built-in Matlab backslash operator or inv function. (b) Test your function by comparing the inverses it computes with the inverses obtained from the built-in inv(A) on a few test matrices. 2.12. If built-in Matlab lu function is called with only two output arguments [L,U] = lu(A) the permutations are incorporated into the output matrix L. The help entry for lu describes L as “psychologically lower triangular.” Modify lutx so that it does the same thing. You can use if nargout == 2, ... to test the number of output arguments. 2.13. (a) Why is M = magic(8) lugui(M) an interesting example? (b) How is the behavior of lugui(M) related to rank(M)? (c) Can you pick a sequence of pivots so that no roundoff error occurs in lugui(M)? 2.14. The pivot selection strategy known as complete pivoting is one of the options available in lugui. It has some slight numerical advantages over partial pivoting. At each state of the elimination the element of largest magnitude in the entire unreduced matrix is selected as pivot. This involves both row and column interchanges and produces two permutation vectors p and q so that L*U = A(p,q) Modify lutx and bslashtx so that they use complete pivoting. 2.15. The function golub in the NCM directory is named after Stanford professor Gene Golub. The function generates test matrices with random integer en- tries. The matrices are very badly conditioned, but Gaussian elimination without pivoting fails to produce the small pivots that would reveal the large condition number.