Computation with Sage

Section 1.3 Computation with Sage

Linear algebra owes its prominence as a powerful scientific tool to the ever-growing power of computers. Carl Cowen, a former president of the Mathematical Association of America, has said, “No serious application of linear algebra happens without a computer.” Indeed, Cowen notes that, in the 1950s, working with a system of 100 equations in 100 variables was difficult. Today, scientists and mathematicians routinely work on problems that are vastly larger. This is only possible because of today’s computing power.

It is therefore important for any student of linear algebra to become comfortable solving linear algebraic problems on a computer. This section will introduce you to a program called Sage that can help. While you may be able to do much of this work on a graphing calculator, you are encouraged to become comfortable with Sage as we will use increasingly powerful features as we encounter their need.

Subsection 1.3.1 Introduction to Sage

There are several ways to access Sage.

If you are reading this book online, there will be embedded “Sage cells” at appropriate places in the text. You have the opportunity to type Sage commands into these cells and execute them, provided you are connected to the Internet. Please be aware that your work will be lost if you reload the page.
Here is a Sage cell containing a command that asks Sage to multiply 5 and 3. You may execute the command by pressing the Evaluate button.
You may also go to cocalc.com, sign up for an account, open a new project, and create a “Sage worksheet.” Once inside the worksheet, you may enter commands as shown here, and evaluate them by pressing Enter on your keyboard while holding down the Shift key.
There is a page of Sage cells at gsvu.edu/s/0Ng. Any results obtained by evaluating one cell are available in other cells. However, your work will be lost when the page is reloaded.

Throughout the text, we will introduce new Sage commands that allow us to explore linear algebraic concepts. These commands are collected and summarized in the reference found in Appendix A.

Activity 1.3.1. Basic Sage commands.

Sage uses the standard operators +, -, *, /, and ^ for the usual arithmetic operations. By entering text in the cell below, ask Sage to evaluate

\begin{equation*} 3 + 4(2^4 - 1) \end{equation*}
Notice that we can create new lines by pressing Enter and entering additional commands on them. What happens when you evaluate this Sage cell?
Notice that we only see the result from the last command. With the print command, we may see earlier results, if we wish.
We may give a name to the result of one command and refer to it in a later command.
Suppose you have three tests in your linear algebra class and your scores are 90, 100, and 98. In the Sage cell below, add your scores together and call the result total. On the next line, find the average of your test scores and print it.
If you are not a programmer, you may ignore this part. If you are an experienced programmer, however, you should know that Sage is written in the Python programming language and that you may enter Python code into a Sage cell.

Subsection 1.3.2 Sage and matrices

When we encounter a matrix, Theorem 1.2.6 tells us that there is exactly one reduced row echelon matrix that is row equivalent to it.

In fact, the uniqueness of this reduced row echelon matrix is what motivates us to define this particular form. When solving a system of linear equations using Gaussian elimination, there are other row equivalent matrices that reveal the structure of the solution space. The reduced row echelon matrix is simply a convenience as it is an agreement we make with one another to seek the same matrix.

An added benefit is that we can ask a computer program, like Sage, to find reduced row echelon matrices for us. We will learn how to do this now that we have a little familiarity with Sage.

First, notice that a matrix has a certain number of rows and columns. For instance, the matrix

\begin{equation*} \left[ \begin{array}{rrrrr} * \amp * \amp * \amp * \amp * \\ * \amp * \amp * \amp * \amp * \\ * \amp * \amp * \amp * \amp * \\ \end{array} \right] \end{equation*}

has three rows and five columns. We consequently refer to this as a \(3\times 5\) matrix.

We may ask Sage to create the \(2\times4\) matrix

\begin{equation*} \left[ \begin{array}{rrrr} -1 \amp 0 \amp 2 \amp 7 \\ 2 \amp 1 \amp -3 \amp -1 \\ \end{array} \right] \end{equation*}

by entering

When evaluated, Sage will confirm the matrix by writing out the rows of the matrix, each inside square brackets.

Notice that there are three separate things (we call them arguments) inside the parentheses: the number of rows, the number of columns, and the entries of the matrix listed by row inside square brackets. These three arguments are separated by commas. Notice that there is no way of specifying whether this is an augmented or coefficient matrix so it will be up to us to interpret our results appropriately.

Sage syntax.

Some common mistakes are

to forget the square brackets around the list of entries,
to omit an entry from the list or to add an extra one,
to forget to separate the rows, columns, and entries by commas, and
to omit the parentheses around the arguments after matrix.

If you see an error message, carefully proofread your input and try again.

Alternatively, you can create a matrix by simply listing its rows, like this

