![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to doc.n CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10dc / cmd / matlab|
Return to doc.n CVS log | Up to [Research Unix] / researchv10dc / cmd / matlab |
1.1 ! root 1: .ds xl MATLAB ! 2: .hy 0 ! 3: .de pp ! 4: .sp ! 5: .ti +5 ! 6: .. ! 7: .de hd ! 8: 'sp 2 ! 9: .if \\n%-1 'tl '\*(xl, page %' ! 10: .if !\\n%-1 'tl '\n(mo/\n(dy/\n(yr' ! 11: 'sp 3 ! 12: .if !\\n%-1 'ns ! 13: .. ! 14: .de fo ! 15: 'bp ! 16: .. ! 17: .wh 0 hd ! 18: .wh -7 fo ! 19: .tr ` ! 20: .de xx ! 21: .sp ! 22: .ti -6 ! 23: .. ! 24: \ ! 25: .sp 7 ! 26: .ce 2 ! 27: MATLAB Users' Guide ! 28: February, 1982 ! 29: .sp 2 ! 30: .ce 3 ! 31: Cleve Moler ! 32: Department of Computer Science ! 33: University of New Mexico ! 34: .sp 2 ! 35: .in +5 ! 36: .ll -5 ! 37: ABSTRACT. ! 38: MATLAB is an interactive computer program that serves as a ! 39: convenient "laboratory" for computations involving matrices. ! 40: It provides easy access to matrix software developed by the ! 41: LINPACK and EISPACK projects. ! 42: The program is written in Fortran and is designed to be readily ! 43: installed under any operating system which permits interactive ! 44: execution of Fortran programs. ! 45: .sp 3 ! 46: .ce ! 47: CONTENTS ! 48: ! 49: 1. Elementary operations page 2 ! 50: 2. MATLAB functions 8 ! 51: 3. Rows, columns and submatrices 9 ! 52: 4. FOR, WHILE and IF 10 ! 53: 5. Commands, text, files and macros 12 ! 54: 6. Census example 13 ! 55: 7. Partial differential equation 19 ! 56: 8. Eigenvalue sensitivity example 23 ! 57: 9. Syntax diagrams 27 ! 58: 10. The parser-interpreter 31 ! 59: 11. The numerical algorithms 34 ! 60: 12. FLOP and CHOP 37 ! 61: 13. Communicating with other programs 41 ! 62: Appendix. The HELP file 46 ! 63: .in -5 ! 64: .ll +5 ! 65: .bp 1 ! 66: \ ! 67: .sp 7 ! 68: .ce 2 ! 69: MATLAB Users' Guide ! 70: February, 1982 ! 71: .sp 2 ! 72: .ce 3 ! 73: Cleve Moler ! 74: Department of Computer Science ! 75: University of New Mexico ! 76: .sp 2 ! 77: .pp ! 78: MATLAB is an interactive computer program that serves as a ! 79: convenient "laboratory" for computations involving matrices. ! 80: It provides easy access to matrix software developed by the ! 81: LINPACK and EISPACK projects [1-3]. The capabilities range from ! 82: standard tasks such as solving simultaneous linear equations and ! 83: inverting matrices, through symmetric and nonsymmetric eigenvalue ! 84: problems, to fairly sophisticated matrix tools such as the singular ! 85: value decomposition. ! 86: .pp ! 87: It is expected that one of MATLAB's primary uses will be in ! 88: the classroom. It should be useful in introductory courses in ! 89: applied linear algebra, as well as more advanced courses in numerical ! 90: analysis, matrix theory, statistics and applications of matrices ! 91: to other disciplines. In nonacademic settings, MATLAB can serve as a ! 92: "desk calculator" for the quick solution of small problems involving ! 93: matrices. ! 94: .pp ! 95: The program is written in Fortran and is designed to be readily ! 96: installed under any operating system which permits interactive ! 97: execution of Fortran programs. The resources required are fairly ! 98: modest. There are less than 7000 lines of Fortran source code, ! 99: including the LINPACK and EISPACK subroutines used. With proper ! 100: use of overlays, it is possible run the system on a minicomputer ! 101: with only 32K bytes of memory. ! 102: .pp ! 103: The size of the matrices that can be handled in MATLAB ! 104: depends upon the amount of storage that is set aside when the ! 105: system is compiled on a particular machine. We have found that ! 106: an allocation of 5000 words for matrix elements is usually quite ! 107: satisfactory. This provides room for several 20 by 20 matrices, ! 108: for example. One implementation on a virtual memory system ! 109: provides 100,000 elements. Since most of the algorithms used ! 110: access memory in a sequential fashion, the large amount of ! 111: allocated storage causes no difficulties. ! 112: .pp ! 113: In some ways, MATLAB resembles SPEAKEASY [4] and, to a ! 114: lesser extent, APL. All are interactive terminal languages ! 115: that ordinarily accept single-line commands or statements, process ! 116: them immediately, and print the results. All have arrays or matrices ! 117: as principal data types. But for MATLAB, the matrix is the only ! 118: data type (although scalars, vectors and text are special cases), ! 119: the underlying system is portable and requires fewer resources, ! 120: and the supporting subroutines are more powerful and, in some ! 121: cases, have better numerical properties. ! 122: .pp ! 123: Together, LINPACK and EISPACK represent the state of the art ! 124: in software for matrix computation. EISPACK is a package of over ! 125: 70 Fortran subroutines for various matrix eigenvalue computations ! 126: that are based for the most part on Algol procedures published by ! 127: Wilkinson, Reinsch and their colleagues [5]. LINPACK is a package ! 128: of 40 Fortran subroutines (in each of four data types) for solving ! 129: and analyzing simultaneous linear equations and related matrix ! 130: problems. Since MATLAB is not primarily concerned with either ! 131: execution time efficiency or storage savings, it ignores most of ! 132: the special matrix properties that LINPACK and EISPACK subroutines ! 133: use to advantage. Consequently, only 8 subroutines from LINPACK ! 134: and 5 from EISPACK are actually involved. ! 135: .pp ! 136: In more advanced applications, MATLAB can be used in conjunction ! 137: with other programs in several ways. ! 138: It is possible to define new MATLAB functions and add them to the system. ! 139: With most operating systems, it is possible to use the local ! 140: file system to pass matrices between MATLAB and other programs. ! 141: MATLAB command and statement input can be obtained from ! 142: a local file instead of from the terminal. ! 143: The most power and flexibility is obtained by using MATLAB ! 144: as a subroutine which is called by other programs. ! 145: .pp ! 146: This document first gives an overview of MATLAB from the user's ! 147: point of view. ! 148: Several extended examples involving data fitting, partial differential ! 149: equations, eigenvalue sensitivity and other topics are included. ! 150: A formal definition of the MATLAB language and ! 151: an brief description of the parser and interpreter are given. ! 152: The system was designed and programmed using techniques described ! 153: by Wirth [6], implemented in nonrecursive, portable Fortran. ! 154: There is a brief discussion of some of the matrix algorithms ! 155: and of their numerical properties. ! 156: The final section describes how MATLAB can be used with other programs. ! 157: The appendix includes the HELP documentation available on-line. ! 158: ! 159: ! 160: 1. Elementary operations ! 161: ! 162: .pp ! 163: MATLAB works with essentially only one kind of object, ! 164: a rectangular matrix with complex elements. If the imaginary parts ! 165: of the elements are all zero, they are not printed, but they still ! 166: occupy storage. In some situations, special meaning is attached ! 167: to 1 by 1 matrices, that is scalars, and to 1 by n and m by 1 ! 168: matrices, that is row and column vectors. ! 169: .pp ! 170: Matrices can be introduced into MATLAB in four different ways: ! 171: -- Explicit list of elements, ! 172: -- Use of FOR and WHILE statements, ! 173: -- Read from an external file, ! 174: -- Execute an external Fortran program. ! 175: .pp ! 176: The explicit list is surrounded by angle brackets, '<' and '>', ! 177: and uses the semicolon ';' to indicate the ends of the rows. ! 178: For example, the input line ! 179: ! 180: A = <1 2 3; 4 5 6; 7 8 9> ! 181: ! 182: will result in the output ! 183: ! 184: A = ! 185: ! 186: 1. 2. 3. ! 187: 4. 5. 6. ! 188: 7. 8. 9. ! 189: ! 190: The matrix A will be saved for later use. The individual elements ! 191: are separated by commas or blanks and can be any MATLAB expressions, ! 192: for example ! 193: ! 194: x = < -1.3, 4/5, 4*atan(1) > ! 195: ! 196: results in ! 197: ! 198: X = ! 199: ! 200: -1.3000 0.8000 3.1416 ! 201: ! 202: The elementary functions available include sqrt, log, exp, sin, cos, ! 203: atan, abs, round, real, imag, and conjg. ! 204: .pp ! 205: Large matrices can be spread across several input lines, with the ! 206: carriage returns replacing the semicolons. The above matrix ! 207: could also have been produced by ! 208: ! 209: A = < 1 2 3 ! 210: 4 5 6 ! 211: 7 8 9 > ! 212: ! 213: .pp ! 214: Matrices can be input from the local file system. ! 215: Say a file named 'xyz' contains five lines of text, ! 216: ! 217: .ne 5 ! 218: A = < ! 219: 1 2 3 ! 220: 4 5 6 ! 221: 7 8 9 ! 222: >; ! 223: ! 224: then the MATLAB statement EXEC('xyz') reads the matrix and ! 225: assigns it to A . ! 226: .pp ! 227: The FOR statement allows the generation of matrices whose ! 228: elements are given by simple formulas. Our example matrix A ! 229: could also have been produced by ! 230: ! 231: for i = 1:3, for j = 1:3, a(i,j) = 3*(i-1)+j; ! 232: ! 233: The semicolon at the end of the line suppresses the printing, ! 234: which in this case would have been nine versions of A with changing ! 235: elements. ! 236: .pp ! 237: Several statements may be given on a line, separated ! 238: by semicolons or commas. ! 239: .pp ! 240: Two consecutive periods anywhere on a line indicate ! 241: continuation. The periods and any following characters are deleted, ! 242: then another line is input and concatenated onto the previous line. ! 243: .pp ! 244: Two consecutive slashes anywhere on a line cause the remainder ! 245: of the line to be ignored. This is useful for inserting comments. ! 246: .pp ! 247: Names of variables are formed by a letter, followed by any ! 248: number of letters and digits, but only the first 4 characters ! 249: are remembered. ! 250: .pp ! 251: The special character prime (') is used to denote the transpose ! 252: of a matrix, so ! 253: ! 254: x = x' ! 255: ! 256: changes the row vector above into the column vector ! 257: ! 258: X = ! 259: ! 260: -1.3000 ! 261: 0.8000 ! 262: 3.1416 ! 263: ! 264: .pp ! 265: Individual matrix elements may be referenced by enclosing their ! 266: subscripts in parentheses. When any element is changed, the entire ! 267: matrix is reprinted. For example, using the above matrix, ! 268: ! 269: a(3,3) = a(1,3) + a(3,1) ! 270: ! 271: results in ! 272: ! 273: A = ! 274: ! 275: 1. 2. 3. ! 276: 4. 5. 6. ! 277: 7. 8. 10. ! 278: ! 279: .pp ! 280: Addition, subtraction and multiplication of matrices are ! 281: denoted by +, -, and * . The operations are performed ! 282: whenever the matrices have the proper dimensions. For example, ! 283: with the above A and x, the expressions A + x and x*A are incorrect ! 284: because A is 3 by 3 and x is now 3 by 1. However, ! 285: ! 286: b = A*x ! 287: ! 288: is correct and results in the output ! 289: ! 290: B = ! 291: ! 292: 9.7248 ! 293: 17.6496 ! 294: 28.7159 ! 295: ! 296: Note that both upper and lower case letters are allowed for input ! 297: (on those systems which have both), but that lower case is converted ! 298: to upper case. ! 299: .pp ! 300: There are two "matrix division" symbols in MATLAB, \\ and / . ! 301: (If your terminal does not have a backslash, use $ instead, or see CHAR.) ! 302: If A and B are matrices, then A\\B and B/A correspond formally ! 303: to left and right multiplication of B by the inverse of A , that is ! 304: inv(A)*B and B*inv(A), but the result is obtained directly ! 305: without the computation of the inverse. In the scalar case, ! 306: 3\\1 and 1/3 have the same value, namely one-third. In general, ! 307: A\\B denotes the solution X to the equation A*X = B and B/A denotes ! 308: the solution to X*A = B. ! 309: .pp ! 310: Left division, A\\B, is defined whenever B has as many rows as A . ! 311: If A is square, it is factored using Gaussian elimination. The ! 312: factors are used to solve the equations A*X(:,j) = B(:,j) ! 313: where B(:,j) denotes the j-th column of B. The result is a ! 314: matrix X with the same dimensions as B. If A is nearly singular ! 315: (according to the LINPACK condition estimator, RCOND), a ! 316: warning message is printed. If A is not square, it is factored ! 317: using Householder orthogonalization with column pivoting. The factors ! 318: are used to solve the under- or overdetermined equations in a least squares ! 319: sense. The result is an m by n matrix X where m is the ! 320: number of columns of A and n is the number of columns of B . ! 321: Each column of X has at most k nonzero components, where k is the ! 322: effective rank of A . ! 323: .pp ! 324: Right division, B/A, can be defined in terms of left division ! 325: by B/A = (A'\\B')'. ! 326: .pp ! 327: For example, since our vector b was computed as A*x, ! 328: the statement ! 329: ! 330: y = A\\b ! 331: ! 332: results in ! 333: ! 334: Y = ! 335: ! 336: -1.3000 ! 337: 0.8000 ! 338: 3.1416 ! 339: ! 340: Of course, y is not exactly equal to x because of the ! 341: roundoff errors involved in both A*x and A\\b , but we are ! 342: not printing enough digits to see the difference. ! 343: The result of the statement ! 344: ! 345: e = x - y ! 346: ! 347: depends upon the particular computer being used. In one case ! 348: it produces ! 349: ! 350: E = ! 351: ! 352: 1.0e-15 * ! 353: ! 354: .3053 ! 355: -.2498 ! 356: .0000 ! 357: ! 358: The quantity 1.0e-15 is a scale factor which multiplies all the ! 359: components which follow. Thus our vectors x and y actually agree to ! 360: about 15 decimal places on this computer. ! 361: .pp ! 362: It is also possible to obtain element-by-element multiplicative ! 363: operations. If A and B have the same dimensions, then ! 364: A .* B denotes the matrix whose elements are simply the products ! 365: of the individual elements of A and B . ! 366: The expressions A ./ B and A .\\ B give the quotients of the ! 367: individual elements. ! 368: .pp ! 369: There are several possible output formats. The statement ! 370: ! 371: long, x ! 372: ! 373: results in ! 374: ! 375: X = ! 376: ! 377: -1.300000000000000 ! 378: .800000000000000 ! 379: 3.141592653589793 ! 380: ! 381: The statement ! 382: ! 383: short ! 384: ! 385: restores the original format. ! 386: .pp ! 387: The expression A**p means A to the p-th power. It is defined ! 388: if A is a square matrix and p is a scalar. If p is an integer ! 389: greater than one, the power is computed by repeated multiplication. ! 390: For other values of p the calculation involves the eigenvalues ! 391: and eigenvectors of A. ! 392: .pp ! 393: Previously defined matrices and matrix expressions can be used ! 394: inside brackets to generate larger matrices, for example ! 395: ! 396: C = <A, b; <4 2 0>*x, x'> ! 397: ! 398: results in ! 399: ! 400: ! 401: .ne 6 ! 402: C = ! 403: ! 404: 1.0000 2.0000 3.0000 9.7248 ! 405: 4.0000 5.0000 6.0000 17.6496 ! 406: 7.0000 8.0000 10.0000 28.7159 ! 407: -3.6000 -1.3000 0.8000 3.1416 ! 408: ! 409: .pp ! 410: There are four predefined variables, EPS, FLOP, RAND and EYE. ! 411: The variable EPS is used as a tolerance is determining such things ! 412: as near singularity and rank. Its initial value is the distance ! 413: from 1.0 to the next largest floating point number on the ! 414: particular computer being used. The user may reset this to ! 415: any other value, including zero. ! 416: EPS is changed by CHOP, which is described in section 12. ! 417: .pp ! 418: The value of RAND is a random variable, with a choice of a uniform or a ! 419: normal distribution. ! 420: .pp ! 421: The name EYE is used in place of I to denote ! 422: identity matrices because I is often used as a subscript or ! 423: as sqrt(-1). The dimensions of EYE are determined by context. ! 424: For example, ! 425: ! 426: B = A + 3*EYE ! 427: ! 428: adds 3 to the diagonal elements of A and ! 429: ! 430: X = EYE/A ! 431: ! 432: is one of several ways in MATLAB to invert a matrix. ! 433: .pp ! 434: FLOP provides a count of the number of floating point ! 435: operations, or "flops", required for each calculation. ! 436: .pp ! 437: A statement may consist of an expression alone, in which case ! 438: a variable named ANS is created and the result stored in ANS for ! 439: possible future use. Thus ! 440: ! 441: A\\A - EYE ! 442: ! 443: is the same as ! 444: ! 445: ANS = A\\A - EYE ! 446: ! 447: (Roundoff error usually causes this result to be a matrix of ! 448: "small" numbers, rather than all zeros.) ! 449: .pp ! 450: All computations are done using either single or double ! 451: precision real arithmetic, whichever is appropriate for the ! 452: particular computer. There is no mixed-precision arithmetic. ! 453: The Fortran COMPLEX data type is not used because many systems ! 454: create unnecessary underflows and overflows with complex operations ! 455: and because some systems do not allow double precision complex ! 456: arithmetic. ! 457: ! 458: ! 459: 2. MATLAB functions ! 460: .pp ! 461: Much of MATLAB's computational power comes from the various ! 462: matrix functions available. The current list includes: ! 463: ! 464: INV(A) - Inverse. ! 465: DET(A) - Determinant. ! 466: COND(A) - Condition number. ! 467: RCOND(A) - A measure of nearness to singularity. ! 468: EIG(A) - Eigenvalues and eigenvectors. ! 469: SCHUR(A) - Schur triangular form. ! 470: HESS(A) - Hessenberg or tridiagonal form. ! 471: POLY(A) - Characteristic polynomial. ! 472: SVD(A) - Singular value decomposition. ! 473: PINV(A,eps) - Pseudoinverse with optional tolerance. ! 474: RANK(A,eps) - Matrix rank with optional tolerance. ! 475: LU(A) - Factors from Gaussian elimination. ! 476: CHOL(A) - Factor from Cholesky factorization. ! 477: QR(A) - Factors from Householder orthogonalization. ! 478: RREF(A) - Reduced row echelon form. ! 479: ORTH(A) - Orthogonal vectors spanning range of A. ! 480: EXP(A) - e to the A. ! 481: LOG(A) - Natural logarithm. ! 482: SQRT(A) - Square root. ! 483: SIN(A) - Trigonometric sine. ! 484: COS(A) - Cosine. ! 485: ATAN(A) - Arctangent. ! 486: ROUND(A) - Round the elements to nearest integers. ! 487: ABS(A) - Absolute value of the elements. ! 488: REAL(A) - Real parts of the elements. ! 489: IMAG(A) - Imaginary parts of the elements. ! 490: CONJG(A) - Complex conjugate. ! 491: SUM(A) - Sum of the elements. ! 492: PROD(A) - Product of the elements. ! 493: DIAG(A) - Extract or create diagonal matrices. ! 494: TRIL(A) - Lower triangular part of A. ! 495: TRIU(A) - Upper triangular part of A. ! 496: NORM(A,p) - Norm with p = 1, 2 or 'Infinity'. ! 497: EYE(m,n) - Portion of identity matrix. ! 498: RAND(m,n) - Matrix with random elements. ! 499: ONES(m,n) - Matrix of all ones. ! 500: MAGIC(n) - Interesting test matrices. ! 501: HILBERT(n) - Inverse Hilbert matrices. ! 502: ROOTS(C) - Roots of polynomial with coefficients C. ! 503: DISPLAY(A,p) - Print base p representation of A. ! 504: KRON(A,B) - Kronecker tensor product of A and B. ! 505: PLOT(X,Y) - Plot Y as a function of X . ! 506: RAT(A) - Find "simple" rational approximation to A. ! 507: USER(A) - Function defined by external program. ! 508: .pp ! 509: Some of these functions have different interpretations when ! 510: the argument is a matrix or a vector and some of them have additional ! 511: optional arguments. Details are given in the HELP document in the ! 512: appendix. ! 513: .pp ! 514: Several of these functions can be used in a generalized ! 515: assignment statement with two or three variables on the left hand ! 516: side. For example ! 517: ! 518: <X,D> = EIG(A) ! 519: ! 520: stores the eigenvectors of A in the matrix X and a diagonal matrix ! 521: containing the eigenvalues in the matrix D. The statement ! 522: ! 523: EIG(A) ! 524: ! 525: simply computes the eigenvalues and stores them in ANS. ! 526: .pp ! 527: Future versions of MATLAB will probably include additional ! 528: functions, since they can easily be added to the system. ! 529: ! 530: ! 531: ! 532: 3. Rows, columns and submatrices ! 533: ! 534: .pp ! 535: Individual elements of a matrix can be accessed by giving ! 536: their subscripts in parentheses, eg. A(1,2), x(i), TAB(ind(k)+1). ! 537: An expression used as a subscript is rounded to the nearest ! 538: integer. ! 539: .pp ! 540: Individual rows and columns can be accessed using a colon ':' ! 541: (or a '|') for the free subscript. ! 542: For example, A(1,:) is the first row of A and A(:,j) is the j-th column. ! 543: Thus ! 544: ! 545: A(i,:) = A(i,:) + c*A(k,:) ! 546: ! 547: adds c times the k-th row of A to the i-th row. ! 548: .pp ! 549: The colon is used in several other ways in MATLAB, but ! 550: all of the uses are based on the following definition. ! 551: ! 552: j:k is the same as <j, j+1, ..., k> ! 553: j:k is empty if j > k . ! 554: j:i:k is the same as <j, j+i, j+2i, ..., k> ! 555: j:i:k is empty if i > 0 and j > k or if i < 0 and j < k . ! 556: ! 557: The colon is usually used with integers, but it is possible to ! 558: use arbitrary real scalars as well. Thus ! 559: ! 560: 1:4 is the same as <1, 2, 3, 4> ! 561: 0: 0.1: 0.5 is the same as <0.0, 0.1, 0.2, 0.3, 0.4, 0.5> ! 562: ! 563: .pp ! 564: In general, a subscript can be a vector. If X and V are ! 565: vectors, then X(V) is <X(V(1)), X(V(2)), ..., X(V(n))> . This ! 566: can also be used with matrices. If V has m components and ! 567: W has n components, then A(V,W) is the m by n matrix formed from ! 568: the elements of A whose subscripts are the elements of V and W. ! 569: Combinations of the colon notation and the indirect subscripting ! 570: allow manipulation of various submatrices. For example, ! 571: ! 572: A(<1,5>,:) = A(<5,1>,:) interchanges rows 1 and 5 of A. ! 573: A(2:k,1:n) is the submatrix formed from rows 2 through k ! 574: and columns 1 through n of A . ! 575: A(:,<3 1 2>) is a permutation of the first three columns. ! 576: ! 577: .pp ! 578: The notation A(:) has a special meaning. On the right hand side ! 579: of an assignment statement, it denotes all the elements of A, ! 580: regarded as a single column. ! 581: When an expression is assigned to A(:), the current dimensions of ! 582: A, rather than of the expression, are used. ! 583: ! 584: ! 585: 4. FOR, WHILE and IF ! 586: ! 587: .pp ! 588: The FOR clause allows statements to be repeated a specific ! 589: number of times. The general form is ! 590: ! 591: FOR variable = expr, statement, ..., statement, END ! 592: ! 593: The END and the comma before it may be omitted. In general, ! 594: the expression may be a matrix, in which case the columns are ! 595: stored one at a time in the variable and the following statements, ! 596: up to the END or the end of the line, are executed. The ! 597: expression is often of the form j:k, and its "columns" are ! 598: simply the scalars from j to k. Some examples (assume n has ! 599: already been assigned a value): ! 600: ! 601: for i = 1:n, for j = 1:n, A(i,j) = 1/(i+j-1); ! 602: ! 603: generates the Hilbert matrix. ! 604: ! 605: for j = 2:n-1, for i = j:n-1, ... ! 606: A(i,j) = 0; end; A(j,j) = j; end; A ! 607: ! 608: changes all but the "outer edge" of the lower triangle and ! 609: then prints the final matrix. ! 610: ! 611: for h = 1.0: -0.1: -1.0, (<h, cos(pi*h)>) ! 612: ! 613: prints a table of cosines. ! 614: ! 615: <X,D> = EIG(A); for v = X, v, A*v ! 616: ! 617: displays eigenvectors, one at a time. ! 618: .pp ! 619: The WHILE clause allows statements to be repeated an ! 620: indefinite number of times. The general form is ! 621: ! 622: WHILE expr relop expr, statement,..., statement, END ! 623: ! 624: where relop is =, <, >, <=, >=, or <> (not equal) . The statements ! 625: are repeatedly executed as long as the indicated comparison ! 626: between the real parts of the first components of the two ! 627: expressions is true. Here are two examples. (Exercise for ! 628: the reader: What do these segments do?) ! 629: ! 630: eps = 1; ! 631: while 1 + eps > 1, eps = eps/2; ! 632: eps = 2*eps ! 633: ! 634: E = 0*A; F = E + EYE; n = 1; ! 635: while NORM(E+F-E,1) > 0, E = E + F; F = A*F/n; n = n + 1; ! 636: E ! 637: ! 638: .pp ! 639: The IF clause allows conditional execution of statements. ! 640: The general form is ! 641: ! 642: IF expr relop expr, statement, ..., statement, ! 643: ELSE statement, ..., statement ! 644: ! 645: The first group of statements are executed if the relation is ! 646: true and the second group are executed if the relation is false. ! 647: The ELSE and the statements following it may be omitted. ! 648: For example, ! 649: ! 650: if abs(i-j) = 2, A(i,j) = 0; ! 651: ! 652: ! 653: 5. Commands, text, files and macros. ! 654: ! 655: .pp ! 656: MATLAB has several commands which control the output format ! 657: and the overall execution of the system. ! 658: .pp ! 659: The HELP command allows on-line access to short portions ! 660: of text describing various operations, functions and special ! 661: characters. The entire HELP document is reproduced in an ! 662: appendix. ! 663: .pp ! 664: Results are usually printed in a scaled fixed point format ! 665: that shows 4 or 5 significant figures. The commands SHORT, LONG, ! 666: SHORT E, LONG E and LONG Z alter the output format, but do not alter ! 667: the precision of the computations or the internal storage. ! 668: .pp ! 669: The WHO, WHAT and WHY commands provide information about ! 670: the functions and variables that are currently defined. ! 671: .pp ! 672: The CLEAR command erases all variables, except EPS, FLOP, RAND and EYE. ! 673: The statement A = <> indicates that a "0 by 0" ! 674: matrix is to be stored in A. This causes A to be erased so ! 675: that its storage can be used for other variables. ! 676: .pp ! 677: The RETURN and EXIT commands cause return to the underlying ! 678: operating system through the Fortran RETURN statement. ! 679: .pp ! 680: MATLAB has a limited facility for handling text. Any string ! 681: of characters delineated by quotes (with two quotes used to ! 682: allow one quote within the string) is saved as a vector of ! 683: integer values with '1' = 1, 'A' = 10, ' ' = 36, etc. ! 684: (The complete list is in the appendix under CHAR.) ! 685: For example ! 686: ! 687: '2*A + 3' is the same as <2 43 10 36 41 36 3> ! 688: ! 689: It is possible, though seldom very meaningful, to use such strings ! 690: in matrix operations. More frequently, the text is used as ! 691: a special argument to various functions. ! 692: ! 693: NORM(A,'inf') computes the infinity norm of A . ! 694: DISPLAY(T) prints the text stored in T . ! 695: EXEC('file') obtains MATLAB input from an external file. ! 696: SAVE('file') stores all the current variables in a file. ! 697: LOAD('file') retrieves all the variables from a file. ! 698: PRINT('file',X) prints X on a file. ! 699: DIARY('file') makes a copy of the complete MATLAB session. ! 700: ! 701: .pp ! 702: The text can also be used in a limited string substitution macro facility. ! 703: If a variable, say T, contains the source text for a MATLAB ! 704: statement or expression, then the construction ! 705: ! 706: > T < ! 707: ! 708: causes T to be executed or evaluated. For example ! 709: ! 710: T = '2*A + 3'; ! 711: S = 'B = >T< + 5' ! 712: A = 4; ! 713: > S < ! 714: ! 715: produces ! 716: ! 717: B = ! 718: ! 719: 16. ! 720: ! 721: Some other examples are given under MACRO in the appendix. ! 722: This facility is useful for fairly short statements and expressions. ! 723: More complicated MATLAB "programs" should use the EXEC facility. ! 724: .pp ! 725: The operations which access external files cannot be handled ! 726: in a completely machine-independent manner by portable Fortran ! 727: code. It is necessary for each particular installation to ! 728: provide a subroutine which associates external text files with ! 729: Fortran logical unit numbers. ! 730: ! 731: ! 732: 6. Census example ! 733: ! 734: .pp ! 735: Our first extended example involves predicting the population of the ! 736: United States in 1980 using extrapolation of ! 737: various fits to the census data from 1900 through 1970. ! 738: There are eight observations, so we begin with the MATLAB statement ! 739: .sp ! 740: n = 8 ! 741: .sp ! 742: The values of the dependent variable, the population in millions, ! 743: can be entered with ! 744: .sp ! 745: .ne 2 ! 746: y = < 75.995 91.972 105.711 123.203 ... ! 747: 131.669 150.697 179.323 203.212>' ! 748: .sp ! 749: In order to produce a reasonably scaled matrix, the independent variable, ! 750: time, is transformed from the interval [1900,1970] to [-1.00,0.75]. ! 751: This can be accomplished directly with ! 752: .sp ! 753: t = -1.0:0.25:0.75 ! 754: .sp ! 755: or in a fancier, but perhaps clearer, way with ! 756: .sp ! 757: t = 1900:10:1970; t = (t - 1940*ones(t))/40 ! 758: .sp ! 759: Either of these is equivalent to ! 760: .sp ! 761: t = <-1 -.75 -.50 -.25 0 .25 .50 .75> ! 762: .pp ! 763: The interpolating polynomial of degree n-1 involves an Vandermonde ! 764: matrix of order n with elements that might be generated by ! 765: .sp ! 766: .ne 2 ! 767: for i = 1:n, for j = 1:n, a(i,j) = t(i)**(j-1); ! 768: .sp ! 769: However, this results in an error caused by 0**0 when i = 5 and j = 1 . ! 770: The preferable approach is ! 771: .sp ! 772: .ne 2 ! 773: A = ones(n,n); ! 774: for i = 1:n, for j = 2:n, a(i,j) = t(i)*a(i,j-1); ! 775: .sp ! 776: Now the statement ! 777: .sp ! 778: cond(A) ! 779: .sp ! 780: produces the output ! 781: .sp ! 782: .ne 3 ! 783: ANS = ! 784: .sp ! 785: 1.1819E+03 ! 786: .sp ! 787: which indicates that transformation of the time variable has resulted ! 788: in a reasonably well conditioned matrix. ! 789: .pp ! 790: The statement ! 791: .sp ! 792: c = A\\y ! 793: .sp ! 794: results in ! 795: .sp ! 796: .ne 10 ! 797: C = ! 798: .sp ! 799: 131.6690 ! 800: 41.0406 ! 801: 103.5396 ! 802: 262.4535 ! 803: -326.0658 ! 804: -662.0814 ! 805: 341.9022 ! 806: 533.6373 ! 807: .sp ! 808: These are the coefficients in the interpolating polynomial ! 809: .sp ! 810: n-1 ! 811: c + c t + ... + c t ! 812: 1 2 n ! 813: .sp ! 814: Our transformation of the time variable has resulted in t = 1 ! 815: corresponding to the year 1980. Consequently, the extrapolated ! 816: population is simply the sum of the coefficients. ! 817: This can be computed by ! 818: .sp ! 819: p = sum(c) ! 820: .sp ! 821: The result is ! 822: .sp ! 823: .ne 3 ! 824: P = ! 825: .sp ! 826: 426.0950 ! 827: .sp ! 828: which indicates a 1980 population of over 426 million. ! 829: Clearly, using the seventh degree interpolating polynomial to ! 830: extrapolate even a fairly short distance beyond the end of the ! 831: data interval is not a good idea. ! 832: .pp ! 833: The coefficients in least squares fits by polynomials of lower degree ! 834: can be computed using fewer than n columns of the matrix. ! 835: .sp ! 836: for k = 1:n, c = A(:,1:k)\\y, p = sum(c) ! 837: .sp ! 838: would produce the coefficients of these fits, as well as the resulting ! 839: extrapolated population. ! 840: If we do not want to print all the coefficients, we can simply ! 841: generate a small table of populations predicted by polynomials ! 842: of degrees zero through seven. ! 843: We also compute the maximum deviation ! 844: between the fitted and observed values. ! 845: .sp ! 846: for k = 1:n, X = A(:,1:k); c = X\\y; ... ! 847: d(k) = k-1; p(k) = sum(c); e(k) = norm(X*c-y,'inf'); ! 848: <d, p, e> ! 849: .sp ! 850: The resulting output is ! 851: .sp ! 852: .ne 8 ! 853: 0 132.7227 70.4892 ! 854: 1 211.5101 9.8079 ! 855: 2 227.7744 5.0354 ! 856: 3 241.9574 3.8941 ! 857: 4 234.2814 4.0643 ! 858: 5 189.7310 2.5066 ! 859: 6 118.3025 1.6741 ! 860: 7 426.0950 0.0000 ! 861: .sp ! 862: The zeroth degree fit, 132.7 million, is the result of fitting a constant ! 863: to the data and is simply the average. ! 864: The results obtained with polynomials of degree one through four ! 865: all appear reasonable. ! 866: The maximum deviation of the degree four fit is slightly greater than the degree ! 867: three, even though the sum of the squares of the deviations is less. ! 868: The coefficients of the highest powers in the fits of degree ! 869: five and six turn out to be negative and the predicted populations ! 870: of less than 200 million are probably unrealistic. ! 871: The hopefully absurd prediction of the interpolating polynomial ! 872: concludes the table. ! 873: .pp ! 874: We wish to emphasize that ! 875: roundoff errors are not significant here. Nearly identical ! 876: results would be obtained on other computers, or with other algorithms. ! 877: The results simply indicate the difficulties associated with ! 878: extrapolation of polynomial fits of even modest degree. ! 879: .pp ! 880: A stabilized fit by a seventh degree polynomial can be obtained ! 881: using the pseudoinverse, but it requires a fairly delicate choice ! 882: of a tolerance. The statement ! 883: .sp ! 884: s = svd(A) ! 885: .sp ! 886: produces the singular values ! 887: .sp ! 888: .ne 10 ! 889: S = ! 890: .sp ! 891: 3.4594 ! 892: 2.2121 ! 893: 1.0915 ! 894: 0.4879 ! 895: 0.1759 ! 896: 0.0617 ! 897: 0.0134 ! 898: 0.0029 ! 899: .sp ! 900: We see that the last three singular values are less than 0.1 , ! 901: consequently, A can be approximately by a matrix of rank five ! 902: with an error less than 0.1 . ! 903: The Moore-Penrose pseudoinverse of this rank five matrix is ! 904: obtained from the singular value decomposition with the following ! 905: statements ! 906: .sp ! 907: c = pinv(A,0.1)*y, p = sum(c), e = norm(a*c-y,'inf') ! 908: .sp ! 909: The output is ! 910: .sp ! 911: .ne 18 ! 912: C = ! 913: .sp ! 914: 134.7972 ! 915: 67.5055 ! 916: 23.5523 ! 917: 9.2834 ! 918: 3.0174 ! 919: 2.6503 ! 920: -2.8808 ! 921: 3.2467 ! 922: .sp ! 923: P = ! 924: .sp ! 925: 241.1720 ! 926: .sp ! 927: E = ! 928: .sp ! 929: 3.9469 ! 930: .sp ! 931: The resulting seventh degree polynomial has coefficients which are much ! 932: smaller than those of the interpolating polynomial given earlier. ! 933: The predicted population and the maximum deviation are reasonable. ! 934: Any choice of the tolerance between the fifth and sixth singular ! 935: values would produce the same results, but choices outside this ! 936: range result in pseudoinverses of different rank and do not work ! 937: as well. ! 938: .pp ! 939: The one term exponential approximation ! 940: .sp ! 941: y(t) = k exp(pt) ! 942: .sp ! 943: can be transformed into a linear approximation by taking logarithms. ! 944: .sp ! 945: .ne 3 ! 946: log(y(t)) = log k + pt ! 947: .sp ! 948: = c + c t ! 949: 1 2 ! 950: .sp ! 951: The following segment makes use of the fact that a function of a vector is ! 952: the function applied to the individual components. ! 953: .sp ! 954: .ne 4 ! 955: X = A(:,1:2); ! 956: c = X\\log(y) ! 957: p = exp(sum(c)) ! 958: e = norm(exp(X*c)-y,'inf') ! 959: .sp ! 960: The resulting output is ! 961: .sp ! 962: .ne 12 ! 963: C = ! 964: .sp ! 965: 4.9083 ! 966: 0.5407 ! 967: .sp ! 968: P = ! 969: .sp ! 970: 232.5134 ! 971: .sp ! 972: E = ! 973: .sp ! 974: 4.9141 ! 975: .sp ! 976: The predicted population and maximum deviation appear satisfactory ! 977: and indicate that the exponential model is a reasonable one to consider. ! 978: .pp ! 979: As a curiousity, we return to the degree six polynomial. ! 980: Since the coefficient of the high order term is negative ! 981: and the value of the polynomial at t = 1 is positive, it must have ! 982: a root at some value of t greater than one. ! 983: The statements ! 984: .sp ! 985: X = A(:,1:7); ! 986: c = X\\y; ! 987: c = c(7:-1:1); //reverse the order of the coefficients ! 988: z = roots(c) ! 989: .sp ! 990: produce ! 991: .sp ! 992: Z = ! 993: .sp ! 994: 1.1023- 0.0000*i ! 995: 0.3021+ 0.7293*i ! 996: -0.8790+ 0.6536*i ! 997: -1.2939- 0.0000*i ! 998: -0.8790- 0.6536*i ! 999: 0.3021- 0.7293*i ! 1000: .sp ! 1001: There is only one real, positive root. ! 1002: The corresponding time on the original scale is ! 1003: .sp ! 1004: 1940 + 40*real(z(1)) ! 1005: .sp ! 1006: = 1984.091 ! 1007: .sp ! 1008: We conclude that ! 1009: the United States population should become zero early in ! 1010: February of 1984. ! 1011: .bp ! 1012: 7. Partial differential equation example ! 1013: ! 1014: .pp ! 1015: Our second extended example is a boundary value problem ! 1016: for Laplace's equation. ! 1017: The underlying physical problem involves the conductivity of a ! 1018: medium with cylindrical inclusions and is considered by ! 1019: Keller and Sachs [7]. ! 1020: .pp ! 1021: Find a function u(x,y) satisfying Laplace's equation ! 1022: .sp ! 1023: u + u = 0 ! 1024: xx yy ! 1025: .sp ! 1026: The domain is a unit square with a quarter circle of radius rho removed ! 1027: from one corner. ! 1028: There are Neumann conditions on the top and bottom edges ! 1029: and Dirichlet conditions on the remainder of the boundary. ! 1030: .sp 2 ! 1031: .ne 22 ! 1032: u = 0 ! 1033: n ! 1034: ! 1035: ------------- ! 1036: | . ! 1037: | . ! 1038: | . ! 1039: | . u = 1 ! 1040: | . ! 1041: | . ! 1042: | . ! 1043: u = 0 | | ! 1044: | | ! 1045: | | ! 1046: | | u = 1 ! 1047: | | ! 1048: | | ! 1049: | | ! 1050: ------------------------ ! 1051: ! 1052: u = 0 ! 1053: n ! 1054: .sp 2 ! 1055: The effective conductivity of an medium is then given by the integral ! 1056: along the left edge, ! 1057: .sp ! 1058: .ne 3 ! 1059: 1 ! 1060: sigma = integral u (0,y) dy ! 1061: 0 n ! 1062: .sp ! 1063: It is of interest to study the relation between the radius rho and ! 1064: the conductivity sigma. ! 1065: In particular, as rho approaches one, sigma becomes infinite. ! 1066: .pp ! 1067: Keller and Sachs use a finite difference approximation. ! 1068: The following technique makes use of the fact that the equation is ! 1069: actually Laplace's equation and leads to a much smaller matrix ! 1070: problem to solve. ! 1071: .pp ! 1072: Consider an approximate solution of the form ! 1073: .sp ! 1074: .ne 3 ! 1075: n 2j-1 ! 1076: u = sum c r cos(2j-1)t ! 1077: j=1 j ! 1078: .sp ! 1079: where r,t are polar coordinates (t is theta). The coefficients ! 1080: are to be determined. ! 1081: For any set of coefficients, this function already satisfies the ! 1082: differential equation because the basis functions are harmonic; ! 1083: it satisfies the normal derivative boundary condition on the bottom ! 1084: edge of the domain because we used cos t in preference to sin t ; ! 1085: and it satisfies the boundary condition on the left edge of the domain ! 1086: because we use only odd multiples of t . ! 1087: .pp ! 1088: The computational task is to find coefficients so that the ! 1089: boundary conditions on the remaining edges are satisfied as well as ! 1090: possible. ! 1091: To accomplish this, pick m points (r,t) on the remaining edges. ! 1092: It is desirable to have m > n and in practice we usually choose ! 1093: m to be two or three times as large as n . ! 1094: Typical values of n are 10 or 20 and of m are 20 to 60. ! 1095: An m by n matrix A is generated. The i,j element is the j-th ! 1096: basis function, or its normal derivative, evaluated at the i-th boundary ! 1097: point. ! 1098: A right hand side with m components is also generated. In this example, ! 1099: the elements of the right hand side are either zero or one. ! 1100: The coefficients are then found by solving the overdetermined set of ! 1101: equations ! 1102: .sp ! 1103: Ac = b ! 1104: .sp ! 1105: in a least squares sense. ! 1106: .pp ! 1107: Once the coefficients have been determined, the approximate solution ! 1108: is defined everywhere on the domain. It is then possible to ! 1109: compute the effective conductivity sigma . In fact, a very simple ! 1110: formula results, ! 1111: .sp ! 1112: .ne 3 ! 1113: n j-1 ! 1114: sigma = sum (-1) c ! 1115: j=1 j ! 1116: .pp ! 1117: To use MATLAB for this problem, the following "program" is first stored in ! 1118: the local computer file system, say under the name "PDE". ! 1119: .bp ! 1120: .nf ! 1121: //Conductivity example. ! 1122: //Parameters --- ! 1123: rho //radius of cylindrical inclusion ! 1124: n //number of terms in solution ! 1125: m //number of boundary points ! 1126: //initialize operation counter ! 1127: flop = <0 0>; ! 1128: //initialize variables ! 1129: m1 = round(m/3); //number of points on each straight edge ! 1130: m2 = m - m1; //number of points with Dirichlet conditions ! 1131: pi = 4*atan(1); ! 1132: //generate points in Cartesian coordinates ! 1133: //right hand edge ! 1134: for i = 1:m1, x(i) = 1; y(i) = (1-rho)*(i-1)/(m1-1); ! 1135: //top edge ! 1136: for i = m2+1:m, x(i) = (1-rho)*(m-i)/(m-m2-1); y(i) = 1; ! 1137: //circular edge ! 1138: for i = m1+1:m2, t = pi/2*(i-m1)/(m2-m1+1); ... ! 1139: x(i) = 1-rho*sin(t); y(i) = 1-rho*cos(t); ! 1140: //convert to polar coordinates ! 1141: for i = 1:m-1, th(i) = atan(y(i)/x(i)); ... ! 1142: r(i) = sqrt(x(i)**2+y(i)**2); ! 1143: th(m) = pi/2; r(m) = 1; ! 1144: //generate matrix ! 1145: //Dirichlet conditions ! 1146: for i = 1:m2, for j = 1:n, k = 2*j-1; ... ! 1147: a(i,j) = r(i)**k*cos(k*th(i)); ! 1148: //Neumann conditions ! 1149: for i = m2+1:m, for j = 1:n, k = 2*j-1; ... ! 1150: a(i,j) = k*r(i)**(k-1)*sin((k-1)*th(i)); ! 1151: //generate right hand side ! 1152: for i = 1:m2, b(i) = 1; ! 1153: for i = m2+1:m, b(i) = 0; ! 1154: //solve for coefficients ! 1155: c = A\\b ! 1156: //compute effective conductivity ! 1157: c(2:2:n) = -c(2:2:n); ! 1158: sigma = sum(c) ! 1159: //output total operation count ! 1160: ops = flop(2) ! 1161: .fi ! 1162: .sp 3 ! 1163: .pp ! 1164: The program can be used within MATLAB by setting the three ! 1165: parameters and then accessing the file. ! 1166: For example, ! 1167: .sp ! 1168: .ne 4 ! 1169: rho = .9; ! 1170: n = 15; ! 1171: m = 30; ! 1172: exec('PDE') ! 1173: .sp ! 1174: The resulting output is ! 1175: .sp ! 1176: RHO = ! 1177: .sp ! 1178: .9000 ! 1179: .sp ! 1180: N = ! 1181: .sp ! 1182: 15. ! 1183: .sp ! 1184: M = ! 1185: .sp ! 1186: 30. ! 1187: .sp ! 1188: .ne 17 ! 1189: C = ! 1190: .sp ! 1191: 2.2275 ! 1192: -2.2724 ! 1193: 1.1448 ! 1194: 0.1455 ! 1195: -0.1678 ! 1196: -0.0005 ! 1197: -0.3785 ! 1198: 0.2299 ! 1199: 0.3228 ! 1200: -0.2242 ! 1201: -0.1311 ! 1202: 0.0924 ! 1203: 0.0310 ! 1204: -0.0154 ! 1205: -0.0038 ! 1206: .sp ! 1207: SIGM = ! 1208: .sp ! 1209: 5.0895 ! 1210: .sp ! 1211: OPS = ! 1212: .sp ! 1213: 16204. ! 1214: .sp ! 1215: .pp ! 1216: A total of 16204 floating point operations were necessary to ! 1217: set up the matrix, solve for the coefficients and compute the ! 1218: conductivity. ! 1219: The operation count is roughly proportional to m*n**2. ! 1220: The results obtained for sigma as a function of rho ! 1221: by this approach are essentially the ! 1222: same as those obtained by the finite difference technique ! 1223: of Keller and Sachs, but the computational effort involved ! 1224: is much less. ! 1225: .bp ! 1226: 8. Eigenvalue sensitivity example ! 1227: ! 1228: .pp ! 1229: In this example, we construct a matrix whose eigenvalues ! 1230: are moderately sensitive to perturbations and then analyze ! 1231: that sensitivity. ! 1232: We begin with the statement ! 1233: ! 1234: B = <3 0 7; 0 2 0; 0 0 1> ! 1235: ! 1236: which produces ! 1237: ! 1238: B = ! 1239: ! 1240: 3. 0. 7. ! 1241: 0. 2. 0. ! 1242: 0. 0. 1. ! 1243: ! 1244: .pp ! 1245: Obviously, the eigenvalues of B are 1, 2 and 3 . ! 1246: Moreover, since B is not symmetric, ! 1247: these eigenvalues are slightly sensitive to perturbation. ! 1248: (The value b(1,3) = 7 was chosen so that the elements ! 1249: of the matrix A below are less than 1000.) ! 1250: .pp ! 1251: We now generate a similarity transformation to ! 1252: disguise the eigenvalues and make them more sensitive. ! 1253: ! 1254: L = <1 0 0; 2 1 0; -3 4 1>, M = L\\L' ! 1255: ! 1256: L = ! 1257: ! 1258: 1. 0. 0. ! 1259: 2. 1. 0. ! 1260: -3. 4. 1. ! 1261: ! 1262: M = ! 1263: ! 1264: 1.0000 2.0000 -3.0000 ! 1265: -2.0000 -3.0000 10.0000 ! 1266: 11.0000 18.0000 -48.0000 ! 1267: ! 1268: The matrix M has determinant equal to 1 and is ! 1269: moderately badly conditioned. ! 1270: The similarity transformation is ! 1271: ! 1272: A = M*B/M ! 1273: ! 1274: A = ! 1275: ! 1276: -64.0000 82.0000 21.0000 ! 1277: 144.0000 -178.0000 -46.0000 ! 1278: -771.0000 962.0000 248.0000 ! 1279: ! 1280: Because det(M) = 1 , the elements of A would be exact ! 1281: integers if there were no roundoff. ! 1282: So, ! 1283: ! 1284: A = round(A) ! 1285: ! 1286: A = ! 1287: ! 1288: -64. 82. 21. ! 1289: 144. -178. -46. ! 1290: -771. 962. 248. ! 1291: ! 1292: .pp ! 1293: This, then, is our test matrix. We can now forget how ! 1294: it was generated and analyze its eigenvalues. ! 1295: ! 1296: <X,D> = eig(A) ! 1297: ! 1298: D = ! 1299: ! 1300: 3.0000 0.0000 0.0000 ! 1301: 0.0000 1.0000 0.0000 ! 1302: 0.0000 0.0000 2.0000 ! 1303: ! 1304: X = ! 1305: ! 1306: -.0891 3.4903 41.8091 ! 1307: .1782 -9.1284 -62.7136 ! 1308: -.9800 46.4473 376.2818 ! 1309: ! 1310: Since A is similar to B, its eigenvalues are also 1, 2 and 3. ! 1311: They happen to be computed in another order by the EISPACK subroutines. ! 1312: The fact that the columns of X, which are the eigenvectors, ! 1313: are so far from being orthonormal is our first indication that ! 1314: the eigenvalues are sensitive. ! 1315: To see this sensitivity, we display more figures of the computed ! 1316: eigenvalues. ! 1317: ! 1318: long, diag(D) ! 1319: ! 1320: ANS = ! 1321: ! 1322: 2.999999999973599 ! 1323: 1.000000000015625 ! 1324: 2.000000000011505 ! 1325: ! 1326: We see that, on this computer, the last five significant figures are ! 1327: contaminated by roundoff error. ! 1328: A somewhat superficial explanation of this is provided by ! 1329: ! 1330: short, cond(X) ! 1331: ! 1332: ANS = ! 1333: ! 1334: 3.2216e+05 ! 1335: ! 1336: The condition number of X gives an upper bound for the relative error in ! 1337: the computed eigenvalues. ! 1338: However, this condition number is affected by scaling. ! 1339: ! 1340: X = X/diag(X(3,:)), cond(X) ! 1341: ! 1342: X = ! 1343: ! 1344: .0909 .0751 .1111 ! 1345: -.1818 -.1965 -.1667 ! 1346: 1.0000 1.0000 1.0000 ! 1347: ! 1348: ANS = ! 1349: ! 1350: 1.7692e+03 ! 1351: ! 1352: .pp ! 1353: Rescaling the eigenvectors so that their last components are ! 1354: all equal to one has two consequences. ! 1355: The condition of X is decreased by over two orders of magnitude. ! 1356: (This is about the minimum condition that can be obtained by ! 1357: such diagonal scaling.) ! 1358: Moreover, it is now apparent that the three eigenvectors are ! 1359: nearly parallel. ! 1360: .pp ! 1361: More detailed information on the sensitivity of the individual ! 1362: eigenvalues involves the left eigenvectors. ! 1363: ! 1364: Y = inv(X'), Y'*A*X ! 1365: ! 1366: Y = ! 1367: ! 1368: -511.5000 259.5000 252.0000 ! 1369: 616.0000 -346.0000 -270.0000 ! 1370: 159.5000 -86.5000 -72.0000 ! 1371: ! 1372: ANS = ! 1373: ! 1374: 3.0000 .0000 .0000 ! 1375: .0000 1.0000 .0000 ! 1376: .0000 .0000 2.0000 ! 1377: ! 1378: We are now in a position to compute the sensitivities of the ! 1379: individual eigenvalues. ! 1380: ! 1381: for j = 1:3, c(j) = norm(Y(:,j))*norm(X(:,j)); end, C ! 1382: ! 1383: C = ! 1384: ! 1385: 833.1092 ! 1386: 450.7228 ! 1387: 383.7564 ! 1388: ! 1389: These three numbers are the reciprocals of ! 1390: the cosines of the angles between the ! 1391: left and right eigenvectors. ! 1392: It can be shown that perturbation of the elements of A ! 1393: can result in a perturbation of the j-th eigenvalue which ! 1394: is c(j) times as large. ! 1395: In this example, the first eigenvalue has the largest sensitivity. ! 1396: .pp ! 1397: We now proceed to show that A is close to a matrix with ! 1398: a double eigenvalue. The direction of the required perturbation ! 1399: is given by ! 1400: ! 1401: E = -1.e-6*Y(:,1)*X(:,1)' ! 1402: ! 1403: E = ! 1404: ! 1405: 1.0e-03 * ! 1406: ! 1407: .0465 -.0930 .5115 ! 1408: -.0560 .1120 -.6160 ! 1409: -.0145 .0290 -.1595 ! 1410: ! 1411: With some trial and error which we do not show, we bracket the point ! 1412: where two eigenvalues of a perturbed A coalesce and then become ! 1413: complex. ! 1414: ! 1415: eig(A + .4*E), eig(A + .5*E) ! 1416: ! 1417: ANS = ! 1418: ! 1419: 1.1500 ! 1420: 2.5996 ! 1421: 2.2504 ! 1422: ! 1423: ANS = ! 1424: ! 1425: 2.4067 + .1753*i ! 1426: 2.4067 - .1753*i ! 1427: 1.1866 + 0.0000*i ! 1428: ! 1429: Now, a bisecting search, driven by the imaginary part of one ! 1430: of the eigenvalues, finds the point where two eigenvalues ! 1431: are nearly equal. ! 1432: ! 1433: r = .4; s = .5; ! 1434: ! 1435: while s-r > 1.e-14, t = (r+s)/2; d = eig(A+t*E); ... ! 1436: if imag(d(1))=0, r = t; else, s = t; ! 1437: ! 1438: long, T ! 1439: ! 1440: T = ! 1441: ! 1442: .450380734134507 ! 1443: ! 1444: .pp ! 1445: Finally, we display the perturbed matrix, which is obviously ! 1446: close to the original, and its pair of nearly equal eigenvalues. ! 1447: (We have dropped a few digits from the long output.) ! 1448: ! 1449: A+t*E, eig(A+t*E) ! 1450: ! 1451: A ! 1452: ! 1453: -63.999979057 81.999958114 21.000230369 ! 1454: 143.999974778 -177.999949557 -46.000277434 ! 1455: -771.000006530 962.000013061 247.999928164 ! 1456: ! 1457: ANS = ! 1458: ! 1459: 2.415741150 ! 1460: 2.415740621 ! 1461: 1.168517777 ! 1462: ! 1463: .pp ! 1464: The first two eigenvectors of A + t*E ! 1465: are almost indistinguishable indicating that ! 1466: the perturbed matrix is almost defective. ! 1467: ! 1468: <X,D> = eig(A+t*E); X = X/diag(X(3,:)) ! 1469: ! 1470: X = ! 1471: ! 1472: .096019578 .096019586 .071608466 ! 1473: -.178329614 -.178329608 -.199190520 ! 1474: 1.000000000 1.000000000 1.000000000 ! 1475: ! 1476: short, cond(X) ! 1477: ! 1478: ANS = ! 1479: ! 1480: 3.3997e+09 ! 1481: ! 1482: ! 1483: 9. Syntax diagrams ! 1484: ! 1485: .pp ! 1486: A formal description of the language acceptable to MATLAB, ! 1487: as well as a flow chart of the MATLAB program, is provided by ! 1488: the syntax diagrams or syntax graphs of Wirth [6]. There are ! 1489: eleven non-terminal symbols in the language: ! 1490: ! 1491: line, statement, clause, expression, term, ! 1492: factor, number, integer, name, command, text . ! 1493: ! 1494: The diagrams define each ! 1495: of the non-terminal symbols using the others and the terminal symbols: ! 1496: ! 1497: letter -- A through Z, ! 1498: digit -- 0 through 9, ! 1499: char -- ( ) ; : + - * / \\ = . , < > ! 1500: quote -- ' ! 1501: ! 1502: ! 1503: .nf ! 1504: .ne 17 ! 1505: line ! 1506: ! 1507: |-----> statement >----| ! 1508: | | ! 1509: |-----> clause >-------| ! 1510: | | ! 1511: -------|-----> expr >---------|------> ! 1512: | | | | ! 1513: | |-----> command >------| | ! 1514: | | | | ! 1515: | |-> > >-> expr >-> < >-| | ! 1516: | | | | ! 1517: | |----------------------| | ! 1518: | | ! 1519: | |-< ; <-| | ! 1520: |--------| |---------| ! 1521: |-< , <-| ! 1522: ! 1523: ! 1524: ! 1525: ! 1526: .ne 13 ! 1527: statement ! 1528: ! 1529: |-> name >--------------------------------| ! 1530: | | | ! 1531: | | |--> : >---| | ! 1532: | | | | | ! 1533: | |-> ( >---|-> expr >-|---> ) >-| ! 1534: | | | | ! 1535: -----| |-----< , <----| |--> = >--> expr >---> ! 1536: | | ! 1537: | |--< , <---| | ! 1538: | | | | ! 1539: |-> < >---> name >---> > >----------------| ! 1540: ! 1541: ! 1542: ! 1543: ! 1544: .ne 14 ! 1545: clause ! 1546: ! 1547: |---> FOR >---> name >---> = >---> expr >--------------| ! 1548: | | ! 1549: | |-> WHILE >-| | ! 1550: |-| |-> expr >---------------------- | ! 1551: | |-> IF >-| | | | | | | | ! 1552: -----| < <= = <> >= > |----> ! 1553: | | | | | | | | ! 1554: | ----------------------> expr >--| ! 1555: | | ! 1556: |---> ELSE >--------------------------------------------| ! 1557: | | ! 1558: |---> END >--------------------------------------------| ! 1559: ! 1560: ! 1561: ! 1562: ! 1563: .ne 11 ! 1564: expr ! 1565: ! 1566: |-> + >-| ! 1567: | | ! 1568: -------|-------|-------> term >----------> ! 1569: | | | | ! 1570: |-> - >-| | |-< + <-| | ! 1571: | | | | ! 1572: |--|-< - <-|--| ! 1573: | | ! 1574: |-< : <-| ! 1575: ! 1576: ! 1577: ! 1578: ! 1579: .ne 9 ! 1580: term ! 1581: ! 1582: ---------------------> factor >----------------------> ! 1583: | | ! 1584: | |-< * <-| | ! 1585: | |-------| | | |-------| | ! 1586: |--| |--|-< / <-|--| |--| ! 1587: |-< . <-| | | |-< . <-| ! 1588: |-< \\ <-| ! 1589: ! 1590: ! 1591: ! 1592: ! 1593: .ne 29 ! 1594: factor ! 1595: ! 1596: |----------------> number >---------------| ! 1597: | | ! 1598: |-> name >--------------------------------| ! 1599: | | | ! 1600: | | |--> : >---| | ! 1601: | | | | | ! 1602: | |-> ( >---|-> expr >-|---> ) >-| ! 1603: | | | | ! 1604: | |-----< , <----| | ! 1605: | | ! 1606: -----|------------> ( >-----> expr >-----> ) >-|-|-------|-----> ! 1607: | | | | | ! 1608: | |--------------| | |-> ' >-| | ! 1609: | | | | | ! 1610: |------------> < >-|---> expr >---|-> > >-| | ! 1611: | | | | | ! 1612: | |--< <---| | | ! 1613: | | | | | ! 1614: | |--< ; <---| | | ! 1615: | | | | | ! 1616: | |--< , <---| | | ! 1617: | | | ! 1618: |------------> > >-----> expr >-----> < >-| | ! 1619: | | | ! 1620: |-----> factor >---> ** >---> factor >----| | ! 1621: | | ! 1622: |------------> ' >-----> text >-----> ' >-------------| ! 1623: ! 1624: ! 1625: ! 1626: ! 1627: .ne 9 ! 1628: number ! 1629: ! 1630: |----------| |-> + >-| ! 1631: | | | | ! 1632: -----> int >-----> . >---> int >-----> E >---------> int >----> ! 1633: | | | | | | ! 1634: | | | |-> - >-| | ! 1635: | | | | ! 1636: |---------------------------------------------| ! 1637: ! 1638: ! 1639: ! 1640: ! 1641: .ne 5 ! 1642: int ! 1643: ! 1644: ------------> digit >-----------> ! 1645: | | ! 1646: |-----------| ! 1647: ! 1648: ! 1649: ! 1650: ! 1651: .ne 7 ! 1652: name ! 1653: ! 1654: |--< letter <--| ! 1655: | | ! 1656: ------> letter >--|--------------|-----> ! 1657: | | ! 1658: |--< digit <--| ! 1659: ! 1660: ! 1661: ! 1662: ! 1663: .ne 9 ! 1664: command ! 1665: ! 1666: |--> name >--| ! 1667: | | ! 1668: --------> name >--------|------------|----> ! 1669: | | ! 1670: |--> char >--| ! 1671: | | ! 1672: |---> ' >----| ! 1673: ! 1674: .ne 11 ! 1675: text ! 1676: ! 1677: |-> letter >--| ! 1678: | | ! 1679: |-> digit >---| ! 1680: ----------------| |--------------> ! 1681: | |-> char >----| | ! 1682: | | | | ! 1683: | |-> ' >-> ' >-| | ! 1684: | | ! 1685: |---------------------| ! 1686: ! 1687: .fi ! 1688: ! 1689: 10. The parser-interpreter ! 1690: ! 1691: .pp ! 1692: The structure of the parser-interpreter is similar to that ! 1693: of Wirth's compiler [6] for his simple language, PL/0 , except that ! 1694: MATLAB is programmed in Fortran, which does not have explicit ! 1695: recursion. The interrelation of the primary subroutines is ! 1696: shown in the following diagram. ! 1697: ! 1698: .ne 48 ! 1699: MAIN ! 1700: | ! 1701: MATLAB |--CLAUSE ! 1702: | | | ! 1703: PARSE-----|--EXPR----TERM----FACTOR ! 1704: | | | | ! 1705: | |-------|-------| ! 1706: | | | | ! 1707: | STACK1 STACK2 STACKG ! 1708: | ! 1709: |--STACKP--PRINT ! 1710: | ! 1711: |--COMAND ! 1712: | ! 1713: | ! 1714: | |--CGECO ! 1715: | | ! 1716: | |--CGEFA ! 1717: | | ! 1718: |--MATFN1--|--CGESL ! 1719: | | ! 1720: | |--CGEDI ! 1721: | | ! 1722: | |--CPOFA ! 1723: | ! 1724: | ! 1725: | |--IMTQL2 ! 1726: | | ! 1727: | |--HTRIDI ! 1728: | | ! 1729: |--MATFN2--|--HTRIBK ! 1730: | | ! 1731: | |--CORTH ! 1732: | | ! 1733: | |--COMQR3 ! 1734: | ! 1735: | ! 1736: |--MATFN3-----CSVDC ! 1737: | ! 1738: | ! 1739: | |--CQRDC ! 1740: |--MATFN4--| ! 1741: | |--CQRSL ! 1742: | ! 1743: | ! 1744: | |--FILES ! 1745: |--MATFN5--| ! 1746: |--SAVLOD ! 1747: .pp ! 1748: Subroutine PARSE controls the interpretation of each statement. ! 1749: It calls subroutines that process the various syntactic quantities ! 1750: such as command, expression, term and factor. A fairly simple ! 1751: program stack mechanism allows these subroutines to recursively ! 1752: "call" each other along the lines allowed by the syntax diagrams. ! 1753: The four STACK subroutines manage the variable memory and perform ! 1754: elementary operations, such as matrix addition and transposition. ! 1755: .pp ! 1756: The four subroutines MATFN1 though MATFN4 are called ! 1757: whenever "serious" matrix computations are required. They ! 1758: are interface routines which call the various LINPACK and ! 1759: EISPACK subroutines. MATFN5 primarily handles the file ! 1760: access tasks. ! 1761: .pp ! 1762: Two large real arrays, STKR and STKI, are used to store all the ! 1763: matrices. Four integer arrays are used to store the names, the row ! 1764: and column dimensions, and the pointers into the real stacks. ! 1765: The following diagram illustrates this storage scheme. ! 1766: .sp 1 ! 1767: .ll 75 ! 1768: .ne 34 ! 1769: .nf ! 1770: TOP IDSTK MSTK NSTK LSTK STKR STKI ! 1771: -- -- -- -- -- -- -- -- -------- -------- ! 1772: | |--->| | | | | | | | | | |----------->| | | | ! 1773: -- -- -- -- -- -- -- -- -------- -------- ! 1774: | | | | | | | | | | | | | | | ! 1775: -- -- -- -- -- -- -- -------- -------- ! 1776: . . . . . . ! 1777: . . . . . . ! 1778: . . . . . . ! 1779: -- -- -- -- -- -- -- -------- -------- ! 1780: BOT | | | | | | | | | | | | | | | ! 1781: -- -- -- -- -- -- -- -- -------- -------- ! 1782: | |--->| X| | | | | 2| | 1| | |----------->| 3.14 | | 0.00 | ! 1783: -- -- -- -- -- -- -- -- -------- -------- ! 1784: | A| | | | | 2| | 2| | |--------- | 0.00 | | 1.00 | ! 1785: -- -- -- -- -- -- -- \\ -------- -------- ! 1786: | E| P| S| | | 1| | 1| | |------- ->| 11.00 | | 0.00 | ! 1787: -- -- -- -- -- -- -- \\ -------- -------- ! 1788: | F| L| O| P| | 1| | 2| | |------ \\ | 21.00 | | 0.00 | ! 1789: -- -- -- -- -- -- -- \\ \\ -------- -------- ! 1790: | E| Y| E| | |-1| |-1| | |--- \\ | | 12.00 | | 0.00 | ! 1791: -- -- -- -- -- -- -- \\ | | -------- -------- ! 1792: | R| A| N| D| | 1| | 1| | |- \\ | | | 22.00 | | 0.00 | ! 1793: -- -- -- -- -- -- -- \\ | \\ \\ -------- -------- ! 1794: | \\ \\ ->| 1.E-15 | | 0.00 | ! 1795: \\ \\ \\ -------- -------- ! 1796: \\ \\ ->| 0.00 | | 0.00 | ! 1797: \\ \\ -------- -------- ! 1798: \\ \\ | 0.00 | | 0.00 | ! 1799: \\ \\ -------- -------- ! 1800: \\ ->| 1.00 | | 0.00 | ! 1801: \\ -------- -------- ! 1802: --->| URAND | | 0.00 | ! 1803: -------- -------- ! 1804: .fi ! 1805: .ll ! 1806: .pp ! 1807: The top portion of the stack is used for temporary variables ! 1808: and the bottom portion for saved variables. The figure shows the ! 1809: situation after the line ! 1810: ! 1811: A = <11,12; 21,22>, x = <3.14, sqrt(-1)>' ! 1812: ! 1813: has been processed. The four permanent names, EPS, FLOP, RAND and EYE, ! 1814: occupy the last four positions of the variable stacks. RAND ! 1815: has dimensions 1 by 1, but whenever its value is requested, ! 1816: a random number generator is used instead. EYE has dimensions ! 1817: -1 by -1 to indicate that the actual dimensions must be ! 1818: determined later by context. The two saved variables have ! 1819: dimensions 2 by 2 and 2 by 1 and so take up a total of 6 ! 1820: locations. ! 1821: .pp ! 1822: Subsequent statements involving A and x will result in ! 1823: temporary copies being made in the top of the stack for use in ! 1824: the actual calculations. Whenever the top of the stack reaches ! 1825: the bottom, a message indicating memory has been ! 1826: exceeded is printed, but the current variables are not ! 1827: affected. ! 1828: .pp ! 1829: This modular structure makes it possible to implement ! 1830: MATLAB on a system with a limited amount of memory. The object ! 1831: code for the MATFN's and the LINPACK-EISPACK subroutines ! 1832: is rarely needed. Although it is not standard, many Fortran ! 1833: operating systems provide some overlay mechanism so that ! 1834: this code is brought into the main memory only when required. ! 1835: The variables, which occupy a relatively small portion of the ! 1836: memory, remain in place, while the subroutines which process ! 1837: them are loaded a few at a time. ! 1838: ! 1839: ! 1840: 11. The numerical algorithms ! 1841: ! 1842: .pp ! 1843: The algorithms underlying the basic MATLAB functions are ! 1844: described in the LINPACK and EISPACK guides [1-3]. ! 1845: The following list gives the subroutines used by these functions. ! 1846: ! 1847: INV(A) - CGECO,CGEDI ! 1848: DET(A) - CGECO,CGEDI ! 1849: LU(A) - CGEFA ! 1850: RCOND(A) - CGECO ! 1851: CHOL(A) - CPOFA ! 1852: SVD(A) - CSVDC ! 1853: COND(A) - CSVDC ! 1854: NORM(A,2) - CSVDC ! 1855: PINV(A,eps) - CSVDC ! 1856: RANK(A,eps) - CSVDC ! 1857: QR(A) - CQRDC,CQRSL ! 1858: ORTH(A) - CQRDC,CSQSL ! 1859: A\\B and B/A - CGECO,CGESL if A is square. ! 1860: - CQRDC,CQRSL if A is not square. ! 1861: EIG(A) - HTRIDI,IMTQL2,HTRIBK if A is Hermitian. ! 1862: - CORTH,COMQR2 if A is not Hermitian. ! 1863: SCHUR(A) - same as EIG. ! 1864: HESS(A) - same as EIG. ! 1865: ! 1866: .pp ! 1867: Minor modifications were made to all these subroutines. The LINPACK ! 1868: routines were changed to replace the Fortran complex arithmetic with ! 1869: explicit references to real and imaginary parts. Since most of the ! 1870: floating point arithmetic is concentrated in a few low-level subroutines ! 1871: which perform vector operations (the Basic Linear Algebra Subprograms), ! 1872: this was not an extensive change. It also facilitated implementation ! 1873: of the FLOP and CHOP features which count and optionally truncate each ! 1874: floating point operation. ! 1875: .pp ! 1876: The EISPACK subroutine COMQR2 was modified to allow access to the ! 1877: Schur triangular form, ordinarily just an intermediate result. IMTQL2 ! 1878: was modified to make computation of the eigenvectors optional. Both ! 1879: subroutines were modified to eliminate the machine-dependent accuracy ! 1880: parameter and all the EISPACK subroutines were changed to include FLOP ! 1881: and CHOP. ! 1882: .pp ! 1883: The algorithms employed for the POLY and ROOTS functions illustrate ! 1884: an interesting aspect of the modern approach to eigenvalue computation. ! 1885: POLY(A) generates the characteristic polynomial of A and ROOTS(POLY(A)) ! 1886: finds the roots of that polynomial, which are, of course, the eigenvalues ! 1887: of A . ! 1888: But both POLY and ROOTS use EISPACK eigenvalues subroutines, ! 1889: which are based on similarity transformations. ! 1890: So the classical approach which characterizes ! 1891: eigenvalues as roots of the characteristic polynomial is actually ! 1892: reversed. ! 1893: .pp ! 1894: If A is an n by n matrix, POLY(A) produces the coefficients ! 1895: C(1) through C(n+1), with C(1) = 1, in ! 1896: ! 1897: DET(z*EYE-A) = C(1)*z**n + ... + C(n)*z + C(n+1) . ! 1898: ! 1899: The algorithm can be expressed compactly using MATLAB: ! 1900: ! 1901: Z = EIG(A); ! 1902: C = 0*ONES(n+1,1); C(1) = 1; ! 1903: for j = 1:n, C(2:j+1) = C(2:j+1) - Z(j)*C(1:j); ! 1904: C ! 1905: ! 1906: This recursion is easily derived by expanding the product ! 1907: ! 1908: (z - z(1))*(z - z(2))* ... * (z-z(n)) . ! 1909: ! 1910: It is possible to prove that POLY(A) produces the coefficients in the ! 1911: characteristic polynomial of a matrix within roundoff error of A . ! 1912: This is true even if the eigenvalues of A are badly conditioned. ! 1913: The traditional algorithms for obtaining the characteristic polynomial ! 1914: which do not use the eigenvalues do not have such satisfactory numerical ! 1915: properties. ! 1916: .pp ! 1917: If C is a vector with n+1 components, ROOTS(C) finds the roots ! 1918: of the polynomial of degree n , ! 1919: ! 1920: p(z) = C(1)*z**n + ... + C(n)*z + C(n+1) . ! 1921: ! 1922: The algorithm simply involves computing the eigenvalues of the companion ! 1923: matrix: ! 1924: ! 1925: A = 0*ONES(n,n) ! 1926: for j = 1:n, A(1,j) = -C(j+1)/C(1); ! 1927: for i = 2:n, A(i,i-1) = 1; ! 1928: EIG(A) ! 1929: ! 1930: It is possible to prove that the results produced are the exact eigenvalues ! 1931: of a matrix within roundoff error of the companion matrix A, but this does ! 1932: not mean that they are the exact roots of a polynomial with coefficients ! 1933: within roundoff error of those in C . There are more accurate, more ! 1934: efficient methods for finding polynomial roots, but this approach has the ! 1935: crucial advantage that it does not require very much additional code. ! 1936: .pp ! 1937: The elementary functions EXP, LOG, SQRT, SIN, COS and ATAN ! 1938: are applied to square matrices by diagonalizing the matrix, applying ! 1939: the functions to the individual eigenvalues and then transforming ! 1940: back. For example, EXP(A) is computed by ! 1941: ! 1942: <X,D> = EIG(A); ! 1943: for j = 1:n, D(j,j) = EXP(D(j,j)); ! 1944: X*D/X ! 1945: ! 1946: This is essentially method number 14 out of the 19 'dubious' ! 1947: possibilities described in [8]. It is dubious because it doesn't ! 1948: always work. The matrix of eigenvectors X can be arbitrarily badly ! 1949: conditioned and all accuracy lost in the computation of X*D/X. ! 1950: A warning message is printed if RCOND(X) is very small, but this ! 1951: only catches the extreme cases. An example of a case not detected ! 1952: is when A has a double eigenvalue, but theoretically only one linearly ! 1953: independent eigenvector associated with it. The computed eigenvalues ! 1954: will be separated by something on the order of the square root of ! 1955: the roundoff level. This separation will be reflected in RCOND(X) ! 1956: which will probably not be small enough to trigger the error message. ! 1957: The computed EXP(A) will be accurate to only half precision. Better ! 1958: methods are known for computing EXP(A), but they do not easily extend ! 1959: to the other five functions and would require a considerable amount ! 1960: of additional code. ! 1961: .pp ! 1962: The expression A**p is evaluated by repeated multiplication if ! 1963: p is an integer greater than 1. Otherwise it is evaluated by ! 1964: ! 1965: <X,D> = EIG(A); ! 1966: for j = 1:n, D(j,j) = EXP(p*LOG(D(j,j))) ! 1967: X*D/X ! 1968: ! 1969: This suffers from the same potential loss of accuracy if X is badly ! 1970: conditioned. It was partly for this reason that the case p = 1 is ! 1971: included in the general case. Comparison of A**1 with A gives some ! 1972: idea of the loss of accuracy for other values of p and for the ! 1973: elementary functions. ! 1974: .pp ! 1975: RREF, the reduced row echelon form, is of some interest in ! 1976: theoretical linear algebra, although it has little computational ! 1977: value. It is included in MATLAB for pedagogical reasons. The ! 1978: algorithm is essentially Gauss-Jordan elimination with detection ! 1979: of negligible columns applied to rectangular matrices. ! 1980: .pp ! 1981: There are three separate places in MATLAB where the rank of a ! 1982: matrix is implicitly computed: ! 1983: in RREF(A), in A\\B for non-square A, and in ! 1984: the pseudoinverse PINV(A). Three different algorithms with three ! 1985: different criteria for negligibility are used and so it is possible ! 1986: that three different values could be produced for the same matrix. ! 1987: With RREF(A), the rank of A is the number of nonzero rows. The ! 1988: elimination algorithm used for RREF is the fastest of the three ! 1989: rank-determining algorithms, but it is the least sophisticated ! 1990: numerically and the least reliable. With A\\B, the algorithm is ! 1991: essentially that used by example subroutine SQRST in chapter 9 of the ! 1992: LINPACK guide. With PINV(A), the algorithm is based on the singular ! 1993: value decomposition and is described in chapter 11 of the LINPACK ! 1994: guide. The SVD algorithm is the most time-consuming, but the most ! 1995: reliable and is therefore also used for RANK(A). ! 1996: .pp ! 1997: The uniformly distributed random numbers in RAND are ! 1998: obtained from the machine-independent random number generator ! 1999: URAND described in [9]. It is possible to switch to normally ! 2000: distributed random numbers, which are obtained using a transformation ! 2001: also described in [9]. ! 2002: .pp ! 2003: The computation of ! 2004: .sp ! 2005: 2 2 ! 2006: sqrt(a + b ) ! 2007: .sp ! 2008: is required in many matrix algorithms, particularly those involving ! 2009: complex arithmetic. ! 2010: A new approach to carrying out this operation is described by ! 2011: Moler and Morrison [10]. ! 2012: It is a cubically convergent algorithm which starts with a and b , ! 2013: rather than with their squares, and thereby avoids destructive arithmetic ! 2014: underflows and overflows. ! 2015: In MATLAB, the algorithm is used for complex modulus, Euclidean vector ! 2016: norm, plane rotations, and the shift calculation in the eigenvalue ! 2017: and singular value iterations. ! 2018: ! 2019: ! 2020: 12. FLOP and CHOP ! 2021: .pp ! 2022: Detailed information about the amount of work involved ! 2023: in matrix calculations and the resulting accuracy ! 2024: is provided by FLOP and CHOP. ! 2025: The basic unit of work is the "flop", or floating point operation. ! 2026: Roughly, one flop is one execution of a Fortran statement like ! 2027: ! 2028: S = S + X(I)*Y(I) ! 2029: ! 2030: or ! 2031: ! 2032: Y(I) = Y(I) + T*X(I) ! 2033: ! 2034: In other words, it consists of one floating point multiplication, ! 2035: together with one floating point addition and the associated indexing and ! 2036: storage reference operations. ! 2037: .pp ! 2038: MATLAB will print the number ! 2039: of flops required for a particular statement when ! 2040: the statement is terminated by an extra comma. For example, the line ! 2041: ! 2042: n = 20; RAND(n)*RAND(n);, ! 2043: ! 2044: ends with an extra comma. Two 20 by 20 random matrices are generated ! 2045: and multiplied together. The result is assigned to ANS, but ! 2046: the semicolon suppresses its printing. ! 2047: The only output is ! 2048: ! 2049: 8800 flops ! 2050: ! 2051: This is n**3 + 2*n**2 flops, n**2 for each random matrix and ! 2052: n**3 for the product. ! 2053: .pp ! 2054: FLOP is a predefined vector with two components. ! 2055: FLOP(1) is the number of flops used by the most recently executed ! 2056: statement, except that statements with zero flops are ignored. ! 2057: For example, after executing the previous statement, ! 2058: ! 2059: flop(1)/n**3 ! 2060: ! 2061: results in ! 2062: ! 2063: ANS = ! 2064: ! 2065: 1.1000 ! 2066: ! 2067: .pp ! 2068: FLOP(2) is the cumulative total of all the flops used since the ! 2069: beginning of the MATLAB session. ! 2070: The statement ! 2071: ! 2072: FLOP = <0 0> ! 2073: ! 2074: resets the total. ! 2075: .pp ! 2076: There are several difficulties associated with keeping a precise ! 2077: count of floating point operations. ! 2078: An addition or subtraction that is not paired with a multiplication ! 2079: is usually counted as a flop. ! 2080: The same is true of an isolated multiplication that is not paired ! 2081: with an addition. ! 2082: Each floating point division counts as a flop. ! 2083: But the number of operations required by system dependent library ! 2084: functions such as square root cannot be counted, so most ! 2085: elementary functions are arbitrarily counted as using only one flop. ! 2086: .pp ! 2087: The biggest difficulty occurs with complex arithmetic. ! 2088: Almost all operations on the real parts of matrices are counted. ! 2089: However, the operations on the complex parts of matrices are counted ! 2090: only when they involve nonzero elements. ! 2091: This means that simple operations on nonreal matrices require ! 2092: only about twice as many flops as the same operations on real matrices. ! 2093: This factor of two is not necessarily an accurate measure of the ! 2094: relative costs of real and complex arithmetic. ! 2095: .pp ! 2096: The result of each floating point operation may also be "chopped" ! 2097: to simulate a computer with a shorter word length. ! 2098: The details of this chopping operation depend upon the format ! 2099: of the floating point word. ! 2100: Usually, the fraction in the floating point word can be ! 2101: regarded as consisting of several octal or hexadecimal digits. ! 2102: The least significant of these digits can be set to zero ! 2103: by a logical masking operation. Thus the statement ! 2104: ! 2105: CHOP(p) ! 2106: ! 2107: causes the p least significant octal or hexadecimal ! 2108: digits in the result of each floating point operation to be set to zero. ! 2109: For example, if the computer being used has an IBM 360 long ! 2110: floating point word with 14 hexadecimal digits in the ! 2111: fraction, then CHOP(8) results in simulation of a computer with only ! 2112: 6 hexadecimal digits in the fraction, i.e. a short floating point word. ! 2113: On a computer such as the CDC 6600 with 16 octal digits, ! 2114: CHOP(8) results in about the same accuracy because the remaining ! 2115: 8 octal digits represent the same number of bits as 6 hexadecimal digits. ! 2116: .pp ! 2117: Some idea of the effect of CHOP on any particular system can be ! 2118: obtained by executing the following statements. ! 2119: ! 2120: long, t = 1/10 ! 2121: long z, t = 1/10 ! 2122: chop(8) ! 2123: long, t = 1/10 ! 2124: long z, t = 1/10 ! 2125: ! 2126: .pp ! 2127: The following Fortran subprograms illustrate more details of ! 2128: FLOP and CHOP. ! 2129: The first subprogram is a simplified example of a system-dependent ! 2130: function used within MATLAB itself. ! 2131: The common variable FLP is essentially the first component ! 2132: of the variable FLOP. ! 2133: The common variable CHP is initially zero, but it is set to p by ! 2134: the statement CHOP(p). ! 2135: To shorten the DATA statement, we assume there are only 6 hexadecimal ! 2136: digits. ! 2137: We also assume an extension of Fortran that allows .AND. to be ! 2138: used as a binary operation between two real variables. ! 2139: ! 2140: .nf ! 2141: .ne 12 ! 2142: REAL FUNCTION FLOP(X) ! 2143: REAL X ! 2144: INTEGER FLP,CHP ! 2145: COMMON FLP,CHP ! 2146: REAL MASK(5) ! 2147: DATA MASK/ZFFFFFFF0,ZFFFFFF00,ZFFFFF000,ZFFFF0000,ZFFF00000/ ! 2148: FLP = FLP + 1 ! 2149: IF (CHP .EQ. 0) FLOP = X ! 2150: IF (CHP .GE. 1 .AND. CHP .LE. 5) FLOP = X .AND. MASK(CHP) ! 2151: IF (CHP .GE. 6) FLOP = 0.0 ! 2152: RETURN ! 2153: END ! 2154: ! 2155: .fi ! 2156: .pp ! 2157: The following subroutine illustrates a typical use of the previous ! 2158: function within MATLAB. It is a simplified version of ! 2159: the Basic Linear Algebra Subprogram that adds a scalar multiple ! 2160: of one vector to another. We assume here that the vectors are ! 2161: stored with a memory increment of one. ! 2162: ! 2163: .nf ! 2164: .ne 11 ! 2165: SUBROUTINE SAXPY(N,TR,TI,XR,XI,YR,YI) ! 2166: REAL TR,TI,XR(N),XI(N),YR(N),YI(N),FLOP ! 2167: IF (N .LE. 0) RETURN ! 2168: IF (TR .EQ. 0.0 .AND. TI .EQ. 0.0) RETURN ! 2169: DO 10 I = 1, N ! 2170: YR(I) = FLOP(YR(I) + TR*XR(I) - TI*XI(I)) ! 2171: YI(I) = YI(I) + TR*XI(I) + TI*XR(I) ! 2172: IF (YI(I) .NE. 0.0D0) YI(I) = FLOP(YI(I)) ! 2173: 10 CONTINUE ! 2174: RETURN ! 2175: END ! 2176: .fi ! 2177: ! 2178: .pp ! 2179: The saxpy operation is perhaps the most fundamental operation within ! 2180: LINPACK. ! 2181: It is used in the computation of the LU, the QR and the SVD factorizations, ! 2182: and in several other places. ! 2183: We see that adding a multiple of one vector with n components to another ! 2184: uses n flops if the vectors are real and between n and 2*n flops if the ! 2185: vectors have nonzero imaginary components. ! 2186: .pp ! 2187: The permanent MATLAB variable EPS is reset by the statement CHOP(p). ! 2188: Its new value is usually the smallest inverse power of two that ! 2189: satisfies the Fortran logical test ! 2190: ! 2191: FLOP(1.0+EPS) .GT. 1.0 ! 2192: ! 2193: However, if EPS had been directly reset to a larger value, the old ! 2194: value is retained. ! 2195: ! 2196: ! 2197: ! 2198: 13. Communicating with other programs ! 2199: .pp ! 2200: There are four different ways MATLAB can be used in conjunction ! 2201: with other programs: ! 2202: -- USER, ! 2203: -- EXEC, ! 2204: -- SAVE and LOAD, ! 2205: -- MATZ, CALL and RETURN . ! 2206: .pp ! 2207: Let us illustrate each of these by the following simple example. ! 2208: ! 2209: n = 6 ! 2210: for i = 1:n, for j = 1:n, a(i,j) = abs(i-j); ! 2211: A ! 2212: X = inv(A) ! 2213: ! 2214: .pp ! 2215: The example A could be introduced into MATLAB by writing ! 2216: the following Fortran subroutine. ! 2217: ! 2218: SUBROUTINE USER(A,M,N,S,T) ! 2219: DOUBLE PRECISION A(1),S,T ! 2220: N = IDINT(A(1)) ! 2221: M = N ! 2222: DO 10 J = 1, N ! 2223: DO 10 I = 1, N ! 2224: K = I + (J-1)*M ! 2225: A(K) = IABS(I-J) ! 2226: 10 CONTINUE ! 2227: RETURN ! 2228: END ! 2229: ! 2230: This subroutine should be compiled and linked into MATLAB in place ! 2231: of the original version of USER. Then the MATLAB statements ! 2232: ! 2233: n = 6 ! 2234: A = user(n) ! 2235: X = inv(A) ! 2236: ! 2237: do the job. ! 2238: .pp ! 2239: The example A could be generated by storing the following ! 2240: text in a file named, say, EXAMPLE . ! 2241: ! 2242: for i = 1:n, for j = 1:n, a(i,j) = abs(i-j); ! 2243: ! 2244: Then the MATLAB statements ! 2245: ! 2246: n = 6 ! 2247: exec('EXAMPLE',0) ! 2248: X = inv(A) ! 2249: ! 2250: have the desired effect. The 0 as the optional second parameter ! 2251: of exec indicates that the text in the file should not be printed ! 2252: on the terminal. ! 2253: .pp ! 2254: The matrices A and X could also be stored in files. ! 2255: Two separate main programs would be involved. ! 2256: The first is: ! 2257: ! 2258: PROGRAM MAINA ! 2259: DOUBLE PRECISION A(10,10) ! 2260: N = 6 ! 2261: DO 10 J = 1, N ! 2262: DO 10 I = 1, N ! 2263: A(I,J) = IABS(I-J) ! 2264: 10 CONTINUE ! 2265: OPEN(UNIT=1,FILE='A') ! 2266: WRITE(1,101) N,N ! 2267: 101 FORMAT('A ',2I4) ! 2268: DO 20 J = 1, N ! 2269: WRITE(1,102) (A(I,J),I=1,N) ! 2270: 20 CONTINUE ! 2271: 102 FORMAT(4Z18) ! 2272: END ! 2273: ! 2274: The OPEN statement may take different forms on different systems. ! 2275: It attaches Fortran logical unit number 1 to the file named A. ! 2276: The FORMAT number 102 may also be system dependent. This particular ! 2277: one is appropriate for hexadecimal computers with an 8 byte ! 2278: double precision floating point word. Check, or modify, MATLAB ! 2279: subroutine SAVLOD. ! 2280: .pp ! 2281: After this program is executed, enter MATLAB and give the following ! 2282: statements: ! 2283: ! 2284: load('A') ! 2285: X = inv(A) ! 2286: save('X',X) ! 2287: ! 2288: If all goes according to plan, this will read the matrix A from the ! 2289: file A, invert it, store the inverse in X and then write the matrix ! 2290: X on the file X . The following program can then access X . ! 2291: ! 2292: PROGRAM MAINX ! 2293: DOUBLE PRECISION X(10,10) ! 2294: OPEN(UNIT=1,FILE='X') ! 2295: REWIND 1 ! 2296: READ (1,101) ID,M,N ! 2297: 101 FORMAT(A4,2I4) ! 2298: DO 10 J = 1, N ! 2299: READ(1,102) (X(I,J),I=1,M) ! 2300: 10 CONTINUE ! 2301: 102 FORMAT(4Z18) ! 2302: ... ! 2303: ... ! 2304: ! 2305: .pp ! 2306: The most elaborate mechanism involves using MATLAB as a ! 2307: subroutine within another program. Communication with the ! 2308: MATLAB stack is accomplished using subroutine MATZ which is ! 2309: distributed with MATLAB, but which is not used by MATLAB itself. ! 2310: The preample of MATZ is: ! 2311: ! 2312: .nf ! 2313: SUBROUTINE MATZ(A,LDA,M,N,IDA,JOB,IERR) ! 2314: INTEGER LDA,M,N,IDA(1),JOB,IERR ! 2315: DOUBLE PRECISION A(LDA,N) ! 2316: C ! 2317: C ACCESS MATLAB VARIABLE STACK ! 2318: C A IS AN M BY N MATRIX, STORED IN AN ARRAY WITH ! 2319: C LEADING DIMENSION LDA. ! 2320: C IDA IS THE NAME OF A. ! 2321: C IF IDA IS AN INTEGER K LESS THAN 10, THEN THE NAME IS 'A'K ! 2322: C OTHERWISE, IDA(1:4) IS FOUR CHARACTERS, FORMAT 4A1. ! 2323: C JOB = 0 GET REAL A FROM MATLAB, ! 2324: C = 1 PUT REAL A INTO MATLAB, ! 2325: C = 10 GET IMAG PART OF A FROM MATLAB, ! 2326: C = 11 PUT IMAG PART OF A INTO MATLAB. ! 2327: C RETURN WITH NONZERO IERR AFTER MATLAB ERROR MESSAGE. ! 2328: C ! 2329: C USES MATLAB ROUTINES STACKG, STACKP AND ERROR ! 2330: .fi ! 2331: ! 2332: .pp ! 2333: The preample of subroutine MATLAB is: ! 2334: ! 2335: ! 2336: .nf ! 2337: SUBROUTINE MATLAB(INIT) ! 2338: C INIT = 0 FOR FIRST ENTRY, NONZERO FOR SUBSEQUENT ENTRIES ! 2339: .fi ! 2340: ! 2341: .pp ! 2342: To do our example, write the following program: ! 2343: ! 2344: DOUBLE PRECISION A(10,10),X(10,10) ! 2345: INTEGER IDA(4),IDX(4) ! 2346: DATA LDA/10/ ! 2347: DATA IDA/'A',' ',' ',' '/, IDX/'X',' ',' ',' '/ ! 2348: CALL MATLAB(0) ! 2349: N = 6 ! 2350: DO 10 J = 1, N ! 2351: DO 10 I = 1, N ! 2352: A(I,J) = IABS(I-J) ! 2353: 10 CONTINUE ! 2354: CALL MATZ(A,LDA,N,N,IDA,1,IERR) ! 2355: IF (IERR .NE. 0) GO TO ... ! 2356: CALL MATLAB(1) ! 2357: CALL MATZ(X,LDA,N,N,IDX,0,IERR) ! 2358: IF (IERR .NE. 0) GO TO ... ! 2359: ... ! 2360: ... ! 2361: ! 2362: When this program is executed, the call to MATLAB(0) ! 2363: produces the MATLAB greeting, then waits for input. ! 2364: The command ! 2365: ! 2366: return ! 2367: ! 2368: sends control back to our example program. The matrix A ! 2369: is generated by the program and sent to the stack by the ! 2370: first call to MATZ. The call to MATLAB(1) produces the ! 2371: MATLAB prompt. Then the statements ! 2372: ! 2373: X = inv(A) ! 2374: return ! 2375: ! 2376: will invert our matrix, put the result on the stack and go ! 2377: back to our program. The second call to MATZ will retrieve X . ! 2378: .pp ! 2379: By the way, this matrix X is interesting. ! 2380: Take a look at round(2*(n-1)*X). ! 2381: ! 2382: ! 2383: ! 2384: ! 2385: Acknowledgement. ! 2386: ! 2387: .pp ! 2388: Most of the ! 2389: work on MATLAB has been carried out at the University of ! 2390: New Mexico, where it is being supported by the National Science Foundation. ! 2391: Additional work has been done ! 2392: during visits to Stanford Linear Accelerator Center, ! 2393: Argonne National Laboratory and Los Alamos Scientific Laboratory, ! 2394: where support has been provided by NSF and ! 2395: the Department of Energy. ! 2396: ! 2397: ! 2398: References ! 2399: .in +5 ! 2400: ! 2401: .ti -5 ! 2402: [1] J. J. Dongarra, J. R. Bunch, C. B. Moler and G. W. Stewart, ! 2403: LINPACK Users' Guide, Society for Industrial and Applied ! 2404: Mathematics, Philadelphia, 1979. ! 2405: ! 2406: .ti -5 ! 2407: [2] B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow, ! 2408: Y. Ikebe, V. C. Klema, C. B. Moler, Matrix Eigensystem ! 2409: Routines -- EISPACK Guide, Lecture Notes in Computer ! 2410: Science, volume 6, second edition, Springer-Verlag, 1976. ! 2411: ! 2412: .ti -5 ! 2413: [3] B. S. Garbow, J. M. Boyle, J. J. Dongarra, C. B. Moler, ! 2414: Matrix Eigensystem Routines -- EISPACK Guide Extension, ! 2415: Lecture Notes in Computer Science, volume 51, Springer-Verlag, ! 2416: 1977. ! 2417: ! 2418: .ti -5 ! 2419: [4] S. Cohen and S. Pieper, SPEAKEASY III Reference Manual, ! 2420: Speakeasy Computing Corp., Chicago, Ill., 1979. ! 2421: ! 2422: .ti -5 ! 2423: [5] J. H. Wilkinson and C. Reinsch, Handbook for Automatic ! 2424: Computation, volume II, Linear Algebra, Springer-Verlag, 1971. ! 2425: ! 2426: .ti -5 ! 2427: [6] Niklaus Wirth, Algorithms + Data Structures = Programs, ! 2428: Prentice-Hall, 1976. ! 2429: ! 2430: .ti -5 ! 2431: [7] H. B. Keller and D. Sachs, ! 2432: "Calculations of the Conductivity of a Medium Containing ! 2433: Cylindrical Inclusions", ! 2434: J. Applied Physics 35, 537-538, 1964. ! 2435: ! 2436: .ti -5 ! 2437: [8] C. B. Moler and C. F. Van Loan, Nineteen Dubious Ways to Compute ! 2438: the Exponential of a Matrix, SIAM Review 20, 801-836, 1979. ! 2439: ! 2440: .ti -5 ! 2441: [9] G. E. Forsythe, M. A. Malcolm and C. B. Moler, Computer Methods ! 2442: for Mathematical Computations, Prentice-Hall, 1977. ! 2443: ! 2444: .ti -5 ! 2445: [10] C. B. Moler and D. R. Morrison, "Replacing square roots by ! 2446: Pythagorean sums", University of New Mexico, ! 2447: Computer Science Department, ! 2448: technical report, submitted for publication, 1980. ! 2449: .in -5 ! 2450: .bp ! 2451: Appendix. The HELP document ! 2452: .in +6 ! 2453: .xx ! 2454: NEWS``MATLAB NEWS dated May, 1981. ! 2455: .br ! 2456: This describes recent or local changes. ! 2457: .br ! 2458: The new features added since the November, 1980, printing of the ! 2459: Users' Guide include DIARY, EDIT, KRON, MACRO, PLOT, RAT, TRIL, TRIU and ! 2460: six element-by-element operations: ! 2461: .br ! 2462: .*```./```.\\```.*.```./.```.\\. ! 2463: .br ! 2464: Some additional capabilities have been added to ! 2465: EXIT, RANDOM, RCOND, SIZE and SVD. ! 2466: .xx ! 2467: INTRO`Welcome to MATLAB. ! 2468: ! 2469: Here are a few sample statements: ! 2470: ! 2471: .nf ! 2472: A = <1 2; 3 4> ! 2473: b = <5 6>' ! 2474: x = A\\b ! 2475: <V,D> = eig(A), norm(A-V*D/V) ! 2476: help \\ , help eig ! 2477: exec('demo',7) ! 2478: .fi ! 2479: ! 2480: For more information, see the MATLAB Users' Guide ! 2481: which is contained in file ... ! 2482: or may be obtained from ... . ! 2483: .xx ! 2484: HELP``HELP gives assistance. ! 2485: .br ! 2486: HELP HELP obviously prints this message. ! 2487: .br ! 2488: To see all the HELP messages, list the file ... . ! 2489: .xx ! 2490: <`````< > Brackets used in forming vectors and matrices. ! 2491: .br ! 2492: <6.9 9.64 SQRT(-1)> is a vector with three elements ! 2493: separated by blanks. <6.9, 9.64, sqrt(-1)> is the same thing. ! 2494: <1+I 2-I 3> and <1 +I 2 -I 3> are not the same. The first ! 2495: has three elements, the second has five. ! 2496: .br ! 2497: <11 12 13; 21 22 23> is a 2 by 3 matrix . The semicolon ! 2498: ends the first row. ! 2499: .sp ! 2500: Vectors and matrices can be used inside < > brackets. ! 2501: .br ! 2502: <A B; C> is allowed if the number of rows of A equals ! 2503: the number of rows of B and the number of columns of A plus ! 2504: the number of columns of B equals the number of columns of C . ! 2505: This rule generalizes in a hopefully obvious way to allow ! 2506: fairly complicated constructions. ! 2507: .sp ! 2508: A = < > stores an empty matrix in A , thereby removing it ! 2509: from the list of current variables. ! 2510: .sp ! 2511: For the use of < and > on the left of the = in multiple ! 2512: assignment statements, see LU, EIG, SVD and so on. ! 2513: .sp ! 2514: In WHILE and IF clauses, <> means less than or greater than, ! 2515: i.e. not equal, < means less than, > means greater than, ! 2516: <= means less than or equal, >= means greater than or equal. ! 2517: .sp ! 2518: For the use of > and < to delineate macros, see MACRO. ! 2519: .xx ! 2520: >`````See < . Also see MACRO. ! 2521: .xx ! 2522: (`````( ) Used to indicate precedence in arithmetic expressions in ! 2523: the usual way. ! 2524: Used to enclose arguments of functions in the usual way. ! 2525: Used to enclose subscripts of vectors and matrices in a manner ! 2526: somewhat more general than the usual way. If X and V ! 2527: are vectors, then X(V) is <X(V(1)), X(V(2)), ..., X(V(N))> . ! 2528: The components of V are rounded to nearest integers and ! 2529: used as subscripts. An error occurs if any such subscript ! 2530: is less than 1 or greater than the dimension of X . ! 2531: Some examples: ! 2532: .br ! 2533: X(3) is the third element of X . ! 2534: .br ! 2535: X(<1 2 3>) is the first three elements of X . So is ! 2536: .br ! 2537: X(<SQRT(2), SQRT(3), 4*ATAN(1)>) . ! 2538: .br ! 2539: If X has N components, X(N:-1:1) reverses them. ! 2540: .br ! 2541: The same indirect subscripting is used in matrices. ! 2542: If V has M components and W has N components, then ! 2543: A(V,W) is the M by N matrix formed from the elements of A ! 2544: whose subscripts are the elements of V and W . ! 2545: For example... ! 2546: .br ! 2547: A(<1,5>,:) = A(<5,1>,:) interchanges rows 1 and 5 of A . ! 2548: .xx ! 2549: )`````See ( . ! 2550: .xx ! 2551: =`````Used in assignment statements and to mean equality in ! 2552: WHILE and IF clauses. ! 2553: .xx ! 2554: .cc - ! 2555: .`````Decimal point. 314/100, 3.14 and .314E1 are all the same. ! 2556: -cc ! 2557: .sp ! 2558: Element-by-element multiplicative operations are obtained ! 2559: using .* , ./ , or .\\ . For example, C = A ./ B is the ! 2560: matrix with elements c(i,j) = a(i,j)/b(i,j) . ! 2561: .sp ! 2562: Kronecker tensor products and quotients are obtained with .*. , ./. ! 2563: and .\\. . See KRON. ! 2564: .sp ! 2565: Two or more points at the end of the line indicate continuation. ! 2566: The total line length limit is 1024 characters. ! 2567: .xx ! 2568: ,`````Used to separate matrix subscripts and function arguments. ! 2569: Used at the end of FOR, WHILE and IF clauses. ! 2570: Used to separate statements in multi-statement lines. In this ! 2571: situation, it may be replaced by semicolon to suppress printing. ! 2572: .xx ! 2573: ;`````Used inside brackets to end rows. ! 2574: .br ! 2575: Used after an expression or statement to suppress printing. ! 2576: .br ! 2577: See SEMI. ! 2578: .xx ! 2579: \\`````Backslash or matrix left division. ! 2580: A\\B is roughly the ! 2581: same as INV(A)*B , except it is computed in a different way. ! 2582: If A is an N by N matrix and B is a column vector with N ! 2583: components, or a matrix with several such columns, then ! 2584: X = A\\B is the solution to the equation A*X = B computed ! 2585: by Gaussian elimination. A warning message is printed if A ! 2586: is badly scaled or nearly singular. ! 2587: .br ! 2588: A\\EYE produces the inverse of A . ! 2589: .sp ! 2590: If A is an M by N matrix with M < or > N and B is ! 2591: a column vector with M components, or a matrix with several ! 2592: such columns, then X = A\\B is the solution in the least ! 2593: squares sense to the under- or overdetermined system of equations ! 2594: A*X = B . ! 2595: The effective rank, K, of A is determined from the QR ! 2596: decomposition with pivoting. A solution X is computed which ! 2597: has at most K nonzero components per column. If K < N ! 2598: this will usually not be the same solution as PINV(A)*B . ! 2599: .br ! 2600: A\\EYE produces a generalized inverse of A . ! 2601: .sp ! 2602: If A and B have the same dimensions, then A .\\ B has elements ! 2603: a(i,j)\\b(i,j) . ! 2604: .sp ! 2605: Also, see EDIT. ! 2606: .xx ! 2607: /`````Slash or matrix right division. B/A is roughly the same ! 2608: as B*INV(A) . More precisely, B/A = (A'\\B')' . See \\ . ! 2609: .sp ! 2610: IF A and B have the same dimensions, then A ./ B has elements ! 2611: a(i,j)/b(i,j) . ! 2612: .sp ! 2613: Two or more slashes together on a line indicate a logical ! 2614: end of line. Any following text is ignored. ! 2615: .xx ! 2616: .c2 - ! 2617: '`````Transpose. X' is the complex conjugate transpose of X . ! 2618: .c2 ! 2619: Quote. 'ANY TEXT' is a vector whose components are the ! 2620: MATLAB internal codes for the characters. ! 2621: '0' = 0, '1' = 1, 'A' = 10, 'Z' = 35, ';' = 53, and so on. ! 2622: A quote within the text is indicated by two quotes. ! 2623: See DISP and FILE . ! 2624: .xx ! 2625: +`````Addition. X + Y . X and Y must have the same dimensions. ! 2626: .xx ! 2627: -`````Subtraction. X - Y . X and Y must have the same dimensions. ! 2628: .xx ! 2629: *`````Matrix multiplication, X*Y . Any scalar (1 by 1 matrix) may ! 2630: multiply anything. Otherwise, the number of columns of X ! 2631: must equal the number of rows of Y . ! 2632: .sp ! 2633: Element-by-element multiplication is obtained with X .* Y . ! 2634: .sp ! 2635: The Kronecker tensor product is denoted by X .*. Y . ! 2636: .sp ! 2637: Powers. X**p is X to the p power. p must be a scalar. ! 2638: If X is a matrix, see FUN . ! 2639: .xx ! 2640: :`````Colon. Used in subscripts, FOR iterations and possibly elsewhere. ! 2641: .br ! 2642: J:K is the same as <J, J+1, ..., K> ! 2643: .br ! 2644: J:K is empty if J > K . ! 2645: .br ! 2646: J:I:K is the same as <J, J+I, J+2I, ..., K> ! 2647: .br ! 2648: J:I:K is empty if I > 0 and J > K or if I < 0 and J < K . ! 2649: .br ! 2650: The colon notation can be used to pick out selected rows, ! 2651: columns and elements of vectors and matrices. ! 2652: .br ! 2653: A(:) is all the elements of A, regarded as a single column. ! 2654: .br ! 2655: A(:,J) is the J-th column of A ! 2656: .br ! 2657: A(J:K) is A(J),A(J+1),...,A(K) ! 2658: .br ! 2659: A(:,J:K) is A(:,J),A(:,J+1),...,A(:,K) and so on. ! 2660: .br ! 2661: For the use of the colon in the FOR statement, See FOR . ! 2662: .xx ! 2663: ABS```ABS(X) is the absolute value, or complex modulus, of ! 2664: the elements of X . ! 2665: .xx ! 2666: ANS```Variable created automatically when expressions are not ! 2667: assigned to anything else. ! 2668: .xx ! 2669: ATAN``ATAN(X) is the arctangent of X . See FUN . ! 2670: .xx ! 2671: BASE``BASE(X,B) is a vector containing the base B representation ! 2672: of X . This is often used in conjunction ! 2673: with DISPLAY. ! 2674: DISPLAY(X,B) is the same as DISPLAY(BASE(X,B)). ! 2675: For example, DISP(4*ATAN(1),16) prints the ! 2676: hexadecimal representation of pi. ! 2677: .xx ! 2678: CHAR``CHAR(K) requests an input line containing a single character ! 2679: to replace MATLAB character number K in the following table. ! 2680: For example, CHAR(45) replaces backslash. ! 2681: CHAR(-K) replaces the alternate character number K. ! 2682: .sp ! 2683: K character alternate name ! 2684: 0 - 9 0 - 9 0 - 9 digits ! 2685: 10 - 35 A - Z a - z letters ! 2686: 36 blank ! 2687: 37 ( ( lparen ! 2688: 38 ) ) rparen ! 2689: 39 ; ; semi ! 2690: 40 : | colon ! 2691: 41 + + plus ! 2692: 42 - - minus ! 2693: 43 * * star ! 2694: 44 / / slash ! 2695: 45 \\ $ backslash ! 2696: 46 = = equal ! 2697: 47 . . dot ! 2698: 48 , , comma ! 2699: 49 ' " quote ! 2700: 50 < [ less ! 2701: 51 > ] great ! 2702: .xx ! 2703: CHOL``Cholesky factorization. CHOL(X) uses only the diagonal and ! 2704: upper triangle of X . The lower triangular is assumed to be ! 2705: the (complex conjugate) transpose of the upper. ! 2706: If X is positive definite, then R = CHOL(X) produces an ! 2707: upper triangular R so that R'*R = X . ! 2708: If X is not positive definite, an error message is printed. ! 2709: .xx ! 2710: CHOP``Truncate arithmetic. CHOP(P) causes P places to be chopped ! 2711: off after each arithmetic operation in subsequent computations. ! 2712: This means P hexadecimal digits on some computers ! 2713: and P octal digits on others. ! 2714: CHOP(0) restores full precision. ! 2715: .xx ! 2716: CLEAR`Erases all variables, except EPS, FLOP, EYE and RAND. ! 2717: .br ! 2718: X = <> erases only variable X . ! 2719: So does CLEAR X . ! 2720: .xx ! 2721: COND``Condition number in 2-norm. COND(X) is the ratio of the ! 2722: largest singular value of X to the smallest. ! 2723: .xx ! 2724: CONJG`CONJG(X) is the complex conjugate of X . ! 2725: .xx ! 2726: COS```COS(X) is the cosine of X . See FUN . ! 2727: .xx ! 2728: DET```DET(X) is the determinant of the square matrix X . ! 2729: .xx ! 2730: DIAG``If V is a row or column vector with N components, ! 2731: DIAG(V,K) is a square matrix of order N+ABS(K) with ! 2732: the elements of V on the K-th diagonal. ! 2733: K = 0 is the main diagonal, K > 0 is above the main diagonal ! 2734: and K < 0 is below the main diagonal. ! 2735: DIAG(V) simply puts V on the main diagonal. ! 2736: .br ! 2737: eg. DIAG(-M:M) + DIAG(ONES(2*M,1),1) + DIAG(ONES(2*M,1),-1) ! 2738: .br ! 2739: produces a tridiagonal matrix of order 2*M+1 . ! 2740: .br ! 2741: IF X is a matrix, DIAG(X,K) is a column vector formed ! 2742: from the elements of the K-th diagonal of X . ! 2743: .br ! 2744: DIAG(X) is the main diagonal of X . ! 2745: .br ! 2746: DIAG(DIAG(X)) is a diagonal matrix . ! 2747: .xx ! 2748: DIARY`DIARY('file') causes a copy of all subsequent terminal input ! 2749: and most of the resulting output to be written on the file. ! 2750: DIARY(0) turns it off. See FILE. ! 2751: .xx ! 2752: DISP``DISPLAY(X) prints X in a compact format. If all the elements ! 2753: of X are integers between 0 and 51, then X is interpreted as MATLAB text ! 2754: and printed accordingly. Otherwise, + , - and blank are printed for ! 2755: positive, negative and zero elements. Imaginary parts are ignored. ! 2756: .br ! 2757: DISP(X,B) is the same as DISP(BASE(X,B)). ! 2758: .xx ! 2759: EDIT``There are no editing features available on most installations and ! 2760: EDIT is not a command. However, on a few systems a command line ! 2761: consisting of a single backslash \\ will cause the local file editor ! 2762: to be called with a copy of the previous input line. When the ! 2763: editor returns control to MATLAB, it will execute the line again. ! 2764: .xx ! 2765: EIG```Eigenvalues and eigenvectors. ! 2766: .br ! 2767: EIG(X) is a vector containing ! 2768: the eigenvalues of a square matrix X . ! 2769: .br ! 2770: <V,D> = EIG(X) produces a diagonal matrix D of eigenvalues ! 2771: and a full matrix V whose columns are the corresponding ! 2772: eigenvectors so that X*V = V*D . ! 2773: .xx ! 2774: ELSE``Used with IF . ! 2775: .xx ! 2776: END```Terminates the scope of FOR, WHILE and IF statements. ! 2777: Without END's, FOR and WHILE repeat all statements up to the ! 2778: end of the line. Each END is paired with the closest previous ! 2779: unpaired FOR or WHILE and serves to terminate its scope. ! 2780: The line ! 2781: .br ! 2782: FOR I=1:N, FOR J=1:N, A(I,J)=1/(I+J-1); A ! 2783: .br ! 2784: would cause A to be printed N**2 times, once for each ! 2785: new element. On the other hand, the line ! 2786: .br ! 2787: FOR I=1:N, FOR J=1:N, A(I,J)=1/(I+J-1); END, END, A ! 2788: .br ! 2789: will lead to only the final printing of A . ! 2790: .br ! 2791: Similar considerations apply to WHILE. ! 2792: .br ! 2793: EXIT terminates execution of loops or of MATLAB itself. ! 2794: .xx ! 2795: EPS```Floating point relative accuracy. A permanent variable whose ! 2796: value is initially the distance from 1.0 to the next largest floating ! 2797: point number. The value is changed by CHOP, and other values ! 2798: may be assigned. EPS is used as a default tolerance by PINV and RANK. ! 2799: .xx ! 2800: EXEC``EXEC('file',k) obtains subsequent MATLAB input from an external ! 2801: file. The printing of input is controlled by the optional parameter k . ! 2802: .br ! 2803: If k = 1 , the input is echoed. ! 2804: .br ! 2805: If k = 2 , the MATLAB prompt <> is printed. ! 2806: .br ! 2807: If k = 4 , MATLAB pauses before each prompt and waits for a null line ! 2808: to continue. ! 2809: .br ! 2810: If k = 0 , there is no echo, prompt or pause. ! 2811: This is the default if the ! 2812: exec command is followed by a semicolon. ! 2813: .br ! 2814: If k = 7 , there will be echos, prompts and pauses. This is useful ! 2815: for demonstrations on video terminals. ! 2816: .br ! 2817: If k = 3 , there will be echos and prompts, but no pauses. ! 2818: This is the the default if the exec command ! 2819: is not followed by a semicolon. ! 2820: .br ! 2821: EXEC(0) causes subsequent input to be obtained from the terminal. ! 2822: An end-of-file has the same effect. ! 2823: .br ! 2824: EXEC's may be nested, i.e. the text in the file may contain EXEC of ! 2825: another file. EXEC's may also be driven by FOR and WHILE loops. ! 2826: .xx ! 2827: EXIT``Causes termination of a FOR or WHILE loop. ! 2828: .br ! 2829: If not in a loop, terminates execution of MATLAB. ! 2830: .xx ! 2831: EXP```EXP(X) is the exponential of X , e to the X . See FUN . ! 2832: .xx ! 2833: EYE```Identity matrix. ! 2834: EYE(N) is the N by N identity matrix. ! 2835: EYE(M,N) is an M by N matrix with 1's on the diagonal and ! 2836: zeros elsewhere. ! 2837: EYE(A) is the same size as A . ! 2838: EYE with no arguments is an identity matrix of whatever ! 2839: order is appropriate in the context. For example, ! 2840: A + 3*EYE adds 3 to each diagonal element of A . ! 2841: .xx ! 2842: FILE``The EXEC, SAVE, LOAD, PRINT and DIARY functions access files. ! 2843: The 'file' parameter takes different forms for different ! 2844: operating systems. ! 2845: On most systems, 'file' may be a string of up to 32 characters ! 2846: in quotes. For example, ! 2847: SAVE('A') or EXEC('matlab/demo.exec') . ! 2848: The string will be used as the name of a file in the local operating system. ! 2849: .br ! 2850: On all systems, 'file' may be a positive integer k less than 10 which ! 2851: will be used as a FORTRAN logical unit number. ! 2852: Some systems then automatically access a file with a name like ! 2853: FORT.k or FORk.DAT. ! 2854: Other systems require a file with a name like FT0kF001 to be assigned to ! 2855: unit k before MATLAB is executed. ! 2856: Check your local installation for details. ! 2857: .xx ! 2858: FLOPS`Count of floating point operations. ! 2859: .br ! 2860: FLOPS is a permanently defined row vector with two elements. ! 2861: FLOPS(1) is the number of floating point operations counted ! 2862: during the previous statement. FLOPS(2) is a cumulative total. ! 2863: FLOPS can be used in the same way as any other vector. ! 2864: FLOPS(2) = 0 resets the cumulative total. ! 2865: In addition, FLOPS(1) will be printed whenever a statement ! 2866: is terminated by an extra comma. For example, ! 2867: .br ! 2868: X = INV(A);, ! 2869: .br ! 2870: or ! 2871: .br ! 2872: COND(A), (as the last statement on the line). ! 2873: .br ! 2874: HELP FLPS gives more details. ! 2875: .xx ! 2876: FLPS``More detail on FLOPS. ! 2877: .br ! 2878: It is not feasible to count absolutely all floating point ! 2879: operations, but most of the important ones are counted. ! 2880: Each multiply and add in a real vector operation such as ! 2881: a dot product or a 'saxpy' counts one flop. Each multiply ! 2882: and add in a complex vector operation counts two flops. ! 2883: Other additions, subtractions and multiplications count ! 2884: one flop each if the result is real and two flops if it is not. ! 2885: Real divisions count one and complex divisions count two. ! 2886: Elementary functions count one if real and two if complex. ! 2887: Some examples. If A and B are real N by N matrices, then ! 2888: .br ! 2889: A + B counts N**2 flops, ! 2890: .br ! 2891: A*B counts N**3 flops, ! 2892: .br ! 2893: A**100 counts 99*N**3 flops, ! 2894: .br ! 2895: LU(A) counts roughly (1/3)*N**3 flops. ! 2896: .xx ! 2897: FOR```Repeat statements a specific number of times. ! 2898: .br ! 2899: FOR variable = expr, statement, ..., statement, END ! 2900: .br ! 2901: The END at the end of a line may be omitted. ! 2902: The comma before the END may also be omitted. ! 2903: The columns of the expression are stored one at a time in ! 2904: the variable and then the following statements, up to the ! 2905: END, are executed. The expression is often of the form ! 2906: X:Y, in which case its columns are simply scalars. ! 2907: Some examples (assume N has already been assigned a value). ! 2908: .nf ! 2909: FOR I = 1:N, FOR J = 1:N, A(I,J) = 1/(I+J-1); ! 2910: FOR J = 2:N-1, A(J,J) = J; END; A ! 2911: FOR S = 1.0: -0.1: 0.0, ... steps S with increments of -0.1 . ! 2912: FOR E = EYE(N), ... sets E to the unit N-vectors. ! 2913: FOR V = A, ... has the same effect as ! 2914: FOR J = 1:N, V = A(:,J); ... except J is also set here. ! 2915: .fi ! 2916: .xx ! 2917: FUN```For matrix arguments X , the functions SIN, COS, ATAN, SQRT, ! 2918: LOG, EXP and X**p are computed using eigenvalues D and ! 2919: eigenvectors V . If <V,D> = EIG(X) then f(X) = V*f(D)/V . ! 2920: This method may give inaccurate results if V is badly ! 2921: conditioned. Some idea of the accuracy can be obtained by ! 2922: comparing X**1 with X . ! 2923: .br ! 2924: For vector arguments, the function is applied to each component. ! 2925: .xx ! 2926: HESS``Hessenberg form. The Hessenberg form of a matrix is zero ! 2927: below the first subdiagonal. If the matrix is symmetric or ! 2928: Hermitian, the form is tridiagonal. <P,H> = HESS(A) produces ! 2929: a unitary matrix P and a Hessenberg matrix H so that A = P*H*P'. ! 2930: By itself, HESS(A) returns H. ! 2931: .xx ! 2932: HILB``Inverse Hilbert matrix. HILB(N) is the inverse of the ! 2933: N by N matrix with elements 1/(i+j-1), which is a famous ! 2934: example of a badly conditioned matrix. The result is exact ! 2935: for N less than about 15, depending upon the computer. ! 2936: .xx ! 2937: IF````Conditionally execute statements. ! 2938: Simple form... ! 2939: .br ! 2940: IF expression rop expression, statements ! 2941: .br ! 2942: where rop is =, <, >, <=, >=, or <> (not equal) . ! 2943: The statements are executed once if the indicated comparison ! 2944: between the real parts of the first components of the two ! 2945: expressions is true, otherwise the statements are skipped. ! 2946: Example. ! 2947: .br ! 2948: IF ABS(I-J) = 1, A(I,J) = -1; ! 2949: .br ! 2950: More complicated forms use END in the same way it is used ! 2951: with FOR and WHILE and use ELSE as an abbreviation for ! 2952: END, IF expression not rop expression . Example ! 2953: .nf ! 2954: FOR I = 1:N, FOR J = 1:N, ... ! 2955: IF I = J, A(I,J) = 2; ELSE IF ABS(I-J) = 1, A(I,J) = -1; ... ! 2956: ELSE A(I,J) = 0; ! 2957: An easier way to accomplish the same thing is ! 2958: A = 2*EYE(N); ! 2959: FOR I = 1:N-1, A(I,I+1) = -1; A(I+1,I) = -1; ! 2960: .fi ! 2961: .xx ! 2962: IMAG``IMAG(X) is the imaginary part of X . ! 2963: .xx ! 2964: INV```INV(X) is the inverse of the square matrix X . A warning ! 2965: message is printed if X is badly scaled or nearly singular. ! 2966: .xx ! 2967: KRON``KRON(X,Y) is the Kronecker tensor product of X and Y . ! 2968: It is also denoted by X`.*.`Y . ! 2969: The result is a large matrix formed by taking all possible products ! 2970: between the elements of X and those of Y . ! 2971: For example, if X is 2 by 3, then X`.*.`Y is ! 2972: ! 2973: < x(1,1)*Y x(1,2)*Y x(1,3)*Y ! 2974: x(2,1)*Y x(2,2)*Y x(2,3)*Y > ! 2975: ! 2976: .br ! 2977: The five-point discrete Laplacian for an n-by-n grid can be ! 2978: generated by ! 2979: ! 2980: T = diag(ones(n-1,1),1); T = T + T'; I = EYE(T); ! 2981: A = T.*.I + I.*.T - 4*EYE; ! 2982: ! 2983: Just in case they might be useful, MATLAB includes ! 2984: constructions called Kronecker tensor quotients, ! 2985: denoted by X`./.`Y and X`.\\.`Y . ! 2986: They are obtained by replacing the elementwise multiplications ! 2987: in X`.*.`Y with divisions. ! 2988: .xx ! 2989: LINES`An internal count is kept of the number of lines of output ! 2990: since the last input. Whenever this count approaches a limit, ! 2991: the user is asked whether or not to suppress printing until ! 2992: the next input. Initially the limit is 25. LINES(N) resets ! 2993: the limit to N . ! 2994: .xx ! 2995: LOAD``LOAD('file') retrieves all the variables from the file . ! 2996: See FILE and SAVE for more details. ! 2997: To prepare your own file for LOADing, change the READs to WRITEs ! 2998: in the code given under SAVE. ! 2999: .xx ! 3000: LOG```LOG(X) is the natural logarithm of X . See FUN . ! 3001: Complex results are produced if X is not positive, or has ! 3002: nonpositive eigenvalues. ! 3003: .xx ! 3004: LONG``Determine output format. All computations are done in ! 3005: complex arithmetic and double precision if it is available. ! 3006: SHORT and LONG merely switch between different output formats. ! 3007: .br ! 3008: SHORT````Scaled fixed point format with about 5 digits. ! 3009: .br ! 3010: LONG`````Scaled fixed point format with about 15 digits. ! 3011: .br ! 3012: SHORT`E``Floating point format with about 5 digits. ! 3013: .br ! 3014: LONG`E```Floating point format with about 15 digits. ! 3015: .br ! 3016: LONG`Z```System dependent format, often hexadecimal. ! 3017: .xx ! 3018: LU````Factors from Gaussian elimination. <L,U> = LU(X) stores a ! 3019: upper triangular matrix in U and a 'psychologically lower ! 3020: triangular matrix', i.e. a product of lower triangular ! 3021: and permutation matrices, in L , so that X = L*U . ! 3022: By itself, LU(X) returns the output from CGEFA . ! 3023: .xx ! 3024: MACRO`The macro facility involves text and inward pointing angle brackets. ! 3025: If STRING is the source text for any MATLAB expression or statement, then ! 3026: .br ! 3027: t = 'STRING'; ! 3028: .br ! 3029: encodes the text as a vector of integers and stores that vector in `t`. ! 3030: DISP(t) will print the text and ! 3031: .br ! 3032: >t< ! 3033: .br ! 3034: causes the text to be interpreted, either as a statement ! 3035: or as a factor in an expression. ! 3036: For example ! 3037: .br ! 3038: t = '1/(i+j-1)'; ! 3039: disp(t) ! 3040: for i = 1:n, for j = 1:n, a(i,j) = >t<; ! 3041: .br ! 3042: generates the Hilbert matrix of order n. ! 3043: .br ! 3044: Another example showing indexed text, ! 3045: .br ! 3046: S = <'x = 3 ' ! 3047: 'y = 4 ' ! 3048: 'z = sqrt(x*x+y*y)'> ! 3049: for k = 1:3, >S(k,:)< ! 3050: .br ! 3051: It is necessary that the strings making up the "rows" of the "matrix" ! 3052: `S` have the same lengths. ! 3053: .xx ! 3054: MAGIC`Magic square. MAGIC(N) is an N by N matrix constructed ! 3055: from the integers 1 through N**2 with equal row and column sums. ! 3056: .xx ! 3057: NORM``For matrices.. ! 3058: .nf ! 3059: NORM(X) is the largest singular value of X . ! 3060: NORM(X,1) is the 1-norm of X . ! 3061: NORM(X,2) is the same as NORM(X) . ! 3062: NORM(X,'INF') is the infinity norm of X . ! 3063: NORM(X,'FRO') is the F-norm, i.e. SQRT(SUM(DIAG(X'*X))) . ! 3064: For vectors.. ! 3065: NORM(V,P) = (SUM(V(I)**P))**(1/P) . ! 3066: NORM(V) = NORM(V,2) . ! 3067: NORM(V,'INF') = MAX(ABS(V(I))) . ! 3068: .fi ! 3069: .xx ! 3070: ONES``All ones. ! 3071: ONES(N) is an N by N matrix of ones. ! 3072: ONES(M,N) is an M by N matrix of ones . ! 3073: ONES(A) is the same size as A and all ones . ! 3074: .xx ! 3075: ORTH``Orthogonalization. Q = ORTH(X) is a matrix with orthonormal ! 3076: columns, i.e. Q'*Q = EYE, which span the same space as the ! 3077: columns of X . ! 3078: .xx ! 3079: PINV``Pseudoinverse. X = PINV(A) produces a matrix X of ! 3080: the same dimensions as A' so that A*X*A = A , ! 3081: X*A*X = X and AX and XA are Hermitian . The ! 3082: computation is based on SVD(A) and any singular values ! 3083: less than ! 3084: a tolerance are treated as zero. The default tolerance is ! 3085: NORM(SIZE(A),'inf')*NORM(A)*EPS. ! 3086: This tolerance may be overridden with X = PINV(A,tol). ! 3087: See RANK. ! 3088: .xx ! 3089: PLOT``PLOT(X,Y) produces a plot of the elements of Y against those ! 3090: of X . ! 3091: PLOT(Y) is the same as PLOT(1:n,Y) where n is the number of elements ! 3092: in Y . ! 3093: PLOT(X,Y,P) or PLOT(X,Y,p1,...,pk) passes the optional parameter vector P ! 3094: or scalars p1 through pk to the plot routine. ! 3095: The default plot routine is a crude printer-plot. ! 3096: It is hoped that an interface to local graphics equipment can be ! 3097: provided. ! 3098: .br ! 3099: An interesting example is ! 3100: t = 0:50; ! 3101: PLOT( t.*cos(t), t.*sin(t) ) ! 3102: .xx ! 3103: POLY``Characteristic polynomial. ! 3104: .br ! 3105: If A is an N by N matrix, POLY(A) is a column vector with N+1 ! 3106: elements which are the coefficients of the characteristic ! 3107: polynomial, DET(lambda*EYE - A) . ! 3108: .br ! 3109: If V is a vector, POLY(V) is a vector whose elements are the ! 3110: coefficients of the polynomial whose roots are the elements of V . ! 3111: For vectors, ROOTS and POLY are inverse functions of each other, ! 3112: up to ordering, scaling, and roundoff error. ! 3113: .br ! 3114: ROOTS(POLY(1:20)) generates Wilkinson's famous example. ! 3115: .xx ! 3116: PRINT`PRINT('file',X) prints X on the file using the current format ! 3117: determined by SHORT, LONG Z, etc. See FILE. ! 3118: .xx ! 3119: PROD``PROD(X) is the product of all the elements of X . ! 3120: .xx ! 3121: QR````Orthogonal-triangular decomposition. ! 3122: .br ! 3123: <Q,R> = QR(X) produces an upper triangular matrix R ! 3124: of the same dimension as X and a unitary matrix Q so ! 3125: that X = Q*R . ! 3126: .br ! 3127: <Q,R,E> = QR(X) produces a permutation matrix E , an ! 3128: upper triangular R with decreasing diagonal elements ! 3129: and a unitary Q so that X*E = Q*R . ! 3130: .br ! 3131: By itself, QR(X) returns the output of CQRDC . ! 3132: TRIU(QR(X)) is R . ! 3133: .xx ! 3134: RAND``Random numbers and matrices. ! 3135: RAND(N) is an N by N matrix with random entries. ! 3136: RAND(M,N) is an M by N matrix with random entries. ! 3137: RAND(A) is the same size as A . ! 3138: RAND with no arguments is a scalar whose value changes each ! 3139: time it is referenced. ! 3140: .br ! 3141: Ordinarily, random numbers are uniformly distributed in ! 3142: the interval (0.0,1.0) . ! 3143: RAND('NORMAL') switches to a normal distribution with mean 0.0 and ! 3144: variance 1.0 . ! 3145: RAND('UNIFORM') switches back to the uniform distribution. ! 3146: .br ! 3147: RAND('SEED') returns the current value of the seed for the generator. ! 3148: RAND('SEED',n) sets the seed to n . ! 3149: RAND('SEED',0) resets the seed to 0, its value when MATLAB is first entered. ! 3150: .xx ! 3151: RANK``Rank. K = RANK(X) is the number of singular values of X ! 3152: that are larger than NORM(SIZE(X),'inf')*NORM(X)*EPS. ! 3153: .br ! 3154: K = RANK(X,tol) is the number of singular values of X ! 3155: that are larger than tol . ! 3156: .xx ! 3157: RCOND`RCOND(X) is an estimate for the reciprocal of the condition ! 3158: of X in the 1-norm obtained by the LINPACK condition estimator. ! 3159: If X is well conditioned, RCOND(X) is near 1.0 . ! 3160: If X is badly conditioned, RCOND(X) is near 0.0 . ! 3161: .br ! 3162: <R, Z> = RCOND(A) sets R to RCOND(A) and also produces a vector Z ! 3163: so that ! 3164: NORM(A*Z,1) = R*NORM(A,1)*NORM(Z,1) ! 3165: .br ! 3166: So, if RCOND(A) is small, then Z is an approximate null vector. ! 3167: .xx ! 3168: RAT```An experimental function which attempts to remove the roundoff error ! 3169: from results that should be "simple" rational numbers. ! 3170: .br ! 3171: RAT(X) approximates each element of X ! 3172: by a continued fraction of the form ! 3173: .sp ! 3174: a/b = d1 + 1/(d2 + 1/(d3 + ... + 1/dk)) ! 3175: .sp ! 3176: with k`<=`len, integer di and abs(di)`<=`max . ! 3177: The default values of the parameters are ! 3178: len = 5 and max = 100. ! 3179: .br ! 3180: RAT(len,max) changes the default values. ! 3181: Increasing either len or max increases the number of possible fractions. ! 3182: .br ! 3183: <A,B> = RAT(X) produces integer matrices A and B so that ! 3184: .sp ! 3185: A ./ B = RAT(X) ! 3186: .sp ! 3187: Some examples: ! 3188: .sp ! 3189: long ! 3190: T = hilb(6), X = inv(T) ! 3191: <A,B> = rat(X) ! 3192: H = A ./ B, S = inv(H) ! 3193: ! 3194: short e ! 3195: d = 1:8, e = ones(d), A = abs(d'*e - e'*d) ! 3196: X = inv(A) ! 3197: rat(X) ! 3198: display(ans) ! 3199: ! 3200: .xx ! 3201: REAL``REAL(X) is the real part of X . ! 3202: .xx ! 3203: RETURN From the terminal, causes return to the operating system or ! 3204: other program which invoked MATLAB. ! 3205: From inside an EXEC, causes return to the invoking EXEC, or to the ! 3206: terminal. ! 3207: .xx ! 3208: RREF``RREF(A) is the reduced row echelon form of the rectangular matrix. ! 3209: RREF(A,B) is the same as RREF(<A,B>) . ! 3210: .xx ! 3211: ROOTS`Find polynomial roots. ! 3212: ROOTS(C) computes the roots of the polynomial whose coefficients ! 3213: are the elements of the vector C . If C has N+1 components, ! 3214: the polynomial is C(1)*X**N + ... + C(N)*X + C(N+1) . ! 3215: See POLY. ! 3216: .xx ! 3217: ROUND`ROUND(X) rounds the elements of X to the nearest integers. ! 3218: .xx ! 3219: SAVE``SAVE('file') stores all the current variables in a file. ! 3220: .br ! 3221: SAVE('file',X) saves only X . See FILE . ! 3222: .br ! 3223: The variables may be retrieved later by LOAD('file') or by ! 3224: your own program using the following code for each matrix. ! 3225: The lines involving XIMAG may be eliminated if everything is known to ! 3226: be real. ! 3227: .sp ! 3228: attach lunit to 'file' ! 3229: REAL or DOUBLE PRECISION XREAL(MMAX,NMAX) ! 3230: REAL or DOUBLE PRECISION XIMAG(MMAX,NMAX) ! 3231: READ(lunit,101) ID,M,N,IMG ! 3232: DO 10 J = 1, N ! 3233: READ(lunit,102) (XREAL(I,J), I=1,M) ! 3234: IF (IMG .NE. 0) READ(lunit,102) (XIMAG(I,J),I=1,M) ! 3235: 10 CONTINUE ! 3236: .sp ! 3237: The formats used are system dependent. The following are typical. ! 3238: See SUBROUTINE SAVLOD in your local implementation of MATLAB. ! 3239: .sp ! 3240: 101 FORMAT(4A1,3I4) ! 3241: 102 FORMAT(4Z18) ! 3242: 102 FORMAT(4O20) ! 3243: 102 FORMAT(4D25.18) ! 3244: .xx ! 3245: SCHUR`Schur decomposition. <U,T> = SCHUR(X) produces an upper ! 3246: triangular matrix T , with the eigenvalues of X on the ! 3247: diagonal, and a unitary matrix U so that X = U*T*U' and ! 3248: U'*U = EYE . By itself, SCHUR(X) returns T . ! 3249: .xx ! 3250: SHORT`See LONG . ! 3251: .xx ! 3252: SEMI``Semicolons at the end of lines will cause, rather than suppress, ! 3253: printing. A second SEMI restores the initial interpretation. ! 3254: .xx ! 3255: SIN```SIN(X) is the sine of X . See FUN . ! 3256: .xx ! 3257: SIZE``If X is an M by N matrix, then SIZE(X) is <M, N> . ! 3258: .br ! 3259: Can also be used with a multiple assignment, ! 3260: <M, N> = SIZE(X) . ! 3261: .xx ! 3262: SQRT``SQRT(X) is the square root of X . See FUN . ! 3263: Complex results are produced if X is not positive, or has ! 3264: nonpositive eigenvalues. ! 3265: .xx ! 3266: STOP``Use EXIT instead. ! 3267: .xx ! 3268: SUM```SUM(X) is the sum of all the elements of X . ! 3269: SUM(DIAG(X)) is the trace of X . ! 3270: .xx ! 3271: SVD```Singular value decomposition. <U,S,V> = SVD(X) produces ! 3272: a diagonal matrix S , of the same dimension as X and with ! 3273: nonnegative diagonal elements in decreasing order, and ! 3274: unitary matrices U and V so that X = U*S*V' . ! 3275: .br ! 3276: By itself, SVD(X) returns a vector containing the singular values. ! 3277: .br ! 3278: <U,S,V> = SVD(X,0) produces the "economy size" decomposition. ! 3279: If X is m by n with m > n, then only the first n columns of U are ! 3280: computed and S is n by n . ! 3281: .xx ! 3282: TRIL``Lower triangle. TRIL(X) is the lower triangular part of X. ! 3283: TRIL(X,K) is the elements on and below the K-th diagonal of X. ! 3284: K = 0 is the main diagonal, K > 0 is above the main diagonal ! 3285: and K < 0 is below the main diagonal. ! 3286: .xx ! 3287: TRIU``Upper triangle. TRIU(X) is the upper triangular part of X. ! 3288: TRIU(X,K) is the elements on and above the K-th diagonal of X. ! 3289: K = 0 is the main diagonal, K > 0 is above the main diagonal ! 3290: and K < 0 is below the main diagonal. ! 3291: .xx ! 3292: USER``Allows personal Fortran subroutines to be linked into MATLAB . ! 3293: The subroutine should have the heading ! 3294: .sp ! 3295: SUBROUTINE USER(A,M,N,S,T) ! 3296: REAL or DOUBLE PRECISION A(M,N),S,T ! 3297: .sp ! 3298: The MATLAB statement Y = USER(X,s,t) results in a call to ! 3299: the subroutine with a copy of the matrix X stored in the ! 3300: argument A , its column and row dimensions in M and N , ! 3301: and the scalar parameters s and t stored in S and T . ! 3302: If s and t are omitted, they are set to 0.0 . ! 3303: After the return, A is stored in Y . ! 3304: The dimensions M and N may be reset within the subroutine. ! 3305: The statement Y = USER(K) results in a call with ! 3306: M = 1, N = 1 and A(1,1) = FLOAT(K) . ! 3307: After the subroutine has been written, it must be compiled ! 3308: and linked to the MATLAB object code within the local operating ! 3309: system. ! 3310: .xx ! 3311: WHAT``Lists commands and functions currently available. ! 3312: .xx ! 3313: WHILE`Repeat statements an indefinite number of times. ! 3314: .br ! 3315: WHILE expr rop expr, statement, ..., statement, END ! 3316: .br ! 3317: where rop is =, <, >, <=, >=, or <> (not equal) . ! 3318: The END at the end of a line may be omitted. ! 3319: The comma before the END may also be omitted. ! 3320: The commas may be replaced by semicolons to avoid printing. ! 3321: The statements are repeatedly executed as long as the ! 3322: indicated comparison between the real parts of the first ! 3323: components of the two expressions is true. ! 3324: Example (assume a matrix A is already defined). ! 3325: .nf ! 3326: E = 0*A; F = E + EYE; N = 1; ! 3327: WHILE NORM(E+F-E,1) > 0, E = E + F; F = A*F/N; N = N + 1; ! 3328: E ! 3329: .fi ! 3330: .xx ! 3331: WHO```Lists current variables. ! 3332: .xx ! 3333: WHY```Provides succinct answers to any questions. ! 3334: .xx ! 3335: //
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.