Chapter 2 - Matrices

Chapter 2.1

How can I compute with matrices?

First, the motivation.  Several of you asked whether Smiley Guy could be transformed by one thing, and then another.  Let's see!

t = var('t') @interact(layout=[['A','B'],['auto_update']]) def _(A=matrix(RDF,[[1,0],[0,1]]),B=matrix(RDF,[[1,0],[0,1]]),auto_update=False): ID = matrix(RDF,[[1,0],[0,1]]) def makeface(M): pll=M*vector((-0.5,0.5)) plr=M*vector((-0.3,0.5)) prl=M*vector((0.3,0.5)) prr=M*vector((0.5,0.5)) left_eye=line([pll,plr])+point(pll,size=5)+point(plr,size=5) right_eye=line([prl,prr],color='green')+point(prl,size=5,color='green')+point(prr,size=5,color='green') mouth=parametric_plot(M*vector([t, -0.15*sin(2*pi*t)-0.5]), (t, -0.5, 0),color='red')+parametric_plot(M*vector([t, -0.15*sin(2*pi*t)-0.5]), (t,0,0.5),color='orange') face=parametric_plot(M*vector([cos(t),sin(t)]), (t,0,pi/2),color='black')+parametric_plot(M*vector([cos(t),sin(t)]), (t,pi/2,pi),color='lavender')+parametric_plot(M*vector([cos(t),sin(t)]), (t,pi,3*pi/2),color='cyan')+parametric_plot(M*vector([cos(t),sin(t)]),(t,3*pi/2,2*pi),color='sienna') return right_eye+left_eye+face+mouth html('smiley guy, then transformed by $A$, and next by $B$') G = graphics_array([[makeface(ID),makeface(A),makeface(B*A)]]) G.show(aspect_ratio=1)

Click to the left again to hide and once more to show the dynamic interactive window

V = vector([1,2]) A = matrix([[1,0],[0,-1]]) B = matrix([[0,1],[1,0]])

How do we apply a transformation to a vector?  We multiply the matrix on the left.

A*V
 (1, -2) (1, -2)

How do we apply the next transformation?  We multiply by the next matrix, on the left again.

B * (A*V)
 (-2, 1) (-2, 1)

Well, it turns out that this is associative; we could just multiply by B times A.

B*A
 [ 0 -1] [ 1 0] [ 0 -1] [ 1 0]
(B*A)*V
 (-2, 1) (-2, 1)

Chapter 2.2

Let's see what happens with Smiley Guy and the tranpose.

t = var('t') @interact(layout=[['A','B'],['auto_update']]) def _(A=matrix(RDF,[[1,0],[0,1]]),B=matrix(RDF,[[1,0],[0,1]]),auto_update=False): ID = matrix(RDF,[[1,0],[0,1]]) def makeface(M): pll=vector((-0.5,0.5))*M plr=vector((-0.3,0.5))*M prl=vector((0.3,0.5))*M prr=vector((0.5,0.5))*M left_eye=line([pll,plr])+point(pll,size=5)+point(plr,size=5) right_eye=line([prl,prr],color='green')+point(prl,size=5,color='green')+point(prr,size=5,color='green') mouth=parametric_plot(vector([t, -0.15*sin(2*pi*t)-0.5])*M, (t, -0.5, 0),color='red')+parametric_plot(vector([t, -0.15*sin(2*pi*t)-0.5])*M, (t,0,0.5),color='orange') face=parametric_plot(vector([cos(t),sin(t)])*M, (t,0,pi/2),color='black')+parametric_plot(vector([cos(t),sin(t)])*M, (t,pi/2,pi),color='lavender')+parametric_plot(vector([cos(t),sin(t)])*M, (t,pi,3*pi/2),color='cyan')+parametric_plot(vector([cos(t),sin(t)])*M,(t,3*pi/2,2*pi),color='sienna') return right_eye+left_eye+face+mouth html('smiley guy, then transformed by $A$, and next by $B$') G = graphics_array([[makeface(ID),makeface(A.transpose()),makeface(A.transpose()*B.transpose())]]) G.show(aspect_ratio=1)

Click to the left again to hide and once more to show the dynamic interactive window

