Practica 1 - Introduccion al Algebra de Matrices con SAGE -- Sage

1 de 14

http://sage.math.canterbury.ac.nz/home/albros/15/print

Practica 1 - Introduccion al Algebra de

Matrices con SAGE
Prctica 1 - Introduccin al lgebra de Matrices con SAGE
DEFINIENDO MATRICES I
Recordamos que cuando definimos una Matriz lo hacemos teniendo en cuenta el conjunto
al que pertenecen sus elementos y hablamos del conjunto M nxm(K), donde K es un
conjunto con estructura algebrica de cuerpo (R o C) o anillo (Z o Q).
En SAGE, las matrices se definen, de forma general, con el comando matrix:
M=matrix(conjunto, filas, columnas, elementos),
donde conjunto es ZZ, QQ, QQbar, RDF o RR, CDF, filas y columnas son las
dimensiones de la matriz y elementos puede se una lista, un vector, un rango o nada. En el
ltimo caso obtenemos una matriz nula.

NOTA: en SAGE solo se muestra la salida de la ltima linea. Si queremos mostrar salidas
intermedias utilizamos los comandos print,html o show, esta ltima muestra la salida
en formato enrriquecido. Tenemos tambin la posibilidad de expresar la salida en LATEX .
IMPORTANTE: en SAGE los ndices empiezan en "cero", es decir, el elemento a 11 de una
matriz es A[0,0]
M=matrix(3,4)
print M
print "--------------------------"
M=matrix(3,4,range(12))
print M
print "--------------------------"
M=matrix(3,4,[1 .. 12])
print M
print "--------------------------"
M=matrix(3,4,[7,6,9,2,3,9,-4,-3,8,6,7,1])
show(M)
print "--------------------------"
print latex(M)
print "--------------------------"
html("el elemento $a_{11}$ es M[0,0]: " + M[0,0].str())
print "--------------------------"

OBJETOS Y METODOS

19/10/2013 17:08

Practica 1 - Introduccion al Algebra de Matrices con SAGE -- Sage

2 de 14

http://sage.math.canterbury.ac.nz/home/albros/15/print

En SAGE todos los elmentos que definimos son "objetos" y como tales, tienen "mtodos".
Por ejemplo, en el caso de las matrices estos metodos son: efectuar operaciones elementales
(add_multiple_of_row, swap_rows y rescale_row), obtener el rango (rank), calcular su
forma reducida por filas (echelon_form), determinar todos los menores de un orden dado
(minors), calcular su determinante (det o determinat), descomposicin LU (LU), etc...
Tenemos que tener presente que algunos de estos mtodos producen cambios en la matriz
original (es una forma de reservar espacio en memoria del procesador).
Si queremos saber todos los mtodos de un objeto determinado, escribimos:
objeto + "." + <TAB>
Si queremos obtener ayuda sobre un mtodo escribimos:
objeto + "." + mtodo + "(" + <TAB>
M.rescale_row(1,1/7) # error debido a que por defecto las matrices
se definen en el conjunto de los enteros

PROBLEMAS AL DEFINIR MATRICES

Los errores que se producen al manipular matrices, en general, vienen determinados por el
hecho de que por defecto las matrices se definen en el conjunto de los nmeros enteros

Tenemos las siguientes posibilidades para cambiar el conjunto de definicin:

1.- Cambiar el anillo/cuerpo de definicicin, mtodo change_ring
2.- Redefinir o definir la matriz indicando directamente el anillo en el que se trabaja.
3.- Especificar elementos del tipo que queremos.
NOTA: Otros errores, evidentemente, se producen al intentar hacer operaciones no
permitidas como sumar matrices de distinto orden o multiplicar matrices incompatibles o
intentar calcular la inversa o el determinante de matrices no cuadradas,... en estos casos
SAGE nos informa del error.
M1=M.change_ring(RDF)
show(M1)
M2=matrix(RDF,3,4,[7,6,9,2,3,9,-4,-3,8,6,7,1])
show(M2)
M3=matrix(3,4,[7.0,6.0,9,2,3,9,-4,-3,8,6,7,1])
show(M3)
#Ver el tipo de matriz
print type(M1)