matrix([ [-1, 0, 2, 7],	
	 [ 2, 1,-3,-1] ])

Activity 1.3.2. Using Sage to find row reduced echelon matrices.

Enter the following matrix into Sage.

\begin{equation*} \left[ \begin{array}{rrrr} -1 \amp -2 \amp 2 \amp -1 \\ 2 \amp 4 \amp -1 \amp 5 \\ 1 \amp 2 \amp 0 \amp 3 \end{array} \right] \end{equation*}
Give the matrix the name \(A\) by entering
```
A = matrix( ..., ..., [ ... ])
	      
```
We may then find its reduced row echelon form by entering
```
A = matrix( ..., ..., [ ... ])
A.rref()
	      
```
A common mistake is to forget the parentheses after rref.
Use Sage to find the reduced row echelon form of the matrix from Item a of this activity.
Use Sage to describe the solution space of the system of linear equations

\begin{equation*} \begin{alignedat}{5} -x_1 \amp \amp \amp \amp \amp {}+{} \amp 2x_4 \amp {}={} \amp 4 \\ \amp \amp 3x_2 \amp {}+{} \amp x_3 \amp {}+{} \amp 2x_4 \amp {}={} \amp 3 \\ 4x_1 \amp {}-{} \amp 3x_2 \amp \amp \amp {}+{} \amp x_4 \amp {}={} \amp 14 \\ \amp \amp 2x_2 \amp {}+{} \amp 2x_3 \amp {}+{} \amp x_4 \amp {}={} \amp 1 \\ \end{alignedat} \end{equation*}
Consider the two matrices:

\begin{equation*} \begin{array}{rcl} A \amp = \amp \left[ \begin{array}{rrrr} 1 \amp -2 \amp 1 \amp -3 \\ -2 \amp 4 \amp 1 \amp 1 \\ -4 \amp 8 \amp -1 \amp 7 \\ \end{array}\right] \\ B \amp = \amp \left[ \begin{array}{rrrrrr} 1 \amp -2 \amp 1 \amp -3 \amp 0 \amp 3 \\ -2 \amp 4 \amp 1 \amp 1 \amp 1 \amp -1 \\ -4 \amp 8 \amp -1 \amp 7 \amp 3 \amp 2 \\ \end{array}\right] \\ \end{array} \end{equation*}

We say that \(B\) is an augmentation of \(A\) because it is obtained from \(A\) by adding some more columns.

Using Sage, define the matrices and compare their reduced row echelon forms. What do you notice about the relationship between the two reduced row echelon forms?
Using the system of equations in Item c, write the augmented matrix corresponding to the system of equations. What did you find for the reduced row echelon form of the augmented matrix?

Now write the coefficient matrix of this system of equations. What does Item d of this activity tell you about its reduced row echelon form?

Sage practices.

Here are some practices that you may find helpful when working with matrices in Sage.

Break the matrix entries across lines, one for each row, for better readability by pressing Enter between rows.
```
A = matrix(2, 4, [ 1, 2, -1, 0,
                  -3, 0,  4, 3 ])
	
```
Print your original matrix to check that you have entered it correctly. You may want to also print a dividing line to separate matrices.
```
A = matrix(2, 2, [ 1, 2,
                   2, 2])
print (A)
print ("---------")
A.rref()
	
```

The last part of the previous activity, Item d, demonstrates something that will be helpful for us in the future. In that activity, we started with a matrix \(A\text{,}\) which we augmented by adding some columns to obtain a matrix \(B\text{.}\) We then noticed that the reduced row echelon form of \(B\) is itself an augmentation of the reduced row echelon form of \(A\text{.}\)

To illustrate, we can consider the reduced row echelon form of the augmented matrix:

\begin{equation*} \left[ \begin{array}{ccc|c} -2 \amp 3 \amp 0 \amp 2 \\ -1 \amp 4 \amp 1 \amp 3 \\ 3 \amp 0 \amp 2 \amp 2 \\ 1 \amp 5 \amp 3 \amp 7 \\ \end{array} \right] \sim \left[ \begin{array}{ccc|c} 1 \amp 0 \amp 0 \amp -4 \\ 0 \amp 1 \amp 0 \amp -2 \\ 0 \amp 0 \amp 1 \amp 7 \\ 0 \amp 0 \amp 0 \amp 0 \\ \end{array} \right] \end{equation*}

We can then determine the reduced row echelon form of the coefficient matrix by looking inside the augmented matrix.

