% 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).