Preface

---

This is a tutorial made solely for the purpose of education and it was designed for students taking Applied Math 0340. It is primarily for students who have some experience using Mathematica. If you have never used Mathematica before and would like to learn more of the basics for this computer algebra system, it is strongly recommended looking at the APMA 0330 tutorial. As a friendly reminder, don'r forget to clear variables in use and/or the kernel.

Finally, the commands in this tutorial are all written in bold black font, while Mathematica output is in normal font. This means that you can copy and paste all commands into Mathematica, change the parameters and run them. You, as the user, are free to use the scripts for your needs to learn the Mathematica program, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately.

Elementary Row Operations

---

Our objective is to find an efficient method of finding the solutions of systems of linear equations. There are three operations that we can perform on the equations of a linear system without altering the set of solutions:

• multiply both sides of an equation by a non-zero constant;
• interchange two equations;
• add a multiple of one equation to another.

These operations, usually called elementary row operations, do not alter the set of solutions of the linear system \( {\bf A}\,{\bf x} = {\bf b} \) since the restrictions on the variables \( x_1 , \, x_2 , \, \ldots , \,x_n \) given by the new equations imply the restrictions given by the old ones. At the same time, these operations really only involve the coefficients of the variables and the right-hand sides of the equations, b.

Since the coefficient matrix A together with the given constant vector b contains all the information we need to use, we unite them into one matrix, called the augmented matrix:

Now rather than manipulating the equations, we can instead manipulate the rows of this augmented matrix. These observations form the motivation behind a method to solve systems of linear equations, known as Gaussian elimination, called after its inventor Carl Friedrich Gauss (1777--1855) from Germany. Various modifications have been made to the method of elimination, and we present one of these, known as the Gauss--Jordan method. It was the German geodesist Wilhelm Jordan (1842-1899) and not a French mathematician Camille Jordan (1838-1922) who introduced the Gauss-Jordan method of solving systems of linear equations. Camille Jordan is credited for Jordan normal form, a well known linear algebra topic.

Elimination produces an upper triangular system, called row echelon form for Gauss elimination and reduced row echelon form for Gauss--Jordan algorithm. The Gauss elimination introduces zeroes below the pivots, while Gauss--Jordan algorithm contains additional phase in which it introduces zeroes above the pivots. Since this phase involves roughly 50% more operations than Gaussian elimination, most computer algorithms are based on the latter method.

Suppose M = [A | b] is the augmented matrix of a linear system A x = b. We can solve the linear system by performing elementary row operations on M. In matlab, these row operations are implemented with the following functions.

Example. Consider the system of linear equations

First, we form the augmented matrix The idea of the elimination procedure is to reduce the augmented matrix to equivalent "upper triangular" matrix. So the first question is how to determine pivots. Depending on this choice, we get the corresponding row echelon form. Therefore, the Gaussian algorithm may lead to different row echelon forms; hence, it is not unique.

From computational point of view, the pivot should be the largest element. If we look at the first column, the largest entry is in the second row. To move it into the first row and choose "4" as a pivot, we need to swap first two rows. It is convenient, but not necessarily, to have 1 in this position, which we call pivot.

Another option is to choose "2" in the first row as pivot because it is not zero. We want to start with a 1 in the first position of the first column. To chieve it, we can either divide each entry in the first row by 2 or switch the first and the third rows because the last row has "1" in the first position already. If we stick with the latter decision, we may want to begin by scaling the first row by 1/2.

The first argument is the index of the row to be replaced, the second argument is the row to form a multiple of, and the final argument is the scale factor. Thus, M.add_multiple_of_row(n,m,A) would replace row n with (row n)+A*(row m).

or without scaling

Since the last entry in the first column is already zero we can move on to the second column. Our first step is to get a 1 in the second position of the second column. Normally we would do this by scaling the row but in this case we can swap the second and third row

Now we want to eliminate the 5 below the 1. This is done by multiplying the second row by −5 and adding it to the third row.

To get a 1 as the last entry in the third column we can scale the third row by −1/2

At this point, the matrix is in echelon form (well, having the 1's down the diagonal of the matrix is not required, but it does make our work easier). All that remains to find the solution is to put the matrix in reduced echelon form which requires that we replace all entries in the first three columns that are not on the main diagonal (where the 1's are) with zeros. We will start with the third column and work our way up and to the left. Remember that this is an augmented matrix and we are going to ignore the right-most column; it just "goes along for the ride."

