Gaussian Elimination

The procedure just described has a major problem. While it may be satisfactory for hand computations with small systems, it is inadequate for a large system. Observe that the transformed coefficients can become very large as we convert to a triangular system.
The method that is called Gaussian elimination avoids this by subtracting $a_{i1}/a_{11}$ times the first equation from the $i^{th}$ equation to make the transformed numbers in the first column equal to zero. We do similarly for the rest of the columns.
We must always guard against dividing by zero. Observe that zeros may be created in the diagonal positions even if they are not present in the original matrix of coefficients. A useful strategy to avoid (if possible) such zero divisors is to rearrange the equations so as to put the coefficient of largest magnitude on the diagonal at each step. This is called pivoting.
Repeat the example of the previous section,

$\begin{displaymath} \left[ \begin{array}{rrrr} 4 & -2 & 1 &15 \\ -3 & -1 & 4 ... ...75 &19.25 \\ 0 & -0.5 & 2.75 & 9.25 \\ \end{array} \right], \end{displaymath}$

$\begin{displaymath} \begin{array}{r} \\ \\ R_3-(-0.5/-2.5)R_2 \rightarrow \... ...& 4.75 & 19.25 \\ 0 & 0 & 1.8 & 5.40 \\ \end{array} \right] \end{displaymath}$
The method we have just illustrated is called Gaussian elimination. (In this example, no pivoting was required to make the largest coefficients be on the diagonal.) Back-substitution, gives us
We shall obtain answers that are just close approximations to the exact answer because of round-off error. When there are many equations, the effects of round-off (the term is applied to the error due to chopping as well as when rounding is used) may cause large effects. In certain cases, the coefficients are such that the results are particularly sensitive to round off; such systems are called ill-conditioned.
if we had stored the ratio of coefficients in place of zero (we show these in parentheses), our final form would have been

$\begin{displaymath} \left[ \begin{array}{rrrr} 4 & -2 & 1 &15 \\ (-0.75) & -2... ...9.25 \\ (0.25) & (0.20) & 1.8 & 5.40 \\ \end{array} \right] \end{displaymath}$
The original matrix can be written as the product:

$\begin{displaymath} A=\underbrace{\left[ \begin{array}{rrr} 1 & 0 & 0 \\ -0.7... ... 0 & -2.5 & 4.75 \\ 0 & 0 & 1.8 \\ \end{array} \right]}_U \end{displaymath}$

This procedure is called a LU decomposition of A. In this case,

$\begin{displaymath} A =L *U \end{displaymath}$

where L is lower-triangular and U is upper-triangular.
Because the determinant of two matrices, $B \times C$ , is the product of each of the determinants, for this example we have

$\begin{displaymath} det(L*U) = det(L) * det(U)= det(U) \end{displaymath}$

because L is triangular and has only ones on its diagonal so that . Thus, for our example, we have

$\begin{displaymath} det(A) = det(U)=(4)*(-2.5)*(1.8)=-18 \end{displaymath}$

because U is upper-triangular and its determinant i is just the product of the diagonal elements.
When there are row interchanges

$\begin{displaymath} det(A) = (-1)^m*u_{11}*\ldots*u_{nn} \end{displaymath}$

where the exponent m represents the number of row interchanges.
Solve the following system of equations using Gaussian elimination. In addition, compute the determinant of the coefficient matrix and the decomposition of this matrix.

$\begin{displaymath} \begin{array}{rrrr} &2x_2& &+x_4=0\\ 2x_1 &+2x_2 &+3x_3 &... ...3 &+2x_4 =-2 \\ 2x_1 &+2x_2 &+3x_3 &+2x_4 =-2 \\ \end{array}\end{displaymath}$

The augmented coefficient matrix is

$\begin{displaymath} \left[ \begin{array}{rrrrr} 0 & 2 & 0 & 1 &0 \\ 2 & 2 & 3... ...& -3 & 0 & 1 &-7 \\ 6 & 1 &-6 &-5 &6 \\ \end{array} \right] \end{displaymath}$

