Просмотреть файл в отдельном окне: f98a3d9dabb8f208cf92c109c34c626c.pdf

1

Matrix Algebra

and Solution

of Matrix Equations

1.1 INTRODUCTION

Computers are best suited for repetitive calculations and for organizing data into

specialized forms. In this chapter, we review the matrix and vector notation and

their manipulations and applications. Vector is a one-dimensional array of numbers

and/or characters arranged as a single column. The number of rows is called the

order of that vector. Matrix is an extension of vector when a set of numbers and/or

characters are arranged in rectangular form. If it has M rows and N column, this

matrix then is said to be of order M by N. When M = N, then we say this square

matrix is of order N (or M). It is obvious that vector is a special case of matrix when

there is only one column. Consequently, a vector is referred to as a column matrix

as opposed to the row matrix which has only one row. Braces are conventionally

used to indicate a vector such as {V} and brackets are for a matrix such as [M].

In writing a computer program, DIMENSION or DIM statements are necessary

to declare that a certain variable is a vector or a matrix. Such statements instruct

the computer to assign multiple memory spaces for keeping the values of that vector

or matrix. When we deal with a large number of different entities in a group, it is

better to arrange these entities in vector or matrix form and refer to a particular

entity by specifying where it is located in that group by pointing to the row (and

column) number(s). Such as in the case of having 100 numbers represented by the

variable names A, B, …, or by A(1) through A(100), the former requires 100 different

characters or combinations of characters and the latter certainly has the advantage

of having only one name. The A(1) through A(100) arrangement is to adopt a vector;

these numbers can also be arranged in a matrix of 10 rows and 10 columns, or 20

rows and five columns depending on the characteristics of these numbers. In the

cases of collecting the engineering data from tests of 20 samples during five different

days, then arranging these 100 data into a matrix of 20 rows and five columns will

be better than of 10 rows and 10 columns because each column contains the data

collected during a particular day.

In the ensuing sections, we shall introduce more definitions related to vector

and matrix such as transpose, inverse, and determinant, and discuss their manipulations such as addition, subtraction, and multiplication, leading to the organizing of

systems of linear algebraic equations into matrix equations and to the methods of

finding their solutions, specifically the Gaussian Elimination method. An apparent

application of the matrix equation is the transformation of the coordinate axes by a

© 2001 by CRC Press LLC

rotation about any one of the three axes. It leads to the derivation of the three basic

transformation matrices and will be elaborated in detail.

Since the interactive operations of modern personal computers are emphasized

in this textbook, how a simple three-dimensional brick can be displayed will be

discussed. As an extended application of the display monitor, the transformation of

coordinate axes will be applied to demonstrate how animation can be designed to

simulate the continuous rotation of the three-dimensional brick. In fact, any threedimensional object could be selected and its motion animated on a display screen.

Programming languages, FORTRAN, QuickBASIC, MATLAB, and Mathematica are to be initiated in this chapter and continuously expanded into higher

levels of sophistication in the later chapters to guide the readers into building a

collection of their own programs while learning the computational methods for

solving engineering problems.

1.2 MANIPULATION OF MATRICES

Two matrices [A] and [B] can be added or subtracted if they are of same order, say

M by N which means both having M rows and N columns. If the sum and difference

matrices are denoted as [S] and [D], respectively, and they are related to [A] and

[B] by the formulas [S] = [A] + [B] and [D] = [A]-[B], and if we denote the elements

in [A], [B], [D], and [S] as aij, bij, dij, and sij for i = 1 to M and j = 1 to N, respectively,

then the elements in [S] and [D] are to be calculated with the equations:

sij = a ij + b ij

(1)

d ij = a ij ? b ij

(2)

and

Equations 1 and 2 indicate that the element in the ith row and jth column of [S]

is the sum of the elements at the same location in [A] and [B], and the one in [D]

is to be calculated by subtracting the one in [B] from that in [A] at the same location.

To obtain all elements in the sum matrix [S] and the difference matrix [D], the index

i runs from 1 to M and the index j runs from 1 to N.

In the case of vector addition and subtraction, only one column is involved (N =

1). As an example of addition and subtraction of two vectors, consider the two

vectors in a two-dimensional space as shown in Figure 1, one vector {V1} is directed

from the origin of the x-y coordinate axes, point O, to the point 1 on the x-axis

which has coordinates (x1,y1) = (4,0) and the other vector {V2} is directed from the

origin O to the point 2 on the y-axis which has coordinates (x2,y2) = (0,3). One may

want to find the resultant of {R} = {V1} + {V2} which is the vector directed from

the origin to the point 3 whose coordinates are (x3,y3) = (4,3), or, one may want to

find the difference vector {D} = {V1} – {V2} which is the vector directed from the

origin O to the point 4 whose coordinates are (x4,y4) = (4,–3). In fact, the vector

{D} can be obtained by adding {V1} to the negative image of {V2}, namely {V2–}

which is a vector directed from the origin O to the point 5 whose coordinates are

(x5,y5). Mathematically, based on Equations 1 and 2, we can have:

© 2001 by CRC Press LLC

? 4 ? ?0 ?

?4?

? 4 ? ?0 ?

?4?

{R} = {V1} + {V2 } = ? ? + ? ? = ? ?

?0 ? ? 3? ? 3 ?

and

{D} = {V1} ? {V2 } = ? ? ? ? ? = ? ?

?0 ? ?3? ??3?

When Equation 1 is applied to two arbitrary two-dimensional vectors which

unlike {V1}, {V2}, and {V2–} but are not on either one of the coordinate axes, such

as {D} and {E} in Figure 1, we then have the sum vector {F} = {D} + {E} which

has components of 1 and –2 units along the x- and y-directions, respectively. Notice

that O467 forms a parallelogram in Figure 1 and the two vectors {D} and {E} are

the two adjacent sides of the parallelogram at O. To find the sum vector {F} of {D}

and {E} graphically, we simply draw a diagonal line from O to the opposite vertex

of the parallelogram — this is the well-known Law of Parallelogram.

It should be evident that to write out a vector which has a large number of rows

will take up a lot of space. If this vector can be rotated to become from one column

to one row, space saving would then be possible. This process is called transposition

as we will be leading to it by first introducing the length of a vector.

For the calculation of the length of a two-dimensional or three-dimensional vector,

such as {V1} and {V2} in Figure 1, it would be a simple matter because they are

oriented along the directions of the coordinate axes. But for the vectors such as {R}

FIGURE 1. Two vectors in a two-dimensional space.

© 2001 by CRC Press LLC

and {D} shown in Figure 1, the calculation of their lengths would need to know the

components of these vectors in the coordinate axes and then apply the Pythagorean

theorem. Since the vector {R} has components equal to rx = 4 and ry = 3 units along

the x- and y-axis, respectively, its length, here denoted with the symbol , is:

{R} = [rx2 + ry2 ] = [42 + 32 ] = 5

0.5

0.5

(3)

To facilitate the calculation of the length of a generalized vector {V} which has

N components, denoted as v1 through vN, its length is to be calculated with the

following formula obtained from extending Equation 3 from two-dimensions to Ndimensions:

{V} = [v12 + v22 + … + v2N ]

0.5

(4)

For example, a three-dimensional vector has components v1 = vx = 4, v2 = vy =

3, and v3 = vz = 12, then the length of this vector is {V} = [42 + 32 + 122]0.5 = 13.

We shall next show that Equation 4 can also be derived through the introduction of

the multiplication rule and transposition of matrices.

1.2 MULTIPLICATION OF MATRICES

A matrix [A] of order L (rows) by M (columns) and a matrix [B] of order M

by N can be multiplied in the order of [A][B] to produce a new matrix [P] of order

L by N. [A][B] is said as [A] post-multiplied by [B], or, [B] pre-multiplied by [A].

The elements in [P] denoted as pij for i = 1 to N and j = 1 to M are to be calculated

by the formula:

p ij =

?

M

k =1

a ik b kj

(5)

Equation 5 indicates that the value of the element pij in the ith row and jth column

of the product matrix [P] is to be calculated by multiplying the elements in the ith

row of the matrix [A] by the corresponding elements in the jth column of the matrix

[B]. It is therefore evident that the number of elements in the ith row of [A] should

be equal to the number of elements in the jth column of [B]. In other words, to

apply Equation 5 for producing a product matrix [P] by multiplying a matrix [A]

on the right by a matrix [B] (or, to say multiplying a matrix [B] on the left by a

matrix [A]), the number of columns of [A] should be equal to the number of row

of [B]. A matrix [A] of order L by M can therefore be post-multiplied by a matrix

[B] of order M by N; but [A] cannot be pre-multiplied by [B] unless L is equal to N!

As a numerical example, consider the case of a square, 3 ? 3 matrix postmultiplied by a rectangular matrix of order 3 by 2. Since L = 3, M = 3, and N = 2,

the product matrix is thus of order 3 by 2.

© 2001 by CRC Press LLC

?1

?4

?

?? 7

3? ? 6

6?? ??5

9?? ??4

2

5

8

?3? ?1(6) + 2(5) + 3( 4) 1( ?3) + 2( ?2) + 3( ?1) ?

?

?

?2 ?? = ?4(6) + 5(5) + 6( 4) 4( ?3) + 5( ?2) + 6( ?1)?

?1?? ?? 7(6) + 8(5) + 9( 4) 7( ?3) + 8( ?2) + 9( ?1)??

? 6 + 10 + 12

= ?? 24 + 25 + 24

??42 + 40 + 32

?3 ? 4 ? 3 ? ? 28

?12 ? 10 ? 5?? = ?? 73

?21 ? 16 ? 9 ?? ??114

?10 ?

?27??

?46??

More exercises are given in the Problems listed at the end of this chapter for

the readers to practice on the matrix multiplications based on Equation 5.

It is of interest to note that the square of the length of a vector {V} which has

N components as defined in Equation 4, {V}2, can be obtained by application of

Equation 5 to {V} and its transpose denoted as {V}T which is a row matrix of order

1 by N (one row and N columns). That is:

{V} = {V}T {V} = v11 + v22 + … + v32

2

(6)

For a L-by-M matrix having elements eij where the row index i ranges from 1

to L and the column index j ranges from 1 to M, the transpose of this matrix when

its elements are designated as trc will have a value equal to ecr where the row index

r ranges from 1 to M and the column index c ranges from 1 to M because this

transpose matrix is of order M by L. As a numerical example, here is a pair of a

3 ? 2 matrix [G] and its 2 ? 3 transpose [H]:

?6

[G] = ??5

3? 2

??4

?3?

?2 ?? and

?1??

?6

??3

[H] = [G]T = ?

2 ?3

5

?2

4?

?1??

If the elements of [G] and [H] are designated respectively as gij and hij, then

hij = gji. For example, from above, we observe that h12 = g21 = 5, h23 = g32 = –1, and

so on. There will be more examples of applications of Equations 5 and 6 in the

ensuing sections and chapters.

Having introduced the transpose of a matrix, we can now conveniently revisit

the addition of {D} and {E} in Figure 1 in algebraic form as {F} = {D} + {E} =

[4 –3]T + [–3 1]T = [4+(–3) –3+1]T = [1 –2]T. The resulting sum vector is indeed

correct as it is graphically verified in Figure 1. The saving of space by use of

transposes of vectors (row matrices) is not evident in this case because all vectors

are two-dimensional; imagine if the vectors are of much higher order.

Another noteworthy application of matrix multiplication and transposition is to

reduce a system of linear algebraic equations into a simple, (or, should we say a

single) matrix equation. For example, if we have three unknowns x, y, and z which

are to be solved from the following three linear algebraic equations:

© 2001 by CRC Press LLC

x + 2 y + 3z = 4

5x + 6y + 7z = 8

?2 x ? 37 = 9

(7)

Let us introduce two vectors, {V} and {R}, which contain the unknown x, y,

and z, and the right-hand-side constants in the above three equations, respectively.

That is:

?x ?

T

V

=

x

y

z

=

{ } [

] ??y??

?? z ??

and

?4?

T

R

=

4

8

9

=

{ } [

] ??8??

??9 ??

(8)

Then, making use of the multiplication rule of matrices, Equation 5, the system

of linear algebraic equations, 7, now can be written simply as:

[C]{V} = {R}

(9)

where the coefficient matrix [C] formed by listing the coefficients of x, y, and z in

first equation in the first row and second equation in the second row and so on. That is,

?1

[C] = ?? 5

???2

2

6

?3

3?

7??

0 ??

There will be more applications of matrix multiplication and transposition in

the ensuing chapters when we discuss how matrix equations, such as [C]{V} = {R},

can be solved by employing the Gaussian Elimination method, and how ordinary

differential equations are approximated by finite differences will lead to the matrix

equations. In the abbreviated matrix form, derivation and explanation of computational methods becomes much simpler.

Also, it can be observed from the expressions in Equation 8 how the transposition

can be conveniently used to define the two vectors not using the column matrices

which take more lines.

FORTRAN VERSION

Since Equations 1 and 2 require repetitive computation of the elements in the

sum matrix [S] and difference matrix [D], machine could certainly help to carry out

this laborous task particularly when matrices of very high order are involved. For

covering all rows and columns of [S] and [D], looping or application of DO statement

of the FORTRAN programming immediately come to mind. The following program

is provided to serve as a first example for generating [S] and [D] of two given

matrices [A] and[B]:

© 2001 by CRC Press LLC

The resulting display on the screen is:

To review FORTRAN briefly, we notice that matrices should be declared as

variables with two subscripts in a DIMENSION statement. The displayed results of

matrices A and B show that the values listed between // in a DATA statment will be

filling into the first column and then second column and so on of a matrix. To instruct

the computer to take the values provided but to fill them into a matrix row-by-row,

a more explicit DATA needs to be given as:

DATA ((A(I,J),J = 1,3),I = 1,3)/1.,4.,7.,2.,5.,8.,3.,6.,9./

When a number needs to be repeated, the * symbol can be conveniently applied

in the DATA statement exemplified by those for the matrix [B].

Some sample WRITE and FORMAT statements are also given in the program.

The first * inside the parentheses of the WRITE statement when replaced by a

number allows a device unit to be specified for saving the message or the values of

the variables listed in the statement. * without being replaced means the monitor

will be the output unit and consequently the message or the value of the variable(s)

will be displayed on screen. The second * inside the parentheses of the WRITE

© 2001 by CRC Press LLC

statement if not replaced by a statement number, in which formats for printing the

listed variables are specified, means “unformatted” and takes whatever the computer

provides. For example, statement number 15 is a FORMAT statement used by the

WRITE statement preceding it. There are 18 variables listed in that WRITE statement

but only six F5.1 codes are specified. F5.1 requests five column spaces and one digit

after the decimal point to be used to print the value of a listed variable. / in a

FORMAT statement causes the print/display to begin at the first column of the next

line. 6F5.1 is, however, enclosed by the inner pair of parentheses that allows it to

be reused and every time it is reused the next six values will be printed or displayed

on next line. The use (*,*) in a WRITE statement has the convenience of viewing

the results and then making a hardcopy on a connected printer by pressing the PrtSc

(Print Screen) key.

INTERACTIVE OPERATION

Program MatxAlgb.1 only allows the two particular matrices having their elements specified in the DATA statement to be added and subtracted. For finding the

sum matrix [S] and difference matrix [D] for any two matrices of same order N, we

ought to upgrade this program to allow the user to enter from keyboard the order

N and then the elements of the two matrices involved. This is interactive operation

of the program and proper messages should be given to instruct the user what to do

which means the program should be user-friendly. The program MatxAlgb.2 listed

below is an attempt to achieve that goal:

© 2001 by CRC Press LLC

The interactive execution of the problem solved by the previous version Matxalgb.1

now can proceed as follows:

© 2001 by CRC Press LLC

The results are identical to those obtained previously. The READ statement

allows the values for the variable(s) to be entered via keyboard. A WRITE statement

has no variable listed serves for need of skipping a line to provide better readability

of the display. Also the I and E format codes are introduced in the statement 10. Iw

where w is an integer in a FORMAT statement requests w columns to be provided

for displaying the value of the integer variable listed in the WRITE statement, in

which the FORMAT statement is utilized. Ew.d where w and d should both be integer

constants requests w columns to be provided for display a real value in the scientific

form and carrying d digits after the decimal point. Ew.d format gives more feasibility

than Fw.d format because the latter may cause an error message of insufficient width

if the value to be displayed becomes too large and/or has a negative sign.

