Consider the following general form of a single linear equation in the single unknown (variable) \(x\text{,}\) here the \(a\) and \(b\) represent any real numbers, \begin{equation*} ax = b. \end{equation*} When given actual values for \(a\) and \(b\text{,}\) does this equation always have a solution? If you think carefully about how we solve such a system then it is easy to see that the answer is “no”. One would solve this by dividing both sides of the equation by \(a\text{,}\) or equivalently multiplying both sides of the equation by \(a^{-1}\) like so: \begin{align*} ax \amp = b\\ a^{-1} ax \amp = a^{-1}b \qquad \text{ multiplicative inverses exist}\\ (a^{-1} a)x \amp = a^{-1}b \qquad \text{ associativity}\\ 1\cdot x \amp = a^{-1}b \qquad \text{ property of inverses}\\ x \amp = a^{-1}b \qquad \text{ 1 is the multiplicative identity} \end{align*} But this is not possible when \(a=0\text{,}\) so for example \(0x=3\) has no solution. However, there is another interesting case, when both \(a=b=0\text{,}\) then we have \(0x=0\text{,}\) and this equation is true for any value assigned to \(x\text{.}\) Thus there are exactly three possibilities:

1. a unique solution,

2. no solution,

3. infinitely many solutions.

What about a system of two linear equations in two unknowns? For example, \begin{equation*} \left\{ \begin{alignedat}{3} 2x \amp {}-{} \amp 2y \amp {}={} \amp {-14} \\ 3x \amp {}+{} \amp 4y \amp {}={} \amp 7 \end{alignedat} \right. \end{equation*} We wish to find the set of all pairs, \((x,y)\) that simultaneously satisfy both equations. To do so, we will use two facts:

• Multiplying both sides of an equation by a nonzero constant does not change the set of points which satisfy it.

• The sum of two equations is again a valid equation, that all solution points must satisfy.

Let's use these two facts to eliminate \(x\) from the second equation by multiplying the first equation by \(-3\text{,}\) the second by \(2\) and then summing them. \begin{equation*} \begin{alignedat}{3} -3[2x \amp {}-{} \amp 2y \amp {}={} \amp {-14}] \\ 2[3x \amp {}+{} \amp 4y \amp {}={} \amp 7] \\ \hline \amp {}~{} \amp 14y \amp {}={} \amp {56~} \\ \end{alignedat} \end{equation*}

The sum is equivalent to the equation \(14y = 56\) which of course has solution \(y=4\text{.}\) So we have determined that all solution points must have \(y\)-coordinate equal to 4. At this point we can plug this new knowledge into either of the two original equations to solve for \(x\text{.}\) Let's plug \(y=4\) into the first equation:

\begin{align*} 2x - 2(4) \amp= {-14}\\ 2x - 8 \amp= -14\\ 2x \amp= -6\\ x \amp= -3 \end{align*} And therefore the solution set is \(\left\{ (-3,4) \right\}\text{.}\)

The method of elimination for finding a solution to a linear system is a recursive algorithm. Suppose you start with a system of \(n\) equations in the \(n\) unknowns, \(x_1, x_2, \ldots, x_n\text{.}\) Pick one equation that contains \(x_1\) and use it to eliminate \(x_1\) from all the other equations. This step will yield a new system of \(n-1\) equations in the \(n-1\) unknowns, \(x_2, x_3, \ldots, x_n\text{.}\) Repeat this process of elimination until you have a single equation in the single unknown \(x_n\text{.}\) This single equation can be solved algebraically. Then we use the previous system of two equations in the two unknowns \(x_{n-1}, x_n\) to solve for \(x_{n-1}\text{.}\) Repeat this process of back substitution until you have solved for each unknown.

We will illustrate the method of elimination with an example system of three equations in three unknowns. \begin{equation*} \left\{ \begin{alignedat}{8} 3x \amp{}-{}\amp 2y \amp{}+{}\amp z \amp{}={}\amp -1 \amp\quad (1) \\ x \amp{}+{}\amp 5y \amp{}+{}\amp 2z \amp{}={}\amp 11 \amp\quad (2) \\ -x \amp{}+{}\amp 2y \amp{}-{}\amp z \amp{}={}\amp 3 \amp\quad (3) \end{alignedat} \right. \end{equation*} The first step is to use equation (1) to eliminate \(x\) from equations (2) and (3). \begin{equation*} \begin{alignedat}{9} 3x \amp{}-{}\amp 2y \amp{}+{}\amp z \amp{}={}\amp -1 \amp\quad\amp (1)\\ -3x \amp{}-{}\amp 15y \amp{}+{}\amp 6z \amp{}={}\amp 33 \amp\quad\amp {-3}(2)\\ \hline \amp{}-{}\amp 17y \amp{}+{}\amp 7z \amp{}={}\amp 32 \amp\quad\amp (1') \end{alignedat} \end{equation*} Next we use equation (1) to eliminate \(x\) from equation (3): \begin{equation*} \begin{alignedat}{9} 3x \amp{}-{}\amp 2y \amp{}+{}\amp z \amp{}={}\amp -1 \amp\quad\amp (1)\\ -3x \amp{}+{}\amp 6y \amp{}-{}\amp 3z \amp{}={}\amp 9 \amp\quad\amp 3(3)\\ \hline \amp{}~{}\amp 4y \amp{}-{}\amp 2z \amp{}={}\amp 8 \amp\quad\amp (2') \end{alignedat} \end{equation*}

We have reduced the \(3\times 3\) system to the \(2 \times 2\) system: \begin{equation*} \left\{ \begin{alignedat}{7} -17y \amp{}+{}\amp 7z \amp{}={}\amp 32 \amp\quad (1') \\ 4y \amp{}-{}\amp 2z \amp{}={}\amp 8 \amp\quad (2') \\ \end{alignedat} \right. \end{equation*} Next we must choose whether to eliminate \(y\) or \(z\text{.}\) Which variable we choose is irrelevant. If we choose to eliminate \(y\text{,}\) then since 17 and 4 are relatively prime, we would have to multiply equation (1') by 4 and equation (2') by 17, but there is a better choice. We will have smaller numbers to deal with if instead we eliminate \(z\text{:}\) \begin{equation*} \begin{alignedat}{6} -34y \amp{}+{}\amp 14z \amp{}={}\amp 64 \amp\quad 2(1') \\ 28y \amp{}-{}\amp 14z \amp{}={}\amp 56 \amp\quad 7(2') \\ \hline -6y \amp{}~{}\amp \amp{}={}\amp 120 \amp\quad (1'') \end{alignedat} \end{equation*} We now have a single equation in the one unknown \(y\) which of course can be solved with basic algebra to yield \(y = -20\text{.}\)

At this point we may substitute \(y=-20\) into either equation 1' or 2' to solve for \(z\text{.}\) Since equation 2' is a little simpler we choose it: \begin{equation*} 4(-20) - 2z = 8 \quad \Rightarrow \quad -2z = 88 \end{equation*} This equation has solution \(z=-44\text{,}\) which we can finally plug into any one of the original three equations to solve for \(x\text{.}\) Here we choose to use equation (2): \begin{align*} x + 5(20) + 2(-44) &= 11\\ x + 12 &= 11\\ x &= -1 \end{align*} Thus the solution to the system is the unique point, \((-1,-20,-44)\text{.}\)

Gaussian elimination is an algorithmic version of elimination that is more amenable to implementation in computer code. A key part of the algorithm is the recognition that many of the symbols in a linear system, are unnecessary and can safely be ignored. This is quite similar to what is done in synthetic division if you are familiar with that way of doing polynomial long division. Before we demonstrate the algorithm we need to define three new terms.

Let's solve the \(3 \times 3\) system from the algorithm statement above using Gaussian elimination.

Recall that an augmented matrix directly corresponds to a linear system. The solution set of a linear system has a geometric interpretation. In the tables which follow, the following symbols will have a special meaning.

\begin{align*} * & = \text{ any nonzero number}\\ \square & = \text{ any number} \end{align*}

The number of \(n\times n\text{,}\) row echelon matrix forms follows a pattern. It follows the pattern found in Pascal's triangle. Once you pick \(n\text{,}\) then just read the \(n^{\text{th}}\) row of Pascal's triangle from left to right to determine the number of \(n\times n\) matrices that have \(0, 1, 2, \ldots, n\) rows of all zeroes. So for example, when \(n=4\) there are exactly 4 unique row echelon forms that correspond to one row of all zeros, 6 unique row echelon forms that correspond to two rows of all zeros, and 4 unique row echelon forms that correspond to three rows of all zeros.

In Sage, a matrix is a special type of data structure. It is more than just a list of rows of numbers. It is an object with many built-in functions. In Python/Sage we construct objects with a special constructor function which in this case is called simply enough, matrix. The first parameter to the matrix constructor is the set of numbers from which the elements of the matrix come from. Here we will assume that all entries will be rationals, which in Sage is denoted QQ, and in math is often written \(\Q\text{.}\) Recall \(\Q\) is defined in terms of the integers, \(\Z\text{:}\) \begin{equation*} \Q = \set{\frac{a}{b}}{a,b \in \Z \text{ and } b \ne 0} . \end{equation*} The following Sage code creates a matrix and then a column vector corresponding to the column of constants from the right side of the system of linear equations.

\begin{equation*} \left\{ \begin{alignedat}{4} 2x \amp {}~{} \amp \amp {}+{} \amp z \amp {}={} \amp 8 \\ 3x \amp {}-{} \amp y \amp {}+{} \amp z \amp {}={} \amp 4 \\ x \amp {}+{} \amp y \amp {}-{} \amp 4z \amp {}={} \amp {-18} \end{alignedat} \right. \end{equation*}

Once we have an augmented matrix, our goal is to apply elementary row operations until we have transformed it into a matrix in row-echelon form often denoted REF. The three elementary row operations are given in the table.

An important thing to know about Python and Sage is that all counting begins with zero, so the first row in a matrix is the zeroth row instead of row 1. Due to this, I find it easiest to use the above commands if you use the terms: zeroth row, oneth row, twoth row, etc. when referring to rows rather than the usual, first row, second row, third row and so forth.

In Python, an object is something that has both state and behavior. State corresponds with the nouns associated with the object and behaviour corresponds with the verbs that are associated with the object. The behaviors of an object often change the object's state.

In the case of a matrix object, its state is the grid of numbers contained within the matrix. Each elementary row operation is a behavior that changes the state of its internal data. For example, each matrix has a swap_rows behavior which changes the order of two rows in the matrix. An object's behaviors/verbs are often called methods. A method can only chage the state of the object to which it belongs..

The cells below, when evaluated in order, correspond with the steps of Gaussian elimination and illustrate all three elmentary row operations. If you evaluate each cell below more than once you will change the matrix \(B\) in ways not intended by the Gaussian elimination algorithm.