19/10/2013 17:08

Practica 1 - Introduccion al Algebra de Matrices con SAGE -- Sage

3 de 14

http://sage.math.canterbury.ac.nz/home/albros/15/print

print type(M2)
print type(M3)

OPERACIONES ELEMENTALES
M=matrix(3,4,[7.0,6,9,2,3,9,-4,-3,8,6,7,1])
show(M)
C=M.copy()
C.rescale_row(0,1/7)
show(C)
M=matrix(3,4,[7.0,6,9,2,3,9,-4,-3,8,6,7,1])
show(M)
C=M.copy()
C.add_multiple_of_row(1,0,-3/7)
show(C)
M=matrix(3,4,[7.0,6,9,2,3,9,-4,-3,8,6,7,1])
show(M)
C=M.copy()
C.swap_rows(1,0)
show(C)

DEFINIENDO MATRICES II
Otra forma de definir matrices con condiciones es utilizando el iterador "lambda"
P=matrix(RDF,3,4, lambda i,j: i+j)
P

DEFINIENDO MATRICES III

O mediante una funcin auxiliar.

breve introduccin a las funciones en SAGE (python)
las funciones en SAGE se definen como en "python", comienzan con la palabra clave "def"
como:
def nombre(argumentos):
cuerpo funcion
return salida
es importante tener en cuenta que el cuerpo de la funcin tiene que ir tabulado
#definimos funcin

19/10/2013 17:08

Practica 1 - Introduccion al Algebra de Matrices con SAGE -- Sage

4 de 14

http://sage.math.canterbury.ac.nz/home/albros/15/print

def A(i,j):
return i+j
#definimos matriz
P=matrix(RDF,3,4, lambda i,j: A(i,j))
P

OPERACIONES CON MATRICES

SUMA (A+B): Para matrices de las mismas dimensiones.
DIFERENCIA (A-B)
PRODUCTO POR ESCALARES (k*A)
PRODUCTO (A*B): el producto se define para matrices compatibles, es decir, nmero de
columna de la matriz que premultiplica igual al nmero de filas de la matriz que
postmultiplica.
M=matrix(RDF,3,4,[7,6,9,2,3,9,-4,-3,8,6,7,1])
show(M)
P=matrix(RDF,3,4, lambda i,j: i+j+1)
show(P)
SUMA=P+M
SUMA
PRODUCTO=P*M
M=matrix(RDF,3,4,[7,6,9,2,3,9,-4,-3,8,6,7,1])
show(M)
P=matrix(RDF,2,3, lambda i,j: i+j+1)
show(P)
PRODUCTO=P*M
PRODUCTO
escalar=7
M_por_escalar= escalar * M
html("M = $" + latex(M) + "$")
print "----------------------------"
html("M_por_escalar = $" + latex(M_por_escalar) + "$")
PRODUCTO=M*M.transpose();
PRODUCTO

DETERMINANTE, TRASPUESTA E INVERSA DE MATRICES CUADRADAS

T=M.transpose()
show(M)
show(T)

19/10/2013 17:08

Practica 1 - Introduccion al Algebra de Matrices con SAGE -- Sage

5 de 14

http://sage.math.canterbury.ac.nz/home/albros/15/print

M.inverse()
M.det()
show((M*T).inverse())
show((M*T).det())

MATRIZ DIAGONAL Y MATRIZ IDENTIDAD

diagonal_matrix(RDF,[4,1,2,3,5,7,8])
identity_matrix(RDF,7)

RANGO DE UNA MATRIZ: Forma escalonada por filas (rref), rango (rank) o estudio
de Menores (minors)

Estudiar el rango de la matriz en funcin de los parmetros m y n:

1
1