Sage interprets a vector as a row or a column, whichever is appropriate.

Also, we can undo Smiley Guy's problems.  We go back to the original setup.

t = var('t') @interact def _(A=matrix(RDF,[[1,0],[0,1]]),auto_update=False): ID = matrix(RDF,[[1,0],[0,1]]) try: B = A^-1 def makeface(M): pll=M*vector((-0.5,0.5)) plr=M*vector((-0.3,0.5)) prl=M*vector((0.3,0.5)) prr=M*vector((0.5,0.5)) left_eye=line([pll,plr])+point(pll,size=5)+point(plr,size=5) right_eye=line([prl,prr],color='green')+point(prl,size=5,color='green')+point(prr,size=5,color='green') mouth=parametric_plot(M*vector([t, -0.15*sin(2*pi*t)-0.5]), (t, -0.5, 0),color='red')+parametric_plot(M*vector([t, -0.15*sin(2*pi*t)-0.5]), (t,0,0.5),color='orange') face=parametric_plot(M*vector([cos(t),sin(t)]), (t,0,pi/2),color='black')+parametric_plot(M*vector([cos(t),sin(t)]), (t,pi/2,pi),color='lavender')+parametric_plot(M*vector([cos(t),sin(t)]), (t,pi,3*pi/2),color='cyan')+parametric_plot(M*vector([cos(t),sin(t)]),(t,3*pi/2,2*pi),color='sienna') return right_eye+left_eye+face+mouth html('smiley guy, then transformed by $A$, and next by $A^{-1}$') G = graphics_array([[makeface(ID),makeface(A),makeface(B*A)]]) G.show(aspect_ratio=1) except: html("Oops! Your matrix $%s$ has no inverse!"%latex(A))

Click to the left again to hide and once more to show the dynamic interactive window

If you ever want to compute an inverse, use Sage.

A = matrix([ [1,2,3] , [4,5,7] , [0,8,2] ])
A^-1
 [-23/17 10/17 -1/34] [ -4/17 1/17 5/34] [ 16/17 -4/17 -3/34] [-23/17 10/17 -1/34] [ -4/17 1/17 5/34] [ 16/17 -4/17 -3/34]

Well, make sure it has one!

B = matrix(3,[1,2,3,4,5,6,7,8,9]); B
 [1 2 3] [4 5 6] [7 8 9] [1 2 3] [4 5 6] [7 8 9]