MORE PROGRAMMING REVIEW

Besides the operation of matrix addition and subtraction, we have also discussed

about the transposition and multiplication of matrices. For further review of computer

programming, it is opportune to incorporate all these matrix algebraic operations

into a single interactive program. In the listing below, three subroutines for matrix

addition and subtraction, transposition, and multiplication named as MatrixSD,

Transpos, and MatxMtpy, respectively, are created to support a program called

MatxAlgb (Matrix Algebra).

© 2001 by CRC Press LLC

© 2001 by CRC Press LLC

The above program shows that Subroutines are independent units all started with

a SUBROUTINE statement which includes a name followed by a pair of parentheses

enclosing a number of arguments. The Subroutines are called in the main program

by specifying which variables or constants should serve as arguments to connect to

the subroutines. Some arguments provide input to the subroutine while other arguments transmit out the results determined by the subroutine. These are referred to

as input arguments and output arguments, respectively. In many instances, an argument may serve a dual role for both input and output purposes. To construct as an

independent unit, a subprogram which can be in the form of a SUBROUTINE, or

FUNCTION (to be elaborated later) must have RETURN and END statements.

It should also be remarked that program MatxAlgb is arranged to handle any

matrix having an order of no higher than 25 by 25. For this restriction and for having

the flexibility of handling any matrices of lesser order, the Lmax, Mmax, and Nmax

arguments are added in all three subroutines in order not to cause any mismatch of

matrix sizes between the main program and the called subroutine when dealing with

any L, M, and N values which are interactively entered via keyboard.

Computed GOTO and arithmetic IF statements are also introduced in the program MatxAlgb. GOTO (i,j,k,…) C will result in going to (execute) the statement

numbered i, j, k, and so on when C has a value equal to 1, 2, 3, and so on, respectively.

IF (Expression) a,b,c will result in going to the statement numbered a, b, or c if the

value calculated by the expression or a single variable is less than, equal to, or,

greater than zero, respectively.

It is important to point out that in describing any derived procedure of numerical

computation, indicial notation such as Equation 5 should always be preferred to

facilitate programming. In that notation, the indices are directly used, or, literally

translated into the index variables for the DO loops as can be seen in Subroutine

MatxMtpy which is developed according to Equation 5. Subroutine MatrixSD is

another example of literally translating Equations 1 and 2. For defining the values

of the element in the following tri-diagonal band matrix:

© 2001 by CRC Press LLC

?1

??3

?

[C ] = ? 0

?

?0

?? 0

2

1

?3

0

0

0

2

1

?3

0

0

0

2

1

?3

0?

0?

?

0?

?

2?

1 ??

we ought not to write 25 separate statements for the 25 elements in this matrix but

derive the indicial formulas for i,j = 1 to 5:

cij = 0,

if j > i + 2, or, j < i ? 2

ci,i +1 = 2,

and

ci,i ?1 = ?3

Then, the matrix [C] can be generated with the DO loops as follows:

The above short program also demonstrates the use of the CONTINUE statement for ending the DO loop(s), and the logical IF statements. The true, or, false

condition of the expression inside the outer pair of parentheses directs the computer

to execute the statement following the parentheses or the next statement immediately

below the current IF statement. Reader may want to practice on deriving indicial

formulas and then write a short program for calculating the elements of the matrix:

?1

?2

?

?3

?

4

[M] = ??

5

?

?6

?7

?

??8

© 2001 by CRC Press LLC

0

1

2

0

0

1

0

0

0

0

0

0

0

0

0

0

0

0

3

4

5

6

7

2

3

4

5

6

1

2

3

4

5

0

1

2

3

4

0

0

1

2

3

0

0

0

1

2

0?

0?

?

0?

?

0?

0?

?

0?

0?

?

1 ??

(10)

As another example of writing a computer program based on indicial notation,

consider the case of calculating ex based on the infinite series:

ex = 1 +

?

=

x1 x 2 x 2

xi

+

+

+…+ +…

1! 2! 3!

i!

?

i=0

xi

i!

(11)

With the understanding that 0! = 1, we have expressed the series as a summation

involving the index i which ranges from zero to infinity. A FUNCTION ExpoFunc

can be developed for calculating ex based on Equation 11 and taking only a finite

number of terms for a partial sum of the series when the contribution of additional

term is less than certain percentage of the sum in magnitude, say 0.001%. This

FUNCTION may be arranged as:

To further show the advantage of adopting vector and matrix notation, here let

us apply FUNCTION ExpoFunc to examine the surface z(x,y) = ex + y above the

rectangular area 0?x?2.0 and 0?y?1.5. The following program, ExpTest, will enable

us to compare the values of ex + y generated by the FUNCTION ExpoFunc and by

the function EXP available in the FORTRAN library (hence called library function).

© 2001 by CRC Press LLC

The resulting printout is:

It is apparent that two approaches produce almost identical results, so the 0.001%

accuracy appears quite adequate for the x and y ranges studied. Also, arranging the

results in vector and matrix forms make the presentation much easy to comprehend.

We have experienced how the summation process for an indicial formula involving a ? should be programmed. Another operation symbol of importance is ? which

is for multiplication of many factors. That is:

N

? a = a a …a

i

1 2

N

(12)

i =1

An obvious application of Equation 12 is for the calculation of factorials. For

example, 5! = ?i for i ranges from 1 to 5. As an exercise, we display the values of

1! through 50! with the following program involving a subroutine IFACTO which

calculates I! for a specified I value:

© 2001 by CRC Press LLC

The resulting print out is (listed in three columns for saving space)

Another application of Equation 12 is for calculation of the binomial coefficients

for a real number r and an integer k defined as:

? r ? r( r ? 1)( r ? 2)…( r ? k + 1)

=

? ?=

? k?

k!

k

?

i =1

r ? i +1

i

(13)

We shall have the occasion of applying Equations 12 and 13 when the finite

differences and Lagrangian interpolation are discussed.

Sample Applications

Program MatxAlgb has been tested interactively, the following are the resulting

displays of four test cases:

© 2001 by CRC Press LLC

© 2001 by CRC Press LLC

QUICKBASIC VERSION

The QuickBASIC language has the advantage over the FORTRAN language

for making quick changes and then running the revised program without compilation.

Furthermore, it offers simple plotting statements. Let us have a QuickBASIC version

of the program MatxAlgb and then discuss its basic features.

© 2001 by CRC Press LLC

© 2001 by CRC Press LLC

Notice that the order limit of 25 needed in the FORTRAN version is removed

in the QuickBASIC version which allows the dim statement to be adjustable. ' is

replacing C in FORTRAN to indicate a comment statement in QuickBASIC. READ

and WRITE in FORTRAN are replaced by INPUT and PRINT in QuickBASIC,

respectively. The DO loop in FORTRAN is replaced by the FOR and NEXT pair

in QuickBASIC.

Sample Applications

When the four cases previously run by the FORTRAN version are executed by

the QuickBASIC version, the screen prompting messages, the interactively entered

data, and the computed results are:

© 2001 by CRC Press LLC

© 2001 by CRC Press LLC

Matrix [P]

Row 1

0.1400E+02

Row 2

0.3200E+02

0.2000E+02

0.4700E+02

MATLAB APPLICATIONS

MATLAB developed by the Mathworks, Inc. offers a quick tool for matrix

manipulations. To load MATLAB after it has been set-up and stored in a subdirectory

of a hard drive, say C, we first switch to this subdirectory by entering (followed by

pressing ENTER)

C:\cd MATLAB

and then switch to its own subdirectory BIN by entering (followed by pressing ENTER)

C:\MATLAB>cd BIN

Next, we type MATLAB to obtain a display of:

C:\MATLAB>BIB>MATLAB

Pressing the ENTER key results in a display of:

>>

which indicates MATLAB is ready to begin. Let us rerun the cases of matrix

subtraction, addition, transposition, and multiplication previously considered in the

FORTRAN and QuickBASIC versions. First, we enter the matrix [A] in the form of:

>> A = [1,2;3,4]

When the ENTER key is pressed, the displayed result is:

© 2001 by CRC Press LLC

A=

1

2

3

4

Notice that the elements of [A] should be entered row by row. While the rows

are separated by ;, in each row elements are separated by comma. After the print

out of the above results, >> sign will again appear. To eliminate the unnecessary line

space (between A = and the first row 1 2), the statement format compact can be entered

as follows (the phrase “pressing ENTER key” will be omitted from now on):

>> format compact, B = [5,6;7,8]

B=

5

6

7

8

Notice that comma is used to separate the statements. To demonstrate matrix subtraction and addition, we can have:

>> A-B

ans =

–4

–4

–4

–4

>> A + B

ans =

6

8

10

12

To apply MATLAB for transposition and multiplication of matrices, we can have:

>> C = [1,2,3;4,5,6]

C=

1

2

3

4

5

6

>> C'

ans =

1

4

2

5

3

6

© 2001 by CRC Press LLC

>> D = [1,2,3;4,5,6]; E = [1,2;2,3;3,4]; P = D*E

P=

14

20

32

47

Notice that MATLAB uses ' (single quote) in place of the superscripted symbol

T for transposition. When ; (semi-colon) follows a statement such as the D statement,

the results will not be displayed. As in FORTRAN and QuickBASIC, * is the

multiplication operator as is used in P = D*E, here involving three matrices not three

single variables. More examples of MATLAB applications including plotting will

ensue. To terminate the MATLAB operation, simply enter quit and then the

RETURN key.

MATHEMATICA APPLICATIONS

To commence the service of Mathematica from Windows setup, simply point

the mouse to it and double click the left button. The Input display bar will appear

on screen, applications of Mathematica can start by entering commands from

keyboard and then press the Shift and Enter keys. To terminate the Mathematica

application, enter Quit[] from keyboard and then press the Shift and Enter keys.

Mathematica commands, statements, and functions are gradually introduced

and applied in increasing degree of difficulty. Its graphic capabilities are also utilized

in presentation of the computed results.

For matrix operations, Mathematica can compute the sum and difference of

two matrices of same order in symbolic forms, such as in the following cases of

involving two matrices, A and B, both of order 2 by 2:

In[1]: = A = {{1,2},{3,4}}; MatrixForm[A]

Out[1]//MatrixForm =

1

2

3

4

In[1]: = is shown on screen by Mathematica while user types in A =

{{1,2},{3,4}}; MatrixForm[A]. Notice that braces are used to enclose the elements

in each row of a matrix, the elments in a same row are separated by commas, and

the rows are also separated by commas. MatrixForm demands that the matrix be

printed in a matrix form. Out[1]//MatrixForm = and the rest are response of Mathematica.

In[2]: = B = {{5,6},{7,8}}; MatrixForm[B]

Out[2]//MatrixForm =

5

7

6

8

© 2001 by CRC Press LLC

In[3]: = MatrixForm[A + B]

Out[3]//MatrixForm =

6

8

10

12

In[4]: = Dif = A-B; MatrixForm[Dif]

Out[4]//MatrixForm =

–4

–4

–4

–4

In[3] and In[4] illustrate how matrices are to be added and subtracted, respectively. Notice that one can either use A + B directly, or, create a variable Dif to

handle the sum and difference matrices.

Also, Mathematica has a function called Transpose for transposition of a

matrix. Let us reuse the matrix A to demonstrate its application:

In[5]: = AT = Transpose[A]; MatrixForm[AT]

Out[5]//MatrixForm =

1

3

2

4

1.3 SOLUTION OF MATRIX EQUATION

Matrix notation offers the convenience of organizing mathematical expression in an

orderly manner and in a form which can be directly followed in coding it into

programming languages, particularly in the case of repetitive computation involving

the looping process. The most notable situation is in the case of solving a system

of linear algebraic equation. For example, if we need to determine a linear equation

y = a1 + a2x which geometrically represents a straight line and it is required to pass

through two specified points (x1,y1) and (x2,y2). To find the values of the coefficients

a1 and a2 in the y equation, two equations can be obtained by substituting the two

given points as:

(1)a1 + (x1 )a 2 = y1

(1)

(1)a1 + (x 2 )a 2 = y2

(2)

and

© 2001 by CRC Press LLC

To facilitate programming, it is advantageous to write the above equations in

matrix form as:

[C]{A} = {Y}

(3)

where:

?1

[C] = ?1

?

x1 ?

? a1 ?

? y1 ?

, {A} = ? ?, and {Y} = ? ?

?

x2 ?

?a 2 ?

?y 2 ?

(4)

The matrix equation 3 in this case is of small order, that is an order of 2. For

small systems, Cramer’s Rule can be conveniently applied which allows the unknown

vector {A} to be obtained by the formula:

[

{A} = [c1 ] [c2 ]

]

T

[C ]

(5)

Equation 5 involves the calculation of three determinants, i.e., , [C1], [C2], and

[C] where [C1] and [C2] are matrices derived from the matrix [C] when the first

and second columns of [C] are replaced by {Y}, respectively. If we denote the

elements of a general matrix [C] of order 2 by cij for i,j = 1,2, the determinant of

[C] by definition is:

[C] = c11c22 ? c12c21

(6)

The general definition of the determinant of a matrix [M] of order N and whose

elements are denoted as mij for i,j = 1,2,…,N is to add all possible product of N

elements selected one from each row but from different column. There are N! such

products and each product carries a positive or negative sign depending on whether

even or odd number of exchanges are necessary for rearranging the N subscripts in

increasing order. For example, in Equation 6, c11 is selected from the first row and

first column of [C] and only c22 can be selected and multiplied by it while the other

possible product is to select c12 from the second row and first column of [C] and

that leaves only c21 from the second row and first column of [C] available as a factor

of the second product. In order to arrange the two subscripts in non-decreasing order,

one exchange is needed and hence the product c12c21 carries a minus sign. We shall

explain this sign convention further when a matrix of order 3 is discussed. However,

it should be evident here that a matrix whose order is large the task of calculating

its determinant would certainly need help from computer. This will be the a topic

discussed in Section 1.5.

Let us demonstrate the application of Cramer’s Rule by having a numerical case.

If the two given points to be passed by the straight line y = a1 + a2x are (x1,y1) =

(1,2) and (x2,y2) = (3,4). Then we can have:

© 2001 by CRC Press LLC

[C ] = 1

1

x1 1

=

x2 1

1

= 1? 3 ?1?1 = 3 ?1 = 2

3

y1

x1 2

=

x2 4

1

= 2 ? 3 ?1? 4 = 6 ? 4 = 2

3

y1 1

=

y2 1

2

= 1? 4 ? 2 ?1 = 4 ? 2 = 2

4

[c ] = y

1

2

and

1

[C ] = 1

2

Consequently, according to Equation 5 we can find the coefficients in the

straight-line equation to be:

a1 = [C1 ]

[C ] = 2 2 = 1

and a 2 = [C2 ] [C] = 2 2 = 1

Hence, the line passing through the points (1,2) and (3,4) is y = a1 + a2x = 1 + x.

Application of Cramer’s Rule can be extended for solving three unknowns from

three linear algebraic equations. Consider the case of finding a plane which passes

three points (xi,yi) for i = 1 to 3. The equation of that plane can first be written as

z = a1 + a2x + a3y. Similar to the derivation of Equation 3, here we substitute the

three given points into the z equation and obtain:

(1)a1 + (x1 )a 2 + (y1 )a 3 = z1

(7)

(1)a1 + (x 2 )a 2 + (y2 )a 3 = z2

(8)

(1)a1 + (x3 )a 2 + (y3 )a 3 = z3

(9)

and

Again, the above three equations can be written in matrix form as:

[C]{A} = {Z}

(10)

where the matrix [C] and the vector {A} previously defined in Equation 4 need to

be reexpanded and redefined as:

1

[C ] = 1

1

© 2001 by CRC Press LLC

x1

x2

x3

y1

a1

? z1 ?

? ?

y 2 , {A} = a 2 , and {Z} = ?z 2 ?

??z3 ??

y3

a3

(11)

And, the Cramer’s Rule for solving Equation 10 can be expressed as:

[

{A} = [C1 ][C2 ] [C3 ]

]

T

[C ]

(12)

where [Ci] for i = 1 to 3 for matrices formed by replacing the ith column of the

matrix [C] by the vector {Z}, respectively. Now, we need the calculation of the

determinant of matrices of order 3. If we denote the element in a matrix [M] as mij

for i,j = 1 to 3, the determinant of [M] can be calculated as:

[M] = m11m 22 m33 ? m11m 23m32 + m12 m 23m31

