MAT 232 Ch. 2

Section 2.1

Sage can help us check all the computations with matrices from this section.

M = matrix(2, [ 1, 2, 3, 4]) N = matrix(2, [-1, 0, 0, 1]) O = matrix(2, [-1, 0, 3, 4, 5, 7])
show(M)
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 1 & 2 \\ 3 & 4 \end{array}\right)
show(3*M)
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 3 & 6 \\ 9 & 12 \end{array}\right)
show(M+N)
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 0 & 2 \\ 3 & 5 \end{array}\right)

But be sure to only add matrices that actually belong together!

show(M+O)
 Traceback (click to the left of this block for traceback) ... TypeError: unsupported operand parent(s) for '+': 'Full MatrixSpace of 2 by 2 dense matrices over Integer Ring' and 'Full MatrixSpace of 2 by 3 dense matrices over Integer Ring' Traceback (most recent call last): File "", line 1, in File "_sage_input_10.py", line 10, in exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("c2hvdyhNK08p"),globals())+"\\n"); execfile(os.path.abspath("___code___.py")) File "", line 1, in File "/tmp/tmpLSeqop/___code___.py", line 2, in exec compile(u'show(M+O) File "", line 1, in File "sage/structure/element.pyx", line 1365, in sage.structure.element.ModuleElement.__add__ (/usr/local/sage-6.9/src/build/cythonized/sage/structure/element.c:11960) File "sage/structure/coerce.pyx", line 1070, in sage.structure.coerce.CoercionModel_cache_maps.bin_op (/usr/local/sage-6.9/src/build/cythonized/sage/structure/coerce.c:9739) TypeError: unsupported operand parent(s) for '+': 'Full MatrixSpace of 2 by 2 dense matrices over Integer Ring' and 'Full MatrixSpace of 2 by 3 dense matrices over Integer Ring'

Let's compare a matrix and its transpose, and then whether a matrix is symmetric.

show(M) show(M.transpose())
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 1 & 2 \\ 3 & 4 \end{array}\right) \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 1 & 3 \\ 2 & 4 \end{array}\right)
M.is_symmetric()
 False False

Section 2.2

The computer can help us to visualize vectors too.  In Sage, there is no distinction between column and row vectors, and indeed they are seen as not intrinsically the same as matrices.  We only write them as the same in this class for convenience, really.

V = vector( [1,3] )

We should often think of vectors as a geometric point.

plot(V, figsize=3)

Here we are multiplying a matrix by a vector.  Is it the same as you had hoped for?

M*V
 (7, 15) (7, 15)
plot(M*V, color='red')+plot(V, figsize=3)

A good way to try to visualize what matrices do, then, is to view multiplying by a matrix as the same thing as moving vectors around.  Rather than just thinking of it as a convenient way to encode a system of linear equations, we can imagine what it does to all vectors at once.  Here is an example.

@interact(layout=[['A','Myvect']]) def _(A = matrix([[1,-1],[0,1]]),Myvect=matrix([[0,-1]])): V = vector([1,3]) W = A*V V1 = vector([-1,2]) W1 = A*V1 L = line([V,V1],linestyle='--') LW = line([W,W1],linestyle='--') Myvect = vector(Myvect[0]) Wyvect = A*Myvect pretty_print(html("Here is the product of the matrix and some vectors")) show(A) show(graphics_array([[plot(V)+plot(V1,color='red')+plot(Myvect,color='green')+plot(L),plot(W)+plot(W1,color='red')+plot(Wyvect,color='green')+plot(LW)]]),aspect_ratio=1)

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

We can see this all a little more systematically by moving a box around.

t = var('t') @interact def _(box=matrix(RDF,[[1,2]]), lin_trans=matrix(RDF,[[1,0],[0,1]]),auto_update=False): a = box[0][0] b = box[0][1] P = line([(0,0),(a,0),(a,b),(0,b),(0,0)],color='red')+point((a,b),color='green',size=20) A,B,C,D = lin_trans[0,0],lin_trans[0,1],lin_trans[1,0],lin_trans[1,1] P += line([(0,0),(A*a,C*a),(A*a+B*b,C*a+D*b),(B*b,D*b),(0,0)],color='blue')+point((a,b),color='green',size=20) pretty_print(html('box to $({},{})$'.format(box[0][0],box[0][1]))) pretty_print(html('transformed by $x={}a+{}b$ and $y={}a+{}b$'.format(lin_trans[0,0],lin_trans[0,1],lin_trans[1,0],lin_trans[1,1]))) P.show(aspect_ratio=1,figsize=4)

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