B^-1
 Traceback (click to the left of this block for traceback) ... ZeroDivisionError: input matrix must be nonsingular Traceback (most recent call last): File "", line 1, in File "_sage_input_19.py", line 10, in exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("Ql4tMQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py")) File "", line 1, in File "/tmp/tmpiVDEES/___code___.py", line 3, in exec compile(u'B**-_sage_const_1 File "", line 1, in File "matrix0.pyx", line 4757, in sage.matrix.matrix0.Matrix.__pow__ (sage/matrix/matrix0.c:26582) File "element.pyx", line 1776, in sage.structure.element.RingElement.__pow__ (sage/structure/element.c:14627) File "element.pyx", line 3371, in sage.structure.element.generic_power_c (sage/structure/element.c:26120) File "matrix0.pyx", line 4657, in sage.matrix.matrix0.Matrix.__invert__ (sage/matrix/matrix0.c:25938) File "matrix_rational_dense.pyx", line 650, in sage.matrix.matrix_rational_dense.Matrix_rational_dense.__invert__ (sage/matrix/matrix_rational_dense.c:9364) File "matrix_rational_dense.pyx", line 756, in sage.matrix.matrix_rational_dense.Matrix_rational_dense.__invert__main (sage/matrix/matrix_rational_dense.c:9962) File "matrix_rational_dense.pyx", line 2655, in sage.matrix.matrix_rational_dense.Matrix_rational_dense._invert_pari (sage/matrix/matrix_rational_dense.c:24339) ZeroDivisionError: input matrix must be nonsingular

Also, be careful with decimals.

A = matrix(RDF,[[1,2,3],[4,5,7],[0,8,2]])
 [1.0 2.0 3.0] [4.0 5.0 7.0] [0.0 8.0 2.0] [1.0 2.0 3.0] [4.0 5.0 7.0] [0.0 8.0 2.0]
A^-1
 [ -1.35294117647 0.588235294118 -0.0294117647059] [ -0.235294117647 0.0588235294118 0.147058823529] [ 0.941176470588 -0.235294117647 -0.0882352941176] [ -1.35294117647 0.588235294118 -0.0294117647059] [ -0.235294117647 0.0588235294118 0.147058823529] [ 0.941176470588 -0.235294117647 -0.0882352941176]
A*A^-1
 [ 1.0 2.77555756156e-17 4.16333634234e-17] [-3.88578058619e-16 1.0 1.11022302463e-16] [ 0.0 0.0 1.0] [ 1.0 2.77555756156e-17 4.16333634234e-17] [-3.88578058619e-16 1.0 1.11022302463e-16] [ 0.0 0.0 1.0]

What the heck????  So ... is this an inverse or not?  What significance do those tiny numbers have?

This is why you need to understand the algorithms and theorems - so that things like this don't happen to you and you don't know what's going on.

Chapter 2.3

This is a fun pattern.

@interact def _(n=[1..10]): M = matrix(n,[1 if i >= j else 0 for i in [1..n] for j in [1..n]]) show(M) show(M^-1)

Chapter 2.4

Block matrices are implemented in Sage.

block_matrix?
 File: /usr/local/sage-5.6-linux-64bit-ubuntu_8.04.4_lts-x86_64-Linux/local/lib/python2.7/site-packages/sage/matrix/constructor.py Type: Definition: block_matrix(*args, **kwds) Docstring: Returns a larger matrix made by concatenating submatrices (rows first, then columns). For example, the matrix [ A B ] [ C D ] is made up of submatrices A, B, C, and D. INPUT: The block_matrix command takes a list of submatrices to add as blocks, optionally preceded by a ring and the number of block rows and block columns, and returns a matrix. The submatrices can be specified as a list of matrices (using nrows and ncols to determine their layout), or a list of lists of matrices, where each list forms a row. ring - the base ring nrows - the number of block rows ncols - the number of block cols sub_matrices - matrices (see below for syntax) subdivide - boolean, whether or not to add subdivision information to the matrix sparse - boolean, whether to make the resulting matrix sparse EXAMPLES: sage: A = matrix(QQ, 2, 2, [3,9,6,10]) sage: block_matrix([ [A, -A], [~A, 100*A] ]) [ 3 9| -3 -9] [ 6 10| -6 -10] [-----------+-----------] [-5/12 3/8| 300 900] [ 1/4 -1/8| 600 1000]  If the number of submatrices in each row is the same, you can specify the submatrices as a single list too: sage: block_matrix(2, 2, [ A, A, A, A ]) [ 3 9| 3 9] [ 6 10| 6 10] [-----+-----] [ 3 9| 3 9] [ 6 10| 6 10]  One can use constant entries: sage: block_matrix([ [1, A], [0, 1] ]) [ 1 0| 3 9] [ 0 1| 6 10] [-----+-----] [ 0 0| 1 0] [ 0 0| 0 1]  A zero entry may represent any square or non-square zero matrix: sage: B = matrix(QQ, 1, 1, [ 1 ] ) sage: C = matrix(QQ, 2, 2, [ 2, 3, 4, 5 ] ) sage: block_matrix([ [B, 0], [0, C] ]) [1|0 0] [-+---] [0|2 3] [0|4 5]  One can specify the number of rows or columns as keywords too: sage: block_matrix([A, -A, ~A, 100*A], ncols=4) [ 3 9| -3 -9|-5/12 3/8| 300 900] [ 6 10| -6 -10| 1/4 -1/8| 600 1000] sage: block_matrix([A, -A, ~A, 100*A], nrows=1) [ 3 9| -3 -9|-5/12 3/8| 300 900] [ 6 10| -6 -10| 1/4 -1/8| 600 1000]  It handles base rings nicely too: sage: R. = ZZ['x'] sage: block_matrix(2, 2, [1/2, A, 0, x-1]) [ 1/2 0| 3 9] [ 0 1/2| 6 10] [-----------+-----------] [ 0 0|x - 1 0] [ 0 0| 0 x - 1] sage: block_matrix(2, 2, [1/2, A, 0, x-1]).parent() Full MatrixSpace of 4 by 4 dense matrices over Univariate Polynomial Ring in x over Rational Field  Subdivisions are optional. If they are disabled, the columns need not line up: sage: B = matrix(QQ, 2, 3, range(6)) sage: block_matrix([ [~A, B], [B, ~A] ], subdivide=False) [-5/12 3/8 0 1 2] [ 1/4 -1/8 3 4 5] [ 0 1 2 -5/12 3/8] [ 3 4 5 1/4 -1/8]  Without subdivisions it also deduces dimensions for scalars if possible: sage: C = matrix(ZZ, 1, 2, range(2)) sage: block_matrix([ [ C, 0 ], [ 3, 4 ], [ 5, 6, C ] ], subdivide=False ) [0 1 0 0] [3 0 4 0] [0 3 0 4] [5 6 0 1]  If all submatrices are sparse (unless there are none at all), the result will be a sparse matrix. Otherwise it will be dense by default. The sparse keyword can be used to override this: sage: A = Matrix(ZZ, 2, 2, [0, 1, 0, 0], sparse=True) sage: block_matrix([ [ A ], [ A ] ]).parent() Full MatrixSpace of 4 by 2 sparse matrices over Integer Ring sage: block_matrix([ [ A ], [ A ] ], sparse=False).parent() Full MatrixSpace of 4 by 2 dense matrices over Integer Ring  Consecutive zero submatrices are consolidated. sage: B = matrix(2, range(4)) sage: C = matrix(2, 8, range(16)) sage: block_matrix(2, [[B,0,0,B],[C]], subdivide=False) [ 0 1 0 0 0 0 0 1] [ 2 3 0 0 0 0 2 3] [ 0 1 2 3 4 5 6 7] [ 8 9 10 11 12 13 14 15]  Ambiguity is not tolerated. sage: B = matrix(2, range(4)) sage: C = matrix(2, 8, range(16)) sage: block_matrix(2, [[B,0,B,0],[C]], subdivide=False) Traceback (click to the left of this block for traceback) ...  File: /usr/local/sage-5.6-linux-64bit-ubuntu_8.04.4_lts-x86_64-Linux/local/lib/python2.7/site-packages/sage/matrix/constructor.py Type: Definition: block_matrix(*args, **kwds) Docstring: Returns a larger matrix made by concatenating submatrices (rows first, then columns). For example, the matrix [ A B ] [ C D ] is made up of submatrices A, B, C, and D. INPUT: The block_matrix command takes a list of submatrices to add as blocks, optionally preceded by a ring and the number of block rows and block columns, and returns a matrix. The submatrices can be specified as a list of matrices (using nrows and ncols to determine their layout), or a list of lists of matrices, where each list forms a row. ring - the base ring nrows - the number of block rows ncols - the number of block cols sub_matrices - matrices (see below for syntax) subdivide - boolean, whether or not to add subdivision information to the matrix sparse - boolean, whether to make the resulting matrix sparse EXAMPLES: sage: A = matrix(QQ, 2, 2, [3,9,6,10]) sage: block_matrix([ [A, -A], [~A, 100*A] ]) [ 3 9| -3 -9] [ 6 10| -6 -10] [-----------+-----------] [-5/12 3/8| 300 900] [ 1/4 -1/8| 600 1000]  If the number of submatrices in each row is the same, you can specify the submatrices as a single list too: sage: block_matrix(2, 2, [ A, A, A, A ]) [ 3 9| 3 9] [ 6 10| 6 10] [-----+-----] [ 3 9| 3 9] [ 6 10| 6 10]  One can use constant entries: sage: block_matrix([ [1, A], [0, 1] ]) [ 1 0| 3 9] [ 0 1| 6 10] [-----+-----] [ 0 0| 1 0] [ 0 0| 0 1]  A zero entry may represent any square or non-square zero matrix: sage: B = matrix(QQ, 1, 1, [ 1 ] ) sage: C = matrix(QQ, 2, 2, [ 2, 3, 4, 5 ] ) sage: block_matrix([ [B, 0], [0, C] ]) [1|0 0] [-+---] [0|2 3] [0|4 5]  One can specify the number of rows or columns as keywords too: sage: block_matrix([A, -A, ~A, 100*A], ncols=4) [ 3 9| -3 -9|-5/12 3/8| 300 900] [ 6 10| -6 -10| 1/4 -1/8| 600 1000] sage: block_matrix([A, -A, ~A, 100*A], nrows=1) [ 3 9| -3 -9|-5/12 3/8| 300 900] [ 6 10| -6 -10| 1/4 -1/8| 600 1000]  It handles base rings nicely too: sage: R. = ZZ['x'] sage: block_matrix(2, 2, [1/2, A, 0, x-1]) [ 1/2 0| 3 9] [ 0 1/2| 6 10] [-----------+-----------] [ 0 0|x - 1 0] [ 0 0| 0 x - 1] sage: block_matrix(2, 2, [1/2, A, 0, x-1]).parent() Full MatrixSpace of 4 by 4 dense matrices over Univariate Polynomial Ring in x over Rational Field  Subdivisions are optional. If they are disabled, the columns need not line up: sage: B = matrix(QQ, 2, 3, range(6)) sage: block_matrix([ [~A, B], [B, ~A] ], subdivide=False) [-5/12 3/8 0 1 2] [ 1/4 -1/8 3 4 5] [ 0 1 2 -5/12 3/8] [ 3 4 5 1/4 -1/8]  Without subdivisions it also deduces dimensions for scalars if possible: sage: C = matrix(ZZ, 1, 2, range(2)) sage: block_matrix([ [ C, 0 ], [ 3, 4 ], [ 5, 6, C ] ], subdivide=False ) [0 1 0 0] [3 0 4 0] [0 3 0 4] [5 6 0 1]  If all submatrices are sparse (unless there are none at all), the result will be a sparse matrix. Otherwise it will be dense by default. The sparse keyword can be used to override this: sage: A = Matrix(ZZ, 2, 2, [0, 1, 0, 0], sparse=True) sage: block_matrix([ [ A ], [ A ] ]).parent() Full MatrixSpace of 4 by 2 sparse matrices over Integer Ring sage: block_matrix([ [ A ], [ A ] ], sparse=False).parent() Full MatrixSpace of 4 by 2 dense matrices over Integer Ring  Consecutive zero submatrices are consolidated. sage: B = matrix(2, range(4)) sage: C = matrix(2, 8, range(16)) sage: block_matrix(2, [[B,0,0,B],[C]], subdivide=False) [ 0 1 0 0 0 0 0 1] [ 2 3 0 0 0 0 2 3] [ 0 1 2 3 4 5 6 7] [ 8 9 10 11 12 13 14 15]  Ambiguity is not tolerated. sage: B = matrix(2, range(4)) sage: C = matrix(2, 8, range(16)) sage: block_matrix(2, [[B,0,B,0],[C]], subdivide=False) Traceback (most recent call last): ... ValueError: insufficient information to determine submatrix widths  Historically, giving only a flat list of submatrices, whose number was a perfect square, would create a new matrix by laying out the submatrices in a square grid. This behavior is now deprecated. sage: A = matrix(2, 3, range(6)) sage: B = matrix(3, 3, range(9)) sage: block_matrix([A, A, B, B]) doctest:...: DeprecationWarning: invocation of block_matrix with just a list whose length is a perfect square is deprecated. See the documentation for details. [0 1 2|0 1 2] [3 4 5|3 4 5] [-----+-----] [0 1 2|0 1 2] [3 4 5|3 4 5] [6 7 8|6 7 8]  Historically, a flat list of matrices whose number is not a perfect square, with no specification of the number of rows or columns, would raise an error. This behavior continues, but could be removed when the deprecation above is completed. sage: A = matrix(2, 3, range(6)) sage: B = matrix(3, 3, range(9)) sage: block_matrix([A, A, A, B, B, B]) Traceback (most recent call last): ... ValueError: must specify nrows or ncols for non-square block matrix. 
A = matrix([[1,2],[4,2]]) B = block_matrix([ [1, A], [0, 1] ]) B
 [1 0|1 2] [0 1|4 2] [---+---] [0 0|1 0] [0 0|0 1] [1 0|1 2] [0 1|4 2] [---+---] [0 0|1 0] [0 0|0 1]

Chapter 2.5

Sage has LU decomposition (and many others).

M = matrix([[1,2],[3,4]]); M
 [1 2] [3 4] [1 2] [3 4]
M.LU()
 ( [0 1] [ 1 0] [ 3 4] [1 0], [1/3 1], [ 0 2/3] ) ( [0 1] [ 1 0] [ 3 4] [1 0], [1/3 1], [ 0 2/3] )

Notice that we get three matrices.  The second and third ones are the L and U matrices; the first one is one that the book refers to about possibly having to permute the rows.

show(M.LU())
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\left(\begin{array}{rr} 0 & 1 \\ 1 & 0 \end{array}\right), \left(\begin{array}{rr} 1 & 0 \\ \frac{1}{3} & 1 \end{array}\right), \left(\begin{array}{rr} 3 & 4 \\ 0 & \frac{2}{3} \end{array}\right)\right)
a,b,c = M.LU()

