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

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

 0:5x

¼ 1:5 (E3

)

We now have an equivalent set of equations that have different coefficients and constants compared to

the original set, but which have the same exact solution:

 4x

þ 3x

¼ 1; (E1)

 4:5x

¼ 3:5; (E2

)

 0:5x

¼ 1:5: (E3

)

Our next aim is to remove the term x

from Equation (E3′) using (E2′). This is done by multiplying (E2′)

by 0.25 and subtracting the resulting equation from (E3′). We obtain a new equation for (E3′):

0:625x

¼ 0:625: (E3

)

Our new set of equivalent linear equations is now in triangular form:

 4x

þ 3x

¼ 1; (E1)

 4:5x

¼ 3:5; (E2

)

0:625x

¼ 0:625: (E3

)

Note that the Gaussian elimination procedure involves modification of only the coefficients of the variables

and the constants of the equation. If the system of linear equations is expressed as an augmented matrix,

similar transformations can be carried out on the row elements of the 3 × 4 augmented matrix.

Back substitution is now employed to solve for the unknowns. Equation (E3″) is the easiest to solve

since the solution is obtained by simply dividing the constant term by the coefficient of x

. Once x

known, it can be plugged into (E2

) to solve for x

. Then x

can be obtained by substituting the values of x

and x

into (E1). Thus,

0:625

¼ 1;

3:5 þ 4:5  1

¼ 2;

1 þ 4  2  3  1

¼ 3:

Next, we generalize this solution method for any syst em of four linear equations as

shown below:

þ a

¼ b

;

þ a

¼ b

;

þ a

¼ b

;

þ a

¼ b

A system of four linear equations can be represented by a 4 × 5 augmented matrix

that contains the coefﬁcients a

of the variables and the constant terms b

on the

right-hand side of the equations:

2.4 Gaussian elimination with backward substitution

aug A ¼

: b

The elements of a matrix can be modiﬁed based on certain rules or allowable

operations called elementary row transformations, which permit the followin g:

(1) interchanging of any two rows;

(2) multiplication of the elements of a row by a non-zero scalar quantity;

(3) addition of a non-zero scalar multiple of the elements of one row to corresponding

elements of another row.

A matrix that is modiﬁed by one or more of these elementary transformations is

called an equivalent matrix since the solution of the corresponding linear equations

remains the same. The purpose of Gaussian elimination when applied to the aug-

mented matrix is to reduce it to an upper triangular matrix. An upper triangular

matrix has non-zero elements on the diagonal and above the diagonal and has only

zeros located below the diagonal. It is the opposite for the case of a lower triangular

matrix, which has non-zero elements only on and below the diagonal.

The elimination process can be broken down into three overall steps:

(1) use a

to zero the elements in column 1 (x

coefﬁcients) in rows 2, 3, and 4;

(2) use a

to zero the elements below it in column 2 (x

coefﬁcients) in rows 3 and 4;

(3) use a

to zero the element below it in column 3 (x

coefﬁcient) in row 4.

Step 1

In step 1, the ﬁrst row is multiplied by an appropriate scalar quantity and then

subtracted from the second row so that a

→ 0. The ﬁrst row is then multiplied by

another pre-determined scalar and subtracted from the third row such that a

→ 0.

This procedure is again repeated for the fourth row so that a

→ 0. The ﬁrst row is

called the pivot row. The pivot row itself remains unchanged, but is used to change

the rows below it. Since a

is the element that is used to zero the other elements

beneath it, it is called the pivot element. It acts as the pivot (a stationary point), about

which the other elements are transformed. The column in which the pivot element

lies, i.e. the column in which the zeroing takes place, is called the pivot column.

To convert a

to zero, the elements of the ﬁrst row are divided by a

, multiplied

by a

, and then subtracted from corresponding elements of the second row. This is

mathematically expressed as

¼ a



1  j  5;

where the prim e indicates the new or modiﬁed value of the element. Hence, when

j = 1 we obtain

¼ a



¼ 0:

To create a zero at the position of a

, the ﬁrst row is multiplied by the quantity

and subtracted from the third row. This is mathematically expressed as

¼ a



1  j  5:

Systems of linear equations

Similarly, for the fourth row, the following elementary row transformation is

carried out:

¼ a



1  j  5:

The algorithm for step 1 is thus given by

for i = 2 to 4

for j = 1 to 5

¼ a



end j

end i

The resulting equivalent augmented matrix is given by

aug A ¼

: b

0 a

: b

0 a

: b

0 a

: b

;

