1439-Section 1.3 Homework Examples.sagews

Section 1.3 Homework Examples

In Section 1.3, we begin to solve systems of equations.  Here are some examples of how to use SAGE to carry out these computations!  Exercise numbers refer to exercises in our textbook.

Let's start with Exercise 3 (a).  The matrix associated with the system in part (a) is $\left[\begin{array}{ccc}2&3&7\\ 1&2&-2\end{array}\right]$To enter this in SAGE, we use the Matrix command.

(If there is a text mess above the matrices: sorry about that.  We are trying to figure out a typesetting problem in the newer version of SAGE.)

A=Matrix(QQ,[[2,3,7],[1,2,-2]])

The "[ ... ]" braces in SAGE are lists.  So, you can think of a row in a matrix as a list of numbers, and you can think of a matrix as a list of rows.  As a side note: we could define a table of numbers without the Matrix command, e.g. just type in

A=[[2,3,7],[1,2,-2]]

But, this just creates a table; matrices are tables for which operations have been defined, such as row reduction and arithmetic.  So, the Matrix command really associates a bunch of useful arithmetic operations with a table of data.  More on QQ later.

Typing A and evaluating now displays the matrix.

A

[ 2 3 7] [ 1 2 -2]

As we did in class, we can put the matrix A in echelon form (using Gauss-Jordan elimination) and read off the solution.  The syntax is as follows:

A.echelon_form()

[ 1 0 20] [ 0 1 -11]