Here's why we need the permutation.  I hope this is clearer than the book's murky allusions.

b*c
 [3 4] [1 2] [3 4] [1 2]
a*b*c
 [1 2] [3 4] [1 2] [3 4]

Why didn't Sage just do what we did?

B = matrix([[1,0],[3,1]]) C = matrix([[1,2],[0,-2]])
show(B);show(C); show(B*C)

The reason is because in numerical routines, we would rather first switch the rows, so that the pivot is 3, not 1 (the biggest element in the column - partial pivoting).  But for our case, we can do it by hand, so we didn't rely on this technique.

LU isn't implemented for all matrices, perhaps since the underlying machinery is not available in arbitrary precision.

N = matrix([[1.,2.],[3.,4.]]); N.LU()
 Traceback (click to the left of this block for traceback) ... TypeError: base ring of the matrix must be exact, not Real Field with 53 bits of precision Traceback (most recent call last): File "", line 1, in File "_sage_input_19.py", line 10, in exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("TiA9IG1hdHJpeChbWzEuLDIuXSxbMy4sNC5dXSk7IE4uTFUoKQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py")) File "", line 1, in File "/tmp/tmp7SemNq/___code___.py", line 3, in exec compile(u'N = matrix([[_sage_const_1p ,_sage_const_2p ],[_sage_const_3p ,_sage_const_4p ]]); N.LU() File "", line 1, in File "matrix2.pyx", line 10302, in sage.matrix.matrix2.Matrix.LU (sage/matrix/matrix2.c:49573) TypeError: base ring of the matrix must be exact, not Real Field with 53 bits of precision

But with the standard decimals in most computers, we do have it.  This uses the core LAPACK library (via Scipy) that we have mentioned before.

N = matrix(RDF,[[1.,2.],[3.,4.]]); show(N.LU())
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\left(\begin{array}{rr} 0.0 & 1.0 \\ 1.0 & 0.0 \end{array}\right), \left(\begin{array}{rr} 1.0 & 0.0 \\ 0.333333333333 & 1.0 \end{array}\right), \left(\begin{array}{rr} 3.0 & 4.0 \\ 0.0 & 0.666666666667 \end{array}\right)\right)