where the primes indicate that the elements have undergone a single modiﬁcation.

Step 2

In step 2, the second row is now the pivot row and a

is the pivot element, which is

used to zero the elements a

and a

below it. Again, no changes occur to the pivot

row elements in this step. The elementary row transformation that results in a

! 0,

and the modiﬁcation of the other elements in row 3, is mathematically expressed as

¼ a



2  j  5;

where the double prime indicates the new or modiﬁed value of the element resulting

from step 2 transformations. These row operations do not change the value of the

zero element located in column 1 in row 3, and therefore performing this calculation

for j = 1 is unnecessary.

To zero the element a

, the second row is multiplied by the quantity a

and is

subtracted from the fourth row. This is mathematically expressed as

¼ a



2  j  5:

The algorithm for step 2 is thus given by

for i = 3 to 4

for j = 2 to 5

¼ a



end j

end i

The resulting equivalent augmented matrix is given by

2.4 Gaussian elimination with backward substitution

aug A ¼

: b

0 a

: b

00a

: b

00a

: b

;

where a single prime indicates that the corresponding elements have had a single

change in value (in step 1), while double primes indicate that the elements have been

modiﬁed tw ice (in steps 1 and 2).

Step 3

The pivot row shifts from the second to the third row in this step. The new pivot

element is a

. Only the last row is modiﬁed in this step in order to zero the a

element. The third row is multiplied by the quantity a

and subtracted from the

fourth row. This is mathematically expressed as

000

¼ a



3  j  5:

The algorithm for step 3 is thus

for i = 4 to 4

for j = 3 to 5

000

¼ a



end j

end i

The resulting equivalent augmented matrix is given by

aug A ¼

: b

0 a

: b

00a

: b

000a

000

: b

000

; (2:17)

where a triple prime indicates that the elements have been modiﬁed three times (in

steps 1, 2, and 3).

The algorithms for the three individual steps discussed above can be consolidated to

arrive at an overall algorithm for the Gaussian elimination method for a 4 × 5

augmented matrix, which can easily be extended to any n × (n + 1) augmented matrix.

In the algorithm, the a

i,n+1

elements correspond to the b

constants for 1 ≤ i ≤ n.

Algorithm for Gaussian elimination withou t pivoting

(1) Loop through k =1ton − 1; k speciﬁes the pivot row (and pivot column).

(2) Assign pivot = a

(3) Loop through i =(k +1)ton; i speciﬁes the row being transformed.

(4) Calculate the multiplier m ¼ a

=pivot:

(5) Loop through j = k to n +1;j speciﬁes the column in which the individual element

is modiﬁed.

(6) Calculate a

¼ a

 ma

Systems of linear equations

(7) End loop j.

(8) End loop i.

(9) End loop k.

Written in MATLAB programming style, the algorithm appears as follows:

for k = 1:n-1

pivot = a(k,k)

for row = k+1:n

m = a(row,k)/pivot;

a(row,k+1:n+1) = a(row,k+1:n+1) – m*a(k,k+1:n+1);

end

Note that the calculation of the zeroed element a(i,k) is excluded in the MATLAB

code above since we already know that its value is zero and is inconsequential

thereafter.

Now that the elimination step is complete, the back substitution step is employed

to obtain the solution. The equivalen t upper triangular 4 × 5 augmented matrix given

by Equation (2.17) produces a set of four linear equations in triangular form.

The solution for x

is most obvious and is given by

000

With x

known, x

can be solved for using the equation speciﬁed by the third row of

aug A in Equation (2.17):

 a

Proceeding with the next higher row, we can solve for x

 a

Finally, x

is obtained from

 a

The algorithm for back substitution can be generalized for any system of n linear

equations once an upper triangular augmented matrix is available, and is detailed below.

Algorithm for back substitution

(1) Calculate x

¼ b

(2) Loop backwards through i = n –1to1.

(3) Calculate x



j¼n

j¼iþ1



(4) End loop i.

Program 2.1 solves a system of n linear equations using the basic Gaussian

elimination with back substitution.

2.4 Gaussian elimination with backward substitution

MATLAB program 2.1

function x = Gaussianbasic(A, b)

% Uses basic Gaussian elimination with backward substitution to solve a set

% of n linear equations in n unknowns.

% Input Variables

% A is the coefﬁcient matrix

% b is the constants vector

% Output Variables

% x is the solution

Aug = [A b];

[n, ncol] = size(Aug);

x = zeros(n,1); %initializing the column vector x

% Elimination steps

for k = 1:n-1