? m12 m 21m 33 + m13m 21m 32 ? m13m 22 m 31

(13)

To give a numerical example, let us consider a plane passing the three points,

(x1,y1,z1) = (1,2,3), (x2,y2,z2) = (–1,0,1), and (x3,y3,z3) = (–4,–2,0). We can then have:

1

[C ] = 1

1

z1

[C1 ] = z2

z3

1

C

=

[ 2] 1

1

x1

x2

x3

y1 1

y2 = 1

y3 1

1

?1

?4

2

0 = ?2

?2

x1

x2

x3

y1 3

y2 = 1

y3 0

1

?1

?4

2

0 =0

?2

z1

z2

z3

y1 1

y2 = 1

y3 1

3

1

0

2

0 =2

?2

1

?1

?4

3

1 = ?4

?2

and

[C ]

3

1

=1

1

x1

x2

x

z1 1

z2 = 1

z 1

According to Equation 13, we find a1 = [C1]/[C] = 0/(–2) = 0, a2 = [C2]/[C] =

2/(–2) = –1, and a3 = [C3]/[C] = –4/(–2) = 2. Thus, the required plane equation is

z = a1 + a2x + a3y = -x + 2y.

QUICKBASIC VERSION

OF THE PROGRAM

CRAMERR

A computer program called CramerR has been developed as a reviewing exercise in programming to solve a matrix equation of order 3 by application of Cramer

© 2001 by CRC Press LLC

Rule and the definition of determinant of a 3 by 3 square matrix according to

Equations 12 and 13, respectively. First, a subroutine called Determ3 is created

explicitly following Equation 13 as listed below:

To interactively enter the elements of the coefficient matrix [C] and also the

elements of the right-hand-side vector {Z} in Equation 12 and to solve for {A}, the

program CramerR can be arranged as:

© 2001 by CRC Press LLC

1.4 PROGRAM GAUSS

Program Gauss is designed for solving N unknowns from N simultaneous, linear

algebraic equations by the Gaussian Elimination method. In matrix notation, the

problem can be described as to solve a vector {X} from the matrix equation:

[C]{X} = {V}

(1)

where [C] is an NxN coefficient matrix and {V} is a Nx1 constant vector, and both

are prescribed. For example, let us consider the following system:

9x1 + x 2 + x3 = 10

(2)

3x1 + 6x 2 + x3 = 14

(3)

2 x1 + 2 x 2 + 3x3 = 3

(4)

If the above three equations are expressed in matrix form as Equation 1, then:

?9

[C] = ??3

??2

1

6

2

1?

1??,

3??

?10 ?

{V} = ??14??,

?? 3 ??

(5,6)

and

? x1 ?

X

=

{ } ??x 2 ?? = x1 x 2 x3

?? x3 ??

[

]

T

(7)

where T designates the transpose of a matrix.

GAUSSIAN ELIMINATION METHOD

A systematic procedure named after Gauss can be employed for solving x1, x2,

and x3 from the above equations. It consists of first dividing Equation 28 by the

leading coefficient, 9, to obtain:

x1 +

1

1

10

x + x =

9 2 9 3 9

(8)

This step is called normalization of the first equation of the system (1). The next

step is to eliminate x1 term from the other (in this case, two) equations. To do that,

we multiply Equation 8 by the coefficients associated with x1 in Equations 3 and 4,

respectively, to obtain:

© 2001 by CRC Press LLC

1

1

10

3x1 + x 2 + x3 =

3

3

3

(9)

2

2

20

x 2 + x3 =

9

9

9

(10)

and

2 x1 +

If we subtract Equation 9 from Equation 3, and subtract Equation 10 from Equation 4, the x1 terms are eliminated. The resulting equations are, respectively:

17

2

32

x 2 + x3 =

3

3

3

(11)

16

25

7

x 2 + x3 =

9

9

9

(12)

and

This completes the first elimination step. The next normalization is applied to

Equation 11, and then the x2 term is to be eliminated from Equation 12. The resulting

equations are:

2

32

(13)

x 2 + x3 =

17

17

and

393

393

(14)

x =?

153 3

153

The last normalization of Equation 14 then gives:

x3 = ?1

(15)

Equations 8, 13, and 15 can be organized in matrix form as:

?1

{V} = ??0

??0

19

1

0

1 9 ? ? x1 ? ? 10 9 ?

? ?

2 17?? ?x 2 ? = ??32 17??

1 ?? ?? x3 ?? ?? ?1 ??

(16)

The coefficient matrix is now a so-called upper triangular matrix since all

elements below the main diagonal are equal to zero.

As x3 is already obtained in Equation 15, the other two unknowns, x2 and x3,

can be obtained by a sequential backward-substitution process. First, Equation 13

can be used to obtain:

x2 =

© 2001 by CRC Press LLC

32 2

32 2

32 + 2

=2

? x3 =

? ( ?1) =

17 17

17 17

17

Once, both x2 and x3 have been calculated, x1 can be obtained from Equation 8 as:

x1 =

10 1

1

10 1

1

10 ? 2 + 1

=1

? x 2 ? x3 =

? (2) ? ( ?1) =

9 9

9

9 9

9

9

To derive a general algorithm for the Gaussian elimination method, let us denote

the elements in [C], {X}, and {V} as ci,j, xi, and vi, respectively. Then the normalization of the first equation can be expressed as:

(c )

1, j new

( ) (c )

= c1, j

old

1,1 old

(17)

and

(v1 )new = (v1 ) old (c1,1 )old

(18)

Equation 17 is to be used for calculating the new coefficient associated with xj

in the first, normalized equation. So, j should be ranged from 2 to N which is the

number of unknowns (equal to 3 in the sample case). The subscripts old and new

are added to indicate the values before and after normalization, respectively. Such

designation is particularly helpful if no separate storage in computer are assigned

for [C] for the values of its elements. Notice that (c1,1)new = 1 is not calculated.

Preserving this diagonal element enables the determinant of [C] to be computed.

(See the topic on matrix inversion and determinant.)

The formulas for the elimination of x1 terms from the second equation are:

(c )

2 , j new

( ) ? (c ) (c )

= c2, j

old

2 ,1 old

1, j old

(19)

for j = 2,3,…,N (there is no need to include j = 1) and

(v2 )new = (v2 )old ? (c2,1 )old (v1 )old

(20)

By changing the subscript 2 in Equations 19 and 20, x1 term in the third equation

can be eliminated. In other words, the general formulas for elimination of x1 terms

from all equation other than the first equation are, for k = 2,3,…,N

(c )

k , j new

for j = 2,3,…,N

© 2001 by CRC Press LLC

( ) ? (c ) (c )

= ck,j

old

k ,l old

1, j old

(21)

(vk )new = (vk )old ? (ck,1 )old (v1 )old

(22)

Instead of normalizing the first equation, we can generalize Equations 17 and

18 for normalization of the ith equation, for i = 1,2,…,N to the expressions:

(c )

i , j new

( ) (c )

= ci. j

old

i ,i old

(23)

for j = i + 1,i + 2,…,N and

(vi )new = (vi ) old (ci,i )old

(24)

Note that (ci,i)new should be equal to 1 but no need to calculate since it is not

involved in later calculation for finding {X}.

Similarly, elimination of xi term from kth equation for k = i + 1,i + 2,…,N

consists of using the general formula:

(c )

k , j new

( ) ? (c ) (c )

= ck,j

k ,i old

old

i , j old

(25)

for j = i + 1,i + 2,…,N and

(vk )new = (vk )old ? (ck,i )old (vi )old

(26)

Backward substitution for finding xi involves the calculation of:

N

x i = vi ?

?c

x

i, j j

(27)

j= i +1

for i = N–1,N–2,…,2,1. Note that xN is already found equal to vN after the Nth

normalization.

Program Gauss listed below in both QuickBASIC and FORTRAN languages

is developed for interactive specification of the number of unknowns, N, and the

values of the elements of [C] and {V}. It proceeds to solve for {X} and prints out

the computed values. Sample applications of both languages are provided immediately following the listed programs.

A subroutine Gauss.Sub is also made available for the frequent need in the

other computer programs which require the solution of matrix equations.

© 2001 by CRC Press LLC

QUICKBASIC VERSION

Sample Application

© 2001 by CRC Press LLC

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Application

GAUSS-JORDAN METHOD

One slight modification of the elimination step will make the backward substitution steps completely unnecessary. That is, during the elimination of the xi terms

from the linear algebraic equations except the ith one, Equations 25 and 26 should

be applied for k equal to 1 through N and excluding k = i. For example, the x3 terms

should be eliminated from the first, second, fourth through Nth equations. In this

manner, after the Nth normalization, [C] becomes an identity matrix and {V} will

have the elements of the required solution {X}. This modified method is called

Gauss-Jordan method.

A subroutine called GauJor is made available based on the above argument. In

this subroutine, a block of statements are also added to include the consideration of

the pivoting technique which is required if ci,i = 0. The normalization steps,

Equations 49 and 50, cannot be implemented if ci,i is equal to zero. For such a

situation, a search for a nonzero ci,k is necessary for i = k + 1,k + 2,…,N. That is,

to find in the kth column of [C] and below the kth row a nonzero element. Once

this nonzero ci,k is found, then we can then interchange the ith and kth rows of [C]

and {V} to allow the normalization steps to be implemented; if no nonzero ci,k can

be found then [C] is singular because the determinant of [C] is equal to zero! This

can be explained by the fact that when ck,k = 0 and no pivoting is possible and the

determinant D of [C] can be calculated by the formula:

N

D = c1,1c2,2 … c k ,k … c N,N =

?c

k =1

where indicates a product of all listed factors.

© 2001 by CRC Press LLC

k ,k

(28)

A subroutine has been written based on the Gauss-Jordan method and called

GauJor.Sub. Both QuickBASIC and FORTRAN versions are made available and

they are listed below.

QUICKBASIC VERSION

© 2001 by CRC Press LLC

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Applications

The same problem previously solved by the program Gauss has been used again

but solved by application of subroutine GauJor. The results obtained with the QuickBASIC and FORTRAN versions are listed, in that order, below:

MATLAB APPLICATIONS

For solving the vector {X} from the matrix equation [C]{X} = {R} when both

the coefficient matrix [C] and the right-hand side vector {R} are specified, MATLAB

simply requires [C] and {R} to be interactively inputted and then uses a statement

X = C\R to obtain the solution vector {X} by multiplying the vector {R} on the left

of the inverse of [C] or dividing {R} on the left by [C]. More details are discussed

in the program MatxAlgb. Here, for providing more examples in MATLAB applications, a m file called GauJor.m is presented below as a companion of the FORTRAN and QuickBASIC versions:

© 2001 by CRC Press LLC

This file GauJor.m should then be added into MATLAB. As an example of

interactive application of this m file, the sample problem used in the FORTRAN

and QuickBASIC versions is again solved by specifying the coefficient matrix [C]

and the right hand side vector {R} to obtain the resulting display as follows:

The results of the vector {X} and determinant D for the coefficient matrix [C]

are same as obtained before.

MATHEMATICA APPLICATIONS

For solving a system of linear algebraic equations which has been arranged in

matrix form as [A]{X} = {R}, Mathematica’s function LinearSolve can be applied

© 2001 by CRC Press LLC

to solve for {X} when the coefficient matrix [A] and the right-hand side vector {R}

are both provided. The following is an example of interactive application:

In[1]: = A = {{3,6,14},{6,14,36},{14,36,98}}

Out[1]: =

{{3, 6, 14}, {6, 14, 36}, {14, 36, 98}}

In[2]: = MatrixForm[A]

Out[2]//MatrixForm: =

3

6

14

6

14

36

14

36

98

In[3]: = R = {9,20,48}

Out[3]: =

{9, 20, 48}

In[4]: = LinearSolve[A,R]

Out[4]: =

{–9,13,–3}

Output[2] and Output[1] demonstrate the difference in display of matrix [A]

when MatrixeForm is requested, or, not requested, respectively. It shows that without

requesting of MatrixForm, some screen space saving can be gained. Output[4] gives

the solution {X} = [–9 13 –3]T for the matrix equation [A]{X} = {R} where the

coefficient matix [A] and vector {R} are provided by Input[1] and Input[3], respectively.

1.5 MATRIX INVERSION, DETERMINANT,

AND PROGRAM MatxInvD

Given a square matrix [C] of order N, its inverse as [C]–1 of the same order is defined

by the equation:

[C][C]?1 = [C]?1[C] = [I]

(1)

where [I] is an identity matrix having elements equal to one along its main diagonal

and equal to zero elsewhere. That is:

.

?1 0

?0 1 0

[I] = ?

? . . .

?

.

?0 0

© 2001 by CRC Press LLC

.

.

.

.

.

.

.

.

.

.

.

.

0?

0?

?

?

?

1?

(2)

To find [C]–1, let cij and dij be the elements at the ith row and jth column of the

matrices [C] and [C]–1, respectively. Both i and j range from 1 to N. Furthermore,

let {Dj} and {Ij} be the jth column of the matrices [C]–1 and [I], respectively. It is

easy to observe that {Ij} has elements all equal to zero except the one in the jth row

which is equal to unity. Also,

{D } = [d d

j

lj 2 j

… d Nj

]

T

(3)

and

[C]?1 = [D1D2 … DN ]

T

(4)

Based on the rules of matrix multiplication, Equation 1 can be interpreted as

[C]{D1} = {I1}, [C]{D2} = {I2}, …, and [C]{DN} = {IN}. This indicates that program

Gauss can be successively employed N times by using the same coefficient matrix

[C] and the vectors {Ii} to find the vectors {Di} for i = 1,2,…,N. Program MatxInvD

is developed with this concept by modifying the program Gauss. It is listed below

along with a sample interactive run.

QUICKBASIC VERSION

© 2001 by CRC Press LLC

Sample Application

FORTRAN VERSION

© 2001 by CRC Press LLC

© 2001 by CRC Press LLC

Sample Applications

MATLAB APPLICATION

MATLAB offers very simple matrix operations. For example, matrix inversion

can be implemented as:

To check if the obtained inversion indeed satisfies the equation [A}[A]–1 = [I]

where [I] is the identity matrix, we enter:

Once [A]–1 becomes available, we can solve the vector {X} in the matrix equation

[A]{X} = {R} if {R} is prescribed, namely {X} = [A]–1{R}. For example, may enter

a {R} vector and find {X} such as:

© 2001 by CRC Press LLC

MATHEMATICA APPLICATIONS

Mathematica has a function called Inverse for inversion of a matrix. Let us

reuse the matrix A that we have entered in earlier examples and continue to demonstrate the application of Inverse:

In[1]: = A = {{1,2},{3,4}}; MatrixForm[A]

Out[1]//MatrixForm =

1

2

3

4

In[2]: = B = {{5,6},{7,8}}; MatrixForm[B]

Out[2]//MatrixForm =

5

6

7

8

In[3]: = MatrixForm[A + B]

Out[3]//MatrixForm =

6

8

10

12

In[4]: = Dif = A-B; MatrixForm[Dif]

Out[4]//MatrixForm =

–4

–4

–4

–4

In[5]: = AT = Transpose[A]; MatrixForm[AT]

Out[5]//MatrixForm =

1

3

2

4

In[6]: = Ainv = Inverse[A]; MatrixForm[Ainv]

Out[6]//MatrixForm =

–2

1

3 ? 1?

?

2 ? 2?

© 2001 by CRC Press LLC

To verify whether or not the inverse matrix Ainv obtained in Output[6] indeed

satisfies the equations [A][A]–1 = [I] which is the identity matrix, we apply Mathematica for matrix multiplication:

In[7]: = Iden = A.Ainv; MatrixForm[Iden]

Out[7]//MatrixForm =

1

0

0

1

A dot is to separate the two matrices A and Ainv which is to be multiplied in that

order. Output[7] proves that the computed matrix, Ainv, is the inverse of A! It should

be noted that D and I are two reserved variables in Mathematica for the determinant

of a matrix and the identity matrix. In their places, here Dif and Iden are adopted,

respectively. For further testing, we show that [A][A]T is a symmetric matrix:

In[8]: = S = A.AT; MatrixForm[S]

Out[8]//MatrixForm =

5

11

11

25

And, the unknown vector {X} in the matrix equation [A]{X} = {R} can be

solved easily if {R} is given and [A]–1 are available:

In[9]: = R = {13,31}; X = Ainv.R

Out[9] = {5, 4}

