page@swan.ulowell.edu (Bob Page) (11/03/88)
Submitted-by: strovink%galaxy-43@afit-ab.arpa (Mark A. Strovink)
Posting-number: Volume 2, Issue 50
Archive-name: applications/matlab/doc.1
# This is a shell archive.
# Remove everything above and including the cut line.
# Then run the rest of the file through sh.
#----cut here-----cut here-----cut here-----cut here----#
#!/bin/sh
# shar: Shell Archiver
# Run the following text with /bin/sh to create:
# doc.1
# This archive created: Wed Nov 2 16:14:31 1988
cat << \SHAR_EOF > doc.1
6/24/81
MATLAB Users' Guide
May, 1981
Cleve Moler
Department of Computer Science
University of New Mexico
ABSTRACT. MATLAB is an interactive computer program
that serves as a convenient "laboratory" for
computations involving matrices. It provides easy
access to matrix software developed by the LINPACK and
EISPACK projects. The program is written in Fortran
and is designed to be readily installed under any
operating system which permits interactive execution of
Fortran programs.
CONTENTS
1. Elementary operations page 2
2. MATLAB functions 8
3. Rows, columns and submatrices 9
4. FOR, WHILE and IF 10
5. Commands, text, files and macros 12
6. Census example 13
7. Partial differential equation 19
8. Eigenvalue sensitivity example 23
9. Syntax diagrams 27
10. The parser-interpreter 31
11. The numerical algorithms 34
12. FLOP and CHOP 37
13. Communicating with other programs 41
Appendix. The HELP file 46
6/24/81
MATLAB Users' Guide
November, 1980
Cleve Moler
Department of Computer Science
University of New Mexico
MATLAB is an interactive computer program that serves as a
convenient "laboratory" for computations involving matrices. It
provides easy access to matrix software developed by the LINPACK
and EISPACK projects [1-3]. The capabilities range from standard
tasks such as solving simultaneous linear equations and inverting
matrices, through symmetric and nonsymmetric eigenvalue problems,
to fairly sophisticated matrix tools such as the singular value
decomposition.
It is expected that one of MATLAB's primary uses will be in
the classroom. It should be useful in introductory courses in
applied linear algebra, as well as more advanced courses in
numerical analysis, matrix theory, statistics and applications of
matrices to other disciplines. In nonacademic settings, MATLAB
can serve as a "desk calculator" for the quick solution of small
problems involving matrices.
The program is written in Fortran and is designed to be
readily installed under any operating system which permits
interactive execution of Fortran programs. The resources
required are fairly modest. There are less than 7000 lines of
Fortran source code, including the LINPACK and EISPACK
subroutines used. With proper use of overlays, it is possible
run the system on a minicomputer with only 32K bytes of memory.
The size of the matrices that can be handled in MATLAB
depends upon the amount of storage that is set aside when the
system is compiled on a particular machine. We have found that
an allocation of 5000 words for matrix elements is usually quite
satisfactory. This provides room for several 20 by 20 matrices,
for example. One implementation on a virtual memory system
provides 100,000 elements. Since most of the algorithms used
access memory in a sequential fashion, the large amount of
allocated storage causes no difficulties.
MATLAB, page 2
In some ways, MATLAB resembles SPEAKEASY [4] and, to a
lesser extent, APL. All are interactive terminal languages that
ordinarily accept single-line commands or statements, process
them immediately, and print the results. All have arrays or
matrices as principal data types. But for MATLAB, the matrix is
the only data type (although scalars, vectors and text are
special cases), the underlying system is portable and requires
fewer resources, and the supporting subroutines are more powerful
and, in some cases, have better numerical properties.
Together, LINPACK and EISPACK represent the state of the art
in software for matrix computation. EISPACK is a package of over
70 Fortran subroutines for various matrix eigenvalue computations
that are based for the most part on Algol procedures published by
Wilkinson, Reinsch and their colleagues [5]. LINPACK is a
package of 40 Fortran subroutines (in each of four data types)
for solving and analyzing simultaneous linear equations and
related matrix problems. Since MATLAB is not primarily concerned
with either execution time efficiency or storage savings, it
ignores most of the special matrix properties that LINPACK and
EISPACK subroutines use to advantage. Consequently, only 8
subroutines from LINPACK and 5 from EISPACK are actually
involved.
In more advanced applications, MATLAB can be used in
conjunction with other programs in several ways. It is possible
to define new MATLAB functions and add them to the system. With
most operating systems, it is possible to use the local file
system to pass matrices between MATLAB and other programs.
MATLAB command and statement input can be obtained from a local
file instead of from the terminal. The most power and
flexibility is obtained by using MATLAB as a subroutine which is
called by other programs.
This document first gives an overview of MATLAB from the
user's point of view. Several extended examples involving data
fitting, partial differential equations, eigenvalue sensitivity
and other topics are included. A formal definition of the MATLAB
language and an brief description of the parser and interpreter
are given. The system was designed and programmed using
techniques described by Wirth [6], implemented in nonrecursive,
portable Fortran. There is a brief discussion of some of the
matrix algorithms and of their numerical properties. The final
section describes how MATLAB can be used with other programs.
The appendix includes the HELP documentation available on-line.
1. Elementary operations
MATLAB works with essentially only one kind of object, a
rectangular matrix with complex elements. If the imaginary parts
of the elements are all zero, they are not printed, but they
MATLAB, page 3
still occupy storage. In some situations, special meaning is
attached to 1 by 1 matrices, that is scalars, and to 1 by n and m
by 1 matrices, that is row and column vectors.
Matrices can be introduced into MATLAB in four different
ways:
-- Explicit list of elements,
-- Use of FOR and WHILE statements,
-- Read from an external file,
-- Execute an external Fortran program.
The explicit list is surrounded by angle brackets, '<' and
'>', and uses the semicolon ';' to indicate the ends of the rows.
For example, the input line
A = <1 2 3; 4 5 6; 7 8 9>
will result in the output
A =
1. 2. 3.
4. 5. 6.
7. 8. 9.
The matrix A will be saved for later use. The individual
elements are separated by commas or blanks and can be any MATLAB
expressions, for example
x = < -1.3, 4/5, 4*atan(1) >
results in
X =
-1.3000 0.8000 3.1416
The elementary functions available include sqrt, log, exp, sin,
cos, atan, abs, round, real, imag, and conjg.
Large matrices can be spread across several input lines,
with the carriage returns replacing the semicolons. The above
matrix could also have been produced by
A = < 1 2 3
4 5 6
7 8 9 >
Matrices can be input from the local file system. Say a
file named 'xyz' contains five lines of text,
MATLAB, page 4
A = <
1 2 3
4 5 6
7 8 9
>;
then the MATLAB statement EXEC('xyz') reads the matrix and
assigns it to A .
The FOR statement allows the generation of matrices whose
elements are given by simple formulas. Our example matrix A
could also have been produced by
for i = 1:3, for j = 1:3, a(i,j) = 3*(i-1)+j;
The semicolon at the end of the line suppresses the printing,
which in this case would have been nine versions of A with
changing elements.
Several statements may be given on a line, separated by
semicolons or commas.
Two consecutive periods anywhere on a line indicate
continuation. The periods and any following characters are
deleted, then another line is input and concatenated onto the
previous line.
Two consecutive slashes anywhere on a line cause the
remainder of the line to be ignored. This is useful for
inserting comments.
Names of variables are formed by a letter, followed by any
number of letters and digits, but only the first 4 characters are
remembered.
The special character prime (') is used to denote the
transpose of a matrix, so
x = x'
changes the row vector above into the column vector
X =
-1.3000
0.8000
3.1416
Individual matrix elements may be referenced by enclosing
their subscripts in parentheses. When any element is changed,
the entire matrix is reprinted. For example, using the above
matrix,
MATLAB, page 5
a(3,3) = a(1,3) + a(3,1)
results in
A =
1. 2. 3.
4. 5. 6.
7. 8. 10.
Addition, subtraction and multiplication of matrices are
denoted by +, -, and * . The operations are performed whenever
the matrices have the proper dimensions. For example, with the
above A and x, the expressions A + x and x*A are incorrect
because A is 3 by 3 and x is now 3 by 1. However,
b = A*x
is correct and results in the output
B =
9.7248
17.6496
28.7159
Note that both upper and lower case letters are allowed for input
(on those systems which have both), but that lower case is
converted to upper case.
There are two "matrix division" symbols in MATLAB, \ and / .
(If your terminal does not have a backslash, use $ instead, or
see CHAR.) If A and B are matrices, then A\B and B/A correspond
formally to left and right multiplication of B by the inverse of
A , that is inv(A)*B and B*inv(A), but the result is obtained
directly without the computation of the inverse. In the scalar
case, 3\1 and 1/3 have the same value, namely one-third. In
general, A\B denotes the solution X to the equation A*X = B and
B/A denotes the solution to X*A = B.
Left division, A\B, is defined whenever B has as many rows
as A . If A is square, it is factored using Gaussian
elimination. The factors are used to solve the equations
A*X(:,j) = B(:,j) where B(:,j) denotes the j-th column of B. The
result is a matrix X with the same dimensions as B. If A is
nearly singular (according to the LINPACK condition estimator,
RCOND), a warning message is printed. If A is not square, it is
factored using Householder orthogonalization with column
pivoting. The factors are used to solve the under- or
overdetermined equations in a least squares sense. The result is
an m by n matrix X where m is the number of columns of A and n is
the number of columns of B . Each column of X has at most k
MATLAB, page 6
nonzero components, where k is the effective rank of A .
Right division, B/A, can be defined in terms of left
division by B/A = (A'\B')'.
For example, since our vector b was computed as A*x, the
statement
y = A\b
results in
Y =
-1.3000
0.8000
3.1416
Of course, y is not exactly equal to x because of the
roundoff errors involved in both A*x and A\b , but we are not
printing enough digits to see the difference. The result of the
statement
e = x - y
depends upon the particular computer being used. In one case it
produces
E =
1.0e-15 *
.3053
-.2498
.0000
The quantity 1.0e-15 is a scale factor which multiplies all the
components which follow. Thus our vectors x and y actually
agree to about 15 decimal places on this computer.
It is also possible to obtain element-by-element
multiplicative operations. If A and B have the same dimensions,
then A .* B denotes the matrix whose elements are simply the
products of the individual elements of A and B . The expressions
A ./ B and A .\ B give the quotients of the individual elements.
There are several possible output formats. The statement
long, x
results in
X =
MATLAB, page 7
-1.300000000000000
.800000000000000
3.141592653589793
The statement
short
restores the original format.
The expression A**p means A to the p-th power. It is
defined if A is a square matrix and p is a scalar. If p is an
integer greater than one, the power is computed by repeated
multiplication. For other values of p the calculation involves
the eigenvalues and eigenvectors of A.
Previously defined matrices and matrix expressions can be
used inside brackets to generate larger matrices, for example
C = <A, b; <4 2 0>*x, x'>
results in
C =
1.0000 2.0000 3.0000 9.7248
4.0000 5.0000 6.0000 17.6496
7.0000 8.0000 10.0000 28.7159
-3.6000 -1.3000 0.8000 3.1416
There are four predefined variables, EPS, FLOP, RAND and
EYE. The variable EPS is used as a tolerance is determining such
things as near singularity and rank. Its initial value is the
distance from 1.0 to the next largest floating point number on
the particular computer being used. The user may reset this to
any other value, including zero. EPS is changed by CHOP, which is
described in section 12.
The value of RAND is a random variable, with a choice of a
uniform or a normal distribution.
The name EYE is used in place of I to denote identity
matrices because I is often used as a subscript or as sqrt(-1).
The dimensions of EYE are determined by context. For example,
B = A + 3*EYE
adds 3 to the diagonal elements of A and
X = EYE/A
MATLAB, page 8
is one of several ways in MATLAB to invert a matrix.
FLOP provides a count of the number of floating point
operations, or "flops", required for each calculation.
A statement may consist of an expression alone, in which
case a variable named ANS is created and the result stored in ANS
for possible future use. Thus
A\A - EYE
is the same as
ANS = A\A - EYE
(Roundoff error usually causes this result to be a matrix of
"small" numbers, rather than all zeros.)
All computations are done using either single or double
precision real arithmetic, whichever is appropriate for the
particular computer. There is no mixed-precision arithmetic.
The Fortran COMPLEX data type is not used because many systems
create unnecessary underflows and overflows with complex
operations and because some systems do not allow double precision
complex arithmetic.
2. MATLAB functions
Much of MATLAB's computational power comes from the various
matrix functions available. The current list includes:
INV(A) - Inverse.
DET(A) - Determinant.
COND(A) - Condition number.
RCOND(A) - A measure of nearness to singularity.
EIG(A) - Eigenvalues and eigenvectors.
SCHUR(A) - Schur triangular form.
HESS(A) - Hessenberg or tridiagonal form.
POLY(A) - Characteristic polynomial.
SVD(A) - Singular value decomposition.
PINV(A,eps) - Pseudoinverse with optional tolerance.
RANK(A,eps) - Matrix rank with optional tolerance.
LU(A) - Factors from Gaussian elimination.
CHOL(A) - Factor from Cholesky factorization.
QR(A) - Factors from Householder orthogonalization.
RREF(A) - Reduced row echelon form.
ORTH(A) - Orthogonal vectors spanning range of A.
EXP(A) - e to the A.
LOG(A) - Natural logarithm.
SQRT(A) - Square root.
SIN(A) - Trigonometric sine.
COS(A) - Cosine.
MATLAB, page 9
ATAN(A) - Arctangent.
ROUND(A) - Round the elements to nearest integers.
ABS(A) - Absolute value of the elements.
REAL(A) - Real parts of the elements.
IMAG(A) - Imaginary parts of the elements.
CONJG(A) - Complex conjugate.
SUM(A) - Sum of the elements.
PROD(A) - Product of the elements.
DIAG(A) - Extract or create diagonal matrices.
TRIL(A) - Lower triangular part of A.
TRIU(A) - Upper triangular part of A.
NORM(A,p) - Norm with p = 1, 2 or 'Infinity'.
EYE(m,n) - Portion of identity matrix.
RAND(m,n) - Matrix with random elements.
ONES(m,n) - Matrix of all ones.
MAGIC(n) - Interesting test matrices.
HILBERT(n) - Inverse Hilbert matrices.
ROOTS(C) - Roots of polynomial with coefficients C.
DISPLAY(A,p) - Print base p representation of A.
KRON(A,B) - Kronecker tensor product of A and B.
PLOT(X,Y) - Plot Y as a function of X .
RAT(A) - Find "simple" rational approximation to A.
USER(A) - Function defined by external program.
Some of these functions have different interpretations when
the argument is a matrix or a vector and some of them have
additional optional arguments. Details are given in the HELP
document in the appendix.
Several of these functions can be used in a generalized
assignment statement with two or three variables on the left hand
side. For example
<X,D> = EIG(A)
stores the eigenvectors of A in the matrix X and a diagonal
matrix containing the eigenvalues in the matrix D. The statement
EIG(A)
simply computes the eigenvalues and stores them in ANS.
Future versions of MATLAB will probably include additional
functions, since they can easily be added to the system.
3. Rows, columns and submatrices
Individual elements of a matrix can be accessed by giving
their subscripts in parentheses, eg. A(1,2), x(i), TAB(ind(k)+1).
An expression used as a subscript is rounded to the nearest
MATLAB, page 10
integer.
Individual rows and columns can be accessed using a colon
':' (or a '|') for the free subscript. For example, A(1,:) is the
first row of A and A(:,j) is the j-th column. Thus
A(i,:) = A(i,:) + c*A(k,:)
adds c times the k-th row of A to the i-th row.
The colon is used in several other ways in MATLAB, but all
of the uses are based on the following definition.
j:k is the same as <j, j+1, ..., k>
j:k is empty if j > k .
j:i:k is the same as <j, j+i, j+2i, ..., k>
j:i:k is empty if i > 0 and j > k or if i < 0 and j < k .
The colon is usually used with integers, but it is possible to
use arbitrary real scalars as well. Thus
1:4 is the same as <1, 2, 3, 4>
0: 0.1: 0.5 is the same as <0.0, 0.1, 0.2, 0.3, 0.4, 0.5>
In general, a subscript can be a vector. If X and V are
vectors, then X(V) is <X(V(1)), X(V(2)), ..., X(V(n))> . This can
also be used with matrices. If V has m components and W has n
components, then A(V,W) is the m by n matrix formed from the
elements of A whose subscripts are the elements of V and W.
Combinations of the colon notation and the indirect subscripting
allow manipulation of various submatrices. For example,
A(<1,5>,:) = A(<5,1>,:) interchanges rows 1 and 5 of A.
A(2:k,1:n) is the submatrix formed from rows 2 through k
and columns 1 through n of A .
A(:,<3 1 2>) is a permutation of the first three columns.
The notation A(:) has a special meaning. On the right hand
side of an assignment statement, it denotes all the elements of
A, regarded as a single column. When an expression is assigned
to A(:), the current dimensions of A, rather than of the
expression, are used.
4. FOR, WHILE and IF
The FOR clause allows statements to be repeated a specific
number of times. The general form is
FOR variable = expr, statement, ..., statement, END
MATLAB, page 11
The END and the comma before it may be omitted. In general, the
expression may be a matrix, in which case the columns are stored
one at a time in the variable and the following statements, up to
the END or the end of the line, are executed. The expression is
often of the form j:k, and its "columns" are simply the scalars
from j to k. Some examples (assume n has already been assigned a
value):
for i = 1:n, for j = 1:n, A(i,j) = 1/(i+j-1);
generates the Hilbert matrix.
for j = 2:n-1, for i = j:n-1, ...
A(i,j) = 0; end; A(j,j) = j; end; A
changes all but the "outer edge" of the lower triangle and then
prints the final matrix.
for h = 1.0: -0.1: -1.0, (<h, cos(pi*h)>)
prints a table of cosines.
<X,D> = EIG(A); for v = X, v, A*v
displays eigenvectors, one at a time.
The WHILE clause allows statements to be repeated an
indefinite number of times. The general form is
WHILE expr relop expr, statement,..., statement, END
where relop is =, <, >, <=, >=, or <> (not equal) . The
statements are repeatedly executed as long as the indicated
comparison between the real parts of the first components of the
two expressions is true. Here are two examples. (Exercise for
the reader: What do these segments do?)
eps = 1;
while 1 + eps > 1, eps = eps/2;
eps = 2*eps
E = 0*A; F = E + EYE; n = 1;
while NORM(E+F-E,1) > 0, E = E + F; F = A*F/n; n = n + 1;
E
The IF clause allows conditional execution of statements.
The general form is
IF expr relop expr, statement, ..., statement,
ELSE statement, ..., statement
The first group of statements are executed if the relation is
MATLAB, page 12
true and the second group are executed if the relation is false.
The ELSE and the statements following it may be omitted. For
example,
if abs(i-j) = 2, A(i,j) = 0;
5. Commands, text, files and macros.
MATLAB has several commands which control the output format
and the overall execution of the system.
The HELP command allows on-line access to short portions of
text describing various operations, functions and special
characters. The entire HELP document is reproduced in an
appendix.
Results are usually printed in a scaled fixed point format
that shows 4 or 5 significant figures. The commands SHORT, LONG,
SHORT E, LONG E and LONG Z alter the output format, but do not
alter the precision of the computations or the internal storage.
The WHO, WHAT and WHY commands provide information about the
functions and variables that are currently defined.
The CLEAR command erases all variables, except EPS, FLOP,
RAND and EYE. The statement A = <> indicates that a "0 by 0"
matrix is to be stored in A. This causes A to be erased so that
its storage can be used for other variables.
The RETURN and EXIT commands cause return to the underlying
operating system through the Fortran RETURN statement.
MATLAB has a limited facility for handling text. Any string
of characters delineated by quotes (with two quotes used to allow
one quote within the string) is saved as a vector of integer
values with '1' = 1, 'A' = 10, ' ' = 36, etc. (The complete list
is in the appendix under CHAR.) For example
'2*A + 3' is the same as <2 43 10 36 41 36 3>
It is possible, though seldom very meaningful, to use such
strings in matrix operations. More frequently, the text is used
as a special argument to various functions.
NORM(A,'inf') computes the infinity norm of A .
DISPLAY(T) prints the text stored in T .
EXEC('file') obtains MATLAB input from an external file.
SAVE('file') stores all the current variables in a file.
LOAD('file') retrieves all the variables from a file.
PRINT('file',X) prints X on a file.
DIARY('file') makes a copy of the complete MATLAB session.
MATLAB, page 13
The text can also be used in a limited string substitution
macro facility. If a variable, say T, contains the source text
for a MATLAB statement or expression, then the construction
> T <
causes T to be executed or evaluated. For example
T = '2*A + 3';
S = 'B = >T< + 5'
A = 4;
> S <
produces
B =
16.
Some other examples are given under MACRO in the appendix. This
facility is useful for fairly short statements and expressions.
More complicated MATLAB "programs" should use the EXEC facility.
The operations which access external files cannot be handled
in a completely machine-independent manner by portable Fortran
code. It is necessary for each particular installation to
provide a subroutine which associates external text files with
Fortran logical unit numbers.
6. Census example
Our first extended example involves predicting the
population of the United States in 1980 using extrapolation of
various fits to the census data from 1900 through 1970. There
are eight observations, so we begin with the MATLAB statement
n = 8
The values of the dependent variable, the population in millions,
can be entered with
y = < 75.995 91.972 105.711 123.203 ...
131.669 150.697 179.323 203.212>'
In order to produce a reasonably scaled matrix, the independent
variable, time, is transformed from the interval [1900,1970] to
[-1.00,0.75]. This can be accomplished directly with
t = -1.0:0.25:0.75
MATLAB, page 14
or in a fancier, but perhaps clearer, way with
t = 1900:10:1970; t = (t - 1940*ones(t))/40
Either of these is equivalent to
t = <-1 -.75 -.50 -.25 0 .25 .50 .75>
The interpolating polynomial of degree n-1 involves an
Vandermonde matrix of order n with elements that might be
generated by
for i = 1:n, for j = 1:n, a(i,j) = t(i)**(j-1);
However, this results in an error caused by 0**0 when i = 5 and
j = 1 . The preferable approach is
A = ones(n,n);
for i = 1:n, for j = 2:n, a(i,j) = t(i)*a(i,j-1);
Now the statement
cond(A)
produces the output
ANS =
1.1819E+03
which indicates that transformation of the time variable has
resulted in a reasonably well conditioned matrix.
The statement
c = A\y
results in
C =
131.6690
41.0406
103.5396
262.4535
-326.0658
-662.0814
341.9022
533.6373
These are the coefficients in the interpolating polynomial
n-1
MATLAB, page 15
c + c t + ... + c t
1 2 n
Our transformation of the time variable has resulted in t = 1
corresponding to the year 1980. Consequently, the extrapolated
population is simply the sum of the coefficients. This can be
computed by
p = sum(c)
The result is
P =
426.0950
which indicates a 1980 population of over 426 million. Clearly,
using the seventh degree interpolating polynomial to extrapolate
even a fairly short distance beyond the end of the data interval
is not a good idea.
The coefficients in least squares fits by polynomials of
lower degree can be computed using fewer than n columns of the
matrix.
for k = 1:n, c = A(:,1:k)\y, p = sum(c)
would produce the coefficients of these fits, as well as the
resulting extrapolated population. If we do not want to print
all the coefficients, we can simply generate a small table of
populations predicted by polynomials of degrees zero through
seven. We also compute the maximum deviation between the fitted
and observed values.
for k = 1:n, X = A(:,1:k); c = X\y; ...
d(k) = k-1; p(k) = sum(c); e(k) = norm(X*c-y,'inf');
<d, p, e>
The resulting output is
0 132.7227 70.4892
1 211.5101 9.8079
2 227.7744 5.0354
3 241.9574 3.8941
4 234.2814 4.0643
5 189.7310 2.5066
6 118.3025 1.6741
7 426.0950 0.0000
The zeroth degree fit, 132.7 million, is the result of fitting a
constant to the data and is simply the average. The results
obtained with polynomials of degree one through four all appear
reasonable. The maximum deviation of the degree four fit is
MATLAB, page 16
slightly greater than the degree three, even though the sum of
the squares of the deviations is less. The coefficients of the
highest powers in the fits of degree five and six turn out to be
negative and the predicted populations of less than 200 million
are probably unrealistic. The hopefully absurd prediction of the
interpolating polynomial concludes the table.
We wish to emphasize that roundoff errors are not
significant here. Nearly identical results would be obtained on
other computers, or with other algorithms. The results simply
indicate the difficulties associated with extrapolation of
polynomial fits of even modest degree.
A stabilized fit by a seventh degree polynomial can be
obtained using the pseudoinverse, but it requires a fairly
delicate choice of a tolerance. The statement
s = svd(A)
produces the singular values
S =
3.4594
2.2121
1.0915
0.4879
0.1759
0.0617
0.0134
0.0029
We see that the last three singular values are less than 0.1 ,
consequently, A can be approximately by a matrix of rank five
with an error less than 0.1 . The Moore-Penrose pseudoinverse of
this rank five matrix is obtained from the singular value
decomposition with the following statements
c = pinv(A,0.1)*y, p = sum(c), e = norm(a*c-y,'inf')
The output is
MATLAB, page 17
C =
134.7972
67.5055
23.5523
9.2834
3.0174
2.6503
-2.8808
3.2467
P =
241.1720
E =
3.9469
The resulting seventh degree polynomial has coefficients which
are much smaller than those of the interpolating polynomial given
earlier. The predicted population and the maximum deviation are
reasonable. Any choice of the tolerance between the fifth and
sixth singular values would produce the same results, but choices
outside this range result in pseudoinverses of different rank and
do not work as well.
The one term exponential approximation
y(t) = k exp(pt)
can be transformed into a linear approximation by taking
logarithms.
log(y(t)) = log k + pt
= c + c t
1 2
The following segment makes use of the fact that a function of a
vector is the function applied to the individual components.
X = A(:,1:2);
c = X\log(y)
p = exp(sum(c))
e = norm(exp(X*c)-y,'inf')
The resulting output is
MATLAB, page 18
C =
4.9083
0.5407
P =
232.5134
E =
4.9141
The predicted population and maximum deviation appear
satisfactory and indicate that the exponential model is a
reasonable one to consider.
As a curiousity, we return to the degree six polynomial.
Since the coefficient of the high order term is negative and the
value of the polynomial at t = 1 is positive, it must have a root
at some value of t greater than one. The statements
X = A(:,1:7);
c = X\y;
c = c(7:-1:1); //reverse the order of the coefficients
z = roots(c)
produce
Z =
1.1023- 0.0000*i
0.3021+ 0.7293*i
-0.8790+ 0.6536*i
-1.2939- 0.0000*i
-0.8790- 0.6536*i
0.3021- 0.7293*i
There is only one real, positive root. The corresponding time on
the original scale is
1940 + 40*real(z(1))
= 1984.091
We conclude that the United States population should become zero
early in February of 1984.
MATLAB, page 19
7. Partial differential equation example
Our second extended example is a boundary value problem for
Laplace's equation. The underlying physical problem involves the
conductivity of a medium with cylindrical inclusions and is
considered by Keller and Sachs [7].
Find a function u(x,y) satisfying Laplace's equation
u + u = 0
xx yy
The domain is a unit square with a quarter circle of radius rho
removed from one corner. There are Neumann conditions on the top
and bottom edges and Dirichlet conditions on the remainder of the
boundary.
u = 0
n
-------------
| .
| .
| .
| . u = 1
| .
| .
| .
u = 0 | |
| |
| |
| | u = 1
| |
| |
| |
------------------------
u = 0
n
The effective conductivity of an medium is then given by the
integral along the left edge,
1
sigma = integral u (0,y) dy
0 n
It is of interest to study the relation between the radius rho
and the conductivity sigma. In particular, as rho approaches
one, sigma becomes infinite.
MATLAB, page 20
Keller and Sachs use a finite difference approximation. The
following technique makes use of the fact that the equation is
actually Laplace's equation and leads to a much smaller matrix
problem to solve.
Consider an approximate solution of the form
n 2j-1
u = sum c r cos(2j-1)t
j=1 j
where r,t are polar coordinates (t is theta). The coefficients
are to be determined. For any set of coefficients, this function
already satisfies the differential equation because the basis
functions are harmonic; it satisfies the normal derivative
boundary condition on the bottom edge of the domain because we
used cos t in preference to sin t ; and it satisfies the
boundary condition on the left edge of the domain because we use
only odd multiples of t .
The computational task is to find coefficients so that the
boundary conditions on the remaining edges are satisfied as well
as possible. To accomplish this, pick m points (r,t) on the
remaining edges. It is desirable to have m > n and in practice
we usually choose m to be two or three times as large as n .
Typical values of n are 10 or 20 and of m are 20 to 60. An
m by n matrix A is generated. The i,j element is the j-th
basis function, or its normal derivative, evaluated at the i-th
boundary point. A right hand side with m components is also
generated. In this example, the elements of the right hand side
are either zero or one. The coefficients are then found by
solving the overdetermined set of equations
Ac = b
in a least squares sense.
Once the coefficients have been determined, the approximate
solution is defined everywhere on the domain. It is then
possible to compute the effective conductivity sigma . In fact,
a very simple formula results,
n j-1
sigma = sum (-1) c
j=1 j
To use MATLAB for this problem, the following "program" is
first stored in the local computer file system, say under the
name "PDE".
MATLAB, page 21
//Conductivity example.
//Parameters ---
rho //radius of cylindrical inclusion
n //number of terms in solution
m //number of boundary points
//initialize operation counter
flop = <0 0>;
//initialize variables
m1 = round(m/3); //number of points on each straight edge
m2 = m - m1; //number of points with Dirichlet conditions
pi = 4*atan(1);
//generate points in Cartesian coordinates
//right hand edge
for i = 1:m1, x(i) = 1; y(i) = (1-rho)*(i-1)/(m1-1);
//top edge
for i = m2+1:m, x(i) = (1-rho)*(m-i)/(m-m2-1); y(i) = 1;
//circular edge
for i = m1+1:m2, t = pi/2*(i-m1)/(m2-m1+1); ...
x(i) = 1-rho*sin(t); y(i) = 1-rho*cos(t);
//convert to polar coordinates
for i = 1:m-1, th(i) = atan(y(i)/x(i)); ...
r(i) = sqrt(x(i)**2+y(i)**2);
th(m) = pi/2; r(m) = 1;
//generate matrix
//Dirichlet conditions
for i = 1:m2, for j = 1:n, k = 2*j-1; ...
a(i,j) = r(i)**k*cos(k*th(i));
//Neumann conditions
for i = m2+1:m, for j = 1:n, k = 2*j-1; ...
a(i,j) = k*r(i)**(k-1)*sin((k-1)*th(i));
//generate right hand side
for i = 1:m2, b(i) = 1;
for i = m2+1:m, b(i) = 0;
//solve for coefficients
c = A\b
//compute effective conductivity
c(2:2:n) = -c(2:2:n);
sigma = sum(c)
//output total operation count
ops = flop(2)
The program can be used within MATLAB by setting the three
parameters and then accessing the file. For example,
rho = .9;
n = 15;
m = 30;
exec('PDE')
The resulting output is
MATLAB, page 22
RHO =
.9000
N =
15.
M =
30.
C =
2.2275
-2.2724
1.1448
0.1455
-0.1678
-0.0005
-0.3785
0.2299
0.3228
-0.2242
-0.1311
0.0924
0.0310
-0.0154
-0.0038
SIGM =
5.0895
OPS =
16204.
A total of 16204 floating point operations were necessary to
set up the matrix, solve for the coefficients and compute the
conductivity. The operation count is roughly proportional to
m*n**2. The results obtained for sigma as a function of rho by
this approach are essentially the same as those obtained by the
finite difference technique of Keller and Sachs, but the
computational effort involved is much less.
MATLAB, page 23
8. Eigenvalue sensitivity example
In this example, we construct a matrix whose eigenvalues are
moderately sensitive to perturbations and then analyze that
sensitivity. We begin with the statement
B = <3 0 7; 0 2 0; 0 0 1>
which produces
B =
3. 0. 7.
0. 2. 0.
0. 0. 1.
Obviously, the eigenvalues of B are 1, 2 and 3 . Moreover,
since B is not symmetric, these eigenvalues are slightly
sensitive to perturbation. (The value b(1,3) = 7 was chosen so
that the elements of the matrix A below are less than 1000.)
We now generate a similarity transformation to disguise the
eigenvalues and make them more sensitive.
L = <1 0 0; 2 1 0; -3 4 1>, M = L\L'
L =
1. 0. 0.
2. 1. 0.
-3. 4. 1.
M =
1.0000 2.0000 -3.0000
-2.0000 -3.0000 10.0000
11.0000 18.0000 -48.0000
The matrix M has determinant equal to 1 and is moderately badly
conditioned. The similarity transformation is
A = M*B/M
A =
-64.0000 82.0000 21.0000
144.0000 -178.0000 -46.0000
-771.0000 962.0000 248.0000
Because det(M) = 1 , the elements of A would be exact integers
if there were no roundoff. So,
MATLAB, page 24
A = round(A)
A =
-64. 82. 21.
144. -178. -46.
-771. 962. 248.
This, then, is our test matrix. We can now forget how it
was generated and analyze its eigenvalues.
<X,D> = eig(A)
D =
3.0000 0.0000 0.0000
0.0000 1.0000 0.0000
0.0000 0.0000 2.0000
X =
-.0891 3.4903 41.8091
.1782 -9.1284 -62.7136
-.9800 46.4473 376.2818
Since A is similar to B, its eigenvalues are also 1, 2 and 3.
They happen to be computed in another order by the EISPACK
subroutines. The fact that the columns of X, which are the
eigenvectors, are so far from being orthonormal is our first
indication that the eigenvalues are sensitive. To see this
sensitivity, we display more figures of the computed eigenvalues.
long, diag(D)
ANS =
2.999999999973599
1.000000000015625
2.000000000011505
We see that, on this computer, the last five significant figures
are contaminated by roundoff error. A somewhat superficial
explanation of this is provided by
short, cond(X)
ANS =
3.2216e+05
The condition number of X gives an upper bound for the relative
error in the computed eigenvalues. However, this condition
MATLAB, page 25
number is affected by scaling.
X = X/diag(X(3,:)), cond(X)
X =
.0909 .0751 .1111
-.1818 -.1965 -.1667
1.0000 1.0000 1.0000
ANS =
1.7692e+03
Rescaling the eigenvectors so that their last components are
all equal to one has two consequences. The condition of X is
decreased by over two orders of magnitude. (This is about the
minimum condition that can be obtained by such diagonal scaling.)
Moreover, it is now apparent that the three eigenvectors are
nearly parallel.
More detailed information on the sensitivity of the
individual eigenvalues involves the left eigenvectors.
Y = inv(X'), Y'*A*X
Y =
-511.5000 259.5000 252.0000
616.0000 -346.0000 -270.0000
159.5000 -86.5000 -72.0000
ANS =
3.0000 .0000 .0000
.0000 1.0000 .0000
.0000 .0000 2.0000
We are now in a position to compute the sensitivities of the
individual eigenvalues.
for j = 1:3, c(j) = norm(Y(:,j))*norm(X(:,j)); end, C
C =
833.1092
450.7228
383.7564
These three numbers are the reciprocals of the cosines of the
angles between the left and right eigenvectors. It can be shown
that perturbation of the elements of A can result in a
MATLAB, page 26
perturbation of the j-th eigenvalue which is c(j) times as large.
In this example, the first eigenvalue has the largest
sensitivity.
We now proceed to show that A is close to a matrix with a
double eigenvalue. The direction of the required perturbation is
given by
E = -1.e-6*Y(:,1)*X(:,1)'
E =
1.0e-03 *
.0465 -.0930 .5115
-.0560 .1120 -.6160
-.0145 .0290 -.1595
With some trial and error which we do not show, we bracket the
point where two eigenvalues of a perturbed A coalesce and then
become complex.
eig(A + .4*E), eig(A + .5*E)
ANS =
1.1500
2.5996
2.2504
ANS =
2.4067 + .1753*i
2.4067 - .1753*i
1.1866 + 0.0000*i
Now, a bisecting search, driven by the imaginary part of one of
the eigenvalues, finds the point where two eigenvalues are nearly
equal.
r = .4; s = .5;
while s-r > 1.e-14, t = (r+s)/2; d = eig(A+t*E); ...
if imag(d(1))=0, r = t; else, s = t;
long, T
T =
.450380734134507
Finally, we display the perturbed matrix, which is obviously
MATLAB, page 27
close to the original, and its pair of nearly equal eigenvalues.
(We have dropped a few digits from the long output.)
A+t*E, eig(A+t*E)
A
-63.999979057 81.999958114 21.000230369
143.999974778 -177.999949557 -46.000277434
-771.000006530 962.000013061 247.999928164
ANS =
2.415741150
2.415740621
1.168517777
The first two eigenvectors of A + t*E are almost
indistinguishable indicating that the perturbed matrix is almost
defective.
<X,D> = eig(A+t*E); X = X/diag(X(3,:))
X =
.096019578 .096019586 .071608466
-.178329614 -.178329608 -.199190520
1.000000000 1.000000000 1.000000000
short, cond(X)
ANS =
3.3997e+09
9. Syntax diagrams
A formal description of the language acceptable to MATLAB,
as well as a flow chart of the MATLAB program, is provided by the
syntax diagrams or syntax graphs of Wirth [6]. There are eleven
non-terminal symbols in the language:
line, statement, clause, expression, term,
factor, number, integer, name, command, text .
The diagrams define each of the non-terminal symbols using the
others and the terminal symbols:
letter -- A through Z,
digit -- 0 through 9,
MATLAB, page 28
char -- ( ) ; : + - * / \ = . , < >
quote -- '
line
|-----> statement >----|
| |
|-----> clause >-------|
| |
-------|-----> expr >---------|------>
| | | |
| |-----> command >------| |
| | | |
| |-> > >-> expr >-> < >-| |
| | | |
| |----------------------| |
| |
| |-< ; <-| |
|--------| |---------|
|-< , <-|
statement
|-> name >--------------------------------|
| | |
| | |--> : >---| |
| | | | |
| |-> ( >---|-> expr >-|---> ) >-|
| | | |
-----| |-----< , <----| |--> = >--> expr >--->
| |
| |--< , <---| |
| | | |
|-> < >---> name >---> > >----------------|
MATLAB, page 29
clause
|---> FOR >---> name >---> = >---> expr >--------------|
| |
| |-> WHILE >-| |
|-| |-> expr >---------------------- |
| |-> IF >-| | | | | | | |
-----| < <= = <> >= > |---->
| | | | | | | |
| ----------------------> expr >--|
| |
|---> ELSE >--------------------------------------------|
| |
|---> END >--------------------------------------------|
expr
|-> + >-|
| |
-------|-------|-------> term >---------->
| | | |
|-> - >-| | |-< + <-| |
| | | |
|--|-< - <-|--|
| |
|-< : <-|
term
---------------------> factor >---------------------->
| |
| |-< * <-| |
| |-------| | | |-------| |
|--| |--|-< / <-|--| |--|
|-< . <-| | | |-< . <-|
|-< \ <-|
MATLAB, page 30
factor
|----------------> number >---------------|
| |
|-> name >--------------------------------|
| | |
| | |--> : >---| |
| | | | |
| |-> ( >---|-> expr >-|---> ) >-|
| | | |
| |-----< , <----| |
| |
-----|------------> ( >-----> expr >-----> ) >-|-|-------|----->
| | | | |
| |--------------| | |-> ' >-| |
| | | | |
|------------> < >-|---> expr >---|-> > >-| |
| | | | |
| |--< <---| | |
| | | | |
| |--< ; <---| | |
| | | | |
| |--< , <---| | |
| | |
|------------> > >-----> expr >-----> < >-| |
| | |
|-----> factor >---> ** >---> factor >----| |
| |
|------------> ' >-----> text >-----> ' >-------------|
number
|----------| |-> + >-|
| | | |
-----> int >-----> . >---> int >-----> E >---------> int >---->
| | | | | |
| | | |-> - >-| |
| | | |
|---------------------------------------------|
int
------------> digit >----------->
| |
|-----------|
MATLAB, page 31
name
|--< letter <--|
| |
------> letter >--|--------------|----->
| |
|--< digit <--|
command
|--> name >--|
| |
--------> name >--------|------------|---->
| |
|--> char >--|
| |
|---> ' >----|
text
|-> letter >--|
| |
|-> digit >---|
----------------| |-------------->
| |-> char >----| |
| | | |
| |-> ' >-> ' >-| |
| |
|---------------------|
10. The parser-interpreter
The structure of the parser-interpreter is similar to that
of Wirth's compiler [6] for his simple language, PL/0 , except
that MATLAB is programmed in Fortran, which does not have
explicit recursion. The interrelation of the primary subroutines
is shown in the following diagram.
MATLAB, page 32
MAIN
|
MATLAB |--CLAUSE
| | |
PARSE-----|--EXPR----TERM----FACTOR
| | | |
| |-------|-------|
| | | |
| STACK1 STACK2 STACKG
|
|--STACKP--PRINT
|
|--COMAND
|
|
| |--CGECO
| |
| |--CGEFA
| |
|--MATFN1--|--CGESL
| |
| |--CGEDI
| |
| |--CPOFA
|
|
| |--IMTQL2
| |
| |--HTRIDI
| |
|--MATFN2--|--HTRIBK
| |
| |--CORTH
| |
| |--COMQR3
|
|
|--MATFN3-----CSVDC
|
|
| |--CQRDC
|--MATFN4--|
| |--CQRSL
|
|
| |--FILES
|--MATFN5--|
|--SAVLOD
Subroutine PARSE controls the interpretation of each
statement. It calls subroutines that process the various
syntactic quantities such as command, expression, term and
factor. A fairly simple program stack mechanism allows these
MATLAB, page 33
subroutines to recursively "call" each other along the lines
allowed by the syntax diagrams. The four STACK subroutines
manage the variable memory and perform elementary operations,
such as matrix addition and transposition.
The four subroutines MATFN1 though MATFN4 are called
whenever "serious" matrix computations are required. They are
interface routines which call the various LINPACK and EISPACK
subroutines. MATFN5 primarily handles the file access tasks.
Two large real arrays, STKR and STKI, are used to store all
the matrices. Four integer arrays are used to store the names,
the row and column dimensions, and the pointers into the real
stacks. The following diagram illustrates this storage scheme.
TOP IDSTK MSTK NSTK LSTK STKR STKI
-- -- -- -- -- -- -- -- -------- --------
| |--->| | | | | | | | | | |----------->| | | |
-- -- -- -- -- -- -- -- -------- --------
| | | | | | | | | | | | | | |
-- -- -- -- -- -- -- -------- --------
. . . . . .
. . . . . .
. . . . . .
-- -- -- -- -- -- -- -------- --------
BOT | | | | | | | | | | | | | | |
-- -- -- -- -- -- -- -- -------- --------
| |--->| X| | | | | 2| | 1| | |----------->| 3.14 | | 0.00 |
-- -- -- -- -- -- -- -- -------- --------
| A| | | | | 2| | 2| | |--------- | 0.00 | | 1.00 |
-- -- -- -- -- -- -- \ -------- --------
| E| P| S| | | 1| | 1| | |------- ->| 11.00 | | 0.00 |
-- -- -- -- -- -- -- \ -------- --------
| F| L| O| P| | 1| | 2| | |------ \ | 21.00 | | 0.00 |
-- -- -- -- -- -- -- \ \ -------- --------
| E| Y| E| | |-1| |-1| | |--- \ | | 12.00 | | 0.00 |
-- -- -- -- -- -- -- \ | | -------- --------
| R| A| N| D| | 1| | 1| | |- \ | | | 22.00 | | 0.00 |
-- -- -- -- -- -- -- \ | \ \ -------- --------
| \ \ ->| 1.E-15 | | 0.00 |
\ \ \ -------- --------
\ \ ->| 0.00 | | 0.00 |
\ \ -------- --------
\ \ | 0.00 | | 0.00 |
\ \ -------- --------
\ ->| 1.00 | | 0.00 |
\ -------- --------
--->| URAND | | 0.00 |
-------- --------
The top portion of the stack is used for temporary variables
and the bottom portion for saved variables. The figure shows the
situation after the line
MATLAB, page 34
A = <11,12; 21,22>, x = <3.14, sqrt(-1)>'
has been processed. The four permanent names, EPS, FLOP, RAND
and EYE, occupy the last four positions of the variable stacks.
RAND has dimensions 1 by 1, but whenever its value is requested,
a random number generator is used instead. EYE has dimensions -1
by -1 to indicate that the actual dimensions must be determined
later by context. The two saved variables have dimensions 2 by 2
and 2 by 1 and so take up a total of 6 locations.
Subsequent statements involving A and x will result in
temporary copies being made in the top of the stack for use in
the actual calculations. Whenever the top of the stack reaches
the bottom, a message indicating memory has been exceeded is
printed, but the current variables are not affected.
This modular structure makes it possible to implement MATLAB
on a system with a limited amount of memory. The object code for
the MATFN's and the LINPACK-EISPACK subroutines is rarely needed.
Although it is not standard, many Fortran operating systems
provide some overlay mechanism so that this code is brought into
the main memory only when required. The variables, which occupy
a relatively small portion of the memory, remain in place, while
the subroutines which process them are loaded a few at a time.
11. The numerical algorithms
The algorithms underlying the basic MATLAB functions are
described in the LINPACK and EISPACK guides [1-3]. The following
list gives the subroutines used by these functions.
INV(A) - CGECO,CGEDI
DET(A) - CGECO,CGEDI
LU(A) - CGEFA
RCOND(A) - CGECO
CHOL(A) - CPOFA
SVD(A) - CSVDC
COND(A) - CSVDC
NORM(A,2) - CSVDC
PINV(A,eps) - CSVDC
RANK(A,eps) - CSVDC
QR(A) - CQRDC,CQRSL
ORTH(A) - CQRDC,CSQSL
A\B and B/A - CGECO,CGESL if A is square.
- CQRDC,CQRSL if A is not square.
EIG(A) - HTRIDI,IMTQL2,HTRIBK if A is Hermitian.
- CORTH,COMQR2 if A is not Hermitian.
SCHUR(A) - same as EIG.
SHAR_EOF
# End of shell archive
exit 0
--
Bob Page, U of Lowell CS Dept. page@swan.ulowell.edu ulowell!page
Have five nice days.