Untitled

895 days ago by pub

Práctica 1 - Introducción 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 $\mathcal{M}_{n x m}(\mathbb{K})$, donde $\mathbb{K}$ es un conjunto con estructura algebráica de cuerpo ($\mathbb{R}$ o $\mathbb{C}$) o anillo ($\mathbb{Z}$ o $\mathbb{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 también 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

En SAGE todos los elmentos que definimos son "objetos" y como tales, tienen "métodos". 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), descomposición $LU$ (LU), etc...

Tenemos que tener presente que algunos de estos métodos producen cambios en la matriz original (es una forma de reservar espacio en memoria del procesador).

Si queremos saber todos los métodos de un objeto determinado, escribimos:

objeto + "." + <TAB>

Si queremos obtener ayuda sobre un método escribimos:

objeto + "." + método + "(" + <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 números enteros

 

Tenemos las siguientes posibilidades para cambiar el conjunto de definición:

1.- Cambiar el anillo/cuerpo de definicición,  método 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) 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 función auxiliar.

breve introducción 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 función tiene que ir tabulado

#definimos función 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, número de columna de la matriz que premultiplica igual al número 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) 
       
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 función de los parámetros $m$ y $n$:

$$\left(\begin{array}{rrr} 1 & 3 & -3 \\ 1 & -1 & 5 \\ 0 & n & m \\ m & -1 & -4 \end{array}\right)$$

Solución: para $m=-2$ y $n=1$ el rango es 2, en cualquier otro caso es 3

 

NOTA: en SAGE el único parámetro (variable) por defecto es $x$, si queremos utilizar otras hay que definirlas antes con "var"

#definimos los parámetros 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 ecuación, en este caso cada #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 of the square matrix: D = 3 var('t,s') example = [[1 if k==j else 0 for k in range(D)] for j in range(D)] example[0][-1] = 2 example[-1][0] = 3 @interact def _(M=input_grid(D,D, default = example, label='Matriz para invertir', to_value=matrix), tt = text_control('bits de precisión usados, ' ' (solo números de coma flotante)'), precision = slider(5,100,5,20,'precisión'), 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) show(M) for n in range(m+1,D): a=M[m,m] if M[n,m]!=0: print "Añadimos %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 "Añadimos %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))) 
       

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

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 
       
M*V 
       
M*V.transpose() 
       

NORMAS VECTORIALES

Definición: Sea $B$ un espacio vectorial  real, $\mathbb{R}^n$. Se define una norma como una función $\| \cdot \|$: $B \longrightarrow \mathbb{R}$ tal que:

  1. $\|x\|\ge 0$, y $\|x\| = 0 \Leftrightarrow x=0$ (definida positiva);
  2. $\|\alpha x\|=|\alpha| \|x\|$ para cualquier escalar (real) $\alpha$ (homogenea); y
  3. $\|x+y\|\le\|x\|+\|y\|$ (desigualdad triangular).

Las normas vectoriales más comunes son:

  1. norma-p definida como $\|x\|_p=\left (\displaystyle\sum_i |x_i|^p \right )^{1/p}$, siendo casos particulares
    1. $\|x\|_1= \displaystyle\sum_i |x_i| $
    2. $\|x\|_2=\left ( \displaystyle\sum_i x_i^2 \right )^{1/2}$
  2. norma-$\infty$ definida como $\|x\|_\infty=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 
       


NORMAS MATRICIALES

Definición: $\| \cdot \|$ es  una norma matricial definida en las matrices de orden $m$ x $n$ si es una norma vectorial en un espacio vectorial de dimensión $n\cdot m$ y

  1. $\|A\|\ge 0$, y $\|A\| = 0 \Leftrightarrow A=0$;
  2. $\|\alpha A\|=|\alpha| \|A\|$ para cualquier escalar (real) $\alpha$; y
  3. $\|A+B\|\le\|A\|+\|B\|$.

Por ejemplo:

  1. norma del máximo: $\|A\|_{max}=max_{ij}|a_{ij}|$.
  2. norma de Frobenius: $\|A\|_F =\left (\displaystyle\sum_i a_{ij}^2 \right )^{1/2}$
  3. norma-$\infty$ (norma fila) definida como $\|A\|_\infty=max_{x\ne0}\frac{\|Ax\|_\infty}{\|x\|_\infty}=max_i\displaystyle\sum_j|a_{ij}|$, máximo de la suma de los valores absolutos de los elementos de cada fila.
  4. norma-1 (norma columna) definida como $\|A\|_1=max_{x\ne0}\frac{\|Ax\|_1}{\|x\|_1}=\|A^t\|_\infty=max_j\displaystyle\sum_i|a_{ij}|$, máximo de la suma de los 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]) 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 máquina al almacenar números finitos (no se puede solventar y aparece siempre), el producido por el método (no siempre usamos métodos exactos para resolver),  y por último errores producidos por la naturaleza del problema. Veamos como cuantificarlos:

