% Refer to Mathematics Matlab tutorial: % The MathWorks Inc, % Matlab Mathematics, % (c) Copyright 1984-2004. % % This file composed of several components of the Matlab % Tutorial provides some useful reminders concerning: % % A - MATRIX POWERS AND EXPONENTIALS % B - EIGENVALUES % C - SINGULAR VALUE DECOMPOSITION % D - VECTOR AND MATRIX NORMS echo on % A - MATRIX POWERS AND EXPONENTIALS % (1) Element-by-Element Powers: % The .^ operator produces element-by-element powers. For example, % X = A.^2 % A = % 1 1 1 % 1 4 9 % 1 9 36 % % (2) Exponentials: % The function % sqrtm(A) % computes A^(1/2) by a more accurate algorithm. The m in sqrtm % distinguishes this function from sqrt(A) which, like A.^(1/2), % does its job element-by-element. % A system of linear, constant coefficient, ordinary differential % equations can be written % % dx/dt = A*x % % where x = x(t) is a vector of functions of t and A is a matrix % independent of t. The solution can be expressed in terms of the % matrix exponential, % % x(t) = e^(t*A)*x(0) % % The function % expm(A) % computes the matrix exponential. An example is provided by the % 3-by-3 coefficient matrix % A = % 0 -6 -1 % 6 2 -16 % -5 20 -10 % and the initial condition, x(0) % x0 = % 1 % 1 % 1 % The matrix exponential is used to compute the solution, x(t), to % the differential equation at 101 points on the interval 0<=t<=1 with % X = []; % for t = 0:.01:1 % X = [X expm(t*A)*x0]; % end % % A three-dimensional phase plane plot obtained with % plot3(X(1,:),X(2,:),X(3,:),'-o') % shows the solution spiraling in towards the origin. This behavior is % related to the eigenvalues of the coefficient matrix, which are % discussed in the next section. echo off disp(' ') disp('Press Any Key to plot ..') disp(' ...the three-dimensional phase plane.') disp(' ') pause h=figure; x0 = [1 1 1]'; A =[0 -6 -1 6 2 -16 -5 20 -10]; X = []; for t = 0:.01:1 X = [X expm(t*A)*x0]; end plot3(X(1,:),X(2,:),X(3,:),'-o'); title(' X = [X expm(t*A)*x0]'); box; disp(' ') disp('Press Any Key to Continue ...') disp(' ') pause close(h) clear all echo on % B - EIGENVALUES % An eigenvalue and eigenvector of a square matrix A are a scalar "lambda" % (denoted here l=lambda) and a nonzero vector v that satisfy: % Av = lv % This section explains: % . Eigenvalue decomposition % . Problems associated with defective (not diagonalizable) matrices % . The use of Schur decomposition to avoid problems associated with % eigenvalue decomposition % % Eigenvalue Decomposition: % With the eigenvalues on the diagonal of a diagonal matrix LAMBDA (denoted % here L=LAMBDA) and the corresponding eigenvectors forming the columns of % a matrix V, we have % AV = VL % if V is nonsingular, this becomes the eigenvalue decomposition % A = VLV^-1 % A good example is provided by the coefficient matrix of the ordinary % differential equation in the previous section % A = [0 -6 -1 % 6 2 -16 % -5 20 -10] % The statement % lambda = eig(A) % produces a column vector containing the eigenvalues. For this matrix, % the eigenvalues are complex. % lambda = % -3.0710 % -2.4645+17.6008i % -2.4645-17.6008i % The real part of each of the eigenvalues is negative, so e^(lambda*t) % approaches zero as t increases. The nonzero imaginary part of two of the % eigenvalues, +/- omega, contributes the oscillatory component, sin(omega*t), % to the solution of the differential equation. % [V,D] = eig(A) % % V = % -0.8326 0.2003 - 0.1394i 0.2003 + 0.1394i % -0.3553 -0.2110 - 0.6447i -0.2110 + 0.6447i % -0.4248 -0.6930 -0.6930 % % D = % -3.0710 0 0 % 0 -2.4645+17.6008i 0 % 0 0 -2.4645-17.6008i % % The first eigenvector is real and the other two vectors are complex % conjugates of each other. All three vectors are normalized to have % Euclidean length, norm(v,2), equal to one. % The matrix V*D*inv(V), which can be written more succinctly as V*D/V, % is within roundoff error of A. And, inv(V)*A*V, or V\A*V, is within % roundoff error of D. % % Defective Matrices: % Some matrices do not have an eigenvector decomposition. These matrices % are defective, or not diagonalizable. For example % A = [ 6 12 19 % -9 -20 -33 % 4 9 15 ] % For this matrix % [V,D] = eig(A) % produces % V = % -0.4741 -0.4082 -0.4082 % 0.8127 0.8165 0.8165 % -0.3386 -0.4082 -0.4082 % % D = % -1.0000 0 0 % 0 1.0000 0 % 0 0 1.0000 % % There is a double eigenvalue at lambda=1. The second and third columns % of V are the same. For this matrix, a full set of linearly independent % eigenvectors does not exist. % The optional "Symbolic Math Toolbox" extends the capabilities of MATLAB by % connecting to Maple, a powerful computer algebra system. One of the % functions provided by the toolbox computes the Jordan Canonical Form. % This is appropriate for matrices like our example, which is 3-by-3 and % has exactly known, integer elements. % [X,J] = jordan(A) % % X = % -1.7500 1.5000 2.7500 % 3.0000 -3.0000 -3.0000 % -1.2500 1.5000 1.2500 % % J = % -1 0 0 % 0 1 1 % 0 0 1 % The Jordan Canonical Form is an important theoretical concept, but it is % not a reliable computational tool for larger matrices, or for matrices % whose elements are subject to roundoff errors and other uncertainties. % % Schur Decomposition in MATLAB Matrix Computations % The MATLAB advanced matrix computations do not require eigenvalue % decompositions. They are based, instead, on the Schur decomposition % % A = U S U^T % % where U is an orthogonal matrix and S is a block upper triangular matrix % with 1-by-1 and 2-by-2 blocks on the diagonal. The eigenvalues are % revealed by the diagonal elements and blocks of S, while the columns of % U provide a basis with much better numerical properties than a set of % eigenvectors. The Schur decomposition of our defective example is % [U,S] = schur(A) % % U = % -0.4741 0.6648 0.5774 % 0.8127 0.0782 0.5774 % -0.3386 -0.7430 0.5774 % % S = % -1.0000 20.7846 -44.6948 % 0 1.0000 -0.6096 % 0 0 1.0000 % The double eigenvalue is contained in the lower 2-by-2 block of S. % If A is complex, "schur" returns the complex Schur form, which is % upper triangular with the eigenvalues of A on the diagonal. echo off pause disp(' ') disp('Press Any Key to Continue ...') disp(' ') echo on % C - SINGULAR VALUE DECOMPOSITION % A singular value and corresponding singular vectors of a rectangular % matrix A are a scalar "sigma" and a pair of vectors u and v that satisfy % A v = sigma u % A^T u = sigma v % With the singular values on the diagonal of a diagonal matrix SIGMA and % the corresponding singular vectors forming the columns of two orthogonal % matrices U and V, we have % A V = U SIGMA % A^T U = V SIGMA % Since U and V are orthogonal, this becomes the singular value % decomposition % A = U SIGMA V^T % The full singular value decomposition of an m-by-n matrix involves an % m-by-m U, an m-by-n SIGMA, and an n-by-n V. In other words, U and V are % both square and SIGMA is the same size as A. % The eigenvalue decomposition is the appropriate tool for analyzing a % matrix when it represents a mapping from a vector space into itself, as % it does for an ordinary differential equation. On the other hand, the % singular value decomposition is the appropriate tool for analyzing a % mapping from one vector space into another vector space, possibly with a % different dimension. Most systems of simultaneous linear equations fall % into this second category. % If A is square, symmetric, and positive definite, then its eigenvalue and % singular value decompositions are the same. But, as A departs from % symmetry and positive definiteness, the difference between the two % decompositions increases. In particular, the singular value decomposition % of a real matrix is always real, but the eigenvalue decomposition of a % real, nonsymmetric matrix might be complex. % For the example matrix % A = % 9 4 % 6 8 % 2 7 % the full singular value decomposition is % [U,S,V] = svd(A) % U = % -0.6105 0.7174 0.3355 % -0.6646 -0.2336 -0.7098 % -0.4308 -0.6563 0.6194 % S = % 14.9359 0 % 0 5.1883 % 0 0 % V = % -0.6925 0.7214 % -0.7214 -0.6925 % You can verify that U*S*V' is equal to A to within roundoff error. % For this small problem, the economy size decomposition is only slightly % smaller. % [U,S,V] = svd(A,0) % U = % -0.6105 0.7174 % -0.6646 -0.2336 % -0.4308 -0.6563 % S = % 14.9359 0 % 0 5.1883 % V = % -0.6925 0.7214 % -0.7214 -0.6925 % Again, U*S*V' is equal to A to within roundoff error. % echo off disp(' ') disp('Press Any Key to Continue ...') disp(' ') pause echo on % D - VECTOR AND MATRIX NORMS % (1) The p-norm of a vector x: % % ||x||_p = [ Sum |x_i|^p ]^1/p % % is computed by norm(x,p). % This is defined by any value of p > 1, but the most common values of p % are 1, 2, and infinity. The default value is p = 2, which corresponds to % Euclidean length. % % (2) The p-norm of a matrix A: % % ||A||_p = max_x [ ||A*x||_p / ||x||_p ] % % can be computed for p = 1, 2, and infinity by norm(A,p). % Again, the default value is p = 2. % echo off F = fix(10*rand(3,2)); [norm(F,1) norm(F) norm(F,inf)] % The MathWorks, Matlab7 - Copyright (c).