ENME462 - Studio 4

Chris Hubley
11/19/2003

Question 1

Given any matrices A and B. Write a MATLAB function to check if the system of equations AX = B is indeterminate or not.

Solution:

syscheck.m 
function [] = syscheck(A,B)
%SYSCHECK denotes whether or not matricies A and B will fit into the
%         equation A x X = B and give solutions for all x-values.
[m1 n1] = size(A);
[m2 n2] = size(B);
if n2 ~= 1
    error('Indeterminant: Matrix B is not a column matrix');
elseif m1 ~= m2
    error('Indeterminant: Matricies A and B must have the same number of rows');
end
[X, pivot] = rref([A,B]);
if length(pivot) ~= m1
    error('Indeterminant: Incorrect number of linearly independant rows');
else
    disp('System is OK');
end
return;

Question 2

Given any matrices A and B. Without using ‘\’ operator and by using rref write a function to compute solutions to AX = B. Make sure it handles all the cases stated above. In the case when rank (A) < n (underdetermined case) just find the particular solution such that the arbitrary variables are assigned value zero. Next using this function find the solutions of the examples 1 – 3 above. Check the correctness of your code by computing the residual vectors in each case when the solution exists. Do your solutions agree with the solutions given by ‘\’ operator for these examples. If not comment why?

Solution:

systemsolver.m 
function [X] = systemsolver(A,B)
K=[A,B];
X=rref(K);
return;

>> A = [1 2 3; 4 5 6; 7 8 10];
>> B = ones(3,1);
>> systemsolver(A,B)

ans =

     1     0     0    -1
     0     1     0     1
     0     0     1     0

>> A = [1 2 3; 4 5 6];
>> B = ones (2,1);
>> systemsolver(A,B)

ans =

     1     0    -1    -1
     0     1     2     1

>> A = [2 -1; 1 10; 1 2];
>> B = ones(3,1);
>> systemsolver(A,B)

ans =

     1     0     0
     0     1     0
     0     0     1
The solution is not necessarily the same as using the '\' operator because A\B requires that rank([A,B])=rank(A). In any other case, it may give an estimation or not be able to complete the operation.

Question 3

Given any square matrix A. Using the MATLAB function rref, Write a function to compute the inverse of A.
Hint: Think what will happen if we compute the Row Echelon matrix of [A, I] where I is the unit matrix of size A.

Solution:

matrixinverse.m 
function [B]=matrixinverse(A)
[a b]=size(A);
K=[A,eye(a,b)];
J=rref(K);
B=J(:,(b+1):2*b);
return;

Question 4

Write a function cofact(A, k, l) to compute the cofactor c[k, l] of the A[k, l] entry of the matrix A.
Hint: First write the minor M[k, l] as a sub matrix.

Solution:

cofact.m 
function [X]=cofact(A,m,n)
[a,b]=size(A);
if m>a | n>b | (a~=b)
    error('Index exceeds matrix dimensions or matrix is not square. Word.')
end
c=1;
for k=1:a
    if k ~= m
        B(c,:)=[A(k,1:(n-1)), A(k,(n+1):b)];
        c=c+1;
    end
end
X=((-1)^(m+n))*det(B);
return;

Question 5

Construct a function d = rDet(A) to compute determinant which, implements the method of cofactor expansion. Make sure that you don’t use the MATLAB function det anywhere in your code.
Hint: Use recursion.
Method of cofactors has a high computational complexity. Therefore it is not recommended for computations with large matrices. It is included here for pedagogical reasons only.

Solution:

rDet.m
Find more information about matrix determinants at:
http://www.richland.cc.il.us/james/lecture/m116/matrices/determinant.html 
function d=rDet(A)

[a,b]=size(A);
if a~=b
    error('Matrix must be square.');
end
if a==1
    d=A; return;
end
m=1; x=0;
for n=1:b
    c=1;
    for k=1:a
        if k ~= m
            B(c,:)=[A(k,1:(n-1)), A(k,(n+1):b)];
            c=c+1;
        end
    end
    if size(B)==[2 2]
        x=A(m,n)*((-1)^(m+n)*(B(1,1)*B(2, 2)-B(2, 1)*B(1, 2)))+x;
    else
       x=A(m,n)*((-1)^(m+n)*rDet(B))+x;
    end