pivot = Aug(k,k);

for row = k + 1:n

m = Aug(row,k)/pivot;

Aug(row,k+1:ncol ) = Aug(row,k+1:nco l) − m*Aug(k,k+1: ncol);

end

% Backward substitution steps

x(n)=Aug(n,ncol)/Aug(n,n);

fork=n-1:-1:1

x(k) = (Aug(k,ncol) - Aug(k,k+1:n)*x(k+1:n))/Aug(k,k);

% Note the inner product

end

Let’s use Gaussianbasic to solve the system of three equations (E1)–(E3)

solved by hand earlier in this section. For this example

A ¼

2 43

123

31 1

; b ¼

The following is typed at the MATLAB prompt:

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

b = [1; 4; 12]; % Column vector

x = Gaussianbasic(A, b)

During the elimination process, if any of the elements along the diagonal becomes

zero, eithe r the elimination method or the backward substitution procedure will fail

due to fatal errors arising from “division by zero.” In such cases, when a diagonal

element becomes zero, the algorithm cannot proceed further. In fact, signiﬁcant

problems arise even if any of the diagonal or pivot elements are not exactly zero but

Systems of linear equations

close to zero. Division by a small number results in magniﬁcation of round-off errors

which are carried through the arithmetic calculations and produce erroneous resul ts.

Example 2.3

Let’s solve the example problem

0:001x

þ 6:4x

¼ 4:801;

22:3x

 8x

¼ 16:3;

using the Gaussian elimination method. As in some of the examples in Chapter 1, we retain only four

significant figures in each number in order to emphasize the limitations of finite-precision arithmetic. The

exact solution to this set of two equations is x

= 1 and x

= 0.75. The augmented matrix is

0:001 6:4 : 4:801

22:3 8 : 16:3



In this example, the pivot element 0.001 is much smaller in magnitude than the other coefficients. A

single elimination step is performed to place a zero in the first column of the second row. The

corresponding matrix is

0:001 6:4 : 4:801

0 142700 : 107000



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

Gaussian elimination 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)

Here, 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



We use Gaussianbasic to solve the system of two equations:

A = [1.1 −0.3; 1 −1];

b = [311.3; 76];

x = Gaussianbasic(A, b)

360.6250

284.6250

The steady state naphthalene concentration in the lung (360.625 μM) is 92.5% that of the feed

concentration. In the liver, the steady state concentration is not as high (284.625 μM; approx. 73% that

of the feed concentration) due to the faster conversion rate in the liver cells of naphthalene to epoxide. If

R is increased, will the steady state concentrations in the liver and lung increase or decrease? Verify your

answer by solving Equations (2.3) and (2.4) for R = 0.8.

2.4 Gaussian elimination with backward substitution

The matrices are not exactly equivalent because the four-digit arithmetic results in a slight loss of

accuracy in the new values. Next, performing back substitution, we obtain

107000

142700

¼ 0:7498

and

4:801  6:4  0:7498

0:001

¼ 2:28:

While the solution for x

may be a sufficiently good approximation, the result for x

is quite inaccurate.

The round-off errors that arise during calculations due to the limits in precision produce a slight error in

the value of x

. This small error is greatly magnified through loss of significant digits during subtractive

cancellation and division by a small number (the pivot element). Suitable modifications of the elimination

technique to avoid such situations are discussed in Section 2.4.2.

2.4.2 Gaussian elimination with pivoting

The procedure of exchanging rows in the augmented matrix either to remove from

the diagonal zeros or small numbers is called partial pivoting or simply pivoting.To

achieve greater accuracy in the ﬁnal result, it is highly recommended that pivoting is

performed for every pivot row. The strategy involves placing the element with the

largest magnitude in the pivot position, such that the pivot element is the largest of

all elements in the pivot column, and using this number to zero the entries lying

below it in the same column. This way the multiplier used to transform elements in

the rows beneath the pivot row will always be less than or equal to one, which will

keep the propagation of round-off errors in check.

Algorithm for partial pivoting

(1) The pivot element is compared with the other elements in the same column. In

MATLAB, the max functi on is used to ﬁnd the largest number in the column.

(2) The row in the augmented matrix that contains the largest element is swapped with

the pivot row.

Using MATLAB

(1) The syntax for the max function is [Y I] = max(X), where X is the vector within

which the maximum element is searche d, Y is the maximum element of X,andI is the

index or location of the element of maximum value.

(2) Array indexing enables the use of an array as an index for another array. For

example, if a is a row vector that contains even numbers in sequential order from 2

