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