Therefore, our solution is complete. We see that the solution to our linear system is
\( x_1 =2, \ x_2 = -2, \ \mbox{and} \ x_3 = 3. \) ■

There is an easy way to check your work, or to carry out these steps in the future. First, let's reload the matrix M.

Theorem: Let A be a \( n \times n \) matrix. Then, if B is a matrix that results from switching one row with another row, then \( \det {\bf A} = - \det {\bf B} . \)

Theorem: Let A be a \( n \times n \) matrix. Then, if B is a matrix that results from adding or subtracting a multiple of one row to another row, then \( \det {\bf A} = \det {\bf B} . \) (The same is true for column operations also.)

Theorem: Let A be a \( n \times n \) that is upper triangular, lower triangular, or diagonal, then its determinant is the product of diagonal entries. ■

Example. We illustrate how to find an inverse matrix using row operations. Consider the matrix

To accomplish this goal, we add a block of the identity matrix and apply row operation to this augmented matrix \( \left[ {\bf A} \,| \,{\bf I} \right] \) until the left side is reduced to I. These operatioons will convert the right side to \( {\bf A}^{-1} ,\) so the final matrix will have the form \( \left[ {\bf I} \,| \,{\bf A}^{-1} \right] .\) Thus,

Naïve Gaussian Elimination Method

The naïve Gaussian Elimination method includes two stages:

• Forward Elimination of Unknowns: the unknown is eliminated in each equation starting with the first equation. This way, the equations are reduced to one equation and one unknown in each equation.
• Back Substitution: starting from the last equation, each of the unknowns is found.

We demonstrate how the method works in the following example.

Example: Consider four equations with four unknowns. First, we enter data.

Then we define the right-hand side vector: Remember that we write vector b in column form, but Mathematica keeps the its enries as a list. Next, we build the augmeted matrix from A and b: Now we perform forward elimination that consists of (n-1) steps. In each step k, the coefficient of the kth unknown will be zeroed from every subsequent equation that follows the kth row. For example, in step 2 (i.e., k = 2), the coefficient of x₂ will be zeroed from rows 3..n. With each step that is conducted, a new matrix is generated until the coefficient matrix is transformed to an upper triangular matrix. The following procedure calculates the upper triangular matrix produced for each step k. The new upper triangular coefficient matrix, A1, can be extracted from the final augmented matrix Aug: Notice that the final row, n, has only one unknown to be solved for. The new right-hand side array/vector, b1, is This is the end of the forward elimination stage. The new upper triangular coefficient matrix A1 and right hand side array b1 permit solving for the solution vector using backward substitution. It begins with solving the last equation as it has only one unknown: The remaining equations can be solved by using the following formula:

First, we define the vector x:

Solving for the nth equation as it has only one unknown: Solving for the remaining n-1 unknowns working backwards from the (n-1)th equation to the first equation: The solution vector X is Using Mathematica's built-in tools, the exact solution is given below. which gives exactly the same solution. ■

Are there any pitfalls of the Naïve Gauss elimination method?

Yes, there are two pitfalls of the Naïve Gauss eimination method.
Division by zero: It is possible for division by zero to occur during the beginning of the n - 1 steps of forward elimination.

Round-off error: The Naïve Gauss elimination method is prone to round-off errors. This is true when there are large numbers of equations as errors propagate. Also, if there is subtraction of numbers from each other, it may create large error. ■

 

 

Mathematica solves linear systems

Elementary Row Operations

Row Echelon Forms

LU-factorization

PLU Factorization

Reflection

Givens rotation

Row Space and Column Space

Null Space or Kernel

Rank and Nullity

Square Matrices

Cramer's Rule

Symmetric and self-adjoint matrices

Cholesky decomposition

Projection Operators

Gram-Schmidt

QR-decomposition

Least Square Approximation

SVD Factorization

Numerical Methods

Markov chains

Inequalities

Miscellaneous

 

Return to Mathematica page

Return to the main page (APMA0330)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Equations
Return to the Part 3 Linear Systems of Ordinary Differential Equations
Return to the Part 4 Non-linear Systems of Ordinary Differential Equations
Return to the Part 5 Numerical Methods
Return to the Part 6 Fourier Series
Return to the Part 7 Partial Differential Equations