LINEAR ALGEBRA

Matrix Manipulation & Solution of Systems of Linear Equations

using

MatLab

By Professor Mohamad H. Hassoun

Department of Electrical & Computer Engineering

Wayne Satate University

Detroit, MI 48202

October 19, 1996

Contents:

I. Vectors & Matrices

II. Matrix Algebra

III. Matrix Determinant and Inverse

IV. Solution of Linear Systems of Equations

V. Simple programs: Creating and Plotting Arrays

I. Vectors & Matrices

Let us start by defining a 1 by 3 vector b (a 3-dimensional row vector)

b = [1 2.5 (-1/3)]; % semicolon after a MatLab statement suppresses display.

and display it

format short % sets numerical display format to short decimal format

b

b =

1.0000 2.5000 -0.3333

Now define a new vector y as the transpose of b (b^T)

y = b'

y =

1.0000

2.5000

-0.3333

This is a quick way of creating a special row vecotr:

z = 1:5

z =

1 2 3 4 5

Let us look at another example:

z=(1:2:10)'

z =

1

3

5

7

9

Note: MatLab allows you to change the number display format

format bank

b

b =

1.00 2.50 -0.33

format rat

b

b =

1 5/2 -1/3

format short e

b

b =

1.0000e+000 2.5000e+000 -3.3333e-001

format long

b

b =

1.00000000000000 2.50000000000000 -0.33333333333333

Now back to the original (default) short format

format short

Let us define a 3 by 3 matrix A

A = [3,2,1; 0,-1/3,1/5; 2,1,-1]

A =

3.0000 2.0000 1.0000

0 -0.3333 0.2000

2.0000 1.0000 -1.0000

and a 3 by 2 matrix B

B = [1 3; -1 1.5; 4 pi]

B =

1.0000 3.0000

-1.0000 1.5000

4.0000 3.1416

Note: pi = 3.1415... is one of a number of predefined Matlab variables. The following are four other predefined variables (eps, inf, NaN, I):

eps

eps =

2.2204e-016

1/0

Warning: Divide by zero

ans =

Inf

0/0

Warning: Divide by zero

ans =

NaN % NaN stands for not a number

i % An imaginary number equals to the square-root of -1

ans =

0 + 1.0000i % See the Module on complex numbers for more details

sqrt(-1)

ans =

0 + 1.0000i

II. Matrix Algebra

The following illustrates various matrix operations (transpose, addition, multiplication)

C = A'

C =

3.0000 0 2.0000

2.0000 -0.3333 1.0000

1.0000 0.2000 -1.0000

A + C

ans =

6.0000 2.0000 3.0000

2.0000 -0.6667 1.2000

3.0000 1.2000 -2.0000

A - A

ans =

0 0 0

0 0 0

0 0 0

3*A

ans =

9.0000 6.0000 3.0000

0 -1.0000 0.6000

6.0000 3.0000 -3.0000

A * B % note that B * A, is not defined (why?)

ans =

5.0000 15.1416

1.1333 0.1283

-3.0000 4.3584

A * b' % b is the vector defined above; A * b is not defined (why?)

ans =

7.6667

-0.9000

4.8333

-A + A'*A

ans =

10.0000 6.0000 0

8.0000 5.4444 0.7333

-1.0000 -0.0667 3.0400

Other means of defining special matrices and vectors (arrays)

a = 1:5

a =

1 2 3 4 5

a1 = 1:2:6, a2 = 10:-1:0

a1 =

1 3 5

a2 =

10 9 8 7 6 5 4 3 2 1 0

E = [1:4; 3:6]

E =

1 2 3 4

3 4 5 6

I = diag(A) % extracts the diagonal elements of matrix A as a column vector

I =

3.0000

-0.3333

-1.0000

diag(I) % creates a diagonal matrix with the vector I on the diagonal

ans =

3.0000 0 0

0 -0.3333 0

0 0 -1.0000

eye(4)

ans =

1 0 0 0

0 1 0 0

0 0 1 0

0 0 0 1

ones(3)

ans =

1 1 1

1 1 1

1 1 1

zeros(2)

ans =

0 0

0 0

R = rand(3:3) % creates a 3 by 3 random matrix with uniform values in [0 1]

R =

0.9304 0.0920 0.7012

0.8462 0.6539 0.9103

0.5269 0.4160 0.7622

Check out these additional operations on the components of a vector or a matrix

b + 1

ans =

2.0000 3.5000 0.6667

which is equivalent to the addition of the two vectors:

b + [1,1,1]

ans =

2.0000 3.5000 0.6667

The following command subtracts 2 from each and every component of the matrix A