(If you are familiar with Java, this probably does not look too weird.  The matrix A is an object, and the function echelon_form() is one of the object's methods.  The use of Matrix when defining A associated a bunch of operations with A; in particular, one of them is the very handy row reduction algorithm.  Some of the syntax in SAGE follows the object-oriented syntax.)

So we can conclude that our system has a single solution $x=20$, $y=-11$.

Next, we turn to Exercise 5 (a).  The matrix for the system is:

B=Matrix(QQ,[[1,1,0,1],[2,2,1,1],[2,2,0,2]])

B.echelon_form()

[ 1 1 0 1] [ 0 0 1 -1] [ 0 0 0 0]

Here, we have a system with one free variable, $x_2$.  The bound variables are$x_1=1-x_2\text{  and  }x_3=-1$

Last, I want to come back to our systems and talk about the QQ in the notation.  Whenever we work with a linear system of equations, we need to specify the set of numbers over which we are allowed to work.  For us, we often specify a "common" set of coefficients, e.g. $\mathbb{Z}$, $\mathbb{Q}$, $\mathbb{R}$, or $\mathbb{C}$.  In row operations, the specified set of numbers is the set from which we draw our scaling constants.  So for example, we can not perform the row operation "multiply row 1 by 1/2 and add to row 2" when the set specified is $\mathbb{Z}$, because 1/2 is not in $\mathbb{Z}$!

So, let's look at an example.

C=Matrix(ZZ,[[2,0,2],[0,1,-1]])

C.echelon_form()

[ 2 0 2] [ 0 1 -1]

So here, we cannot legally use row operations involving only multiplication by constants from $\mathbb{Z}$ to change the (1,1) entry to a 1.  But, if we change the specified set of numbers to $\mathbb{Q}$, we have:

D=Matrix(QQ,[[2,0,2],[0,1,-1]])

D.echelon_form()

[ 1 0 1] [ 0 1 -1]

For almost all of the computations we perform in this class, we will want to use the reals $\mathbb{R}$ (denoted by RR in SAGE) or the complexes $\mathbb{C}$ (denoted by CC in SAGE).  For row reduction, working with $\mathbb{Q}$ will work whenever only rational numbers appear in the matrix, and makes the results appear much cleaner than they do when we specify $\mathbb{R}$ as the set of constants; compare the following two calculations.

G=Matrix(QQ,[[2,0,1],[0,1,-1]])

G.echelon_form()

[ 1 0 1/2] [ 0 1 -1]

H=Matrix(RR,[[2,0,1],[0,1,-1]])

H.echelon_form()

[ 1.00000000000000 0.000000000000000 0.500000000000000] [0.000000000000000 1.00000000000000 -1.00000000000000]

Chapter 2: Computations and Examples

In chapter 2, we look at matrix arithmetic.  A computer algebra system can be extremely helpful in dealing with problems requiring a lot of matrix arithmetic.  Just as with the integer/rational number arithmetic you learn in elementary school, we can learn some very important ideas from performing basic computations by hand.  And, just as in elementary school, there are some application that require hundreds, thousands, (or even trillions) of basic calculations; it is not just tedious to perform these operations by hand, it is silly.

Let's see how SAGE handles the basic arithmetic operations for $2\times 2$ matrices.

A=Matrix(QQ,[[2,1],[0,1]])
B=Matrix(QQ,[[-1,3],[2,-2]])
print A
print B

[2 1] [0 1] [-1 3] [ 2 -2]

A+B

[ 1 4] [ 2 -1]

3*A

[6 3] [0 3]

3*A-4*B

[-10 -5] [ -8 -5]

A*B

[ 0 4] [ 2 -2]

B*A

[-2 2] [ 4 0]

Notice that A and B are not commutative!

What about matrix powers?  We saw this come up when we investigated Markov chains.

A^2

[4 3] [0 1]

(Side note: This may sound a little weird, but it is important to check, in any CAS, what the notation "A^2" is actually doing.  There are at least two other matrix products that I can think of that are heavily used, and so the CAS may be using a different product when it computes "A^2".  An easy check is to just run a few simple examples.)

A*A

[4 3] [0 1]

A^2==A*A

True

One of the things that CAS's are good at is keeping track of lots of significant digits, as in the following example.

C=A^200

Let's find out what the first row, first column entry is.  SAGE does C-style indexing of arrays; in practical terms, this means that the first item in the array A is accessed via A[0] (not A[1]).  The second entry is A[1], the third is A[2], and so forth.

The first entry in the first array in an array of arrays is A[0,0], the second entry in the first array is A[0,1], and so forth.

How do we convert back and forth from linear algebra numbering to CAS numbering?  Well, $A_{i,j}$=A[i-1,j-1].  For my matrix C=A^200, if I wish to access the entry $C_{11}$, I type C[0,0] and hit .

C[0,0]

1606938044258990275541962092341162602522202993782792835301376

That is a long integer.  (Note for the programming geeks out there: SAGE and Python perform integer arithmetic as follows.  A certain amount of memory is devoted to your running program, and as soon as it starts running, it starts eating up that memory.  Storing a long integer like this takes up more memory than rounding it and storing it with only a few significant figures and an exponent for scientific notation; but, it is more precise.  In any arithmetic operation involving integers, Python will just keep increasing the length of any integer when necessary, until it runs out of memory.  On the plus side, this allows for great precision in computations, which can be very important if we have to work with large numbers of significant figures.  On the negative side, this practice requires dynamic memory allocation for the integer type (so computations might run really slowly, compared to C), and some care is required to make sure that your program does not run out of memory.)

Let's take a look at a few other entries of the matrix C:

C[0,1]

1606938044258990275541962092341162602522202993782792835301375

C[0,0]-C[0,1]

1

The second row of the matrix C is

C[1]

(0, 1)

The transpose of a matrix A can be obtained in either of two ways.

Recalling that A is

A

[2 1] [0 1]

we compute the transpose as

A.transpose()

[2 0] [1 1]

or as

transpose(A)

[2 0] [1 1]

B=transpose(A)*A
B

[4 2] [2 2]

Notice that B is symmetric!

transpose(B)==B

True

Let's suppose we wish to deal with vectors, and their outer and inner products.

We can try to define vectors as

u=Matrix(QQ,[2,-1,1])
v=Matrix(QQ,[3,4,1])
print u
print v

[ 2 -1 1] [3 4 1]

but this creates row vectors instead of column vectors.  This is fine if you want row vectors, but most of the time we deal with column vectors.  So, try defining u and v as

u=transpose(Matrix(QQ,[2,-1,1]))
v=transpose(Matrix(QQ,[3,4,1]))

u

[ 2] [-1] [ 1]

The inner product of $u$ and $v$ is $u^Tv$, which we compute as

transpose(u)*v

[3]

Similarly, the outer product of $u$ and $v$ is the matrix $uv^T$.

u*transpose(v)

[ 6 8 2] [-3 -4 -1] [ 3 4 1]

Examples from Section 2.3

Exercise 9 page 85:  We are given the transition matrix for a Markov chain,

A=Matrix(QQ,[[1/10,3/10,0],[0,4/10,1],[9/10,3/10,0]])
A

[1/10 3/10 0] [ 0 2/5 1] [9/10 3/10 0]

x=transpose(Matrix(QQ,[5/10,5/10,0]))

x1=A*x
x1

[1/5] [1/5] [3/5]

x2=A*x1
x2

[ 2/25] [17/25] [ 6/25]

We can show that the matrix is a Markov chain by summing the columns.  It is easier to sum the rows of the transpose:

[sum(r) for r in transpose(A)]

[1, 1, 1]

Since all of the columns sum to 1, we have a Markov chain.

Next, we can go a little beyond the exercise.  What is the steady state?

v=(A^1000)*x
v

[431034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379402818298732410068704274813173273449057439083568517230153482346337773599426099247535338599716659213757160557775238507997271486212181922774454850824583038213565828836692377534349581962958893256195868746109938147758135478192530120997502569834560290405674759141610661852121638545982972451167084943615804018596950641986860317863394594497297106351644490650937884704112876535920342593546914172850150194735065143095488303473391589999276316153345694708918780756326175892324285452805310227567546316847380018403590754825858664825528704048027590633850296944254279753575545467809675306327449841178411416189771774710614515140905833601277028159097669543400483614235991283946370355407152196034994684635617/2500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000] [646551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068844143430405446392562202275344430805381566614945606846880443346190342926314873948380569969291289580856396023694775984724727831980358444941264774475292558750643079843508148569163423951434361153881408229300968124766608608560213698400858323065152017709397384535985912755314216608317133377826991680593901296839497144148060505836470213882829173454542663327977370130432975333729005768037851120024346527258467225504313305822905942305755906609824853023577598118302529960141368263848500825643146288785835892063852478214584973341019232722310262500654162247128523268250057972539629470769790957690125101263628758451456196663842095192239534168711355486377945271460983280469362691436289420865025494720691/1250000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000] [775862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482758620689655172413793103448275862068965517241379310344827586206896551724137931034482908894840456697146171320636137864940179427686540269076085630961281540547944152855703521461700761624530047394835209522553272849827101187343015600224831844285148011476291325327323570134172384436041314795288125602708647304687042482200780784035135674175530471786417512637249928237382760793178931695196393387724055069717018670463664977737044546739270182693107375035021172796621645870377383587101156750748000405895885084880796525389211870627004599243926023007068764187392978019497688121146161105580948197468704288744971388492432830507351884364841378561488673709924338587111065752132968243441338381282970708386473091531409976014243903503479619483843625842842042155114904261720268962234954325923001/2500000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000]

Gak!  Those are some long numbers.  Can we simplify them?

[[float(r) for r in s] for s in v]

[[0.17241379310344826], [0.51724137931034475], [0.31034482758620685]]

Verifying, we have

A*v==v

False

Hmmm, what is going on?  Actually, we are looking at a "situation", numerically speaking.  The 1000th, or 10,000th, or the millionth  iterate need not be the fixed point.  But the size of the difference between iterates must be getting smaller (in analysis, we have a Cauchy sequence).  Let's see what the size of the iterates is doing as time goes on.

def dist(u,v):
    return sqrt(float((transpose(u-v)*(u-v))[0,0]))
data=[]
u=(A^1000)*x
for j in range(20):
    v=A*u
    data.append(dist(u,v))
    u=v

data

[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

Interesting!  When we truncate the result using float, the computer can not tell the difference between the 1000-1020th iterates.

We can find the fixed point in a different way as well.  We need to solve the equation $Ax=x$ for the vector $x$.  Equivalently, we need to solve $(A-I)x=0$, together with the constraints that all $x_i\geq 0$ and $\sum_ix_i=1$.

A=Matrix(QQ,[[1/10-1,3/10,0,0],[0,4/10-1,1,0],[9/10,3/10,0-1,0],[1,1,1,1]])
A

[-9/10 3/10 0 0] [ 0 -3/5 1 0] [ 9/10 3/10 -1 0] [ 1 1 1 1]

A.echelon_form()

[ 1 0 0 5/29] [ 0 1 0 15/29] [ 0 0 1 9/29] [ 0 0 0 0]

We can read our solution directly off of the last column (excluding the last 0).  How does this stack up against our numerical answer?

u=transpose(Matrix(QQ,[5/29,15/29,9/29]))
A=Matrix(QQ,[[1/10,3/10,0],[0,4/10,1],[9/10,3/10,0]])
x=transpose(Matrix(QQ,[1/2,1/2,0]))
v=(A^1000)*x
dist(u,v)

0.0

As far as the computer's floating point calculations can tell, the vectors are the same.  What can the rational arithmetic tell us?

First, let's get the square of the distance between the two (the computed fixed point and the one we estimated by powers of $A$).

d=(transpose(u-v)*(u-v))[0,0]

Now, $d$ is a rational number but the numerator and denominator are large integers.  How large?  Well, let's count the digits in each.  We can do this by converting to a string and finding the length of the string, as in

print len(d.numerator().digits())
print len(d.denominator().digits())

1382 2002

This is a very, very small number, on the order of $10^{-620}$.  So, pretty small.

From a pure mathematics standpoint, we still have not proved that the iterates converge to the fixed point.  There is a theorem that does that, though (the Perron-Frobenius Theorem in matrix analysis leads us to a theorem about Markov chains.  PF is one of my top ten favorite theorems!).

Examples from Section 2.5 (Matrix Inverses)

In section 2.5 we begin to compute matrix inverses and study their properties.  In this part, we'll see how SAGE computes matrix inverses, and some of the pitfalls of computing inverses with computers.

So, let's start with the matrix in exercise 1 (a) page 111.

A=Matrix(QQ,[[1,-2,1],[0,2,0],[-1,0,1]])

We can compute the inverse pretty easily...

A.inverse()

[ 1/2 1/2 -1/2] [ 0 1/2 0] [ 1/2 1/2 1/2]

t=var('t')

Let's look at the rotation matrix from 1 (e).

B=Matrix([[cos(t),-sin(t)],[sin(t),cos(t)]]);
B

[ cos(t) -sin(t)] [ sin(t) cos(t)]

Inverting, we find...

Binv=B.inverse();
Binv

[-sin(t)^2/((sin(t)^2/cos(t) + cos(t))*cos(t)^2) + 1/cos(t) sin(t)/((sin(t)^2/cos(t) + cos(t))*cos(t))] [ -sin(t)/((sin(t)^2/cos(t) + cos(t))*cos(t)) 1/(sin(t)^2/cos(t) + cos(t))]

Ummm... that's not quite what we expected from our calculations by hand.  Well, whatever!  Let's try plugging in some angles of rotation and see what happens...

Binv(t=pi/3)

[ 1/2 1/2*sqrt(3)] [-1/2*sqrt(3) 1/2]

Looks great.  What happens if we try to rotate by a right angle, $\pi/2$ radians?

Binv(t=pi/2)

Traceback (most recent call last): File "<stdin>", line 1, in <module> File "_sage_input_21.py", line 9, in <module> exec compile(ur'open("___code___.py","w").write("# -*- coding: utf-8 -*-\n" + _support_.preparse_worksheet_cell(base64.b64decode("Qmludih0PXBpLzIp"),globals())+"\n"); execfile(os.path.abspath("___code___.py"))' + '\n', '', 'single') File "", line 1, in <module> File "/tmp/tmpU9J44S/___code___.py", line 3, in <module> exec compile(ur'Binv(t=pi/_sage_const_2 )' + '\n', '', 'single') File "", line 1, in <module> File "matrix_symbolic_dense.pyx", line 572, in sage.matrix.matrix_symbolic_dense.Matrix_symbolic_dense.__call__ (sage/matrix/matrix_symbolic_dense.c:3802) File "expression.pyx", line 3370, in sage.symbolic.expression.Expression.__call__ (sage/symbolic/expression.cpp:15477) File "ring.pyx", line 624, in sage.symbolic.ring.SymbolicRing._call_element_ (sage/symbolic/ring.cpp:6563) File "expression.pyx", line 3221, in sage.symbolic.expression.Expression.substitute (sage/symbolic/expression.cpp:14867) RuntimeError: power::eval(): division by zero

Looking back at the expression for $B^{-1}$, we see that ther are some $\cos(t)$ terms in the denominators of some expressions, and without simplifying first, we do indeed have division by zero.

So here is one pitfall in inverting symbolic matrices: your CAS might not be very good at simplifying expressions in the way that we usually would.

Let's try exercise 13 (b).  Here, we have a symbolic expression in our matrix, and we would like to know which values of that symbol lead to an invertible matrix.

k=var('k')
C=Matrix([[1,0,1],[0,1,1],[k,0,1]])

C.echelon_form()

[1 0 0] [0 1 0] [0 0 1]

So, according to SAGE, the value of $k$ does not matter here.  We will always be able to row-reduce to the identity, and thus our matrix $C$ will always be invertible.  But wait... just look at the matrix $C$:

C

[1 0 1] [0 1 1] [k 0 1]

If $k=1$ then the first and last rows are the same.  Thus, the first row reduction step can zero out the third row, leaving us with a matrix that has less than full rank, and is thus not invertible.

So why did SAGE tell us that $C$ row reduces to the identity?  Most computer algebra systems assume that algebraic variables are generic, meaning that the operations we employ only have to be valid almost all of the time.  (Technically, this operation is valid generically in $\mathbb{R}^n$ means that the operation works for all possible values of the variables involved except for values lying on an algebraic variety, which is an intersection of a bunch of curves defined by polynomials.  So, mathematicians are pretty precise about what we mean when we say generically.)

But this can mean that determining whether a matrix with variable quantities in it is invertible can be pretty difficult with a CAS (unless you get used to playing around with it).  Using determinants from section 2.6 can help clear this up:

C.det()

-k + 1

So, this says to us that the determinant is nonzero as long as $k\neq 1$.  That's exactly what we needed to know!

As it turns out, I disagree just a little with what the book says about determinants in the beginning of section 2.6.  Determinants are pretty important for applied mathematics, in the sense that they provide fast and unambiguous characterizations of the invertibility of matrices, which can be especially useful when dealing with computing devices.  Looking back at our rotation matrix,

B

[ cos(t) -sin(t)] [ sin(t) cos(t)]

we have a determinant of

B.det()

sin(t)^2 + cos(t)^2