Finally, we should try to visualize the solutions to a non-homogeneous system.

This is what we mean by the solution to a homogeneous system.  In this case, it's just the solutions to the "system" $2x+y=0$, which is multiples of the vector $\left(\begin{array}{r}1 \\-2\end{array}\right)$.

var('s,t') V = vector([1,-2]) P = parametric_plot(t*V,(t,-2,2)) P += V.plot(width=5,color='red') show(P,figsize=4)

But what about solutions to the system $2x+y=3$?  We can see by "inspection" (i.e. guess and check) that $\left(\begin{array}{r}1 \\1\end{array}\right)$ is a solution, so we add all vectors from the homogeneous system solutions to that vector, to obtain all solutions to the new, non-homogeneous system.

var('s,t') V = vector([1,-2]) P = parametric_plot(t*V,(t,-2,2)) P += V.plot(width=5,color='red') P += parametric_plot(vector([1,1])+t*V,(t,-2,2)) P += vector([1,1]).plot(color='green') show(P)
var('s,t') V = vector([1,-2]) P = parametric_plot(t*V,(t,-2,2)) P += V.plot(width=5,color='red') P += parametric_plot(vector([0,3])+t*V,(t,-2,2)) P += vector([0,3]).plot(color='green') show(P)

The same thing works out great for systems in as many dimensions as you want.

var('s,t') V = vector([1/2,0,1]) W = vector([-1/2,1,0]) P = parametric_plot3d(s*W+t*V,(s,-5,5),(t,-5,5)) P += V.plot(thickness=5,color='red') P += W.plot(thickness=5,color='red') P += parametric_plot3d(vector([2,1,0])+s*W+t*V,(s,-5,5),(t,-5,5)) P += parametric_plot3d(vector([3,0,1])+s*W+t*V,(s,-5,5),(t,-5,5)) P += vector([3,0,1]).plot(thickness=5,color='green') P += vector([2,1,0]).plot(thickness=5,color='green') show(P)

Section 2.3

Let's think about how to instantiate matrix multiplication.

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)

We can verify that multiplication doesn't always work as with numbers - not only is it not commutative, it sometimes doesn't even make sense the other direction!