We cannot permit a zero in the $a_{11}$ position because that element is the pivot in reducing the first column. We could interchange the first row with any of the other rows to avoid a zero divisor, but interchanging the first and fourth rows is our best choice. This gives

$\begin{displaymath} \left[ \begin{array}{rrrrr} 6 & 1 &-6 &-5 &6 \\ 2 & 2 & 3... ...& -3 & 0 & 1 &-7 \\ 0 & 2 & 0 & 1 &0 \\ \end{array} \right] \end{displaymath}$

$\begin{displaymath} \left[ \begin{array}{rrrrr} 6 & 1 &-6 &-5 &6 \\ 0 &1.6667... ... & 4 &4.3333&-11 \\ 0 & 2 & 0 & 1 &0 \\ \end{array} \right] \end{displaymath}$

We again interchange before reducing the second column, not because we have a zero divisor, but because we want to preserve accuracy. Interchanging the second and third rows puts the element of largest magnitude on the diagonal.

$\begin{displaymath} \left[ \begin{array}{rrrrr} 6 & 1 &-6 &-5 &6 \\ 0 &-3.666... ...67& 5 &3.6667&-4 \\ 0 & 2 & 0 & 1 &0 \\ \end{array} \right] \end{displaymath}$

Now we reduce in the second column

$\begin{displaymath} \left[ \begin{array}{rrrrr} 6 & 1 &-6 &-5 &6 \\ 0 &-3.666... ...1\\ 0 & 0 & 2.1818 & 3.3636 &-5.9999 \\ \end{array} \right] \end{displaymath}$

No interchange is indicated in the third column. Reducing, we get

$\begin{displaymath} \left[ \begin{array}{rrrrr} 6 & 1 &-6 &-5 &6 \\ 0 &-3.666... ...9.0001\\ 0 & 0 & 0 & 1.5600 &-3.1199 \\ \end{array} \right] \end{displaymath}$

Back-substitution gives . The correct answers are . In this calculation we have carried five significant figures and rounded each calculation. Even so, we do not have five-digit accuracy in the answers. The discrepancy is due to round off.
In this example, if we had replaced the zeros below the main diagonal with the ratio of coefficients at each step, the resulting augmented matrix would be

$\begin{displaymath} \left[ \begin{array}{rrrrr} 6 & 1 &-6 &-5 &6 \\ (0.66667)... ... (-0.54545) & (0.32) & 1.5600 &-3.1199 \\ \end{array} \right] \end{displaymath}$

This gives a LU decomposition as

$\begin{displaymath} \left[ \begin{array}{rrrr} 1&0&0&0\\ 0.66667 & 1 & 0 & 0 ... ... & 6.8182 &5.6364\\ 0& 0 &0 & 1.5600 \\ \end{array} \right] \end{displaymath}$

It should be noted that the product or the matrices produces a permutation of the original matrix, call it , where

$\begin{displaymath} A'=\left[ \begin{array}{rrrr} 6 & 1 &-6 &-5 \\ 4 & -3 & 0... ... \\ 2 & 2 & 3 & 2 \\ 0 & 2 & 0 & 1 \\ \end{array} \right] \end{displaymath}$

The determinant of the original matrix of coefficients can be easily computed according to the formula

$\begin{displaymath} det(A)=(-1)^2*(6)*(-3.6667)*(6.8182)*(1.5600)=-234.0028 \end{displaymath}$

which is close to the exact solution: -234. The exponent 2 is required, because there were two row interchanges in solving this system. To summarize
1. The solution to the four equations
2. The determinant of the coefficient matrix
3. A decomposition of the matrix, , which is just the original matrix, , after we have interchanged its rows in the process.
```
>> A = [0 2 0 1 0 ;  2 2 3 2  -2 ;  4 -3 0 1 -7 ;  6 1 -6 -5 6]
>> [L, U, P]=lu(A)
>> b = [ 0 -2 -7 6]
>> A\b
ans =
-0.5000
 1.0000
 0.3333
-2.0000
```

2004-11-03