end
d=x;
return;

Question 6

The adjoint matrix adj(A) of the matrix A is also of interest in linear algebra. Write a function adj(A) to compute adjoint of matrix.
Hint: Just use the cofact function and the MATLAB function reshape.

Solution:

adjoint.m
Find more information about adjoint matrices at:
http://www.bath.ac.uk/~cs1spw/notes/Maths4Apps/nav.htm 
function Z=adjoint(A)
[a,b]=size(A);
for m=1:a
    for n=1:b
        c=1;
        for k=1:a
            if k ~= m
                B(c,:)=[A(k,1:(n-1)), A(k,(n+1):b)];
                c=c+1;
            end
        end
        X(m,n)=((-1)^(m+n))*det(B);
    end
end
Z=transpose(X);
return;

Question 7

Show that the adjoint matrix and the inverse matrix satisfy the equation inv(A) = adj(A)/det(A)

Solution:

Please click here for the proof.
A =

     4     2     5
     5     3     7
     1     9     2

>> Q=adjoint(A)

Q =

   -57    41    -1
    -3     3    -3
    42   -34     2

>> R=det(A)

R =

   -24

>> X=Q/R

X =

    2.3750   -1.7083    0.0417
    0.1250   -0.1250    0.1250
   -1.7500    1.4167   -0.0833

>> Y=inv(A)

Y =

    2.3750   -1.7083    0.0417
    0.1250   -0.1250    0.1250
   -1.7500    1.4167   -0.0833
Since X = Y, it follows that inv(A) = adj(A)/det(A)

Question 8

The Caley-Hamilton Theorem states that each matrix satisfies its characteristic equation, i.e., chpol(A) = 0, where the last zero stands for the matrix of zeros of the appropriate dimension.

Write a MATLAB code to verify the above theorem.

Solution:

Please click here for the proof.
Here is some information on eigenvectors/eigenvalues:
http://www.mathphysics.com/calc/eigen.html
>> A=magic(3)

A =

     8     1     6
     3     5     7
     4     9     2

>> C=poly(A)

C =

    1.0000  -15.0000  -24.0000  360.0000

>> X=polyvalm(C,A)

X =

  1.0e-012 *

    0.4547   -0.5116   -0.0853
   -0.3340    0.1705    0.0568
   -0.2700    0.1705         0
X is approximately a matrix of zeros. The residual values are due to roundoff errors.

Question 9

Find all the false statements from the list given below.
Note: To show a statement is false you have to find a counterexample (use MATLAB for that, if necessary).
Let A be a square matrix of size n.
  1. Number of distinct eigenvalues of A = n.
    FALSE
    A matrix "A" of size n × n has "n" eigenvalues, however they are not necessisarily distinct. 
    E=eig(eye(3))

E =

     1
     1
     1
  2. Distinct eigenvalues will have linearly independent (we can also say distinct) eigenvectors.
    TRUE
    (ref: http://oregonstate.edu/dept/math/CalculusQuestStudyGuides/vcalc/eigen/eigen.html)
  3. For an eigenvalue of multiplicity = 1, there is one and only one linearly independent eigenvector.
    TRUE
  4. For an eigenvalue of multiplicity > 1, there is one and only one linearly independent eigenvector.
    FALSE
    >> A=eye(3);
>> [V D]=eig(A)

V =

     1     0     0
     0     1     0
     0     0     1


D =

     1     0     0
     0     1     0
     0     0     1
  5. If all the eigenvalues are distinct then the matrix is diagonalizable.
    TRUE
    (ref: http://isolatium.uhh.hawaii.edu/linear/ch6/ch62.htm)
  6. If all the eigenvalues are not distinct then the matrix is not diagonalizable.
    FALSE
    >> A=eye(3);
>> [E, S]=eig(A);
>> S*A*inv(S)

ans =

     1     0     0
     0     1     0
     0     0     1

>> E

E =

     1     0     0
     0     1     0
     0     0     1
  7. Trace of A = sum of the eigenvalues.
    TRUE
    (ref: http://www.quantlet.com/mdstat/scripts/mva/htmlbook/mvahtmlnode14.html)
  8. Determinant of A = product of eigenvalues.
    TRUE
    (ref: http://turnbull.mcs.st-and.ac.uk/~john/geometry/Lectures/L3.html)