A = matrix([[0,-1],[1,0]]) B = matrix([[1,2,3],[4,5,6]]) show([A,B])
 \newcommand{\Bold}[1]{\mathbf{#1}}\left[\left(\begin{array}{rr} 0 & -1 \\ 1 & 0 \end{array}\right), \left(\begin{array}{rrr} 1 & 2 & 3 \\ 4 & 5 & 6 \end{array}\right)\right]
show(A*B)
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} -4 & -5 & -6 \\ 1 & 2 & 3 \end{array}\right)
show(B*A)
 Traceback (click to the left of this block for traceback) ... TypeError: unsupported operand parent(s) for '*': 'Full MatrixSpace of 2 by 3 dense matrices over Integer Ring' and 'Full MatrixSpace of 2 by 2 dense matrices over Integer Ring' Traceback (most recent call last): File "", line 1, in File "_sage_input_14.py", line 10, in exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("c2hvdyhCKkEp"),globals())+"\\n"); execfile(os.path.abspath("___code___.py")) File "", line 1, in File "/tmp/tmpZdBPFw/___code___.py", line 2, in exec compile(u'show(B*A) File "", line 1, in File "sage/structure/element.pyx", line 2949, in sage.structure.element.Matrix.__mul__ (/usr/local/sage-6.9/src/build/cythonized/sage/structure/element.c:23250) File "sage/structure/coerce.pyx", line 1070, in sage.structure.coerce.CoercionModel_cache_maps.bin_op (/usr/local/sage-6.9/src/build/cythonized/sage/structure/coerce.c:9739) TypeError: unsupported operand parent(s) for '*': 'Full MatrixSpace of 2 by 3 dense matrices over Integer Ring' and 'Full MatrixSpace of 2 by 2 dense matrices over Integer Ring'

Here is a visual version of an example from class.   Note that we again get things in the opposite order, so to speak.

t = var('t') @interact def _(box=matrix(RDF,[[1,2]]), lin_trans=matrix(RDF,[[1,0],[0,1]]), lin_trans2=matrix(RDF,[[1,0],[0,1]]),auto_update=False): a = box[0][0] b = box[0][1] P = line([(0,0),(a,0),(a,b),(0,b),(0,0)],color='red')+point((a,b),color='green',size=20) A,B,C,D = lin_trans[0,0],lin_trans[0,1],lin_trans[1,0],lin_trans[1,1] P += line([(0,0),(A*a,C*a),(A*a+B*b,C*a+D*b),(B*b,D*b),(0,0)],color='blue') LT = lin_trans2*lin_trans P2 = line([(0,0),(a,0),(a,b),(0,b),(0,0)],color='red')+point((a,b),color='green',size=20) A2,B2,C2,D2 = LT[0,0],LT[0,1],LT[1,0],LT[1,1] P2 += line([(0,0),(A2*a,C2*a),(A2*a+B2*b,C2*a+D2*b),(B2*b,D2*b),(0,0)],color='blue') graphics_array([P,P2]).show(figsize=4)

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

And of course, here is Smiley Guy!

t = var('t') @interact def _(auto_update=False,A=matrix(RDF,[[1,0],[0,1]])): pll=A*vector((-0.5,0.5)) plr=A*vector((-0.3,0.5)) prl=A*vector((0.3,0.5)) prr=A*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(A*vector([t, -0.15*sin(2*pi*t)-0.5]), (t, -0.5, 0),color='red')+parametric_plot(A*vector([t, -0.15*sin(2*pi*t)-0.5]), (t,0,0.5),color='orange') face=parametric_plot(A*vector([cos(t),sin(t)]), (t,0,pi/2),color='black')+parametric_plot(A*vector([cos(t),sin(t)]), (t,pi/2,pi),color='lavender')+parametric_plot(A*vector([cos(t),sin(t)]), (t,pi,3*pi/2),color='cyan')+parametric_plot(A*vector([cos(t),sin(t)]),(t,3*pi/2,2*pi),color='sienna') P=right_eye+left_eye+face+mouth pretty_print(html('smiley guy transformed by $A$')) P.show(aspect_ratio=1, figsize=4)

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

Here is Smiley Guy transformed by two different matrices.

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 pretty_print(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

Just because he is so fun, here are a few more variations I've used in the past.

Rotation Smiley Guy:

t = var('t') @interact def _(theta=pi/2,phi=pi/3,auto_update=False): A=matrix(RDF,[[cos(theta),-sin(theta)],[sin(theta),cos(theta)]]) B=matrix(RDF,[[cos(phi),-sin(phi)],[sin(phi),cos(phi)]]) 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 pretty_print(html('Matrices are ${}$ and ${}$'.format(latex(A),latex(B)))) pretty_print(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

Not so Smiley Guy, courtesy of Natou Rudenberg:

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.4,0.47)) prr=M*vector((0.-0.4,0.43)) left_eye=parametric_plot(M*vector([t, 0.3+(0.04-(t+0.4)^2)^(1/2)]), (t, -0.5, -0.3))+parametric_plot(M*vector([t, 0.6-(0.04-(t+0.4)^2)^(1/2)]), (t, -0.5,-0.3))+point(prr,size=8) right_eye=parametric_plot(M*vector([t, 0.3+(0.04-(t-0.4)^2)^(1/2)]), (t, 0.3, 0.5),color='green')+parametric_plot(M*vector([t, 0.6-(0.04-(t-0.4)^2)^(1/2)]), (t, 0.3, 0.5),color='green')+point(prl,size=8) mouth=parametric_plot(M*vector([t, -1.5+(1-t^2)^(1/2)]), (t, -0.5, 0.5),color='orange') nose1 = M*vector((0.2,-0.1)) nose2 = M*vector((0.0,-0.1)) nose3 = M*vector((0,0.3)) nose_bottom = line([nose1,nose2]) nose = line([nose1,nose3]) 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='red')+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+nose_bottom+nose pretty_print(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

Smiley Guy of the transpose.  Notice that although it is using the transpose matrix, we also multiply all the matrices on the right, so there no difference in behavior!  This is because Sage interprets all vectors as columns or vectors, as needed for things to be defined properly.

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 pretty_print(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)

Section 2.4

Can we undo what we have done to Smiley Guy?

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 pretty_print(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

In this cell, I try to do it automatically.

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 pretty_print(html('smiley guy, then transformed by $A$, and next by $A^{-1}$, which is')) show(latex(B)) G = graphics_array([[makeface(ID),makeface(A),makeface(B*A)]]) G.show(aspect_ratio=1) except: pretty_print(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: Matrix is singular Traceback (most recent call last): File "", line 1, in File "_sage_input_63.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/tmpQutDPS/___code___.py", line 3, in exec compile(u'B**-_sage_const_1 File "", line 1, in File "sage/matrix/matrix_integer_dense.pyx", line 1042, in sage.matrix.matrix_integer_dense.Matrix_integer_dense.__pow__ (/usr/local/sage-6.9/src/build/cythonized/sage/matrix/matrix_integer_dense.c:11506) File "sage/matrix/matrix_integer_dense.pyx", line 3851, in sage.matrix.matrix_integer_dense.Matrix_integer_dense.__invert__ (/usr/local/sage-6.9/src/build/cythonized/sage/matrix/matrix_integer_dense.c:33250) File "sage/matrix/matrix_integer_dense.pyx", line 3812, in sage.matrix.matrix_integer_dense.Matrix_integer_dense._invert_flint (/usr/local/sage-6.9/src/build/cythonized/sage/matrix/matrix_integer_dense.c:33084) ZeroDivisionError: Matrix is singular

Here we confirm Example 2.4.6 from the book.

M = matrix(3, [2, 7, 1, 1, 4, -1, 1, 3, 0]) show(M^-1)
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} -\frac{3}{2} & -\frac{3}{2} & \frac{11}{2} \\ \frac{1}{2} & \frac{1}{2} & -\frac{3}{2} \\ \frac{1}{2} & -\frac{1}{2} & -\frac{1}{2} \end{array}\right)
M * M^-1
 [1 0 0] [0 1 0] [0 0 1] [1 0 0] [0 1 0] [0 0 1]

We can solve the linear system $$\begin{array}{lr}2x+7y+z & =1 \\x+4y-z & =2\\ x+3y+0z & =3\end{array}$$ by using an inverse matrix.

M^-1 * vector([1,2,3])
 (12, -3, -2) (12, -3, -2)

We check our work:

M * vector([12, -3, -2])
 (1, 2, 3) (1, 2, 3)

Chapter 2.7

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.

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 might do on our own?

B = matrix([[1,0],[3,1]]) C = matrix([[1,2],[0,-2]])
show(B);show(C); show(B*C)
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 1 & 0 \\ 3 & 1 \end{array}\right) \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 1 & 2 \\ 0 & -2 \end{array}\right) \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rr} 1 & 2 \\ 3 & 4 \end{array}\right)

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 wouldn'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_5.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/tmp6HQH6y/___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 "sage/matrix/matrix2.pyx", line 11584, in sage.matrix.matrix2.Matrix.LU (/usr/local/sage-6.9/src/build/cythonized/sage/matrix/matrix2.c:79561) 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)

For another example, let's see this one from another text.  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)
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrrr} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -\frac{1}{4} & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ -\frac{1}{4} & -\frac{1}{15} & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & -\frac{4}{15} & -\frac{2}{7} & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & -\frac{15}{56} & -\frac{1}{12} & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & -\frac{7}{24} & -\frac{26}{89} & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & -\frac{24}{89} & -\frac{208}{2415} & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & -\frac{712}{2415} & -\frac{2623}{8948} & 1 \end{array}\right) \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrrrr} 4 & -1 & -1 & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{15}{4} & -\frac{1}{4} & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{56}{15} & -\frac{16}{15} & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{24}{7} & -\frac{2}{7} & -1 & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{89}{24} & -\frac{13}{12} & -1 & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{2415}{712} & -\frac{26}{89} & -1 \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{8948}{2415} & -\frac{2623}{2415} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{30305}{8948} \end{array}\right)

Here is how this would be used in practice.

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)