Let's look at Exercise 31.  I'll use exact matrices, not approximations, though in practice this would be done with decimals.

A = matrix([[4, -1, -1, 0,0,0,0,0], [-1,4,0,-1,0,0,0,0], [-1,0,4,-1,-1,0,0,0], [0,-1,-1,4,0,-1,0,0], [0,0,-1,0,4,-1,-1,0], [0,0,0,-1,-1,4,0,-1], [0,0,0,0,-1,0,4,-1], [0,0,0,0,0,-1,-1,4]]); show(A)
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrrr} 4 & -1 & -1 & 0 & 0 & 0 & 0 & 0 \\ -1 & 4 & 0 & -1 & 0 & 0 & 0 & 0 \\ -1 & 0 & 4 & -1 & -1 & 0 & 0 & 0 \\ 0 & -1 & -1 & 4 & 0 & -1 & 0 & 0 \\ 0 & 0 & -1 & 0 & 4 & -1 & -1 & 0 \\ 0 & 0 & 0 & -1 & -1 & 4 & 0 & -1 \\ 0 & 0 & 0 & 0 & -1 & 0 & 4 & -1 \\ 0 & 0 & 0 & 0 & 0 & -1 & -1 & 4 \end{array}\right)

In this case, there is no permutation needed, so we get the LU decomposition without further ado.

