King M.R., Mody N.A. Numerical and Statistical Methods for Bioengineering: Applications in MATLAB

Подождите немного. Документ загружается.

For large matr ices, the number of multiplication/division operations and addi-

tion/subtraction operations required to carry out the eliminat ion routine is rather

formidable. For large sparse matrices , i.e. matrices that consist of many zeros,

traditional Gaussian elimination results in many wasteful ﬂoating-point operations

since addition or subtraction of any numb er with zero or multiplication of two zeros

are unnecessary steps. Tridiagonal matrices are sparse matrices that have non-zero

elements only along the main diagonal, and one diagonal above and below the main,

and are often encountered in engineering problems (see Problems 2.7– 2.9). Such

matrices that have non-zero elements only along the main diagonal and one or more

other smaller diagonals to the left/right of the main diagonal are called banded

matrices. Special elimination methods, such as the Thomas method, which are faster

than O(n

) are employed to reduce tridiagonal (or banded) matrices to the appro-

priate triangular form for performing back substitution.

2.5 LU factorization

Analysis of the number of ﬂoating-point operations (ﬂops) carried out by a numerical

procedure allows one to compare the efﬁciency of different computational methods

in terms of the computational time and effort expended to achieve the same goal. As

mentioned earlier, Gaussian elimination requires O(n

/3) multiplication/division oper-

ations and O(n

/3) addition/subtraction operations to reduce an augmented matrix

that represents n equations with n unknowns to a upper triangular matrix (see

Section 1.6.1 for a discussion on “order of magnitude”). Back substitution is an O

) process and is thus achieved much more quickly than the elimination method for

large matrices. If we wish to solve systems of equations whose coefﬁcient matrix A

remains the same but have different b vectors using the Gaussian elimination method,

we would be forced to carry out the same operations to reduce the coefﬁcients in A in

the augmented matrix every time. This is a wasteful process considering that most of

the computational work is performed during the elimination process.

This problem can be addressed by decoupling the elimination steps carried out on

A from that on b . A can be converted into a product of two matrices: a lower

triangular matrix L and an upper triangular matrix U, such that

A ¼ LU:

The matrix equation Ax = b now becomes

LUx ¼ b: (2:18)

Once A has been factored into a product of two triangular matrices, the solution is

straightforward. First, we substitute

y ¼ Ux (2:19)

into Equation (2.18) to obtain

Ly ¼ b : (2:20)

In Equation (2.20), L is a lower triangular matrix, and y can be easily solved for by the

forward substitution method, which is similar to the backward substitution process

except that the sequence of equations used to obtain the values of the unknowns starts

with the ﬁrst equation and proceeds sequentially downwards to the next equation.

Once y is known, Equation (2.19) can easily be solved for x using the backward

substitution method. This procedure used to solve linear systems of equations is called

2.5 LU factorization

the LU factorization method or LU decomposition method, and is well-suited to solving

problems with similar coefﬁcient matrices. There are three steps involved in the LU

factorization method as compared to the two-step Gauss elimination method.

(1) Factorization of A into two triangular matrices L (lower triangular) and U

(upper triangular).

(2) Forward substitution method to solve Ly = b.

(3) Backward substitution method to solve Ux = y.

There are two important questions that arise here: ﬁrst, how can A be factored

into L and U, and second how computationally efﬁcient is the factorization step?

The answer to both questions lies with the Gaussian elimination method, for it is

this very procedure that points us in the direction of how to obtain the matrices

L and U. As discussed in Section 2.4, Gaussian elimination converts A into an

upper triangular matrix U. The technique used to facilitate this conversion process,

i.e. the use of multipliers by which the pivot row is added to another row, produces

the L matrix. Thus, factoring A into L and U is as computationally efﬁcient as

performing the Gaussian elimination procedure on A. The advantage of using LU

factorization stems from the fact that once A is factorized, the solution of any

system of equations for which A is a coefﬁcient m atrix becomes a two-step process

