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