a,L,U = A.LU(); show(L); show(U)
b = vector([5,15,0,10,0,10,20,30])
x = U\(L\b); x
 (827/209, 1377/209, 886/209, 1546/209, 1171/209, 1831/209, 1967/209, 2517/209) (827/209, 1377/209, 886/209, 1546/209, 1171/209, 1831/209, 1967/209, 2517/209)
A*x
 (5, 15, 0, 10, 0, 10, 20, 30) (5, 15, 0, 10, 0, 10, 20, 30)

Yippee, it worked!

Notice how ugly the inverse is.  And imagine this for a much more realistic situation - not eight nodes, but tens of thousands, perhaps while you are trying to model heat changes in an airplane battery to make sure it doesn't catch fire...

show(A^-1)
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrrr} \frac{8948}{30305} & \frac{2623}{30305} & \frac{2864}{30305} & \frac{1544}{30305} & \frac{964}{30305} & \frac{689}{30305} & \frac{303}{30305} & \frac{248}{30305} \\ \frac{2623}{30305} & \frac{8948}{30305} & \frac{1544}{30305} & \frac{2864}{30305} & \frac{689}{30305} & \frac{964}{30305} & \frac{248}{30305} & \frac{303}{30305} \\ \frac{2864}{30305} & \frac{1544}{30305} & \frac{9912}{30305} & \frac{3312}{30305} & \frac{3167}{30305} & \frac{1792}{30305} & \frac{964}{30305} & \frac{689}{30305} \\ \frac{1544}{30305} & \frac{2864}{30305} & \frac{3312}{30305} & \frac{9912}{30305} & \frac{1792}{30305} & \frac{3167}{30305} & \frac{689}{30305} & \frac{964}{30305} \\ \frac{964}{30305} & \frac{689}{30305} & \frac{3167}{30305} & \frac{1792}{30305} & \frac{9912}{30305} & \frac{3312}{30305} & \frac{2864}{30305} & \frac{1544}{30305} \\ \frac{689}{30305} & \frac{964}{30305} & \frac{1792}{30305} & \frac{3167}{30305} & \frac{3312}{30305} & \frac{9912}{30305} & \frac{1544}{30305} & \frac{2864}{30305} \\ \frac{303}{30305} & \frac{248}{30305} & \frac{964}{30305} & \frac{689}{30305} & \frac{2864}{30305} & \frac{1544}{30305} & \frac{8948}{30305} & \frac{2623}{30305} \\ \frac{248}{30305} & \frac{303}{30305} & \frac{689}{30305} & \frac{964}{30305} & \frac{1544}{30305} & \frac{2864}{30305} & \frac{2623}{30305} & \frac{8948}{30305} \end{array}\right)