of O(n

There are several different methods that have been developed to perform the LU

factorization of A. Here, we discuss the Doolittle method, a technique that is easily

derived from the Gaussian elimination process. In this factorization method, U is

similar to the upper triangular matrix that is obtained as the end result of Gaussian

elimination, while L has ones along its diagonal, and the multipliers used to carry out

the elimination are placed as entries below the diagonal. Before we learn how to

implement the Doolittle method, let’s familiarize ourselves with how this method is

derived from the Gaussian elimination method. We will ﬁrst consider the case in

which pivoting of the rows is not performed.

2.5.1 LU factorization without pivoting

Recall that the Gaussian elimination method involves zeroing all the elements below

the diagonal, starting from the ﬁrst or leftmost to the last or rightmost column.

When the ﬁrst row is the pivot row, the elements in the ﬁrst column below the

topmost entry are zeroed, and the other elements in the second to nth rows are

transformed using the following formula:

¼ a



1  j  n; 2  i  n;

where the prime indicates the new or modiﬁed value of the element. In this formula

is called the multiplier since its products with the elements of the pivot row are

correspondingly subtracted from the elements of the other rows located below. The

multiplier has two subscripts that denote the row and the column location of the

element that is converted to zero after applying the transformation.

Systems of linear equations

The transformation operations on the elements of the ﬁrst column of the coef-

ﬁcient matrix can be compactly and efﬁciently represented by a single matrix multi-

plication operation. Let A be a 4 × 4 non-singular matrix. We create a matrix M

that

contains ones along the diagonal (as in an identity matrix) and negatives of the

multipliers situated in the appropriate places along the ﬁrst column of M

. Observe

that the product of M

and A produces the desired transformati on:

A ¼

1 000

m

100

m

010

m

001

0 a

Similarly, matrices M

and M

can be created to zero the elements of the coef-

ﬁcient matrix in the second and third columns. The matrix multiplications are

sequentially performed to arrive at our end prod uct – an upper triangular matrix:

ðM

AÞ¼

1000

0100

0 m

0 a

00a

;

ðM

AÞ¼

10 0 0

01 0 0

00 1 0

00m

0 a

00a

0 a

00a

000a

000

¼ U:

Thus, we have

U ¼ M

Left-multiplying with the inverse of M

on both sides of the equation, we get

ðM

1

ÞM

A ¼ M

1

Again, left-multiplying with the inverse of M

on both sides of the equation

1

¼ IÞ (see Equation (2.15)), we get

1

A ¼ M

1

Since IM

¼ M

(see Equation (2.15))

A ¼ M

1

Finally, left-multip lying with the inverse of M

on both sides of the equation, we get

1

A ¼ A ¼ M

1

U ¼ LU;

where

L ¼ M

1

2.5 LU factorization

Now let’s evaluate the product of the inverses of M

to verify if L is indeed a lower

triangular matrix. Determining the product M

1

will yield clues on how to

construct L from the multipliers used to perform the elimination process.

The inverse of

1000

m

100

m

010

m

001

is simply (calculate this yourself using Equation (2.16))

1

1 000

100

010

001

Similarly,

1

1000

0100

0 m

; M

1

10 0 0

01 0 0

00 1 0

00m

Now,

1

1000

100

010

001

1000

0100

0 m

1000

100

and

ðM

1

ÞM

1

1000

100

10 0 0

01 0 0

00 1 0

00m

1000

100

¼ L:

Thus, L is simply a lower triangular matrix in which the elements below the

diagonal are the multipliers whose positions coincide with their subscript nota-

tions. A multiplier used t o zero the entry in the second row, ﬁrst column is placed

in the second row and ﬁrst column in L. While carrying out the transformation

A → U, no extra effort is expended in producing L since the multipliers are

calculated along the way.

As mentioned earlier, the determinant of any triangular matrix is easily obtained

from the product of all the elements along the diagonal. Since L has only ones on its

diagonal, its determinant, L

¼ 1, and U will have a non-zero determinant only if all

Systems of linear equations

the elements along the diagonal are non-zero. If even one element on the diagonal is

zero, then U

¼ 0. Since A

¼ L

, the determinant of A can be obtained by

simply multiplying the elements along the diagonal of U.If A

¼ 0, a unique

solution of the existing set of equations cannot be determined.

Example 2.4

We use the LU factorization method to solve the system of equations given by Equations (E1)–(E3). The

coefficient matrix is

A ¼

2 43

123

311

LU factorization/decomposition

We begin the procedure by initializing the L and U matrices:

L ¼

100

010

001

; U ¼

2 43

123

311

We divide the factorization procedure into two steps. In the first step we zero the elements in the first

column below the pivot row (first row), and in the second step we zero the element in the second column

below the pivot row (second row).

Step 1: The multipliers m

and m

are calculated to be

; m

The transformation

¼ u

 m

; 1  j  3; 2  i  3;

produces the following revised matrices:

L ¼

100

and U ¼

2 43

04

07

Step 2: The multiplier m

is calculated as

The transformation

¼ u

 m

; 2  j  3; i ¼ 3;

produces the following revised matrices:

L ¼

100

and U ¼

2 43

04

The LU factorization of A is complete. If the matrix product LU is computed, one obtains the original

matrix A. Next, we use the forward and backward substitution methods to obtain the solution.

2.5 LU factorization

Forward substitution

First, we solve the equation Ly = b:

100

This matrix equation represents the following three equations:

¼ 1; (E4)

þ y

¼ 4; (E5)

þ y

¼ 12: (E6)

Equation (E4) gives us y

directly. Substituting the value of y

into Equation (E5), we get y

= 3.5. Equation

(E6) then gives us y

= 35/8.

Backward substitution

Next we solve the equation Ux = y:

2 43

04

This matrix equation represents the following three equations:

 4x

þ 3x

¼ 1; (E7)



; (E8)

: (E9)

Equation (E9) directly yields x

= 1. Working in the upwards direction, we use Equation (E8) to solve for

= 2, and then use Equation (E7) to obtain x

=3.

In the above example, L and U were created as separate matrices. The storage and

manipulation of the elements of L and U can be achieved most efﬁci ently by housing

the elements in a single matrix A. Matrix elements of A can be transformed in situ as

in Gaussian elimination to yield U. The multipliers of L can be stored beneath the

diagonal of A, with the ones along the diagonal being implied. For a 4 × 4 non-

singular matrix A, the factorized matr ices L and U can be stored in a single matrix

given by

000

Systems of linear equations

2.5.2 LU factorization with pivoting

Situations can arise in which row interchanges are necessary either to limit the round-off

error generated during the factorization or substitution process, or to remove a zero from

the diagonal of U. If pivoting is performed during the factorization process, information

on the row swapping performed mustbestoredsothatelementsoftheb vector can be

permuted accordingly. Permutation matrices, which are identity matrices with permuted

rows (or columns), are used for this purpose. A permutation matrix has only one non-

zero entry equal to one in any column and in any row, and has a determinant equal to

one. Pre-multiplication of A or b by an appropriate permutation matrix P effectively

swaps rows in accordance with the structure of P. (Post-multiplication of any matrix with

a permutation matrix will result in swapping of columns.)

For example, if permutation matrix P is an identity matrix in which the ﬁrst two

rows have been interchanged, then

P ¼

0100

1000

0010

0001

and

PA ¼

0100

1000

0010

0001

Here, the pre-multiplication operation with P interchanged the ﬁrst two rows of A.

The permutation matrix P is generated during the factorization process. Prior to

any row swapping operations, P is initialized as an identity matrix, and any row

swaps executed while creating L and U are also performed on P.

Example 2.5

Here Equations (E1)–(E3) are solved again using LU factorization with pivoting. Matrices P and A are

initialized as follows:

P ¼

100

010

001

; A ¼

2 43

123

311

The largest entry is searched for in the first column or pivot column of A. Since a

, rows 1 and

3 are swapped. We also interchange the first and third rows of P to obtain

P ¼

001

010

100

; A ¼

311

123

2 43

Elimination is carried out using the first row, which is the pivot row. The multipliers are stored in

the matrix A itself below the diagonal in place of the zeroed elements. The multipliers used to zero the

two entries below the diagonal in the first column are m

¼ 1=3 and m

¼ 2=3.

We obtain

A ¼

31 1



2.5 LU factorization

Next, the pivot column shifts to the second column and the largest entry is searched for in this

column. Since a

, rows 2 and 3 are interchanged in both A and P to obtain

P ¼

001

100

010

A ¼

31 1



Note that m

is calculated as 5=14. Transformation of the third row using the pivot (second) row yields

A ¼

311



The LU factorization process is complete and the matrices L and U can be extracted from A:

L ¼

100



and U ¼

311

0 

00

Since A is permuted during the factorization process, the LU factorization obtained does not correspond

to A but to the product PA. In other words,

LU ¼ PA

and

LUx ¼ PAx ¼ Pb

Here,

Pb ¼

We use the forward substitution process to solve Ly = Pb, which represents the following set of equations:

¼ 12;

þ y

¼ 1;



þ y

¼ 4;

which has the solution y

= 12, y

= −7, and y

= −5/2.

The equation set Ux = y is solved via backward substitution:

þ x

¼ 12;



¼7;



¼

;

which has the solution x

=3,x

=2,x

=1.

Systems of linear equations

2.5.3 The MATLAB lu function

MATLAB provides a built-in function to resolve a matrix A into two triangular

matrices: a lower triangular matrix L and an upper triangular matrix U. The

function lu is versatile; it can factorize both square and rectangular matrices, and

can accept a variety of parameters. A commonly used syntax for the lu function is

discussed here. For more information on the lu function, type help lu in the

Command Window. The lu functi on uses partial pivoting and therefore outputs the

permutation matrix P along with the factorized matrices L and U.

One syntax form is

[L,U,P] = lu(A)

where LU = PA. Let

A=[2−43;12−3;311]

2 −43

12−3

311

Applying the lu function to factorize A,

[L, U, P] = lu(A)

gives the following result:

1.0000 0 0

0.6667 1.0000 0

0.3333 −0.3571 1.0000

3.0000 1.0000 1.0000

0 −4.6667 2.3333

00−2.5000

001

100

010

Another syntax form for the lu function is

[L,U] = lu(A)

where A = L

U, and L

= P

−1

[L,U] = lu(A)

gives the following result:

0.6667 1.0000 0

0.3333 −0.3571 1.0000

1.0000 0 0

3.0000 1.0000 1.0000

0 −4.6667 2.3333

00−2.5000

2.5 LU factorization

2.6 The MATLAB backslash (\) operator

Solutions to systems of linear equations can be obtained in MATLAB with the use

of the backslash operator. If we desire the solution of the matrix equation Ax = b,

where A is an n × n matrix, we can instruct MATLAB to solve the equation by

typing

x = A\b.

The backslash operator, also known as the “matrix left-division operator,” can be

interpreted in several ways, such as “b divided by A”or“b multiplied on the left by

the inverse of A.” However, MATLAB does not explicitly compute A

−1

when solving

a linea r system of equations, but instead uses efﬁcient algorithms such as Gaussian

Box 2.1C Drug development and toxicity studies: animal-on-a-chip

LU decomposition is used to solve the system of equations specified by Equations (2.3) and (2.4),

which are reproduced below:

2  1:5RðÞC

lung

 0 :5RC

liver

¼ 779:3  780R; (2:3)

lung

 C

liver

¼ 76: (2:4)

We choose R=0.6, i.e. 60% of the exit stream is recycled back into the chip.

Equation (2.3) becomes

1:1C

lung

 0: 3C

liver

¼ 311:3:

For this problem,

A ¼

1:1 0:3

1 1



; b ¼

311:3



A = [1.1 −0.3; 1 −1]; b= [311.3; 76];

[L, U, P] = lu(A)

1.0000 0

0.9091 1.0000

1.1000 −0.3000

0 −0.7273

On performing forward (Ly = b) and backward (Ux = y) substitutions by hand (verify this yourself),

the solution

lung

¼ 360:6uM;

liver

¼ 284:6uM

is obtained.

Systems of linear equations