to 40, and b is a vector containing integ ers from 1 to 10, and if we use b as the index

into a, the numeric values of the individual elements of b specify the position of the

desired elements of a. In our example the index b chooses the ﬁrst ten elements of a.

This is demonstrated below.

a = 2:2:30

24681012141618202224262830

b = 1:10

12345678910

Systems of linear equations

c = a(b)

2468101214161820

Program 2.1 is rewritten to incorporate partial pivoting. The program uses the max

function to determine the largest element in the pivot column, and also uses array

indexing to swap rows of the matrix.

MATLAB program 2.2

function x = Gaussianpivoting(A, b)

% Uses Gaussian elimination with pivoting to solve a set

% of n linear equations in n unknowns.

% Input Variables

% A is the coefﬁcient matrix

% b is the constants vector

% Output Variables

% x is the solution

Aug = [A b];

[n, ncol] = size(Aug);

x = zeros(n,1); % initializing the column vector x

% Elimination with Pivoting

for k = 1:n−1

[maxele, l] = max(abs(Aug(k:n,k)));

% l is the location of max value between k and n

p=k+l− 1; % row to be interchanged with pivot row

if (p ~= k) % if p is not equal to k

Aug([k p],:) = Aug([p k],:); % array index used for efﬁcient swap

end

pivot = Aug(k,k);

for row = k + 1:n

m = Aug(row,k)/pivot;

Aug(row,k+1:ncol ) = Aug(row,k+1:ncol) - m*Aug (k,k+1:ncol);

end

% Backward substitution steps

x(n)=Aug(n,ncol)/Aug(n,n);

fork=n-1:-1:1

x(k) = (Aug(k,ncol) - Aug(k,k+1:n)*x(k+1:n))/ Aug(k,k);

% Note the inner product

end

Example 2.3 (continued)

Let’s now swap the first and second rows so that the pivot element is now larger than all the other elements

in the pivot column. The new augmented matrix is

22:3 8 : 16:3

0:001 6:4 : 4:801



2.4 Gaussian elimination with backward substitution

After using the pivot row to transform the second row, we obtain

22:3 8 : 16:3

06:400 : 4:800



Again, due to finite precision, the original augmented matrix and the transformed matrix are not exactly

equivalent. Performing back substitution, we obtain

4:800

6:400

¼ 0:75

and

16:3 þ 8  0:75

22:3

¼ 1

We have arrived at the exact solution!

Although partial pivoting is useful in safeguarding against corruption of the results by

round-off error is some cases, it is not a foolproof strategy. If elements in a coefﬁcient

matrix are of disparate magnitudes, then multiplication operations with the relatively

larger numbers can also be problematic. The number of multiplication/division oper-

ations (ﬂoating-point operations) carried out during Gaussian elimination of an n × n

matrix A is O(n

/3), and the number of addition/subtraction operations carried out is

also O(n

/3). This can be compared to the number of ﬂoating-point operations that are

required to carry out the backward substitution procedure, which is O(n

), or the

number of comparisons that are needed to perform partial pivoting, which is also O

). Gaussian elimination entails many more arithmetic operations than that performed

by either backward substitution or the row comparisons and swapping steps. Gaussian

elimination is therefore the most computationally intensive step of the entire solution

method, posing a bottleneck in terms of computational time and resources. For

large matrices, the number of arithmetic operations that are carried out during the

elimination process is considerable and can result in sizeable accumulation of round-off

errors.

Several methods have been devised to resolve such situations. One strategy is

maximal or complete pivoting, in which the maximum value of the pivot element is

searched amongst all rows and columns (except rows above the pivot row and

columns to the left of the pivot column). The row and column containing the largest

element are then swapped with the pivot row and pivot column. This entails a more

complicated procedure since column interchanges require that the order of the vari-

ables be interchanged too. Permutation matrices, which are identity matrices that

contain the same row and column interchanges as the augmented matrix, must be

created to keep track of the interchanges performed during the elimination operations.

Another method that is used to reduce the magnitude of the round-off errors is

called scaled partial pivoting. In this method every element in the pivot column is ﬁrst

scaled by the largest element in its row. The element of largest absolute magnitude in

the pivot row p and in each row below the pivot row serves as the scale factor for that

row, and will be different for different rows. If s

¼ max

1jn



for p ≤ i ≤ n are the

scale factors, then the the pivot element, and therefore the new pivot row, is chosen

based on the following criteria:

¼ max

pin

The kth row is then interchanged with the current pivot row.

Systems of linear equations