Chapter 2.6

Just in case you don't believe that the series works...

I = identity_matrix(2) C = matrix([[0.,.5],[.6,.2]]) show((I-C)^-1)
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 1.60000000000000 & 1.00000000000000 \\ 1.20000000000000 & 2.00000000000000 \end{array}\right)
show(sum(C^n for n in [0..10]))
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 1.58835092800000 & 0.987112896000000 \\ 1.18453547520000 & 1.98319608640000 \end{array}\right)

Chapter 2.7

Computer graphics are everywhere.  How did we make Smiley Guy?

t = var('t') @interact() def _(A=matrix(RDF,[[1,0],[0,1]])): def makeface(M): pll=M*vector((-0.5,0.5)) plr=M*vector((-0.3,0.5)) prl=M*vector((0.3,0.5)) prr=M*vector((0.5,0.5)) mouth=parametric_plot(M*vector([t, -0.15*sin(2*pi*t)-0.5]), (t, -0.5, 0),color='red')+parametric_plot(M*vector([t, -0.15*sin(2*pi*t)-0.5]), (t,0,0.5),color='orange') face=parametric_plot(M*vector([cos(t),sin(t)]), (t,0,pi/2),color='black')+parametric_plot(M*vector([cos(t),sin(t)]), (t,pi/2,pi),color='lavender')+parametric_plot(M*vector([cos(t),sin(t)]), (t,pi,3*pi/2),color='cyan')+parametric_plot(M*vector([cos(t),sin(t)]),(t,3*pi/2,2*pi),color='sienna') left_eye=line([pll,plr])+point(pll,size=5)+point(plr,size=5) right_eye=line([prl,prr],color='green')+point(prl,size=5,color='green')+point(prr,size=5,color='green') return right_eye+left_eye+face+mouth html('smiley guy, transformed by $A$') G = makeface(A) G.show(aspect_ratio=1)