A - 2

ans =

1.0000 0 -1.0000

-2.0000 -2.3333 -1.8000

0 -1.0000 -3.0000

b .* b % component-wise multiplication of two vectors: no space between '.' and '*'.

ans =

1.0000 6.2500 0.1111

A .^ 2 % component-wise squaring of matrix elements

ans =

9.0000 4.0000 1.0000

0 0.1111 0.0400

4.0000 1.0000 1.0000

Compare this last result to A*A = A^2

A^2

ans =

11.0000 6.3333 2.4000

0.4000 0.3111 -0.2667

4.0000 2.6667 3.2000

A*A

ans =

11.0000 6.3333 2.4000

0.4000 0.3111 -0.2667

4.0000 2.6667 3.2000

Some logical operations

M = A > 0

M =

1 1 1

0 0 1

1 1 0

N = A <= 1

N =

0 0 1

1 1 1

0 1 1

min(b)

ans =

-0.3333

max(A)

ans =

3 2 1

More operations: Eigenvalues of a matrix

e = eig(A)

e =

3.4970

-1.4662

-0.3641

we can display the individual values as follows

e(1)

ans =

3.4970

e(2)

ans =

-1.4662

e(3)

ans =

-0.3641

Eigenvectors: [V,D] = EIG(A) produces a diagonal matrix D of eigenvalues and a full matrix V whose columns are the corresponding eigenvectors so that A*V = V*D.

[V,D] = eig(A)

V =

0.9117 -0.1412 -0.4768

0.0214 -0.1721 0.8688

0.4102 0.9749 -0.1335

D =

3.4970 0 0

0 -1.4662 0

0 0 -0.3641

We may now verify the second eigenvector-eigenvalue pair (v₂ , l₂) to see that Av₂ = l₂v₂

A*V(:,2)

ans =

0.2071

0.2524

-1.4294

D(2,2)*V(:,2)

ans =

0.2071

0.2524

-1.4294

Characteristic polynomial associated with a square matrix (|lI - A| = 0):

p = poly(A)

p =

1.0000 -1.6667 -5.8667 -1.8667

which represents the polynomial l³ - 1.67 l² - 5.87 l - 1.87 = 0. One may also synthesize a matrix (companion matrix) from a given characteristic polynomial

compan(p)

ans =

1.6667 5.8667 1.8667

1.0000 0 0

0 1.0000 0

Verification of result

poly(ans)

ans =

1.0000 -1.6667 -5.8667 -1.8667

Rank of a matrix

r = rank(A) % Matrix rank

r =

3

Here is a rank computation of another 3 by 3 matrix

rank([1, 1, 1; 2, 2, 2; 1, 2, 1]) % note that the second row is a linear scaling of the first

ans =

2

trace(A) % Matrix trace

ans =

1.6667

Compare 'ans' to the sum of the diagonal elements of A

A(1,1) + A(2,2) + A(3,3)

ans =

1.6667

Size of a matrix/vector:

size(A)

ans =

3 3

size(b)

ans =

1 3

k = length(b)

k =

3

abs(A)

ans =

3.0000 2.0000 1.0000

0 0.3333 0.2000

2.0000 1.0000 1.0000

III. Matrix Determinant and Inverse

Let us define a new matrix A1

A1 = [1/2, 1/2,1 ,-1; 1, -1, -1, -1; -1, -1/4, 1, 1;1,-1,-8,1]

A1 =

0.5000 0.5000 1.0000 -1.0000

1.0000 -1.0000 -1.0000 -1.0000

-1.0000 -0.2500 1.0000 1.0000

1.0000 -1.0000 -8.0000 1.0000

The determinant of A1 is

D = det(A1)

D =

0.6250