\begin{equation*} \left[ \begin{array}{ccc} -2 \amp 3 \amp 0 \\ -1 \amp 4 \amp 1 \\ 3 \amp 0 \amp 2 \\ 1 \amp 5 \amp 3 \\ \end{array} \right] \sim \left[ \begin{array}{ccc} 1 \amp 0 \amp 0 \\ 0 \amp 1 \amp 0 \\ 0 \amp 0 \amp 1 \\ 0 \amp 0 \amp 0 \\ \end{array} \right] \end{equation*}

If we trace through the steps in the Gaussian elimination algorithm carefully, we see that this is a general principle, which we now state.

Proposition 1.3.1. Augmentation Principle.

If matrix \(B\) is an augmentation of matrix \(A\text{,}\) then the reduced row echelon form of \(B\) is an augmentation of the reduced row echelon form of \(A\text{.}\)

Subsection 1.3.3 Computational effort

At the beginning of this section, we indicated that linear algebra has become more prominent as computers have grown more powerful. Computers, however, still have limits. Let’s consider how much effort is expended when we ask to find the reduced row echelon form of a matrix. We will measure, very roughly, the effort by the number of times the algorithm requires us to multiply or add two numbers.

We will assume that our matrix has the same number of rows as columns, which we call \(n\text{.}\) We are mainly interested in the case when \(n\) is very large, which is when we need to worry about how much effort is required.

Let’s first consider the effort required for each of our row operations.

Scaling a row multiplies each of the \(n\) entries in a row by some number, which requires \(n\) operations.
Interchanging two rows requires no multiplications or additions so we won’t worry about the effort required by an interchange.
A replacement requires us to multiply each entry in a row by some number, which takes \(n\) operations, and then add the resulting entries to another row, which requires another \(n\) operations. The total number of operations is \(2n\text{.}\)

Our goal is to transform a matrix to its reduced row echelon form, which looks something like this:

\begin{equation*} \left[ \begin{array}{cccc} 1 \amp 0 \amp \ldots \amp 0 \\ 0 \amp 1 \amp \ldots \amp 0 \\ \vdots \amp \vdots \amp \ddots \amp 0 \\ 0 \amp 0 \amp \ldots \amp 1 \\ \end{array} \right]\text{.} \end{equation*}

We roughly perform one replacement operation for every 0 entry in the reduced row echelon matrix. When \(n\) is very large, most of the \(n^2\) entries in the reduced row echelon form are 0 so we need roughly \(n^2\) replacements. Since each replacement operation requires \(2n\) operations, the number of operations resulting from the needed replacements is roughly \(n^2(2n) = 2n^3\text{.}\)

Each row is scaled roughly one time so there are roughly \(n\) scaling operations, each of which requires \(n\) operations. The number of operations due to scaling is roughly \(n^2\text{.}\)

Therefore, the total number of operations is roughly

\begin{equation*} 2n^3 + n^2\text{.} \end{equation*}

When \(n\) is very large, the \(n^2\) term is much smaller than the \(n^3\) term. We therefore state that

Observation 1.3.2.

The number of operations required to find the reduced row echelon form of an \(n\times n\) matrix is roughly proportional to \(n^3\text{.}\)

This is a very rough measure of the effort required to find the reduced row echelon form; a more careful accounting shows that the number of arithmetic operations is roughly \(\frac23 n^3\text{.}\) As we have seen, some matrices require more effort than others, but the upshot of this observation is that the effort is proportional to \(n^3\text{.}\) We can think of this in the following way: If the size of the matrix grows by a factor of 10, then the effort required grows by a factor of \(10^3 = 1000\text{.}\)

While today’s computers are powerful, they cannot handle every problem we might ask of them. Eventually, we would like to be able to consider matrices that have \(n=10^{12}\) (a trillion) rows and columns. In very broad terms, the effort required to find the reduced row echelon matrix will require roughly \((10^{12})^3 = 10^{36}\) operations.

To put this into context, imagine we need to solve a linear system with a trillion equations and a trillion variables and that we have a computer that can perform a trillion, \(10^{12}\text{,}\) operations every second. Finding the reduced row echelon form would take about \(10^{16}\) years. At this time, the universe is estimated to be approximately \(10^{10}\) years old. If we started the calculation when the universe was born, we’d be about one-millionth of the way through.

This may seem like an absurd situation, but we’ll see in Subsection 4.5.3 how we use the results of such a computation every day. Clearly, we will need some better tools to deal with really big problems like this one.

Subsection 1.3.4 Summary

We learned some basic features of Sage with an emphasis on finding the reduced row echelon form of a matrix.

