381 - A17 - Linear Systems with Gaussian Elimination

John Travis

Mississippi College

Solving Linear Systems

Linear systems of equations arise from many applications...often with a huge number of equations and a corresponding huge number of unknowns. In general, a system of m linear equations in the n unknowns $x_1, x_2, ... , x_n$ is of the form

$ a_{11}x_1 + a_{12}x_2 + a_{13}x_3 + ... +  a_{1n}x_n  = b_1$ 

$ a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + ... +  a_{2n}x_n  = b_2$ 

$ a_{31}x_1 + a_{32}x_2 + a_{33}x_3 + ... +  a_{3n}x_n  = b_3$ 

...

$ a_{m1}x_1 + a_{m2}x_2 + a_{m3}x_3 + ... +  a_{mn}x_n  = b_m$ 

 and by collecting the given coefficients $ a_{kj} $ into a mxn matrix A

$ A = \begin{bmatrix} a_{11} & a_{12} & ... & a_{1n}\\ a_{21} & a_{22} & ... & a_{2n}\\ ... & ... & ... & ... \\ a_{m1} & a_{m2} & ... & a_{mn} \end{bmatrix} $

and the given numbers $b_k$ into a mx1 vector b gives the matrix equation

$Ax=b$.

Since solving a linear equation is a finite process, there should be no "trucation error" in general when devising algorithms to determine a solution.  So, the major source of error will arise up due to the use of floating point arithmetic...that is "roundoff error".  

Note, in this worksheet we will be running in Matlab/Octave mode so that all calculations are done using floating point arithmetic and so that any resulting errors in the answers generated will be due to the use of floating point arithmetic.  We will also concernn ourselves with only solving a square linear system written in matrix form as $Ax = b$  where $A$ is a nxn matrix that is given and $b$ is a similar nx1 vector that has also been provided.  Our goal will be to solve for the unknown vector $x$. 

Notice, for this diagonal matrix the first row can be written out as the linear equation

$2x_1 + 0x_2 + 0x_3 + 0x_4 + 0x_5 = 12$ or $2x_1 = 12$

and since there is only one variable to deal with then solving for $x_1 = 6$ is trivial. The second row can be written out as the linear equation 

$0x_1 - 3x_2 + 0x_3 + 0x_4 + 0x_5 = 11$ and thus $x_2 = -11/3$. 

Write out the actual equations for the other rows. 

Therefore, for diagonal matrices it is easy to develop an algorithm for generating the solution of $Ax=b$...

Go ahead and create your own diagonal system of size 8x8 and solve using the code above.

Ok. So, diagonal systems are trivial to solve.  But what about more complex systems?  

Consider a "lower triangular" system which might be something of the form

$2x_1                        = 4$

$-x_1 + 3x_2             = 9$

$ 5x_1 - 4x_2 + 6x_3 = 0$

This gives a matrix form $Ax=b$ with

$ A = \begin{bmatrix}2 & 0 & 0 \\-1 & 3 & 0 \\5 & -4 & 6\end{bmatrix}, x = \begin{bmatrix}x_1 \\x_2 \\x_3\end{bmatrix}, b=\begin{bmatrix}4 \\9 \\0\end{bmatrix}$

Notice that when written out, the first equation has only one variable, the second has only two, and the third has only three.  For larger lower triangular matrices this pattern would continue.  This leads us to succesively solve such a system by getting the value for $x_1$ and then using it to reduce the number of variables left in all the other equations by one.  Indeed, the second equation would now only have one unknown variable.  Find it and broadcast to the remaining equations and repeat giving only one new variable found for each new equation considered. One can easily follow this patter to show that the final solution to this lower triangular system is $x_1 = 2, x_2 = 11/3, x_3 = 7/9$.

This is known as "forward substitution".

For the following matrix, write out the actual equations proscribed using variables $x_1, x_2, x_3, x_4$.

Go ahead and create your own lower triangular system of size 6x6 and solve by modifying the code above.

Similarly, one can have a system with an "upper triangular matrix" that only has non-zero entries on and above the main diagonal.  We can solve such a system now by starting with the last equation and working our way back up through the other equations.  The last equation will only have one variable, the second to last will only have two, etc.  Doing this approach is called "backward solving" and is very similar to forward solving.

Go ahead and create your own upper triangular system of size 7x7 and solve by modifying the code above.

Next, let's consider a general nxn system which doesn not have any special structure (such as diagonal or triangular).  For these, we use a process known as Gaussian Elimination that utilizes the following general linear algebra properties:

1. Any equation can be multiplied by a non-zero number and the solution to the resulting system does not change.
2. Any two equations can be swapped and the solution of the resulting system does not change.
3. Any multiple of one equation can be added to any other equation and then replace the second equation by this sum and the solution of the resulting system does not change.

Using these repeatedly in a manner that aims for creating zero entries below the diagonal is the essence of Gaussian Elimination.  Once this is done, one will have a new upper triangular system that has the same solution as the original system.  Finish by backsolving.

First, consider a situation in which there is no need to swap any rows...

For the fun of it let's create a bigger A and use a given x to determine the vector b.  Then, pretend you don't know x and solve $Ax=b$.

Notice that we haven't used one of the allowable row operations...that is, we haven't ever swapped rows.  Sometimes this might be necessary or at least "adivsable".  See what happens when using the matrix provided below.  What is the problem and why does this fail?

This matrix did not perform so well due to the fact that one of the pivot elements turned out to be zero. Since we have to divide by that element to proceed, the algorithm bombs.

If we proceed blindly through Gaussian Elimination without considering swapping rows then it is possible that we can end up with a pivot element that equals 0 for which dividing in the algorithm is not possible.  However, if we can swap the row with a zero in the desired pivot spot with another row below that has a non-zero then we can be good to go.

Using the problem presented above, let's swap rows 2 and 3 (from the beginning) and see if that works.

To see what makes this different, write out the equations from the matrix problem given above and then write out the equations for the matrix problem given below.

And so the Gaussian Elimination algorithm worked if we were to know ahead of time where the zeros might pop up.

In reality and especially for larger systems, there is no way that we might know when this zero pivot event might occur.  Indeed, the algorithm itself should determine by itself when such a swap is necessary.  To do this will in theory require rows to be swapped dynamically as the algorithm progresses and doing such a swap is tedious.  An easier way to do swaps is by using pointers to the rows and simply swapping pointers instead of entire rows.

Initially, set up pointers so that p(1) = 1, p(2) = 2, p(3) = 3, ... , p(n) = n so that the kth row points to the physical kth row.  Then, to swap (say) rows 4 and 7 would require setting r(4)=7 and r(7)=4.  Watch how this works in the Gaussian Elimination with Pivoting algorithm below.

Note, the final output for the pointers shows where the kth row ended up logically.

Notice that when a zero is found in the (logically) next diagonal position, we swap with the first row below that row that has a non-zero in that same column.  

Gaussian Elimination with Scaled Partial Pivoting looks in the step-th column at ALL the rows below the step-th and to minimize roundoff selects the swap row to be the one with the largest magnitude term in that column.

Revise the algorithm to find the largest in absolute value element in the step-th column below the step-th and swap regardless if one with larger magnitude is found.  No need then to look for zeros.