3
1
n
1

3
5
m

Solucin: para m = 2 y n = 1 el rango es 2, en cualquier otro caso es 3

NOTA: en SAGE el nico parmetro (variable) por defecto es x, si queremos utilizar otras
hay que definirlas antes con "var"
#definimos los parmetros
var('n,m')
#definimos la matriz
C=matrix([[1,3,-3],[1,-1,5],[0,n,m],[m,-1,-4]])
C.rref() #
C.rank() #
C.minors(3)
# el comando solve nos permite calcular los valores para los que se
anula una ecuacin, en este caso cada

19/10/2013 17:08

Practica 1 - Introduccion al Algebra de Matrices con SAGE -- Sage

6 de 14

http://sage.math.canterbury.ac.nz/home/albros/15/print

#uno de los menores

solve(C.minors(3),m,n)
C(n=1).minors(3)
C(m=-2).minors(3)
C(n=1,m=-2).minors(3)
C(n=1,m=-2).minors(2)
C(n=1,m=-2).rref()
C(n=1,m=-2).rank()

CALCULO DE LA INVERSA POR EL METODO DE GAUSS

(ejemplo interactivo)
#Choose the size
D = 3
var('t,s')
example = [[1 if
example[0][-1] =
example[-1][0] =

D of the square matrix:

k==j else 0 for k in range(D)] for j in range(D)]

2
3

@interact
def _(M=input_grid(D,D, default = example,
label='Matriz para invertir', to_value=matrix),
tt = text_control('bits de precisin usados, '
' (solo nmeros de coma flotante)'),
precision = slider(5,100,5,20,'precisin'),
auto_update=False):
if det(M)==0:
print 'Fallo: La Matriz no es Invertible'
return
if M.base_ring() == RR:
M = M.apply_map(RealField(precision))
N=M
M=M.augment(identity_matrix(D),subdivide=true)
print 'Construimos la matriz aumentada'
show(M)
for m in range(0,D-1):
if M[m,m] == 0:
lista = [(M[j,m],j) for j in range(m,D)]
maxi, c = max(lista)
M[c,:],M[m,:]=M[m,:],M[c,:]
print 'Permutamos filas %d y %d'%(m+1,c+1)

19/10/2013 17:08

Practica 1 - Introduccion al Algebra de Matrices con SAGE -- Sage

7 de 14

http://sage.math.canterbury.ac.nz/home/albros/15/print

show(M)
for n in range(m+1,D):
a=M[m,m]
if M[n,m]!=0:
print "Aadimos %s veces la fila %d a la fila
%d"%(-M[n,m]/a, m+1, n+1)
M=M.with_added_multiple_of_row(n,m,-M[n,m]/a)
show(M)
for m in range(D-1,-1,-1):
for n in range(m-1,-1,-1):
a=M[m,m]
if M[n,m]!=0:
print "Aadimos %s veces la fila %d a la fila
%d"%(-M[n,m]/a, m+1, n+1)
M=M.with_added_multiple_of_row(n,m,-M[n,m]/a)
show(M)
for m in range(0,D):
if M[m,m]!=1:
print 'Dividimos la fila %d por %s'%(m+1,M[m,m])
M = M.with_row_set_to_multiple_of_row(m,m,1/M[m,m])
show(M)
M=M.submatrix(0,D,D)
print 'Conservamos la submatriz de la derecha, que contiene la
inversa'
html('$$M^{-1}=%s$$'%latex(M))
print 'Comprobamos que, efectivamente es la inversa'
html('$$M^{-1}*M=%s*%s=%s$$'%(latex(M),latex(N),latex(M*N)))

VECTORES O MATRICES FILA/COLUMNA

Si definimos una matriz fila o columna como un vector, SAGE entiende perfectamente
cuando necesita utilizarlo como fila o como columna.
V=vector(RDF,4)
V
V=vector(RDF,(i*3 for i in [1 .. 4]))
V
V=vector(RDF,range(4))
V
M*V
V*V
V=matrix(RDF,1,4,range(4))
V