Sage can perform basic arithmetic using standard operators. Sage can also save results from one command to be reused in a later command.
We may define matrices in Sage and find the reduced row echelon form using the rref command.
We saw an example of the Augmentation Principle, which we then stated as a general principle.
We saw that the computational effort required to find the reduced row echelon form of an \(n\times n\) matrix is proportional to \(n^3\text{.}\)

Appendix A contains a reference outlining the Sage commands that we have encountered.

Exercises 1.3.5 Exercises

1.

Consider the linear system

\begin{equation*} \begin{alignedat}{4} x \amp {}+{} \amp 2y \amp {}-{} \amp z \amp {}={} \amp 1 \\ 3x \amp {}+{} \amp 2y \amp {}+{} \amp 2z \amp {}={} \amp 7 \\ -x \amp \amp \amp {}+{} \amp 4z \amp {}={} \amp -3 \\ \end{alignedat} \end{equation*}

Write this system as an augmented matrix and use Sage to find a description of the solution space.

2.

Shown below are some traffic patterns in the downtown area of a large city. The figures give the number of cars per hour traveling along each road. Any car that drives into an intersection must also leave the intersection. This means that the number of cars entering an intersection in an hour is equal to the number of cars leaving the intersection.

Let’s begin with the following traffic pattern.
1. How many cars per hour enter the upper left intersection? How many cars per hour leave this intersection? Use this to form a linear equation in the variables \(x\text{,}\) \(y\text{,}\) \(z\text{,}\) and \(w\text{.}\)
2. Form three more linear equations from the other three intersections to form a linear system having four equations in four variables. Then use Sage to find the solution space to this system.
3. Is there exactly one solution or infinitely many solutions? Explain why you would expect this given the information provided.
Another traffic pattern is shown below.
1. Once again, write a linear system for the quantities \(x\text{,}\) \(y\text{,}\) \(z\text{,}\) and \(w\) and solve the system using the Sage cell below.
2. What can you say about the solution of this linear system? Is there exactly one solution or infinitely many solutions? Explain why you would expect this given the information provided.
3. What is the smallest possible amount of traffic flowing through \(x\text{?}\)

3.

A typical problem in thermodynamics is to find the steady-state temperature distribution inside a thin plate if we know the temperature around the boundary. Let \(T_1, T_2, \ldots, T_6\) be the temperatures at the six nodes inside the plate as shown below.

The temperature at a node is approximately the average of the four nearest nodes: for instance,

\begin{equation*} T_1 = (10 + 15 + T_2 + T_4)/4\text{,} \end{equation*}

which we may rewrite as

\begin{equation*} 4T_1 - T_2 - T_4 = 25\text{.} \end{equation*}

Set up a linear system to find the temperature at these six points inside the plate. Then use Sage to solve the linear system. If all the entries of the matrix are integers, Sage will compute the reduced row echelon form using rational numbers. To view a decimal approximation of the results, you may use

A.rref().numerical_approx(digits=4)

In the real world, the approximation becomes better as we add more and more points into the grid. This is a situation where we may want to solve a linear system having millions of equations and millions of variables.

4.

The fuel inside model rocket motors is a black powder mixture that ideally consists of 60% charcoal, 30% potassium nitrate, and 10% sulfur by weight.

Suppose you work at a company that makes model rocket motors. When you come into work one morning, you learn that yesterday’s first shift made a perfect batch of fuel. The second shift, however, misread the recipe and used 50% charcoal, 20% potassium nitrate and 30% sulfur. Then the two batches were mixed together. A chemical analysis shows that there are 100.3 pounds of charcoal in the mixture and 46.9 pounds of potassium nitrate.

Assuming the first shift produced \(x\) pounds of fuel and the second \(y\) pounds, set up a linear system in terms of \(x\) and \(y\text{.}\) How many pounds of fuel did the first shift produce and how many did the second shift produce?
How much sulfur would you expect to find in the mixture?

5.

This exercise is about balancing chemical reactions.

Chemists denote a molecule of water as \(\text{H}_2\text{O}\text{,}\) which means it is composed of two atoms of hydrogen (H) and one atom of oxygen (O). The process by which hydrogen burns is described by the chemical reaction

\begin{equation*} x\thinspace \text{H}_2 + y\thinspace\text{O}_2 \to z\thinspace \text{H}_2\text{O} \end{equation*}