det(A1')

ans =

0.6250

compare this result to the product of all eigenvalues of A1

e = eig(A1)

e =

1.4307 + 3.0894i

1.4307 - 3.0894i

-1.3205

-0.0408

Note that eigenvalues can be complex (here i² = -1). Check the module on complex numbers for details.

e(1)*e(2)*e(3)*e(4) % this is the product of the four eigenvalues of A1

ans =

0.6250

The following is the inverse of matrix A1

inv(A1)

ans =

-18.0000 -1.2000 -15.2000 -4.0000

0 -0.8000 -0.8000 0

-4.0000 -0.2000 -3.2000 -1.0000

-14.0000 -1.2000 -11.2000 -3.0000

Check A1^-1 A1

inv(A1)*A1

ans =

1.0000 0.0000 0.0000 0

0 1.0000 0 0

0 0.0000 1.0000 0

0 0 0.0000 1.0000

A1*inv(A1)

ans =

1.0000 0.0000 0.0000 0

0 1.0000 0.0000 0

0 0.0000 1.0000 0

0.0000 0.0000 0.0000 1.0000

IV. Solution of Linear Systems of Equations

Consider the linear system of three unknowns given by

a₁₁x₁ + a₁₂ x₂ + a₁₃ x₃ = y₁

a₂₁x₁ + a₂₂ x₂ + a₂₃ x₃ = y₂

a₃₁x₁ + a₃₂ x₂ + a₃₃ x₃ = y₃

which can be expressed in matrix form as Ax = y. We may now solve for x in Ax = y as x = A^-1 y:

First, let A be as defined above and let y be

y = [1; 2; 3]

y =

1

2

3

x = inv(A) * y

x =

4.4643

-6.1071

-0.1786

A preferrable solution is one based on LU factorization approach (modification of Gaussian elimination) since it is faster and generally more accurate; here it is

x = A\y

x =

4.4643

-6.1071

-0.1786

Another Example with symbolic display of results:

G = sym('[1/2, 1/2,1 ,-1; 1, -1, -1, -1; -1, -1/4, 1, 1;1,-1,-8,1]')

G =

[1/2, 1/2, 1 , -1]

[ 1, -1, -1, -1]

[ -1, -1/4, 1, 1]

[ 1, -1, -8, 1]

d = sym('[0;-10;0;-1]')

d =

[ 0]

[-10]

[ 0]

[ -1]

x = linsolve(G,d) % linsolve is equivalent to the backslash operator used above

x =

[16]

[ 8]

[ 3]

[15]

Linsolve also allows for symbolic A2 and/or b2. Let us consider an example:

A2 = sym('1,2,t;0,1,0;t,1,0')

b2 = sym('2;1/t;3')

x = linsolve(A2,b2)

A2 =

[1,2,t]

[0,1,0]

[t,1,0]

b2 =

[ 2]

[1/t]

[ 3]

x =

[ (3*t-1)/t^2]

[ 1/t]

[(2*t^2-5*t+1)/t^3]

One can simplify or alter the above symbolic result in several ways by using symbolic manipulation commands such as 'simplify', 'collect', 'expand', and 'factor'. The above result for x may be simplified to

simplify

ans =

[(2*t+2+3*t^2)/t]

[ 1/t]

[ (2*t^2+1)/t]

The same solution for x may be obtained by using direct 'symbolic' inversion of matrix A2, as follows:

inverse(A2)

ans =

[ 0, -1/t, 1/t]

[ 0, 1, 0]

[1/t, -(-1+2*t)/t^2, -1/t^2]

Next, the symbolic maltiplication 'symmul' command is used to solve for x (check also 'symadd', 'symdiv', and 'symsub')

x = symmul(ans,b2)

x =

[ (3*t-1)/t^2]

[ 1/t]

[(2*t^2-5*t+1)/t^3]

Here are some simplifications of the last expression:

collect

ans =

[ -1/t^2+3/t]

[ 1/t]

[2/t-5/t^2+1/t^3]

factor

ans =

[ (3*t-1)/t^2]

[ 1/t]

[(2*t^2-5*t+1)/t^3]

expand

ans =

[ -1/t^2+3/t]

[ 1/t]

[2/t-5/t^2+1/t^3]

V. Simple programs: Creating and Plotting Arrays

A "for loop" may be used to create a 10-dimensional row vector/array, as follows:

for n=1:10

z(n)=sin(n*pi/10);

end

z

z =

Columns 1 through 7

0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090

Columns 8 through 10

0.5878 0.3090 0.0000

A more efficient way of creating the above vector is to avoid the for loop

n=1:10;

z=sin(n*pi/10)

z =

Columns 1 through 7

0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090

Columns 8 through 10

0.5878 0.3090 0.0000

Plotting Arrays

Let us now generate a simple plot for the array z

plot(z)

For a better plot, use

plot(z, '*')

title('My First Plot')

xlabel('Position Number in the Array'), ylabel('Value')

Now, take a look at this 2-D plot

for n = 1:10

for m=1:10

O(n,m) = sin((n+m)*pi/10);

end

end

mesh(O)

A better view can be arrived at by rotating the plot by 90 degrees:

mesh(rot90(O))

Yet, an alternative way to visualize the array O is possible by using the pseudocolor plot command 'pcolor'

pcolor(O)

Here is a bipolar version of the above plot

pcolor(sign(O))

colormap(gray(2))

We conclude this module with the following random pattern:

pcolor(randn(10,10))

colormap(gray(2))