Dada una norma $\|\cdot\|$ y dos vectores $u$ y $v$ llamamos distancia a $\|u-v\|$.

Si el vector $v$ es una aproximación del vector $u$, es decir, $v=u+\delta u$. Llamaremos error absoluto ($\epsilon_a$) a la distancia entre $u$ y $v$, es decir, $\|u-v\|$ o $\|\delta u\|$, y llamaremos error relativo ($\epsilon_r$) a $\frac{\|u-v\|}{\|u\|}$ o $\frac{\|\delta u\|}{\|u\|}$. El error relativo nos da una mejor visión de la aproximación y se suele expresar como $\epsilon_r \cdot 100 \%$.

Veamos esto en $\mathbb{R}$:

Supongamos, por un lado, $x=35287$ y su aproximación $x_a = 35287.1$, y por otro lado $y = 0.35287$ y su aproximación $y_a = 0.352871$ entonces:

$\epsilon_a(x)=|x-x_a|=0.1$ y $\epsilon_r(x)=\frac{|x-x_a|}{|x|}=2.8339\cdot10^{-6}$

$\epsilon_a(y)=|y-y_a|=1\cdot 10^{-6}$ y $\epsilon_r(y)=\frac{|y-y_a|}{|y|}=2.8339\cdot10^{-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() +"%" 
       
#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: 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

Número de condición:

Se define para una matriz asociada a un sistema, y nos indica en que medida varía su solución, en relación a la variación de los datos. Si la variación es grande diremos que el sistema está mal condicionado. Esto detecta uno de los errores numéricos, el dado por el planteamiento del problema.


El número de condición es relativo a la norma matricial utilizada y se calcula $cond(A)=\|A\|\|A^{-1}\|$.


Cuanto menor sea el número de condición menor será la sensibilidad a que cambien los resultados con una pequeña variación. Si el número de condición es cercano a 1 se dice que la matriz está bien condicionada.

Ejemplos:

1.- Consideremos el sistema $\left(\begin{array}{rr} 1 & 2 \\ 1.00001 & 2 \end{array}\right) \cdot \left(\begin{array}{r}x \\ y \end{array}\right)= \left(\begin{array}{rr} 3 \\ 3.00001 \end{array}\right)$ y el sistema perturbado $\left(\begin{array}{rr} 1 & 2 \\ 0.99999 & 2 \end{array}\right) \cdot \left(\begin{array}{r}x \\ y \end{array}\right)= \left(\begin{array}{rr} 3 \\ 3.00001 \end{array}\right)$ ¿que ocurre al resolverlos?

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 
       
print "error absoluto = " + ((sol_1-sol_2).norm(oo)).str() print "error relativo = " + (((sol_1-sol_2).norm(oo)/sol_1.norm(oo))*100).str() + "%" 
       

Calculemos en número de condición de la Matriz

condM=M.norm(oo)*M.inverse().norm(oo) print "Número de Condición = " + condM.str() 
       

2.- Consideremos el sistema $\left(\begin{array}{rr} 3 & 4 \\ 4 & -3 \end{array}\right) \cdot \left(\begin{array}{r}x \\ y \end{array}\right)= \left(\begin{array}{rr} 7 \\ 1 \end{array}\right)$ y el sistema perturbado $\left(\begin{array}{rr} 3 & 4 \\ 4 & -2.9 \end{array}\right) \cdot \left(\begin{array}{r}x \\ y \end{array}\right)= \left(\begin{array}{rr} 7 \\ 1 \end{array}\right)$ ¿que ocurre al resolverlos?

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_1-sol_2).norm(oo)/sol_1.norm(oo))*100).str() + "%" 
       

Observamos ahora que la diferencia es pequeña, veamos ahora que el número de condición debe ser próximo a uno

condM=M.norm(oo)*M.inverse().norm(oo) print "Número de Condición = " + condM.str() 
       

3.- Consideremos un sistema cuya matriz asociada es la de Vandermonde de orden 20 x 20 y el término 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 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()