19/10/2013 17:08

Practica 1 - Introduccion al Algebra de Matrices con SAGE -- Sage

8 de 14

http://sage.math.canterbury.ac.nz/home/albros/15/print

M*V
M*V.transpose()

NORMAS VECTORIALES
Definicin: Sea B un espacio vectorial real, R n . Se define una norma como una funcin
: B R tal que:
1.
2.
3.

x 0 , y x = 0 x = 0 (definida positiva);
x = ||x para cualquier escalar (real) (homogenea); y
x + y x + y (desigualdad triangular).

Las normas vectoriales ms comunes son:

a.

norma-p definida como x p = (|x i | p )

1/p

, siendo casos particulares

x 1 = |x i |
i

x 2 = (x 2i )

1/2

b.

norma- definida como x = max i |x i |

V=vector(RDF,(i*3 for i in [1 .. 4]))

V
norma_2=sqrt(V[0]**2+V[1]**2+V[2]**2+V[3]**2)
norma_2
norma_2=sqrt(V*V)
norma_2
norma_2=sqrt(sum(V[i]**2 for i in [0 .. 3]))
norma_2
V.norm(2)
norma_inf=max(abs(V[0]),abs(V[1]),abs(V[2]),abs(V[3]))
norma_inf
norma_inf=max(abs(V[i]) for i in range(4))
norma_inf
norma_inf=V.norm(Infinity)
norma_inf

19/10/2013 17:08

Practica 1 - Introduccion al Algebra de Matrices con SAGE -- Sage

9 de 14

http://sage.math.canterbury.ac.nz/home/albros/15/print

NORMAS MATRICIALES
Definicin: es una norma matricial definida en las matrices de orden m x n si es una
norma vectorial en un espacio vectorial de dimensin n m y
1.
2.
3.

A 0 , y A = 0 A = 0 ;
A = ||A para cualquier escalar (real) ; y
A + B A + B .

Por ejemplo:
a.

norma del mximo: A max = max ij|a ij| .

(a 2ij)
i

1/2

b.

norma de Frobenius: A F =

c.

norma- (norma fila) definida como A = max x0

Ax
x

= max i |a ij| ,
j

mximo de la suma de los valores absolutos de los elementos de cada fila.

d. norma-1
(norma
columna)
definida
como
Ax 1
t
A 1 = max x0 x = A = max j |a ij| , mximo de la suma de los
1

valores absolutos de los elementos de cada columna.

M=matrix(RDF,3,4,[7,6,9,2,3,9,-4,-3,8,6,7,1])
show(M)
P=matrix(RDF,2,3, lambda i,j: i+j+1)
show(P)
#determinamos las dimensiones de la matriz (filas,columnas)
print M.dimensions()
print "filas
= " + str(M.dimensions()[0])
print "columnas = " + str(M.dimensions()[1])
#distintas formas de establecer rangos o conjuntos de ndices
print range(M.dimensions()[0])
print [0..M.dimensions()[0]]
print [0..M.dimensions()[0]-1]
norma_1=max(sum(abs(M[i,j]) for i in range(M.dimensions()[0])) for j
in range(M.dimensions()[1]))
norma_1
norma_1=M.norm(1)
norma_1
norma_inf=max(sum(abs(M[i,j]) for j in [0 .. M.dimensions()[1]-1])
for i in [0 .. M.dimensions()[0]-1])

19/10/2013 17:08

Practica 1 - Introduccion al Algebra de Matrices con SAGE -- Sage

10 de 14

http://sage.math.canterbury.ac.nz/home/albros/15/print

norma_inf
norma_inf=M.norm(infinity)
print norma_inf
#otra forma
norma_inf=M.norm(oo)
norma_inf

ERRORES ABSOLUTOS Y RELATIVOS