The solution of x1 = 5 and x2 = 4 do satisfy the equations x1 + 2x2 = 13 and 3x1

+ 4x2 = 31.

TRANSFORMATION

OF

COORDINATE SYSTEMS, ROTATION,

AND

ANIMATION

Matrix algebra can be effectively applied for transformation of coordinate systems. When the cartesian coordinate system, x-y-z, is rotated by an angle z about

the z-axis to arrive at the system x-y-z as shown in Figure 2, where z and z axes

coincide and directed outward normal to the plane of paper, the new coordinates of

a typical point P whose coordinates are (xP,yP,zP) can be easily obtained as follows:

x ?P = OP cos(? P ? ? z ) = (OP cos ? P ) cos ? z + (OP sin ? P ) sin ? z

= x P cos ? z + y p sin ? z

y ?P = OP sin(? P ? ? z ) = (OP sin ? P ) cos ? z ? (OP cos ? P ) sin ? z

= x p sin ? z + y p sin ? z

© 2001 by CRC Press LLC

FIGURE 2. The cartesian coordinate system, x-y-z, is rotated by an angle z about the zaxis to arrive at the system x-y-z.

and

z?P = z P

In matrix notation, we may define {P} = [xP yP zP]T and {P'} = [xP' yP' zP']T and

write the above equations as {P'} = [Tz]{P} where the transformation matrix for a

rotation of z-axis by z is:

? cos ?z

[Tz ] = ??? sin ?z

?? 0

sin ?z

cos ?z

0

0?

0 ??

1 ??

(5)

In a similar manner, it can be shown that the transformation matrices for rotating

about the x- and y-axes by angles x and y, respectively, are:

?1

[Tx ] = ??0

??0

0

cos ?x

? sin ?x

0 ?

?

sin ?x ?

cos ?x ??

(6)

? sin ? y ?

?

0 ?

cos ? y ??

(7)

and

?cos ? y

?

Ty = ? 0

? sin ?

y

?

[ ]

© 2001 by CRC Press LLC

0

1

0

FIGURE 3. Point P whose coordinates are (xP,yP,zP) is rotated to the point P' by a rotation

of z.

It is interesting to note that if a point P whose coordinates are (xP,yP,zP) is rotated

to the point P' by a rotation of z as shown in Figure 3, the coordinates of P' can

be easily obtained by the formula {P'} = [Rz]{P} where [Rz] = [Tz]T. If the rotation

is by an angle x or y, then {P'} = [Rx]{P} or {P'} = [Ry]{P} where [Rx] = [Tx]T

and [Ry] = [Ty]T.

Having discussed about transformations and rotations of coordinate systems,

we are ready to utilize the derived formulas to demonstrate the concept of animation. Motion can be simulated by first generating a series of rotated views of

a three-dimensional object, and showing them one at a time. By erasing each

displayed view and then showing the next one at an adequate speed, a smooth

motion of the object is achievable to produce the desired animation. Program

Animate1.m is developed to demonstrate this concept of animation by using a

4 2 3 brick and rotating it about the x-axis by an angle of 25° and then

rotating about the y-axis as many revolutions as desired. The front side of the

block (x-y plane) is marked with a character F, and the right side (y-z plane) is

marked with a character R, and the top side (x-z plane) is marked with a character

T for helping the viewer to have a better three-dimensional perspective of the

rotated brick (Figure 4). The x-rotation prior to y-rotation is needed to tilt the top

side of the brick toward the front. The speed of animation is controlled by a

parameter Damping. This parameter and the desired number of y-revolutions,

Ncycle, are both to be interactively specified by the viewer (Figure 5).

© 2001 by CRC Press LLC

FIGURE 4. The characters F, R, and T help the viewer to have a better three-dimensional

perspective of the rotated brick.

FIGURE 5. The speed of animation is controlled by a parameter Damping. This parameter and

the desired number of y-revolutions, Ncycle, are both to be interactively specified by the viewer.

© 2001 by CRC Press LLC

FUNCTION ANIMATE1(NCYCLE,DAMPING)

Notice that the coordinates for the corners of the brick are defined in arrays xb,

yb, and zb. The coordinates of the points to be connected by linear segments for

drawing the characters F, R, ant T are defined in arrays xf, yf, and zf, and xr, yr,

and zr, and xt, yt, and zt, respectively.

The equations in deriving [Rx] ( = [Tx]T) and [Ry] ( = [Ty]T) are applied for xand y- rotations in the above program. Angle increments of 5 and 10° are arranged

for the x- and y-rotations, respectively. The rotated views are plotted using the new

coordinates of the points, (xbn,ybn,zbn), (xfn,yfn,zfn), etc. Not all of these new

arrays but only those needed in subsequent plot are calculated in this m file.

MATLAB command clg is used to erase the graphic window before a new

rotated view the brick is displayed. The speed of animation is retarded by the “hold”

loops in both x- and y-rotations involving the interactively entered value of the

parameter Damping. The MATLAB command pause enables Figure 4 to be read

and requires the viewer to press any key on the keyboard to commence the animation.

Notice that a statement begins with a % character making that a comment statement,

and that % can also be utilized for spacing purpose.

The xs and ys arrays allow the graphic window to be scaled by plotting them

and then held (by command hold) so that all subsequent plots are using the same

© 2001 by CRC Press LLC

FIGURE 6. Animation of a rotating brick.

scales in both x- and y-directions. The values in xs and ys arrays also control where

to properly place the texts in Figure 4 as indicated in the text statements.

QUICKBASIC VERSION

A QuickBASIC version of the program Animate1.m called Animate1.QB also

is provided. It uses commands GET and PUT to animate the rotation of the 4 3

2 brick. More features have been added to show the three principal views of the

brick and also the rotated view at the northeast corner of screen, as illustrated in

Figure 6.

The window-viewport transformation of the rotated brick for displaying on the

screen is implemented through the functions FNTX and FNTY. The actual ranges

of the x and y measurements of the points used for drawing the brick are described

by the values of V1 and V2, and V3 and V4, respectively. These ranges are mapped

onto the screen matching the ranges of W1 and W2, and W4 and W3, respectively.

The rotated views of the brick are stored in arrays S1 through S10 using the

GET command. Animation retrieves these views by application of the PUT command. Presently, animation is set for 10 y-swings (Ncycle = 10 in the program

Animate1.m, arranged in Line 600). The parameter Damping described in the

program Animate1.m here is set equal to 1500 (in Line 695).

© 2001 by CRC Press LLC

© 2001 by CRC Press LLC

1.6 PROBLEMS

MATRIX ALGEBRA

1. Calculate the product [A][B][C] by (1) finding [T] = [A][B] and then

[T][C], and (2) finding [T] = [B][C] and then [A][T] where:

?1

[A] = ?

?4

2

5

3?

6??

?6

[B] = ??4

??2

5?

3??

1??

? ?1

??3

[C ] = ?

?2 ?

?4??

2. Calculate [A][B] of the two matrices given above and then take the

transpose of product matrix. Is it equal to the product of [B]T[A]T?

3. Are ([A][B][C])T and the product [C]T[B]T[A]T identical to each other?

© 2001 by CRC Press LLC

4. Apply the QuickBASIC and FORTRAN versions of the program MatxAlgb to verify the results of Problems 1, 2, and 3.

5. Repeat Problem 4 but use MATLAB.

6. Apply the program MatxInvD to find [C]–1 of the matrix [C] given in

Problem 1 and also to ([C]T)–1. Is ([C]–1)T equal to ([C]T)–1?

7. Repeat Problem 6 but use MATLAB.

8. For statistical analysis of a set of N given data X1, X2, …, XN, it is often

necessary to calculate the mean, m, and standard deviation, 5, by use of

the formulas:

m=

1

(X + X2 + … + X N )

N 1

and

[

]

2

2

2

1

? = ?? ( X1 ? m ) + ( X 2 ? m ) + … + ( X N ? m ) ??

?

?N

0.5

Use indicial notation to express the above two equations and then develop

a subroutine meanSD(X,N,RM,SD) for taking the N values of X to

compute the real value of mean, RM, and standard deviation, SD.

9. Express the ith term in the following series in indicial notation and then

write an interactive program SinePgrm allowing input of the x value to

calculate sin(x) by terminating the series when additional term contributes

less than 0.001% of the partial sum of series in magnitude:

Sin x =

x1 x3 x 5

? + ?…

1! 3! 5!

Notice that Sin(x) is an odd function so the series contains only terms of

odd powers of x and the series carries alternating signs. Compare the

result of the program SinePgrm with those obtained by application of the

library function Sin available in FORTRAN and QuickBASIC.

10. Same as Problem 9, but for the cosine series:

Cos x = 1 ?

x2 x 4 x6

+

?

+…

2! 4! 6!

Notice that Cos(x) is an even function so the series contains only terms

of even powers of x and the series also carries alternating signs.

11. Repeat Problem 4 but use Mathematica.

12. Repeat Problem 6 but use Mathematica.

© 2001 by CRC Press LLC

GAUSS

1. Run the program GAUSS to solve the problem:

?1

?4

?

?? 7

2

5

8

3 ? ? x1 ? ? 2 ?

? ?

6 ?? ?x 2 ? = ?? 8 ??

10 ?? ?? x3 ?? ??14??

2. Run the program GAUSS to solve the problem:

?0

?4

?

?? 7

2

5

8

3? ? x1 ? ??1?

? ?

6?? ?x 2 ? = ?? 8 ??

9?? ?? x3 ?? ??14 ??

What kind of problem do you encounter? “Divided by zero” is the message! This happens because the coefficient associated with x1 in the first

equation is equal to zero and the normalization in the program GAUSS

cannot be implemented. In this case, the order of the given equations

needs to be interchanged. That is to put the second equation on top or

find below the first equation an equation which has a coefficient associated

with x1 not equal to zero and is to be interchanged with the first equation.

This procedure is called “pivoting.” Subroutine GauJor has such a feature

incorporated, apply it for solving the given matrix equation.

3. Modify the program GAUSS by following the Gauss-Jordan elimination

procedure and excluding the back-substitution steps. Name this new program GauJor and test it by solving the matrix equations given in Problems

1 and 2.

4. Show all details of the normalization, elimination, and backward substitution steps involved in solving the following equations by application of

Gaussian Elimination method:

4x1 + 2x2 – 3x3 = 8

5x1 – 3x2 + 7x3 = 26

–x1 + 9x2 – 8x3 = –10

5. Present every normalization and elimination steps involved in solving the

following system of linear algebraic equations by the Gaussian Elimination Method:

5x1 – 2x2 + 2x3 = 9, –2x1 + 7x2 – 2x3 = 9, and 2x1 – 2x2 + 9x3 = 41

© 2001 by CRC Press LLC

6. Apply the Gauss-Jordan elimination method to solve for x1, x2, and x3

from the following equations:

?0

?2

?

??4

?1? ? x1 ? ?1?

? ?

3 ?? ?x 2 ? = ??1??

7 ?? ?? x3 ?? ??1??

1

9

24

Show every normalization, elimination, and pivoting (if necessary) steps

of your calculation.

7. Solve the matrix equation [A]{X} = {C} by Gauss-Jordan method

where:

?3

?2

?

??4

2

5

1

1 ? ? x1 ? ??2 ?

? ?

?1?? ?x 2 ? = ?? ?3??

7 ?? ?? x3 ?? ?? 3 ??

Show every interchange of rows (if you are required to do pivoting before

normalization), normalization, and elimination steps by indicating the

changes in [A] and {C}.

8. Apply the program GauJor to solve Problem 7.

9. Present every normalization and elimination steps involved in solving the

following system of linear algebraic equations by the Gauss-Jordan Elimination Method:

5x1 – 2x2 + x3 = 4

–2x1 + 7x2 – 2x3 = 9

x1 – 2x2 + 9x3 = 40

10.

11.

12.

13.

14.

Apply the program Gauss to solve Problem 9 described above.

Use MATLAB to solve the matrix equation given in Problem 7.

Use MATLAB to solve the matrix equation given in Problem 9.

Use Mathematica to solve the matrix equation given in Problem 7.

Use Mathematica to solve the matrix equation given in Problem 9.

MATRIX INVERSION

1. Run the program MatxInvD for finding the inverse of the matrix:

?3

[A] = ??0

??2

© 2001 by CRC Press LLC

0

5

0

2?

0 ??

3??

2. Write a program Invert3 which inverts a given 3 ? 3 matrix [A] by using

the cofactor method. A subroutine COFAC should be developed for calculating the cofactor of the element at Ith row and Jth column of [A] in

term of the elements of [A] and the user-specified values of I and J. Let

the inverse of [A] be designated as [AI] and the determinant of [A] be

designated as D. Apply the developed program Invert3 to generate all

elements of [AI] by calling the subroutine COFAC and by using D.

3. Write a QuickBASIC or FORTRAN program MatxSorD which will

perform the addition and subtraction of two matrices of same order.

4. Write a QuickBASIC or FORTRAN program MxTransp which will

perform the transposition of a given matrix.

5. Translate the FORTRAN subroutine MatxMtpy into a MATLAB m file

so that by entering the matrices [A] and [B] of order L by M and M by

N, respectively, it will produce a product matrix [P] of order L by N.

6. Enter MATLAB commands interactively first a square matrix [A] and

then calculate its trace.

7. Use MATLAB commands to first define the elements in its upper right

corner including the diagonal, and then use the symmetric properties to

define those in the lower left corner.

8. Convert either QuickBasic or FORTRAN version of the program MatxInvD into a MATLAB function file MatxInvD.m with a leading statement

function [Cinv,D] = MatxInvD(C,N)

9. Apply the program MatxInvD to invert the matrix:

?1

[A] = ??5

??8

3

6

9

4?

7 ??

10 ??

Verify the answer by using Equation 1.

10. Repeat Problem 9 but by MATLAB operation.

11. Apply the program MatxInvD to invert the matrix:

??9

[A] = ?? ?3

???6

?1

?4

?7

?2 ?

?5??

?8??

Verify the answer by using Equation 1.

12. Repeat Problem 11 but by MATLAB operations.

13. Derive [Rx] and verify that it is indeed equal to [Tx]T. Repeat for [Ry] and

[Rz].

14. Apply MATLAB to generate a matrix [Rz] for ?z = 45° and then to use

[Rz] to find the rotated coordinates of a point P whose coordinates before

rotation are (1,–2,5).

© 2001 by CRC Press LLC

FIGURE 7. Problem 18.

15. What will be the coordinates for the point P mentioned in Problem 14 if

the coordinate axes are rotated counterclockwise about the z-axis by 45°?

Use MATLAB to find your answer.

16. Apply MATLAB to find the location of a point whose coordinates are

(1,2,3) after three rotations in succession: (1) about y-axis by 30°, (2)

about z-axis by 45° and then (3) about x-axis by –60°.

17. Change m file Animate1.m to animate just the rotation of the front (F)

side of the 4 2 3 brick in the graphic window.

18. Write a MATLAB m file for animation of pendulum swing1 as shown in

Figure 7.

19. Write a MATLAB m file for animation of a bouncing ball1 using an

equation of y = 3e–0.1xsin(2x + 1.5708) as shown in Figure 8.

20. Write a MATLAB m file for animation of the motion of crank-piston

system as shown in Figure 9.

21. Write a MATLAB m file to animate the vibrating system of a mass

attached to a spring as shown in Figure 10.

© 2001 by CRC Press LLC

FIGURE 8. Problem 19.

22. Write a MATLAB m file to animate the motion of a cam-follower system

as shown in Figure 11.

23. Write a MATLAB m file to animate the rotary motion of a wankel cam

as shown in Figure 12.

24. Repeat Problem 9 but by Mathematica operation.

25. Repeat Problem 11 but by Mathematica operation.

26. Repeat Problem 14 but by Mathematica operation.

27. Repeat Problem 15 but by Mathematica operation.

28. Repeat Problem 16 but by Mathematica operation.

© 2001 by CRC Press LLC

FIGURE 9. Problem 20.

© 2001 by CRC Press LLC

FIGURE 10. Problem 21.

FIGURE 11. Problem 22.

© 2001 by CRC Press LLC

FIGURE 12. Problem 23.

1.7 REFERENCE

1. Y. C. Pao, “On Development of Engineering Animation Software,” in Computers in

Engineering, edited by K. Ishii, ASME Publications, New York, 1994, pp. 851–855.

© 2001 by CRC Press LLC

Matrix Algebra

and Solution

of Matrix Equations

1.1 INTRODUCTION

Computers are best suited for repetitive calculations and for organizing data into

specialized forms. In this chapter, we review the matrix and vector notation and

their manipulations and applications. Vector is a one-dimensional array of numbers

and/or characters arranged as a single column. The number of rows is called the

order of that vector. Matrix is an extension of vector when a set of numbers and/or

characters are arranged in rectangular form. If it has M rows and N column, this

matrix then is said to be of order M by N. When M = N, then we say this square

matrix is of order N (or M). It is obvious that vector is a special case of matrix when

there is only one column. Consequently, a vector is referred to as a column matrix

as opposed to the row matrix which has only one row. Braces are conventionally

used to indicate a vector such as {V} and brackets are for a matrix such as [M].

In writing a computer program, DIMENSION or DIM statements are necessary

to declare that a certain variable is a vector or a matrix. Such statements instruct

the computer to assign multiple memory spaces for keeping the values of that vector

or matrix. When we deal with a large number of different entities in a group, it is

better to arrange these entities in vector or matrix form and refer to a particular

entity by specifying where it is located in that group by pointing to the row (and

column) number(s). Such as in the case of having 100 numbers represented by the

variable names A, B, …, or by A(1) through A(100), the former requires 100 different

characters or combinations of characters and the latter certainly has the advantage

of having only one name. The A(1) through A(100) arrangement is to adopt a vector;

these numbers can also be arranged in a matrix of 10 rows and 10 columns, or 20

rows and five columns depending on the characteristics of these numbers. In the

cases of collecting the engineering data from tests of 20 samples during five different

days, then arranging these 100 data into a matrix of 20 rows and five columns will

be better than of 10 rows and 10 columns because each column contains the data

collected during a particular day.

In the ensuing sections, we shall introduce more definitions related to vector

and matrix such as transpose, inverse, and determinant, and discuss their manipulations such as addition, subtraction, and multiplication, leading to the organizing of

systems of linear algebraic equations into matrix equations and to the methods of

finding their solutions, specifically the Gaussian Elimination method. An apparent

application of the matrix equation is the transformation of the coordinate axes by a

© 2001 by CRC Press LLC

rotation about any one of the three axes. It leads to the derivation of the three basic

transformation matrices and will be elaborated in detail.

Since the interactive operations of modern personal computers are emphasized

in this textbook, how a simple three-dimensional brick can be displayed will be

discussed. As an extended application of the display monitor, the transformation of

coordinate axes will be applied to demonstrate how animation can be designed to

simulate the continuous rotation of the three-dimensional brick. In fact, any threedimensional object could be selected and its motion animated on a display screen.

Programming languages, FORTRAN, QuickBASIC, MATLAB, and Mathematica are to be initiated in this chapter and continuously expanded into higher

levels of sophistication in the later chapters to guide the readers into building a

collection of their own programs while learning the computational methods for

solving engineering problems.

1.2 MANIPULATION OF MATRICES

Two matrices [A] and [B] can be added or subtracted if they are of same order, say

M by N which means both having M rows and N columns. If the sum and difference

matrices are denoted as [S] and [D], respectively, and they are related to [A] and

[B] by the formulas [S] = [A] + [B] and [D] = [A]-[B], and if we denote the elements

in [A], [B], [D], and [S] as aij, bij, dij, and sij for i = 1 to M and j = 1 to N, respectively,

then the elements in [S] and [D] are to be calculated with the equations:

sij = a ij + b ij

(1)

d ij = a ij ? b ij

(2)

and

Equations 1 and 2 indicate that the element in the ith row and jth column of [S]

is the sum of the elements at the same location in [A] and [B], and the one in [D]

is to be calculated by subtracting the one in [B] from that in [A] at the same location.

To obtain all elements in the sum matrix [S] and the difference matrix [D], the index

i runs from 1 to M and the index j runs from 1 to N.

In the case of vector addition and subtraction, only one column is involved (N =

1). As an example of addition and subtraction of two vectors, consider the two

vectors in a two-dimensional space as shown in Figure 1, one vector {V1} is directed

from the origin of the x-y coordinate axes, point O, to the point 1 on the x-axis

which has coordinates (x1,y1) = (4,0) and the other vector {V2} is directed from the

origin O to the point 2 on the y-axis which has coordinates (x2,y2) = (0,3). One may

want to find the resultant of {R} = {V1} + {V2} which is the vector directed from

the origin to the point 3 whose coordinates are (x3,y3) = (4,3), or, one may want to

find the difference vector {D} = {V1} – {V2} which is the vector directed from the

origin O to the point 4 whose coordinates are (x4,y4) = (4,–3). In fact, the vector

{D} can be obtained by adding {V1} to the negative image of {V2}, namely {V2–}

which is a vector directed from the origin O to the point 5 whose coordinates are

(x5,y5). Mathematically, based on Equations 1 and 2, we can have:

© 2001 by CRC Press LLC

? 4 ? ?0 ?

?4?

? 4 ? ?0 ?

?4?

{R} = {V1} + {V2 } = ? ? + ? ? = ? ?

?0 ? ? 3? ? 3 ?

and

{D} = {V1} ? {V2 } = ? ? ? ? ? = ? ?

?0 ? ?3? ??3?

When Equation 1 is applied to two arbitrary two-dimensional vectors which

unlike {V1}, {V2}, and {V2–} but are not on either one of the coordinate axes, such

as {D} and {E} in Figure 1, we then have the sum vector {F} = {D} + {E} which

has components of 1 and –2 units along the x- and y-directions, respectively. Notice

that O467 forms a parallelogram in Figure 1 and the two vectors {D} and {E} are

the two adjacent sides of the parallelogram at O. To find the sum vector {F} of {D}

and {E} graphically, we simply draw a diagonal line from O to the opposite vertex

of the parallelogram — this is the well-known Law of Parallelogram.

It should be evident that to write out a vector which has a large number of rows

will take up a lot of space. If this vector can be rotated to become from one column

to one row, space saving would then be possible. This process is called transposition

as we will be leading to it by first introducing the length of a vector.

For the calculation of the length of a two-dimensional or three-dimensional vector,

such as {V1} and {V2} in Figure 1, it would be a simple matter because they are

oriented along the directions of the coordinate axes. But for the vectors such as {R}

FIGURE 1. Two vectors in a two-dimensional space.

© 2001 by CRC Press LLC

and {D} shown in Figure 1, the calculation of their lengths would need to know the

components of these vectors in the coordinate axes and then apply the Pythagorean

theorem. Since the vector {R} has components equal to rx = 4 and ry = 3 units along

the x- and y-axis, respectively, its length, here denoted with the symbol , is:

{R} = [rx2 + ry2 ] = [42 + 32 ] = 5

0.5

0.5

(3)

To facilitate the calculation of the length of a generalized vector {V} which has

N components, denoted as v1 through vN, its length is to be calculated with the

following formula obtained from extending Equation 3 from two-dimensions to Ndimensions:

{V} = [v12 + v22 + … + v2N ]

0.5

(4)

For example, a three-dimensional vector has components v1 = vx = 4, v2 = vy =

3, and v3 = vz = 12, then the length of this vector is {V} = [42 + 32 + 122]0.5 = 13.

We shall next show that Equation 4 can also be derived through the introduction of

the multiplication rule and transposition of matrices.

1.2 MULTIPLICATION OF MATRICES

A matrix [A] of order L (rows) by M (columns) and a matrix [B] of order M

by N can be multiplied in the order of [A][B] to produce a new matrix [P] of order

L by N. [A][B] is said as [A] post-multiplied by [B], or, [B] pre-multiplied by [A].

The elements in [P] denoted as pij for i = 1 to N and j = 1 to M are to be calculated

by the formula:

p ij =

?

M

k =1

a ik b kj

(5)

Equation 5 indicates that the value of the element pij in the ith row and jth column

of the product matrix [P] is to be calculated by multiplying the elements in the ith

row of the matrix [A] by the corresponding elements in the jth column of the matrix

[B]. It is therefore evident that the number of elements in the ith row of [A] should

be equal to the number of elements in the jth column of [B]. In other words, to

apply Equation 5 for producing a product matrix [P] by multiplying a matrix [A]

on the right by a matrix [B] (or, to say multiplying a matrix [B] on the left by a

matrix [A]), the number of columns of [A] should be equal to the number of row

of [B]. A matrix [A] of order L by M can therefore be post-multiplied by a matrix

[B] of order M by N; but [A] cannot be pre-multiplied by [B] unless L is equal to N!

As a numerical example, consider the case of a square, 3 ? 3 matrix postmultiplied by a rectangular matrix of order 3 by 2. Since L = 3, M = 3, and N = 2,

the product matrix is thus of order 3 by 2.

© 2001 by CRC Press LLC

?1

?4

?

?? 7

3? ? 6

6?? ??5

9?? ??4

2

5

8

?3? ?1(6) + 2(5) + 3( 4) 1( ?3) + 2( ?2) + 3( ?1) ?

?

?

?2 ?? = ?4(6) + 5(5) + 6( 4) 4( ?3) + 5( ?2) + 6( ?1)?

?1?? ?? 7(6) + 8(5) + 9( 4) 7( ?3) + 8( ?2) + 9( ?1)??

? 6 + 10 + 12

= ?? 24 + 25 + 24

??42 + 40 + 32

?3 ? 4 ? 3 ? ? 28

?12 ? 10 ? 5?? = ?? 73

?21 ? 16 ? 9 ?? ??114

?10 ?

?27??

?46??

More exercises are given in the Problems listed at the end of this chapter for

the readers to practice on the matrix multiplications based on Equation 5.

It is of interest to note that the square of the length of a vector {V} which has

N components as defined in Equation 4, {V}2, can be obtained by application of

Equation 5 to {V} and its transpose denoted as {V}T which is a row matrix of order

1 by N (one row and N columns). That is:

{V} = {V}T {V} = v11 + v22 + … + v32

2

(6)

For a L-by-M matrix having elements eij where the row index i ranges from 1

to L and the column index j ranges from 1 to M, the transpose of this matrix when

its elements are designated as trc will have a value equal to ecr where the row index

r ranges from 1 to M and the column index c ranges from 1 to M because this

transpose matrix is of order M by L. As a numerical example, here is a pair of a

3 ? 2 matrix [G] and its 2 ? 3 transpose [H]:

?6

[G] = ??5

3? 2

??4

?3?

?2 ?? and

?1??

?6

??3

[H] = [G]T = ?

2 ?3

5

?2

4?

?1??

If the elements of [G] and [H] are designated respectively as gij and hij, then

hij = gji. For example, from above, we observe that h12 = g21 = 5, h23 = g32 = –1, and

so on. There will be more examples of applications of Equations 5 and 6 in the

ensuing sections and chapters.

Having introduced the transpose of a matrix, we can now conveniently revisit

the addition of {D} and {E} in Figure 1 in algebraic form as {F} = {D} + {E} =

[4 –3]T + [–3 1]T = [4+(–3) –3+1]T = [1 –2]T. The resulting sum vector is indeed

correct as it is graphically verified in Figure 1. The saving of space by use of

transposes of vectors (row matrices) is not evident in this case because all vectors

are two-dimensional; imagine if the vectors are of much higher order.

Another noteworthy application of matrix multiplication and transposition is to

reduce a system of linear algebraic equations into a simple, (or, should we say a

single) matrix equation. For example, if we have three unknowns x, y, and z which

are to be solved from the following three linear algebraic equations:

© 2001 by CRC Press LLC

x + 2 y + 3z = 4

5x + 6y + 7z = 8

?2 x ? 37 = 9

(7)

Let us introduce two vectors, {V} and {R}, which contain the unknown x, y,

and z, and the right-hand-side constants in the above three equations, respectively.

That is:

?x ?

T

V

=

x

y

z

=

{ } [

] ??y??

?? z ??

and

?4?

T

R

=

4

8

9

=

{ } [

] ??8??

??9 ??

(8)

Then, making use of the multiplication rule of matrices, Equation 5, the system

of linear algebraic equations, 7, now can be written simply as:

[C]{V} = {R}

(9)

where the coefficient matrix [C] formed by listing the coefficients of x, y, and z in

first equation in the first row and second equation in the second row and so on. That is,

?1

[C] = ?? 5

???2

2

6

?3

3?

7??

0 ??

There will be more applications of matrix multiplication and transposition in

the ensuing chapters when we discuss how matrix equations, such as [C]{V} = {R},

can be solved by employing the Gaussian Elimination method, and how ordinary

differential equations are approximated by finite differences will lead to the matrix

equations. In the abbreviated matrix form, derivation and explanation of computational methods becomes much simpler.

Also, it can be observed from the expressions in Equation 8 how the transposition

can be conveniently used to define the two vectors not using the column matrices

which take more lines.

FORTRAN VERSION

Since Equations 1 and 2 require repetitive computation of the elements in the

sum matrix [S] and difference matrix [D], machine could certainly help to carry out

this laborous task particularly when matrices of very high order are involved. For

covering all rows and columns of [S] and [D], looping or application of DO statement

of the FORTRAN programming immediately come to mind. The following program

is provided to serve as a first example for generating [S] and [D] of two given

matrices [A] and[B]:

© 2001 by CRC Press LLC

The resulting display on the screen is:

To review FORTRAN briefly, we notice that matrices should be declared as

variables with two subscripts in a DIMENSION statement. The displayed results of

matrices A and B show that the values listed between // in a DATA statment will be

filling into the first column and then second column and so on of a matrix. To instruct

the computer to take the values provided but to fill them into a matrix row-by-row,

a more explicit DATA needs to be given as:

DATA ((A(I,J),J = 1,3),I = 1,3)/1.,4.,7.,2.,5.,8.,3.,6.,9./

When a number needs to be repeated, the * symbol can be conveniently applied

in the DATA statement exemplified by those for the matrix [B].

Some sample WRITE and FORMAT statements are also given in the program.

The first * inside the parentheses of the WRITE statement when replaced by a

number allows a device unit to be specified for saving the message or the values of

the variables listed in the statement. * without being replaced means the monitor

will be the output unit and consequently the message or the value of the variable(s)

will be displayed on screen. The second * inside the parentheses of the WRITE

© 2001 by CRC Press LLC

statement if not replaced by a statement number, in which formats for printing the

listed variables are specified, means “unformatted” and takes whatever the computer

provides. For example, statement number 15 is a FORMAT statement used by the

WRITE statement preceding it. There are 18 variables listed in that WRITE statement

but only six F5.1 codes are specified. F5.1 requests five column spaces and one digit

after the decimal point to be used to print the value of a listed variable. / in a

FORMAT statement causes the print/display to begin at the first column of the next

line. 6F5.1 is, however, enclosed by the inner pair of parentheses that allows it to

be reused and every time it is reused the next six values will be printed or displayed

on next line. The use (*,*) in a WRITE statement has the convenience of viewing

the results and then making a hardcopy on a connected printer by pressing the PrtSc

(Print Screen) key.

INTERACTIVE OPERATION

Program MatxAlgb.1 only allows the two particular matrices having their elements specified in the DATA statement to be added and subtracted. For finding the

sum matrix [S] and difference matrix [D] for any two matrices of same order N, we

ought to upgrade this program to allow the user to enter from keyboard the order

N and then the elements of the two matrices involved. This is interactive operation

of the program and proper messages should be given to instruct the user what to do

which means the program should be user-friendly. The program MatxAlgb.2 listed

below is an attempt to achieve that goal:

© 2001 by CRC Press LLC

The interactive execution of the problem solved by the previous version Matxalgb.1

now can proceed as follows:

© 2001 by CRC Press LLC

The results are identical to those obtained previously. The READ statement

allows the values for the variable(s) to be entered via keyboard. A WRITE statement

has no variable listed serves for need of skipping a line to provide better readability

of the display. Also the I and E format codes are introduced in the statement 10. Iw

where w is an integer in a FORMAT statement requests w columns to be provided

for displaying the value of the integer variable listed in the WRITE statement, in

which the FORMAT statement is utilized. Ew.d where w and d should both be integer

constants requests w columns to be provided for display a real value in the scientific

form and carrying d digits after the decimal point. Ew.d format gives more feasibility

than Fw.d format because the latter may cause an error message of insufficient width

if the value to be displayed becomes too large and/or has a negative sign.

MORE PROGRAMMING REVIEW

Besides the operation of matrix addition and subtraction, we have also discussed

about the transposition and multiplication of matrices. For further review of computer

programming, it is opportune to incorporate all these matrix algebraic operations

into a single interactive program. In the listing below, three subroutines for matrix

addition and subtraction, transposition, and multiplication named as MatrixSD,

Transpos, and MatxMtpy, respectively, are created to support a program called

MatxAlgb (Matrix Algebra).

© 2001 by CRC Press LLC

© 2001 by CRC Press LLC

The above program shows that Subroutines are independent units all started with

a SUBROUTINE statement which includes a name followed by a pair of parentheses

enclosing a number of arguments. The Subroutines are called in the main program

by specifying which variables or constants should serve as arguments to connect to

the subroutines. Some arguments provide input to the subroutine while other arguments transmit out the results determined by the subroutine. These are referred to

as input arguments and output arguments, respectively. In many instances, an argument may serve a dual role for both input and output purposes. To construct as an

independent unit, a subprogram which can be in the form of a SUBROUTINE, or

FUNCTION (to be elaborated later) must have RETURN and END statements.

It should also be remarked that program MatxAlgb is arranged to handle any

matrix having an order of no higher than 25 by 25. For this restriction and for having

the flexibility of handling any matrices of lesser order, the Lmax, Mmax, and Nmax

arguments are added in all three subroutines in order not to cause any mismatch of

matrix sizes between the main program and the called subroutine when dealing with

any L, M, and N values which are interactively entered via keyboard.

Computed GOTO and arithmetic IF statements are also introduced in the program MatxAlgb. GOTO (i,j,k,…) C will result in going to (execute) the statement

numbered i, j, k, and so on when C has a value equal to 1, 2, 3, and so on, respectively.

IF (Expression) a,b,c will result in going to the statement numbered a, b, or c if the

value calculated by the expression or a single variable is less than, equal to, or,

greater than zero, respectively.

It is important to point out that in describing any derived procedure of numerical

computation, indicial notation such as Equation 5 should always be preferred to

facilitate programming. In that notation, the indices are directly used, or, literally

translated into the index variables for the DO loops as can be seen in Subroutine

MatxMtpy which is developed according to Equation 5. Subroutine MatrixSD is

another example of literally translating Equations 1 and 2. For defining the values

of the element in the following tri-diagonal band matrix:

© 2001 by CRC Press LLC

?1

??3

?

[C ] = ? 0

?

?0

?? 0

2

1

?3

0

0

0

2

1

?3

0

0

0

2

1

?3

0?

0?

?

0?

?

2?

1 ??

we ought not to write 25 separate statements for the 25 elements in this matrix but

derive the indicial formulas for i,j = 1 to 5:

cij = 0,

if j > i + 2, or, j < i ? 2

ci,i +1 = 2,

and

ci,i ?1 = ?3

Then, the matrix [C] can be generated with the DO loops as follows:

The above short program also demonstrates the use of the CONTINUE statement for ending the DO loop(s), and the logical IF statements. The true, or, false

condition of the expression inside the outer pair of parentheses directs the computer

to execute the statement following the parentheses or the next statement immediately

below the current IF statement. Reader may want to practice on deriving indicial

formulas and then write a short program for calculating the elements of the matrix:

?1

?2

?

?3

?

4

[M] = ??

5

?

?6

?7

?

??8

© 2001 by CRC Press LLC

0

1

2

0

0

1

0

0

0

0

0

0

0

0

0

0

0

0

3

4

5

6

7

2

3

4

5

6

1

2

3

4

5

0

1

2

3

4

0

0

1

2

3

0

0

0

1

2

0?

0?

?

0?

?

0?

0?

?

0?

0?

?

1 ??

(10)

As another example of writing a computer program based on indicial notation,

consider the case of calculating ex based on the infinite series:

ex = 1 +

?

=

x1 x 2 x 2

xi

+

+

+…+ +…

1! 2! 3!

i!

?

i=0

xi

i!

(11)

With the understanding that 0! = 1, we have expressed the series as a summation

involving the index i which ranges from zero to infinity. A FUNCTION ExpoFunc

can be developed for calculating ex based on Equation 11 and taking only a finite

number of terms for a partial sum of the series when the contribution of additional

term is less than certain percentage of the sum in magnitude, say 0.001%. This

FUNCTION may be arranged as:

To further show the advantage of adopting vector and matrix notation, here let

us apply FUNCTION ExpoFunc to examine the surface z(x,y) = ex + y above the

rectangular area 0?x?2.0 and 0?y?1.5. The following program, ExpTest, will enable

us to compare the values of ex + y generated by the FUNCTION ExpoFunc and by

the function EXP available in the FORTRAN library (hence called library function).

© 2001 by CRC Press LLC

The resulting printout is:

It is apparent that two approaches produce almost identical results, so the 0.001%

accuracy appears quite adequate for the x and y ranges studied. Also, arranging the

results in vector and matrix forms make the presentation much easy to comprehend.

We have experienced how the summation process for an indicial formula involving a ? should be programmed. Another operation symbol of importance is ? which

is for multiplication of many factors. That is:

N

? a = a a …a

i

1 2

N

(12)

i =1

An obvious application of Equation 12 is for the calculation of factorials. For

example, 5! = ?i for i ranges from 1 to 5. As an exercise, we display the values of

1! through 50! with the following program involving a subroutine IFACTO which

calculates I! for a specified I value:

© 2001 by CRC Press LLC

The resulting print out is (listed in three columns for saving space)

Another application of Equation 12 is for calculation of the binomial coefficients

for a real number r and an integer k defined as:

? r ? r( r ? 1)( r ? 2)…( r ? k + 1)

=

? ?=

? k?

k!

k

?

i =1

r ? i +1

i

(13)

We shall have the occasion of applying Equations 12 and 13 when the finite

differences and Lagrangian interpolation are discussed.

Sample Applications

Program MatxAlgb has been tested interactively, the following are the resulting

displays of four test cases:

© 2001 by CRC Press LLC

© 2001 by CRC Press LLC

QUICKBASIC VERSION

The QuickBASIC language has the advantage over the FORTRAN language

for making quick changes and then running the revised program without compilation.

Furthermore, it offers simple plotting statements. Let us have a QuickBASIC version

of the program MatxAlgb and then discuss its basic features.

© 2001 by CRC Press LLC

© 2001 by CRC Press LLC

Notice that the order limit of 25 needed in the FORTRAN version is removed

in the QuickBASIC version which allows the dim statement to be adjustable. ' is

replacing C in FORTRAN to indicate a comment statement in QuickBASIC. READ

and WRITE in FORTRAN are replaced by INPUT and PRINT in QuickBASIC,

respectively. The DO loop in FORTRAN is replaced by the FOR and NEXT pair

in QuickBASIC.

Sample Applications

When the four cases previously run by the FORTRAN version are executed by

the QuickBASIC version, the screen prompting messages, the interactively entered

data, and the computed results are:

© 2001 by CRC Press LLC

© 2001 by CRC Press LLC

Matrix [P]

Row 1

0.1400E+02

Row 2

0.3200E+02

0.2000E+02

0.4700E+02

MATLAB APPLICATIONS

MATLAB developed by the Mathworks, Inc. offers a quick tool for matrix

manipulations. To load MATLAB after it has been set-up and stored in a subdirectory

of a hard drive, say C, we first switch to this subdirectory by entering (followed by

pressing ENTER)

C:\cd MATLAB

and then switch to its own subdirectory BIN by entering (followed by pressing ENTER)

C:\MATLAB>cd BIN

Next, we type MATLAB to obtain a display of:

C:\MATLAB>BIB>MATLAB

Pressing the ENTER key results in a display of:

>>

which indicates MATLAB is ready to begin. Let us rerun the cases of matrix

subtraction, addition, transposition, and multiplication previously considered in the

FORTRAN and QuickBASIC versions. First, we enter the matrix [A] in the form of:

>> A = [1,2;3,4]

When the ENTER key is pressed, the displayed result is:

© 2001 by CRC Press LLC

A=

1

2

3

4

Notice that the elements of [A] should be entered row by row. While the rows

are separated by ;, in each row elements are separated by comma. After the print

out of the above results, >> sign will again appear. To eliminate the unnecessary line

space (between A = and the first row 1 2), the statement format compact can be entered

as follows (the phrase “pressing ENTER key” will be omitted from now on):

>> format compact, B = [5,6;7,8]

B=

5

6

7

8

Notice that comma is used to separate the statements. To demonstrate matrix subtraction and addition, we can have:

>> A-B

ans =

–4

–4

–4

–4

>> A + B

ans =

6

8

10

12

To apply MATLAB for transposition and multiplication of matrices, we can have:

>> C = [1,2,3;4,5,6]

C=

1

2

3

4

5

6

>> C'

ans =

1

4

2

5

3

6

© 2001 by CRC Press LLC

>> D = [1,2,3;4,5,6]; E = [1,2;2,3;3,4]; P = D*E

P=

14

20

32

47

Notice that MATLAB uses ' (single quote) in place of the superscripted symbol

T for transposition. When ; (semi-colon) follows a statement such as the D statement,

the results will not be displayed. As in FORTRAN and QuickBASIC, * is the

multiplication operator as is used in P = D*E, here involving three matrices not three

single variables. More examples of MATLAB applications including plotting will

ensue. To terminate the MATLAB operation, simply enter quit and then the

RETURN key.

MATHEMATICA APPLICATIONS

To commence the service of Mathematica from Windows setup, simply point

the mouse to it and double click the left button. The Input display bar will appear

on screen, applications of Mathematica can start by entering commands from

keyboard and then press the Shift and Enter keys. To terminate the Mathematica

application, enter Quit[] from keyboard and then press the Shift and Enter keys.

Mathematica commands, statements, and functions are gradually introduced

and applied in increasing degree of difficulty. Its graphic capabilities are also utilized

in presentation of the computed results.

For matrix operations, Mathematica can compute the sum and difference of

two matrices of same order in symbolic forms, such as in the following cases of

involving two matrices, A and B, both of order 2 by 2:

In[1]: = A = {{1,2},{3,4}}; MatrixForm[A]

Out[1]//MatrixForm =

1

2

3

4

In[1]: = is shown on screen by Mathematica while user types in A =

{{1,2},{3,4}}; MatrixForm[A]. Notice that braces are used to enclose the elements

in each row of a matrix, the elments in a same row are separated by commas, and

the rows are also separated by commas. MatrixForm demands that the matrix be

printed in a matrix form. Out[1]//MatrixForm = and the rest are response of Mathematica.

In[2]: = B = {{5,6},{7,8}}; MatrixForm[B]

Out[2]//MatrixForm =

5

7

6

8

© 2001 by CRC Press LLC

In[3]: = MatrixForm[A + B]

Out[3]//MatrixForm =

6

8

10

12

In[4]: = Dif = A-B; MatrixForm[Dif]

Out[4]//MatrixForm =

–4

–4

–4

–4

In[3] and In[4] illustrate how matrices are to be added and subtracted, respectively. Notice that one can either use A + B directly, or, create a variable Dif to

handle the sum and difference matrices.

Also, Mathematica has a function called Transpose for transposition of a

matrix. Let us reuse the matrix A to demonstrate its application:

In[5]: = AT = Transpose[A]; MatrixForm[AT]

Out[5]//MatrixForm =

1

3

2

4

1.3 SOLUTION OF MATRIX EQUATION

Matrix notation offers the convenience of organizing mathematical expression in an

orderly manner and in a form which can be directly followed in coding it into

programming languages, particularly in the case of repetitive computation involving

the looping process. The most notable situation is in the case of solving a system

of linear algebraic equation. For example, if we need to determine a linear equation

y = a1 + a2x which geometrically represents a straight line and it is required to pass

through two specified points (x1,y1) and (x2,y2). To find the values of the coefficients

a1 and a2 in the y equation, two equations can be obtained by substituting the two

given points as:

(1)a1 + (x1 )a 2 = y1

(1)

(1)a1 + (x 2 )a 2 = y2

(2)

and

© 2001 by CRC Press LLC

To facilitate programming, it is advantageous to write the above equations in

matrix form as:

[C]{A} = {Y}

(3)

where:

?1

[C] = ?1

?

x1 ?

? a1 ?

? y1 ?

, {A} = ? ?, and {Y} = ? ?

?

x2 ?

?a 2 ?

?y 2 ?

(4)

The matrix equation 3 in this case is of small order, that is an order of 2. For

small systems, Cramer’s Rule can be conveniently applied which allows the unknown

vector {A} to be obtained by the formula:

[

{A} = [c1 ] [c2 ]

]

T

[C ]

(5)

Equation 5 involves the calculation of three determinants, i.e., , [C1], [C2], and

[C] where [C1] and [C2] are matrices derived from the matrix [C] when the first

and second columns of [C] are replaced by {Y}, respectively. If we denote the

elements of a general matrix [C] of order 2 by cij for i,j = 1,2, the determinant of

[C] by definition is:

[C] = c11c22 ? c12c21

(6)

The general definition of the determinant of a matrix [M] of order N and whose

elements are denoted as mij for i,j = 1,2,…,N is to add all possible product of N

elements selected one from each row but from different column. There are N! such

products and each product carries a positive or negative sign depending on whether

even or odd number of exchanges are necessary for rearranging the N subscripts in

increasing order. For example, in Equation 6, c11 is selected from the first row and

first column of [C] and only c22 can be selected and multiplied by it while the other

possible product is to select c12 from the second row and first column of [C] and

that leaves only c21 from the second row and first column of [C] available as a factor

of the second product. In order to arrange the two subscripts in non-decreasing order,

one exchange is needed and hence the product c12c21 carries a minus sign. We shall

explain this sign convention further when a matrix of order 3 is discussed. However,

it should be evident here that a matrix whose order is large the task of calculating

its determinant would certainly need help from computer. This will be the a topic

discussed in Section 1.5.

Let us demonstrate the application of Cramer’s Rule by having a numerical case.

If the two given points to be passed by the straight line y = a1 + a2x are (x1,y1) =

(1,2) and (x2,y2) = (3,4). Then we can have:

© 2001 by CRC Press LLC

[C ] = 1

1

x1 1

=

x2 1

1

= 1? 3 ?1?1 = 3 ?1 = 2

3

y1

x1 2

=

x2 4

1

= 2 ? 3 ?1? 4 = 6 ? 4 = 2

3

y1 1

=

y2 1

2

= 1? 4 ? 2 ?1 = 4 ? 2 = 2

4

[c ] = y

1

2

and

1

[C ] = 1

2

Consequently, according to Equation 5 we can find the coefficients in the

straight-line equation to be:

a1 = [C1 ]

[C ] = 2 2 = 1

and a 2 = [C2 ] [C] = 2 2 = 1

Hence, the line passing through the points (1,2) and (3,4) is y = a1 + a2x = 1 + x.

Application of Cramer’s Rule can be extended for solving three unknowns from

three linear algebraic equations. Consider the case of finding a plane which passes

three points (xi,yi) for i = 1 to 3. The equation of that plane can first be written as

z = a1 + a2x + a3y. Similar to the derivation of Equation 3, here we substitute the

three given points into the z equation and obtain:

(1)a1 + (x1 )a 2 + (y1 )a 3 = z1

(7)

(1)a1 + (x 2 )a 2 + (y2 )a 3 = z2

(8)

(1)a1 + (x3 )a 2 + (y3 )a 3 = z3

(9)

and

Again, the above three equations can be written in matrix form as:

[C]{A} = {Z}

(10)

where the matrix [C] and the vector {A} previously defined in Equation 4 need to

be reexpanded and redefined as:

1

[C ] = 1

1

© 2001 by CRC Press LLC

x1

x2

x3

y1

a1

? z1 ?

? ?

y 2 , {A} = a 2 , and {Z} = ?z 2 ?

??z3 ??

y3

a3

(11)

And, the Cramer’s Rule for solving Equation 10 can be expressed as:

[

{A} = [C1 ][C2 ] [C3 ]

]

T

[C ]

(12)

where [Ci] for i = 1 to 3 for matrices formed by replacing the ith column of the

matrix [C] by the vector {Z}, respectively. Now, we need the calculation of the

determinant of matrices of order 3. If we denote the element in a matrix [M] as mij

for i,j = 1 to 3, the determinant of [M] can be calculated as:

[M] = m11m 22 m33 ? m11m 23m32 + m12 m 23m31

? m12 m 21m 33 + m13m 21m 32 ? m13m 22 m 31

(13)

To give a numerical example, let us consider a plane passing the three points,

(x1,y1,z1) = (1,2,3), (x2,y2,z2) = (–1,0,1), and (x3,y3,z3) = (–4,–2,0). We can then have:

1

[C ] = 1

1

z1

[C1 ] = z2

z3

1

C

=

[ 2] 1

1

x1

x2

x3

y1 1

y2 = 1

y3 1

1

?1

?4

2

0 = ?2

?2

x1

x2

x3

y1 3

y2 = 1

y3 0

1

?1

?4

2

0 =0

?2

z1

z2

z3

y1 1

y2 = 1

y3 1

3

1

0

2

0 =2

?2

1

?1

?4

3

1 = ?4

?2

and

[C ]

3

1

=1

1

x1

x2

x

z1 1

z2 = 1

z 1

According to Equation 13, we find a1 = [C1]/[C] = 0/(–2) = 0, a2 = [C2]/[C] =

2/(–2) = –1, and a3 = [C3]/[C] = –4/(–2) = 2. Thus, the required plane equation is

z = a1 + a2x + a3y = -x + 2y.

QUICKBASIC VERSION

OF THE PROGRAM

CRAMERR

A computer program called CramerR has been developed as a reviewing exercise in programming to solve a matrix equation of order 3 by application of Cramer

© 2001 by CRC Press LLC

Rule and the definition of determinant of a 3 by 3 square matrix according to

Equations 12 and 13, respectively. First, a subroutine called Determ3 is created

explicitly following Equation 13 as listed below:

To interactively enter the elements of the coefficient matrix [C] and also the

elements of the right-hand-side vector {Z} in Equation 12 and to solve for {A}, the

program CramerR can be arranged as:

© 2001 by CRC Press LLC

1.4 PROGRAM GAUSS

Program Gauss is designed for solving N unknowns from N simultaneous, linear

algebraic equations by the Gaussian Elimination method. In matrix notation, the

problem can be described as to solve a vector {X} from the matrix equation:

[C]{X} = {V}

(1)

where [C] is an NxN coefficient matrix and {V} is a Nx1 constant vector, and both

are prescribed. For example, let us consider the following system:

9x1 + x 2 + x3 = 10

(2)

3x1 + 6x 2 + x3 = 14

(3)

2 x1 + 2 x 2 + 3x3 = 3

(4)

If the above three equations are expressed in matrix form as Equation 1, then:

?9

[C] = ??3

??2

1

6

2

1?

1??,

3??

?10 ?

{V} = ??14??,

?? 3 ??

(5,6)

and

? x1 ?

X

=

{ } ??x 2 ?? = x1 x 2 x3

?? x3 ??

[

]

T

(7)

where T designates the transpose of a matrix.

GAUSSIAN ELIMINATION METHOD

A systematic procedure named after Gauss can be employed for solving x1, x2,

and x3 from the above equations. It consists of first dividing Equation 28 by the

leading coefficient, 9, to obtain:

x1 +

1

1

10

x + x =

9 2 9 3 9

(8)

This step is called normalization of the first equation of the system (1). The next

step is to eliminate x1 term from the other (in this case, two) equations. To do that,

we multiply Equation 8 by the coefficients associated with x1 in Equations 3 and 4,

respectively, to obtain:

© 2001 by CRC Press LLC

1

1

10

3x1 + x 2 + x3 =

3

3

3

(9)

2

2

20

x 2 + x3 =

9

9

9

(10)

and

2 x1 +

If we subtract Equation 9 from Equation 3, and subtract Equation 10 from Equation 4, the x1 terms are eliminated. The resulting equations are, respectively:

17

2

32

x 2 + x3 =

3

3

3

(11)

16

25

7

x 2 + x3 =

9

9

9

(12)

and

This completes the first elimination step. The next normalization is applied to

Equation 11, and then the x2 term is to be eliminated from Equation 12. The resulting

equations are:

2

32

(13)

x 2 + x3 =

17

17

and

393

393

(14)

x =?

153 3

153

The last normalization of Equation 14 then gives:

x3 = ?1

(15)

Equations 8, 13, and 15 can be organized in matrix form as:

?1

{V} = ??0

??0

19

1

0

1 9 ? ? x1 ? ? 10 9 ?

? ?

2 17?? ?x 2 ? = ??32 17??

1 ?? ?? x3 ?? ?? ?1 ??

(16)

The coefficient matrix is now a so-called upper triangular matrix since all

elements below the main diagonal are equal to zero.

As x3 is already obtained in Equation 15, the other two unknowns, x2 and x3,

can be obtained by a sequential backward-substitution process. First, Equation 13

can be used to obtain:

x2 =

© 2001 by CRC Press LLC

32 2

32 2

32 + 2

=2

? x3 =

? ( ?1) =

17 17

17 17

17

Once, both x2 and x3 have been calculated, x1 can be obtained from Equation 8 as:

x1 =

10 1

1

10 1

1

10 ? 2 + 1

=1

? x 2 ? x3 =

? (2) ? ( ?1) =

9 9

9

9 9

9

9

To derive a general algorithm for the Gaussian elimination method, let us denote

the elements in [C], {X}, and {V} as ci,j, xi, and vi, respectively. Then the normalization of the first equation can be expressed as:

(c )

1, j new

( ) (c )

= c1, j

old

1,1 old

(17)

and

(v1 )new = (v1 ) old (c1,1 )old

(18)

Equation 17 is to be used for calculating the new coefficient associated with xj

in the first, normalized equation. So, j should be ranged from 2 to N which is the

number of unknowns (equal to 3 in the sample case). The subscripts old and new

are added to indicate the values before and after normalization, respectively. Such

designation is particularly helpful if no separate storage in computer are assigned

for [C] for the values of its elements. Notice that (c1,1)new = 1 is not calculated.

Preserving this diagonal element enables the determinant of [C] to be computed.

(See the topic on matrix inversion and determinant.)

The formulas for the elimination of x1 terms from the second equation are:

(c )

2 , j new

( ) ? (c ) (c )

= c2, j

old

2 ,1 old

1, j old

(19)

for j = 2,3,…,N (there is no need to include j = 1) and

(v2 )new = (v2 )old ? (c2,1 )old (v1 )old

(20)

By changing the subscript 2 in Equations 19 and 20, x1 term in the third equation

can be eliminated. In other words, the general formulas for elimination of x1 terms

from all equation other than the first equation are, for k = 2,3,…,N

(c )

k , j new

for j = 2,3,…,N

© 2001 by CRC Press LLC

( ) ? (c ) (c )

= ck,j

old

k ,l old

1, j old

(21)

(vk )new = (vk )old ? (ck,1 )old (v1 )old

(22)

Instead of normalizing the first equation, we can generalize Equations 17 and

18 for normalization of the ith equation, for i = 1,2,…,N to the expressions:

(c )

i , j new

( ) (c )

= ci. j

old

i ,i old

(23)

for j = i + 1,i + 2,…,N and

(vi )new = (vi ) old (ci,i )old

(24)

Note that (ci,i)new should be equal to 1 but no need to calculate since it is not

involved in later calculation for finding {X}.

Similarly, elimination of xi term from kth equation for k = i + 1,i + 2,…,N

consists of using the general formula:

(c )

k , j new

( ) ? (c ) (c )

= ck,j

k ,i old

old

i , j old

(25)

for j = i + 1,i + 2,…,N and

(vk )new = (vk )old ? (ck,i )old (vi )old

(26)

Backward substitution for finding xi involves the calculation of:

N

x i = vi ?

?c

x

i, j j

(27)

j= i +1

for i = N–1,N–2,…,2,1. Note that xN is already found equal to vN after the Nth

normalization.

Program Gauss listed below in both QuickBASIC and FORTRAN languages

is developed for interactive specification of the number of unknowns, N, and the

values of the elements of [C] and {V}. It proceeds to solve for {X} and prints out

the computed values. Sample applications of both languages are provided immediately following the listed programs.

A subroutine Gauss.Sub is also made available for the frequent need in the

other computer programs which require the solution of matrix equations.

© 2001 by CRC Press LLC

QUICKBASIC VERSION

Sample Application

© 2001 by CRC Press LLC

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Application

GAUSS-JORDAN METHOD

One slight modification of the elimination step will make the backward substitution steps completely unnecessary. That is, during the elimination of the xi terms

from the linear algebraic equations except the ith one, Equations 25 and 26 should

be applied for k equal to 1 through N and excluding k = i. For example, the x3 terms

should be eliminated from the first, second, fourth through Nth equations. In this

manner, after the Nth normalization, [C] becomes an identity matrix and {V} will

have the elements of the required solution {X}. This modified method is called

Gauss-Jordan method.

A subroutine called GauJor is made available based on the above argument. In

this subroutine, a block of statements are also added to include the consideration of

the pivoting technique which is required if ci,i = 0. The normalization steps,

Equations 49 and 50, cannot be implemented if ci,i is equal to zero. For such a

situation, a search for a nonzero ci,k is necessary for i = k + 1,k + 2,…,N. That is,

to find in the kth column of [C] and below the kth row a nonzero element. Once

this nonzero ci,k is found, then we can then interchange the ith and kth rows of [C]

and {V} to allow the normalization steps to be implemented; if no nonzero ci,k can

be found then [C] is singular because the determinant of [C] is equal to zero! This

can be explained by the fact that when ck,k = 0 and no pivoting is possible and the

determinant D of [C] can be calculated by the formula:

N

D = c1,1c2,2 … c k ,k … c N,N =

?c

k =1

where indicates a product of all listed factors.

© 2001 by CRC Press LLC

k ,k

(28)

A subroutine has been written based on the Gauss-Jordan method and called

GauJor.Sub. Both QuickBASIC and FORTRAN versions are made available and

they are listed below.

QUICKBASIC VERSION

© 2001 by CRC Press LLC

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Applications

The same problem previously solved by the program Gauss has been used again

but solved by application of subroutine GauJor. The results obtained with the QuickBASIC and FORTRAN versions are listed, in that order, below:

MATLAB APPLICATIONS

For solving the vector {X} from the matrix equation [C]{X} = {R} when both

the coefficient matrix [C] and the right-hand side vector {R} are specified, MATLAB

simply requires [C] and {R} to be interactively inputted and then uses a statement

X = C\R to obtain the solution vector {X} by multiplying the vector {R} on the left

of the inverse of [C] or dividing {R} on the left by [C]. More details are discussed

in the program MatxAlgb. Here, for providing more examples in MATLAB applications, a m file called GauJor.m is presented below as a companion of the FORTRAN and QuickBASIC versions:

© 2001 by CRC Press LLC

This file GauJor.m should then be added into MATLAB. As an example of

interactive application of this m file, the sample problem used in the FORTRAN

and QuickBASIC versions is again solved by specifying the coefficient matrix [C]

and the right hand side vector {R} to obtain the resulting display as follows:

The results of the vector {X} and determinant D for the coefficient matrix [C]

are same as obtained before.

MATHEMATICA APPLICATIONS

For solving a system of linear algebraic equations which has been arranged in

matrix form as [A]{X} = {R}, Mathematica’s function LinearSolve can be applied

© 2001 by CRC Press LLC

to solve for {X} when the coefficient matrix [A] and the right-hand side vector {R}

are both provided. The following is an example of interactive application:

In[1]: = A = {{3,6,14},{6,14,36},{14,36,98}}

Out[1]: =

{{3, 6, 14}, {6, 14, 36}, {14, 36, 98}}

In[2]: = MatrixForm[A]

Out[2]//MatrixForm: =

3

6

14

6

14

36

14

36

98

In[3]: = R = {9,20,48}

Out[3]: =

{9, 20, 48}

In[4]: = LinearSolve[A,R]

Out[4]: =

{–9,13,–3}

Output[2] and Output[1] demonstrate the difference in display of matrix [A]

when MatrixeForm is requested, or, not requested, respectively. It shows that without

requesting of MatrixForm, some screen space saving can be gained. Output[4] gives

the solution {X} = [–9 13 –3]T for the matrix equation [A]{X} = {R} where the

coefficient matix [A] and vector {R} are provided by Input[1] and Input[3], respectively.

1.5 MATRIX INVERSION, DETERMINANT,

AND PROGRAM MatxInvD

Given a square matrix [C] of order N, its inverse as [C]–1 of the same order is defined

by the equation:

[C][C]?1 = [C]?1[C] = [I]

(1)

where [I] is an identity matrix having elements equal to one along its main diagonal

and equal to zero elsewhere. That is:

.

?1 0

?0 1 0

[I] = ?

? . . .

?

.

?0 0

© 2001 by CRC Press LLC

.

.

.

.

.

.

.

.

.

.

.

.

0?

0?

?

?

?

1?

(2)

To find [C]–1, let cij and dij be the elements at the ith row and jth column of the

matrices [C] and [C]–1, respectively. Both i and j range from 1 to N. Furthermore,

let {Dj} and {Ij} be the jth column of the matrices [C]–1 and [I], respectively. It is

easy to observe that {Ij} has elements all equal to zero except the one in the jth row

which is equal to unity. Also,

{D } = [d d

j

lj 2 j

… d Nj

]

T

(3)

and

[C]?1 = [D1D2 … DN ]

T

(4)

Based on the rules of matrix multiplication, Equation 1 can be interpreted as

[C]{D1} = {I1}, [C]{D2} = {I2}, …, and [C]{DN} = {IN}. This indicates that program

Gauss can be successively employed N times by using the same coefficient matrix

[C] and the vectors {Ii} to find the vectors {Di} for i = 1,2,…,N. Program MatxInvD

is developed with this concept by modifying the program Gauss. It is listed below

along with a sample interactive run.

QUICKBASIC VERSION

© 2001 by CRC Press LLC

Sample Application

FORTRAN VERSION

© 2001 by CRC Press LLC

© 2001 by CRC Press LLC

Sample Applications

MATLAB APPLICATION

MATLAB offers very simple matrix operations. For example, matrix inversion

can be implemented as:

To check if the obtained inversion indeed satisfies the equation [A}[A]–1 = [I]

where [I] is the identity matrix, we enter:

Once [A]–1 becomes available, we can solve the vector {X} in the matrix equation

[A]{X} = {R} if {R} is prescribed, namely {X} = [A]–1{R}. For example, may enter

a {R} vector and find {X} such as:

© 2001 by CRC Press LLC

MATHEMATICA APPLICATIONS

Mathematica has a function called Inverse for inversion of a matrix. Let us

reuse the matrix A that we have entered in earlier examples and continue to demonstrate the application of Inverse:

In[1]: = A = {{1,2},{3,4}}; MatrixForm[A]

Out[1]//MatrixForm =

1

2

3

4

In[2]: = B = {{5,6},{7,8}}; MatrixForm[B]

Out[2]//MatrixForm =

5

6

7

8

In[3]: = MatrixForm[A + B]

Out[3]//MatrixForm =

6

8

10

12

In[4]: = Dif = A-B; MatrixForm[Dif]

Out[4]//MatrixForm =

–4

–4

–4

–4

In[5]: = AT = Transpose[A]; MatrixForm[AT]

Out[5]//MatrixForm =

1

3

2

4

In[6]: = Ainv = Inverse[A]; MatrixForm[Ainv]

Out[6]//MatrixForm =

–2

1

3 ? 1?

?

2 ? 2?

© 2001 by CRC Press LLC

To verify whether or not the inverse matrix Ainv obtained in Output[6] indeed

satisfies the equations [A][A]–1 = [I] which is the identity matrix, we apply Mathematica for matrix multiplication:

In[7]: = Iden = A.Ainv; MatrixForm[Iden]

Out[7]//MatrixForm =

1

0

0

1

A dot is to separate the two matrices A and Ainv which is to be multiplied in that

order. Output[7] proves that the computed matrix, Ainv, is the inverse of A! It should

be noted that D and I are two reserved variables in Mathematica for the determinant

of a matrix and the identity matrix. In their places, here Dif and Iden are adopted,

respectively. For further testing, we show that [A][A]T is a symmetric matrix:

In[8]: = S = A.AT; MatrixForm[S]

Out[8]//MatrixForm =

5

11

11

25

And, the unknown vector {X} in the matrix equation [A]{X} = {R} can be

solved easily if {R} is given and [A]–1 are available:

In[9]: = R = {13,31}; X = Ainv.R

Out[9] = {5, 4}

The solution of x1 = 5 and x2 = 4 do satisfy the equations x1 + 2x2 = 13 and 3x1

+ 4x2 = 31.

TRANSFORMATION

OF

COORDINATE SYSTEMS, ROTATION,

AND

ANIMATION

Matrix algebra can be effectively applied for transformation of coordinate systems. When the cartesian coordinate system, x-y-z, is rotated by an angle z about

the z-axis to arrive at the system x-y-z as shown in Figure 2, where z and z axes

coincide and directed outward normal to the plane of paper, the new coordinates of

a typical point P whose coordinates are (xP,yP,zP) can be easily obtained as follows:

x ?P = OP cos(? P ? ? z ) = (OP cos ? P ) cos ? z + (OP sin ? P ) sin ? z

= x P cos ? z + y p sin ? z

y ?P = OP sin(? P ? ? z ) = (OP sin ? P ) cos ? z ? (OP cos ? P ) sin ? z

= x p sin ? z + y p sin ? z

© 2001 by CRC Press LLC

FIGURE 2. The cartesian coordinate system, x-y-z, is rotated by an angle z about the zaxis to arrive at the system x-y-z.

and

z?P = z P

In matrix notation, we may define {P} = [xP yP zP]T and {P'} = [xP' yP' zP']T and

write the above equations as {P'} = [Tz]{P} where the transformation matrix for a

rotation of z-axis by z is:

? cos ?z

[Tz ] = ??? sin ?z

?? 0

sin ?z

cos ?z

0

0?

0 ??

1 ??

(5)

In a similar manner, it can be shown that the transformation matrices for rotating

about the x- and y-axes by angles x and y, respectively, are:

?1

[Tx ] = ??0

??0

0

cos ?x

? sin ?x

0 ?

?

sin ?x ?

cos ?x ??

(6)

? sin ? y ?

?

0 ?

cos ? y ??

(7)

and

?cos ? y

?

Ty = ? 0

? sin ?

y

?

[ ]

© 2001 by CRC Press LLC

0

1

0

FIGURE 3. Point P whose coordinates are (xP,yP,zP) is rotated to the point P' by a rotation

of z.

It is interesting to note that if a point P whose coordinates are (xP,yP,zP) is rotated

to the point P' by a rotation of z as shown in Figure 3, the coordinates of P' can

be easily obtained by the formula {P'} = [Rz]{P} where [Rz] = [Tz]T. If the rotation

is by an angle x or y, then {P'} = [Rx]{P} or {P'} = [Ry]{P} where [Rx] = [Tx]T

and [Ry] = [Ty]T.

Having discussed about transformations and rotations of coordinate systems,

we are ready to utilize the derived formulas to demonstrate the concept of animation. Motion can be simulated by first generating a series of rotated views of

a three-dimensional object, and showing them one at a time. By erasing each

displayed view and then showing the next one at an adequate speed, a smooth

motion of the object is achievable to produce the desired animation. Program

Animate1.m is developed to demonstrate this concept of animation by using a

4 2 3 brick and rotating it about the x-axis by an angle of 25° and then

rotating about the y-axis as many revolutions as desired. The front side of the

block (x-y plane) is marked with a character F, and the right side (y-z plane) is

marked with a character R, and the top side (x-z plane) is marked with a character

T for helping the viewer to have a better three-dimensional perspective of the

rotated brick (Figure 4). The x-rotation prior to y-rotation is needed to tilt the top

side of the brick toward the front. The speed of animation is controlled by a

parameter Damping. This parameter and the desired number of y-revolutions,

Ncycle, are both to be interactively specified by the viewer (Figure 5).

© 2001 by CRC Press LLC

FIGURE 4. The characters F, R, and T help the viewer to have a better three-dimensional

perspective of the rotated brick.

FIGURE 5. The speed of animation is controlled by a parameter Damping. This parameter and

the desired number of y-revolutions, Ncycle, are both to be interactively specified by the viewer.

© 2001 by CRC Press LLC

FUNCTION ANIMATE1(NCYCLE,DAMPING)

Notice that the coordinates for the corners of the brick are defined in arrays xb,

yb, and zb. The coordinates of the points to be connected by linear segments for

drawing the characters F, R, ant T are defined in arrays xf, yf, and zf, and xr, yr,

and zr, and xt, yt, and zt, respectively.

The equations in deriving [Rx] ( = [Tx]T) and [Ry] ( = [Ty]T) are applied for xand y- rotations in the above program. Angle increments of 5 and 10° are arranged

for the x- and y-rotations, respectively. The rotated views are plotted using the new

coordinates of the points, (xbn,ybn,zbn), (xfn,yfn,zfn), etc. Not all of these new

arrays but only those needed in subsequent plot are calculated in this m file.

MATLAB command clg is used to erase the graphic window before a new

rotated view the brick is displayed. The speed of animation is retarded by the “hold”

loops in both x- and y-rotations involving the interactively entered value of the

parameter Damping. The MATLAB command pause enables Figure 4 to be read

and requires the viewer to press any key on the keyboard to commence the animation.

Notice that a statement begins with a % character making that a comment statement,

and that % can also be utilized for spacing purpose.

The xs and ys arrays allow the graphic window to be scaled by plotting them

and then held (by command hold) so that all subsequent plots are using the same

© 2001 by CRC Press LLC

FIGURE 6. Animation of a rotating brick.

scales in both x- and y-directions. The values in xs and ys arrays also control where

to properly place the texts in Figure 4 as indicated in the text statements.

QUICKBASIC VERSION

A QuickBASIC version of the program Animate1.m called Animate1.QB also

is provided. It uses commands GET and PUT to animate the rotation of the 4 3

2 brick. More features have been added to show the three principal views of the

brick and also the rotated view at the northeast corner of screen, as illustrated in

Figure 6.

The window-viewport transformation of the rotated brick for displaying on the

screen is implemented through the functions FNTX and FNTY. The actual ranges

of the x and y measurements of the points used for drawing the brick are described

by the values of V1 and V2, and V3 and V4, respectively. These ranges are mapped

onto the screen matching the ranges of W1 and W2, and W4 and W3, respectively.

The rotated views of the brick are stored in arrays S1 through S10 using the

GET command. Animation retrieves these views by application of the PUT command. Presently, animation is set for 10 y-swings (Ncycle = 10 in the program

Animate1.m, arranged in Line 600). The parameter Damping described in the

program Animate1.m here is set equal to 1500 (in Line 695).

© 2001 by CRC Press LLC

© 2001 by CRC Press LLC

1.6 PROBLEMS

MATRIX ALGEBRA

1. Calculate the product [A][B][C] by (1) finding [T] = [A][B] and then

[T][C], and (2) finding [T] = [B][C] and then [A][T] where:

?1

[A] = ?

?4

2

5

3?

6??

?6

[B] = ??4

??2

5?

3??

1??

? ?1

??3

[C ] = ?

?2 ?

?4??

2. Calculate [A][B] of the two matrices given above and then take the

transpose of product matrix. Is it equal to the product of [B]T[A]T?

3. Are ([A][B][C])T and the product [C]T[B]T[A]T identical to each other?

© 2001 by CRC Press LLC

4. Apply the QuickBASIC and FORTRAN versions of the program MatxAlgb to verify the results of Problems 1, 2, and 3.

5. Repeat Problem 4 but use MATLAB.

6. Apply the program MatxInvD to find [C]–1 of the matrix [C] given in

Problem 1 and also to ([C]T)–1. Is ([C]–1)T equal to ([C]T)–1?

7. Repeat Problem 6 but use MATLAB.

8. For statistical analysis of a set of N given data X1, X2, …, XN, it is often

necessary to calculate the mean, m, and standard deviation, 5, by use of

the formulas:

m=

1

(X + X2 + … + X N )

N 1

and

[

]

2

2

2

1

? = ?? ( X1 ? m ) + ( X 2 ? m ) + … + ( X N ? m ) ??

?

?N

0.5

Use indicial notation to express the above two equations and then develop

a subroutine meanSD(X,N,RM,SD) for taking the N values of X to

compute the real value of mean, RM, and standard deviation, SD.

9. Express the ith term in the following series in indicial notation and then

write an interactive program SinePgrm allowing input of the x value to

calculate sin(x) by terminating the series when additional term contributes

less than 0.001% of the partial sum of series in magnitude:

Sin x =

x1 x3 x 5

? + ?…

1! 3! 5!

Notice that Sin(x) is an odd function so the series contains only terms of

odd powers of x and the series carries alternating signs. Compare the

result of the program SinePgrm with those obtained by application of the

library function Sin available in FORTRAN and QuickBASIC.

10. Same as Problem 9, but for the cosine series:

Cos x = 1 ?

x2 x 4 x6

+

?

+…

2! 4! 6!

Notice that Cos(x) is an even function so the series contains only terms

of even powers of x and the series also carries alternating signs.

11. Repeat Problem 4 but use Mathematica.

12. Repeat Problem 6 but use Mathematica.

© 2001 by CRC Press LLC

GAUSS

1. Run the program GAUSS to solve the problem:

?1

?4

?

?? 7

2

5

8

3 ? ? x1 ? ? 2 ?

? ?

6 ?? ?x 2 ? = ?? 8 ??

10 ?? ?? x3 ?? ??14??

2. Run the program GAUSS to solve the problem:

?0

?4

?

?? 7

2

5

8

3? ? x1 ? ??1?

? ?

6?? ?x 2 ? = ?? 8 ??

9?? ?? x3 ?? ??14 ??

What kind of problem do you encounter? “Divided by zero” is the message! This happens because the coefficient associated with x1 in the first

equation is equal to zero and the normalization in the program GAUSS

cannot be implemented. In this case, the order of the given equations

needs to be interchanged. That is to put the second equation on top or

find below the first equation an equation which has a coefficient associated

with x1 not equal to zero and is to be interchanged with the first equation.

This procedure is called “pivoting.” Subroutine GauJor has such a feature

incorporated, apply it for solving the given matrix equation.

3. Modify the program GAUSS by following the Gauss-Jordan elimination

procedure and excluding the back-substitution steps. Name this new program GauJor and test it by solving the matrix equations given in Problems

1 and 2.

4. Show all details of the normalization, elimination, and backward substitution steps involved in solving the following equations by application of

Gaussian Elimination method:

4x1 + 2x2 – 3x3 = 8

5x1 – 3x2 + 7x3 = 26

–x1 + 9x2 – 8x3 = –10

5. Present every normalization and elimination steps involved in solving the

following system of linear algebraic equations by the Gaussian Elimination Method:

5x1 – 2x2 + 2x3 = 9, –2x1 + 7x2 – 2x3 = 9, and 2x1 – 2x2 + 9x3 = 41

© 2001 by CRC Press LLC

6. Apply the Gauss-Jordan elimination method to solve for x1, x2, and x3

from the following equations:

?0

?2

?

??4

?1? ? x1 ? ?1?

? ?

3 ?? ?x 2 ? = ??1??

7 ?? ?? x3 ?? ??1??

1

9

24

Show every normalization, elimination, and pivoting (if necessary) steps

of your calculation.

7. Solve the matrix equation [A]{X} = {C} by Gauss-Jordan method

where:

?3

?2

?

??4

2

5

1

1 ? ? x1 ? ??2 ?

? ?

?1?? ?x 2 ? = ?? ?3??

7 ?? ?? x3 ?? ?? 3 ??

Show every interchange of rows (if you are required to do pivoting before

normalization), normalization, and elimination steps by indicating the

changes in [A] and {C}.

8. Apply the program GauJor to solve Problem 7.

9. Present every normalization and elimination steps involved in solving the

following system of linear algebraic equations by the Gauss-Jordan Elimination Method:

5x1 – 2x2 + x3 = 4

–2x1 + 7x2 – 2x3 = 9

x1 – 2x2 + 9x3 = 40

10.

11.

12.

13.

14.

Apply the program Gauss to solve Problem 9 described above.

Use MATLAB to solve the matrix equation given in Problem 7.

Use MATLAB to solve the matrix equation given in Problem 9.

Use Mathematica to solve the matrix equation given in Problem 7.

Use Mathematica to solve the matrix equation given in Problem 9.

MATRIX INVERSION

1. Run the program MatxInvD for finding the inverse of the matrix:

?3

[A] = ??0

??2

© 2001 by CRC Press LLC

0

5

0

2?

0 ??

3??

2. Write a program Invert3 which inverts a given 3 ? 3 matrix [A] by using

the cofactor method. A subroutine COFAC should be developed for calculating the cofactor of the element at Ith row and Jth column of [A] in

term of the elements of [A] and the user-specified values of I and J. Let

the inverse of [A] be designated as [AI] and the determinant of [A] be

designated as D. Apply the developed program Invert3 to generate all

elements of [AI] by calling the subroutine COFAC and by using D.

3. Write a QuickBASIC or FORTRAN program MatxSorD which will

perform the addition and subtraction of two matrices of same order.

4. Write a QuickBASIC or FORTRAN program MxTransp which will

perform the transposition of a given matrix.

5. Translate the FORTRAN subroutine MatxMtpy into a MATLAB m file

so that by entering the matrices [A] and [B] of order L by M and M by

N, respectively, it will produce a product matrix [P] of order L by N.

6. Enter MATLAB commands interactively first a square matrix [A] and

then calculate its trace.

7. Use MATLAB commands to first define the elements in its upper right

corner including the diagonal, and then use the symmetric properties to

define those in the lower left corner.

8. Convert either QuickBasic or FORTRAN version of the program MatxInvD into a MATLAB function file MatxInvD.m with a leading statement

function [Cinv,D] = MatxInvD(C,N)

9. Apply the program MatxInvD to invert the matrix:

?1

[A] = ??5

??8

3

6

9

4?

7 ??

10 ??

Verify the answer by using Equation 1.

10. Repeat Problem 9 but by MATLAB operation.

11. Apply the program MatxInvD to invert the matrix:

??9

[A] = ?? ?3

???6

?1

?4

?7

?2 ?

?5??

?8??

Verify the answer by using Equation 1.

12. Repeat Problem 11 but by MATLAB operations.

13. Derive [Rx] and verify that it is indeed equal to [Tx]T. Repeat for [Ry] and

[Rz].

14. Apply MATLAB to generate a matrix [Rz] for ?z = 45° and then to use

[Rz] to find the rotated coordinates of a point P whose coordinates before

rotation are (1,–2,5).

© 2001 by CRC Press LLC

FIGURE 7. Problem 18.

15. What will be the coordinates for the point P mentioned in Problem 14 if

the coordinate axes are rotated counterclockwise about the z-axis by 45°?

Use MATLAB to find your answer.

16. Apply MATLAB to find the location of a point whose coordinates are

(1,2,3) after three rotations in succession: (1) about y-axis by 30°, (2)

about z-axis by 45° and then (3) about x-axis by –60°.

17. Change m file Animate1.m to animate just the rotation of the front (F)

side of the 4 2 3 brick in the graphic window.

18. Write a MATLAB m file for animation of pendulum swing1 as shown in

Figure 7.

19. Write a MATLAB m file for animation of a bouncing ball1 using an

equation of y = 3e–0.1xsin(2x + 1.5708) as shown in Figure 8.

20. Write a MATLAB m file for animation of the motion of crank-piston

system as shown in Figure 9.

21. Write a MATLAB m file to animate the vibrating system of a mass

attached to a spring as shown in Figure 10.

© 2001 by CRC Press LLC

FIGURE 8. Problem 19.

22. Write a MATLAB m file to animate the motion of a cam-follower system

as shown in Figure 11.

23. Write a MATLAB m file to animate the rotary motion of a wankel cam

as shown in Figure 12.

24. Repeat Problem 9 but by Mathematica operation.

25. Repeat Problem 11 but by Mathematica operation.

26. Repeat Problem 14 but by Mathematica operation.

27. Repeat Problem 15 but by Mathematica operation.

28. Repeat Problem 16 but by Mathematica operation.

© 2001 by CRC Press LLC

FIGURE 9. Problem 20.

© 2001 by CRC Press LLC

FIGURE 10. Problem 21.

FIGURE 11. Problem 22.

© 2001 by CRC Press LLC

FIGURE 12. Problem 23.

1.7 REFERENCE

1. Y. C. Pao, “On Development of Engineering Animation Software,” in Computers in

Engineering, edited by K. Ishii, ASME Publications, New York, 1994, pp. 851–855.

© 2001 by CRC Press LLC

Чтобы не видеть никакую рекламу на сайте, нужно стать

**VIP**-пользователем.

Это можно сделать совершенно бесплатно. Читайте подробности тут.