Click to the left again to hide and once more to show the dynamic interactive window

Note how nearly everything involves a vector times a matrix.  We are moving points around with a vector.

Let's try an easier example.

@interact() def _(A=matrix(RDF,[[1,0],[0,1]])): def maketriangle(M): vertex1 = M*vector((0,0)) vertex2 = M*vector((1,3)) vertex3 = M*vector((-1,2)) edges = line([vertex1,vertex2,vertex3,vertex1])+point(vertex1,size=40,color='black')+point(vertex2,size=40,color='red')+point(vertex3,size=40,color='green')+point((0,0),size=0) return edges html('triangle, transformed by $%s$'%latex(A)) G = maketriangle(A) + maketriangle(identity_matrix(2)) G.show(aspect_ratio=1)

Click to the left again to hide and once more to show the dynamic interactive window

We specify the vertices, and then we move them with a transformation.  THEN, and only then, do we connect the lines; this way we don't have to think about how to move the lines, just the vertices.

Problem: why can't we move it away from the origin?

Answer: we can - using block matrices and homogeneous coordinates!

@interact() def _(v = matrix([[0],[0]])): A=matrix(RDF,[[1,0],[0,1]]) def maketriangle(M,tr): M = block_matrix([[M,tr],[0,identity_matrix(1)]]) vertex1 = vector((M*vector((0,0,1)) )[:2] ) vertex2 = vector((M*vector((1,3,1)) )[:2] ) vertex3 = vector((M*vector((-1,2,1)) )[:2] ) edges = line([vertex1,vertex2,vertex3,vertex1])+point(vertex1,size=40,color='black')+point(vertex2,size=40,color='red')+point(vertex3,size=40,color='green')+point((0,0),size=0) return edges html('triangle, translated by $%s$'%latex(v)) G = maketriangle(A,v) + maketriangle(identity_matrix(2),matrix([[0],[0]])) G.show(aspect_ratio=1)

Click to the left again to hide and once more to show the dynamic interactive window

We can combine them.  Recall that the order of translation and rotation matters!

@interact def _(v = matrix([[0],[0]]) , A = matrix(RDF,[[1,0],[0,1]]) ): def maketriangle(M,tr): M = block_matrix([[M,tr],[0,identity_matrix(1)]]) vertex1 = vector((M*vector((0,0,1)) )[:2] ) vertex2 = vector((M*vector((1,3,1)) )[:2] ) vertex3 = vector((M*vector((-1,2,1)) )[:2] ) edges = line([vertex1,vertex2,vertex3,vertex1])+point(vertex1,size=40,color='black')+point(vertex2,size=40,color='red')+point(vertex3,size=40,color='green')+point((0,0),size=0) return edges html('triangle, translated by $%s$ then transformed by $%s$'%(latex(v),latex(A))) html('then vice versa') G = maketriangle(A,v) + maketriangle(identity_matrix(2),matrix([[0],[0]])) H = maketriangle(A,A*v) + maketriangle(identity_matrix(2),matrix([[0],[0]])) show(graphics_array([[H,G]]),aspect_ratio=1)