Cuando trabajamos con ordenadores irremediablemente se producen errores en los calculos,
estos pueden ser de tres tipos; el que se produce por la mquina al almacenar nmeros
finitos (no se puede solventar y aparece siempre), el producido por el mtodo (no siempre
usamos mtodos exactos para resolver), y por ltimo errores producidos por la naturaleza
del problema. Veamos como cuantificarlos:
Dada una norma y dos vectores u y v llamamos distancia a u v .
Si el vector v es una aproximacin del vector u, es decir, v = u + u . Llamaremos error
absoluto ( a ) a la distancia entre u y v, es decir, u v o u , y llamaremos error
uv
u
relativo ( r ) a u o u . El error relativo nos da una mejor visin de la aproximacin y
se suele expresar como r 100% .
Veamos esto en R:
Supongamos, por un lado, x = 35287 y su aproximacin x a = 35287.1 , y por otro lado
y = 0.35287 y su aproximacin y a = 0.352871 entonces:

a (x) = |x x a | = 0.1 y r (x) =

xxa
x

a (y) = |y y a | = 1 10 6 y r (y) =

= 2.8339 10 6

yya
y

= 2.8339 10 6

Y con vectores
V=vector(RDF,4,[1,2,4,5])
show(V)
V_prima=vector(RDF,4,[V[i]+(0.5-random()) for i in [0 .. 3]])
show(V_prima)
deltav=V-V_prima;
deltav
#norma 2
error_2=deltav.norm(2)
print "error absoluto = " + error_2.str()
print "error realtivo = " +
((error_2/V.norm(2))*100).n(digits=4).str() +"%"

19/10/2013 17:08

Practica 1 - Introduccion al Algebra de Matrices con SAGE -- Sage

11 de 14

http://sage.math.canterbury.ac.nz/home/albros/15/print

#norma infinito
error_inf=deltav.norm(oo)
print "error absoluto = " + error_inf.str()
print "error realtivo = " +
((error_inf/V.norm(oo))*100).n(digits=4).str() +"%"
U=vector(RDF,4,[V[i]*1e-6 for i in [0 .. 3]])
U_prima=vector(RDF,4,[V_prima[i] *1e-6 for i in [0 .. 3]])
show(U)
show(U_prima)
deltau=U-U_prima;
show(deltau)
print "Norma 2"
#norma 2
error_2=deltau.norm(2)
print " error absoluto = " + error_2.str()
print " error relativo = " +
((error_2/U.norm(2))*100).n(digits=4).str() +"%"
print "Norma Infinito"
#norma infinito
error_inf=deltau.norm(oo)
print " error absoluto = " + error_inf.str()
print " error realtivo = " +
((error_inf/U.norm(oo))*100).n(digits=4).str() +"%"

DEFINIENDO MATRICES GRANDES

def a(i,j):
if i==j:
return i+1;
elif i>j:
return 0;
else:
return (i+1)*(j+1);
A=matrix(10,10,lambda i,j: a(i,j));A
A=matrix(10,10,lambda i,j: i+1 if i==j else 0 if i>j else (i+1)*
(j+1));A
def a(i,j):
if i==j:
return i+1;
elif abs(i-j)==1:
return 2;
elif abs(i-j)==3:
return 4;
else:

19/10/2013 17:08

Practica 1 - Introduccion al Algebra de Matrices con SAGE -- Sage

12 de 14

http://sage.math.canterbury.ac.nz/home/albros/15/print

return 0;
A=matrix(10,10,lambda i,j: a(i,j))
show(A)
A=matrix(10,10,lambda i,j: i+1 if i==j else 2 if abs(i-j)==1 else 4
if abs(i-j)== 3 else 0)
show(A)
#Matriz de Hilbert
A=matrix(RR,10,10,lambda i,j: 1/(i+j+1));show(A)
A=matrix(QQ,10,10,lambda i,j: 1/(i+j+1));show(A)
A=matrix(QQbar,10,10,lambda i,j: 1/(i+j+1));show(A)

CONDICIONAMIENTO DE UNA MATRIZ

Nmero de condicin:
Se define para una matriz asociada a un sistema, y nos indica en que medida vara su
solucin, en relacin a la variacin de los datos. Si la variacin es grande diremos que el
sistema est mal condicionado. Esto detecta uno de los errores numricos, el dado por el
planteamiento del problema.
El nmero de condicin es relativo a la norma matricial utilizada y se calcula
cond(A) = AA 1 .
Cuanto menor sea el nmero de condicin menor ser la sensibilidad a que cambien los
resultados con una pequea variacin. Si el nmero de condicin es cercano a 1 se dice que
la matriz est bien condicionada.
Ejemplos:
1.- Consideremos el sistema (

perturbado (

1
0.99999

2
x
3
)( )=(
) y el sistema
y
2
3.00001
x
2
3
)( ) =(
) que ocurre al resolverlos?
2
y
3.00001
1
1.00001

M=matrix(RR,2,2,[1,2,1.00001,2])
b=matrix(RR,2,1,[3,3.00001])
sol_1=M.solve_right(b)
sol_1
M=matrix(RR,2,2,[1,2,0.99999,2])
b=matrix(RR,2,1,[3,3.00001])
sol_2=M.solve_right(b)
sol_2

19/10/2013 17:08

Practica 1 - Introduccion al Algebra de Matrices con SAGE -- Sage

13 de 14

http://sage.math.canterbury.ac.nz/home/albros/15/print

print "error absoluto = " + ((sol_1-sol_2).norm(oo)).str()

print "error relativo = " + (((sol_1sol_2).norm(oo)/sol_1.norm(oo))*100).str() + "%"

Calculemos en nmero de condicin de la Matriz

condM=M.norm(oo)*M.inverse().norm(oo)
print "Nmero de Condicin = " + condM.str()

2.-

Consideremos

perturbado (

3
4

el

sistema

3
4

4
x
7
)( )=( )
3
y
1

el

sistema

4
x
7
) ( ) = ( ) que ocurre al resolverlos?
y
2.9
1

M=matrix(RR,2,2,[3,4,4,-3])
b=matrix(RR,2,1,[7,1])
sol_1=M.solve_right(b)
sol_1
M=matrix(RR,2,2,[3,4,4,-2.9])
b=matrix(RR,2,1,[7,1.09757])
sol_2=M.solve_right(b)
sol_2
print "error absoluto = " + ((sol_1-sol_2).norm(oo)).str()
print "error relativo = " + (((sol_1sol_2).norm(oo)/sol_1.norm(oo))*100).str() + "%"

Observamos ahora que la diferencia es pequea, veamos ahora que el nmero de condicin
debe ser prximo a uno
condM=M.norm(oo)*M.inverse().norm(oo)
print "Nmero de Condicin = " + condM.str()

3.- Consideremos un sistema cuya matriz asociada es la de Vandermonde de orden 20 x 20 y

el trmino independiente es [7, -2, 3, 4, 5, 2, 1, 5, 7, 8, 8, 9, 0, 0, 0, 0, 0, 0, 0, 0] t
A=matrix(RDF,20,20, lambda i,j: (i+1)^j)
b=vector(RDF,[7,-2,3,4,5,2,1,5,7,8,8,9,0,0,0,0,0,0,0,0])
solucion=A.solve_right(b)
solucion.n(digits=2)
A*solucion
deltae=solucion-A*solucion

19/10/2013 17:08

Practica 1 - Introduccion al Algebra de Matrices con SAGE -- Sage

14 de 14

http://sage.math.canterbury.ac.nz/home/albros/15/print

print "error absoluto = " + (deltae.norm(oo)).str()

print "error relativo = " +
((deltae.norm(oo)/solucion.norm(oo))*100).str() + "%"
print "numero de condicion = " +
(A.norm(oo)*A.inverse().norm(oo)).str()

19/10/2013 17:08