A Brief Introduction to Sage for Math 300

By Brian Gill

First, some credit. I borrowed some of the ideas included in this worksheet from two sources:

One is a Sage Tutorial by Jason Grout and Leslie Hogben titled "AIM Basic Sage Linear Algebra" (which is a published notebook freely available on this site). The second source is Robert Beezer's open source linear algebra textbook...

Sage notebooks are based on cells, which are boxes one the page like the ones below which contain code. You will rarely see or need text like this in a notebook - I've just put it here to make this notebook and the instructions in it a little easier to read. To execute or evaluate a cell, you can either click in the cell and press the SHIFT and ENTER keys at the same time, or you can just click on the "execute" link that appears just below the cell once you've clicked in the cell. For example, click in the cell below and execute it to perform a basic calculation.

7*5-3

Other basic mathematical operations are performed in the standard ways as well. For example:

5^3

Now let's define our first matrix. We provide Sage with a list of the rows of the matrix, with the lists enclosed in square brackets, as follows:

A = matrix([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]])

Note that Sage creates the matrix and stores it as the object "A", but does not display any output. If you want to display the matrix, simply enter A alone and evaluate the expression:

Sage only cares about the order in which you actually execute things, NOT the order in which they appear on the screen. You can scroll up and down through a worksheet changing things, and Sage does not care about the order that they appear on the screen. For example, later in this worksheet, we'll redefine the matrix A to be something else. But if we scroll back up to the lines above and re-evaluate them, the matrix A will revert to the one defined above rather than one of the ones defined later in the worksheet. Sage does not care about the order that things appear in the worksheet--it only cares about the chronological order in which they are actually evaluated.

We've defined a matrix A above, but unfortunately there is a problem with that matrix for many purposes. To illustrate the problem, let's try to row reduce the matrix into reduced echelon form. To do this, we enter the following command:

A.echelon_form()

Something is wrong! This matrix is clearly not in echelon form at all. The leading entry in the second row is not a one. What happened?

Whenever we enter a matrix, Sage attaches a number system (a "ring") to that matrix. Since our matrix had only integers as entries, Sage guessed that we were working in a universe that consisted ONLY of integers, and it attached the "integer ring" to the matrix. When it tried to perform row operations on the matrix, the next step would be to multiply the second row by 1/4. But 1/4 does not exist in a universe of only integers, so it was impossible to for Sage to do that. So Sage quit and could not do anything more to reduce the matrix since it was working only with integers.

To fix the problem, we need to explicitly tell Sage what number system to work with. We have several choices:

ZZ: the set of integers

QQ: the rational numbers

RR: the real numbers

CC: the complex numbers

For purposes of this course, it will usually be best to use the rational numbers QQ. This will allow Sage to use fractions, but still force Sage to do exact integer arithmetic and report answers as fractions (like 137/100 rather than 1.37). Allowing Sage to work with real or complex numbers lets Sage use decimal approximations and potentially introduces problems with round-off errors. Sometimes an entry in a matrix which is supposed to be exactly 0 ends up becoming something like 1.392x10^(-19). This is then followed by a series of other erroneous row operations that try to eliminate that nonzero entry in the matrix. Whenever possible, we'll avoid this problem by using rationals and exact arithmetic.

Let's re-enter our matrix, this time telling it that our universe is the rational numbers QQ. Note that we can enter more than one command in a single cell in Sage.

B = matrix(QQ,[[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]])
B

Let's try again to get the matrix into reduced echelon form:

B.echelon_form()

That's much better. Working in the rationals, Sage correctly does the computations as expected.

Let's create a larger matrix with random entries and then reduce it to reduced echelon form:

D=random_matrix(QQ,20,21,2^5)
D

D.echelon_form()

Note that Sage has performed exact arithmetic and given the solution to the system as very complicated exact rational expressions. In this case, it might be nicer to have decimal approximations of these solutions (because they would be much easier to read and interpret). However, it is much better to first SOLVE the system using exact arithmetic (as we have done), and then get a decimal approximation of the solution rather than using real numbers throughout and allowing rounding errors at every stage of the computations. We can get a decimal approximation of the solution as follows. The command "D.echelon_form()" asks Sage to reduce D to echelon form in its universe (QQ), then the matrix command converts that final row reduced matrix to a matrix in the new universe RR. This gives us the solution in decimal form.

E=matrix(RR,D.echelon_form())
E

Matrix Operations

Basic matrix operations in Sage are quite simple. Let's start by defining a few matrices:

A=matrix(QQ,[[1,2,0],[3,0,1]])
B=matrix(QQ,[[1,2],[0,1],[-2,4]])
C=matrix(QQ,[[3,-2,5],[0,-1,2]])
D=matrix(QQ,[[1,0],[3,2]])
A, B, C, D

Matrix addition is quite simple:

A+C

Scalar multiplication:

5*A

Matrix multiplication:

A*B

We can also find the transpose of a matrix:

transpose(A)

Note that we attempt an operation that is undefined, such as multiplying matrices when the number of rows of the first matrix does not match the number of columns of the second, Sage will return an error:

A*C

We get a similar error if we try to add matrices of different sizes:

A+B

We can also multiply a matrix by itself several times:

D*D*D*D*D

Or we can do the same thing by using an exponent for shorthand notation:

D^5

Vectors can be entered in the same way as matrices:

u = matrix(QQ,[[1],[2],[3]])
u

A*u