This means that \(x\) molecules of hydrogen \(\text{H}_2\) combine with \(y\) molecules of oxygen \(\text{O}_2\) to produce \(z\) water molecules. The number of hydrogen atoms is the same before and after the reaction; the same is true of the oxygen atoms.
1. In terms of \(x\text{,}\) \(y\text{,}\) and \(z\text{,}\) how many hydrogen atoms are there before the reaction? How many hydrogen atoms are there after the reaction? Find a linear equation in \(x\text{,}\) \(y\text{,}\) and \(z\) by equating these quantities.
2. Find a second linear equation in \(x\text{,}\) \(y\text{,}\) and \(z\) by equating the number of oxygen atoms before and after the reaction.
3. Find the solutions of this linear system. Why are there infinitely many solutions?
4. In this chemical setting, \(x\text{,}\) \(y\text{,}\) and \(z\) should be positive integers. Find the solution where \(x\text{,}\) \(y\text{,}\) and \(z\) are the smallest possible positive integers.
Now consider the reaction where potassium permanganate and manganese sulfate combine with water to produce manganese dioxide, potassium sulfate, and sulfuric acid:

\begin{equation*} x_1\thinspace \text{K}\text{Mn}\text{O}_4 + x_2\thinspace \text{Mn}\text{S}\text{O}_4 + x_3\thinspace \text{H}_2\text{O} \to x_4\thinspace \text{Mn}\text{O}_2 + x_5\thinspace \text{K}_2\text{S}\text{O}_4 + x_6\thinspace \text{H}_2\text{S}\text{O}_4. \end{equation*}
As in the previous exercise, find the appropriate values for \(x_1, x_2, \ldots, x_6\) to balance the chemical reaction.

6.

We began this section by stating that increasing computational power has helped linear algebra assume a prominent role as a scientific tool. Later, we looked at one computational limitation: once a matrix gets to be too big, it is not reasonable to apply Gaussian elimination to find its reduced row echelon form.

In this exercise, we will see another limitation: computer arithmetic with real numbers is only an approximation because computers represent real numbers with only a finite number of bits. For instance, the number pi

\begin{equation*} \pi = 3.141592653589793238462643383279502884197169399\ldots \end{equation*}

would be approximated inside a computer by, say,

\begin{equation*} \pi\approx 3.141592653589793 \end{equation*}

Most of the time, this is not a problem. However, when we perform millions or even billions of arithmetic operations, the error in these approximations starts to accumulate and can lead to results that are wildly inaccurate. Here are two examples demonstrating this.

Let’s first see an example showing that computer arithmetic really is an approximation. First, consider the linear system

\begin{equation*} \begin{alignedat}{4} x \amp {}+{} \amp \frac12y \amp {}+{} \amp \frac13z \amp {}={} \amp 1 \\ \frac12x \amp {}+{} \amp \frac13y \amp {}+{} \amp \frac14z \amp {}={} \amp 0 \\ \frac13x \amp {}+{} \amp \frac14y \amp {}+{} \amp \frac15z \amp {}={} \amp 0 \\ \end{alignedat} \end{equation*}
If the coefficients are entered into Sage as fractions, Sage will find the exact reduced row echelon form. Find the exact solution to this linear system.
Now let’s ask Sage to compute with real numbers. We can do this by representing one of the coefficients as a decimal. For instance, the same linear system can be represented as

\begin{equation*} \begin{alignedat}{4} x \amp {}+{} \amp 0.5y \amp {}+{} \amp \frac13z \amp {}={} \amp 1 \\ \frac12x \amp {}+{} \amp \frac13y \amp {}+{} \amp \frac14z \amp {}={} \amp 0 \\ \frac13x \amp {}+{} \amp \frac14y \amp {}+{} \amp \frac15z \amp {}={} \amp 0 \\ \end{alignedat} \end{equation*}
Most computers do arithmetic using either 32 or 64 bits. To magnify the problem so that we can see it better, we will ask Sage to do arithmetic using only 10 bits as follows.
What does Sage give for the solution now? Compare this to the exact solution that you found previously.
Some types of linear systems are particularly sensitive to errors resulting from computers’ approximate arithmetic. For instance, suppose we are interested in the linear system

\begin{equation*} \begin{alignedat}{3} x \amp {}+{} \amp y \amp {}={} \amp 2 \\ x \amp {}+{} 1.001\amp y \amp {}={} \amp 2 \\ \end{alignedat} \end{equation*}
Find the solution to this linear system.
Suppose now that the computer has accumulated some error in one of the entries of this system so that it incorrectly stores the system as

\begin{equation*} \begin{alignedat}{3} x \amp {}+{} \amp y \amp {}={} \amp 2 \\ x \amp {}+{} 1.001\amp y \amp {}={} \amp 2.001 \\ \end{alignedat} \end{equation*}
Find the solution to this linear system.
Notice how a small error in one of the entries in the linear system leads to a solution that has a dramatically large error. Fortunately, this is an issue that has been well studied, and there are techniques that mitigate this type of behavior.

Understanding Linear Algebra

Search Results: