GIIT_AMA_1718_T3_Tutorial_2

478 days ago by GIIT_AMA_1718

Manejo de matrices con Sage

Comenzamos definiendo una matriz. Por defecto, pondremos RDF en la primera línea de la matriz. Esto nos dice que la matriz es de números reales y por tanto podemos utilizar las operaciones habituales (Sage puede trabajar también con matrices mucho más complejas). La matriz se crea a partir de una lista con sus filas.

A = matrix(RDF,[[4,3,5,6],[6,-1,2,-4],[-3,0,-2,-8]]) show(A) 
       

Análogamente, definimos un vector.

v = vector(RDF,[1,2,3,4]) v 
       

Podemos multiplicar matrices por un escalar y sumar o multiplicar matrices, siempre que las dimensiones cuadren. En particular, podemos multiplicar una matriz por un vector.

show(4*A) 
       
A*v 
       

Nota: No se debe asignar una matriz a otra, es decir operaciones del tipo:

B=A

Se debe hacer usando copy.

B=copy(A) show(B) 
       
show(A+3*B) 
       

Para hacer referencia a un elemento de la matriz, se indica dicho elemento entre corchetes. Sage comienza numerando por 0. Así, $(0,0)$ será el elemento de la esquina superior izquierda de $A$, $(0,3)$ el de la esquina superior derecha, etc.

A[0,0] 
       
A[0,3] 
       

Podemos utilizar las operaciones habituales sobre un elemento. Podemos redefinir un elemento de la matriz, con lo que se modifica la misma.

A[0,0]=A[2,2]+3^2 show(A) 
       

Trabajar con elementos de matrices es un proceso muy lento. Normalmente queremos trabajar con filas o columnas. Para ello se utiliza : para indicar que se toman todos los elementos (de la fila o de la columna). De nuevo recordemos que comenzamos contando por $0$.

A[1,:] 
       
A[:,2] 
       

Podemos modificar una columna entera:

A[:,2]=A[:,0]+A[:,1] show(A) 
       

E intercambiar filas con swap.rows pasando como parámetro los números de filas que queramos intercambiar.

A.swap_rows(0,1) show(A) 
       

Practica: Consideramos la matriz

 

$$A=\begin{pmatrix} 5.0 & 1.0 & 2.0 & 1.0\\ 6.0 & 4.0 & -1.0 & 2.0 \\ -3.0 & 6.0 & -17.0 & 3.0 \end{pmatrix} $$

  1. Crea la matriz A.
  2. Guarda en la segunda fila (A[1,:]) el resultado de restarle a la segunda fila la primera (A[0,:]) multiplicada por 6/5. Prueba evaluar varias veces para ver cómo esta operación crea y elimina un cero en la posición $(1,2)$.
  3. Guarda en la tercera fila, el resultado de restarle a la tercera fila la primera fila multiplicada por el primer elemento de la tercera fila (A[2,0]) dividido entre el primer elemento de la primera fila (A[0,0]). Prueba evaluar varias veces para ver que esta operación mantiene el cero.
  4. Para seguir aplicando eliminación gaussiana con pivote, ahora debemos intercambiar las filas 2 y 3 de la matriz, de modo que el elemento mayor esté en la diagonal. (Nota: si repetimos esta operación se deshace el cambio; podemos evitar eso escribiendo "if abs(A[2,1])>abs(A[1,1]):" antes de la operación.
  5. Modifica la tercera fila para que aparezca un cero debajo de la diagonal.
  6. (Opcional) Crea un programa que aplique el método de Gauss con pivote a matrices de 3x4.
 
       

 

Resolución de sistemas lineales

Eliminación gaussiana con pivote parcial

En el apartado anterior, ya hemos visto cómo aplicar la eliminación gaussiana con pivote parcial para obtener una matriz triangular superior. Ahora sólo tenemos que resolver el sistema. Volvemos a escribir la matriz.

M=matrix(RDF,[[ 5, 1, 2, 1], [ 0, 33/5, -79/5, 18/5], [ 0, 0, 109/33, -8/11]]) show(M) 
       

Para obtener la solución, despejamos "de abajo a arriba". Guardaremos las soluciones en un vector.

x=vector(RDF,3) 
       

Despejamos la tercera incógnita de la tercera línea de la matriz ampliada.

x[2]=M[2,3]/M[2,2];x 
       

La segunda incógnita (x[1]) de la segunda línea (usando que ya conocemos la tercera incógnita).

x[1]=(M[1,3]-M[1,2]*x[2])/M[1,1];x 
       

Y finalmente la primera.

x[0]=(M[0,3]-M[0,1]*x[1]-M[0,2]*x[2])/M[0,0];x 
       

Podemos automatizar la sustitución regresiva con la siguiente función:

%auto def sustitucion_regresiva(U_aumentada): n=U_aumentada.nrows() x=vector(RDF,n) for i in [(n-1),(n-2)..0]: x[i]=( U_aumentada[i,n] - sum([U_aumentada[i,j]*x[j] for j in [i+1..n-1] ]) )/U_aumentada[i,i] return x 
       
sustitucion_regresiva(M) 
       

Automatizando el proceso

En Sage podemos resolver diréctamente el sistema $Ax=b$ el comando solve_right. El método utilizado dependerá del tipo de matriz.

A=M[:,:3] b=vector([M[i,3] for i in [0..2]]) show(A,b) 
       
A.solve_right(b) 
       

Practica:

Resuelve el siguiente sistema mediante eliminación gaussiana con pivote. Compara la solución obtenida con la que da solve_right.

$$ \begin{array}{r} x+z=1  \\ 2x+y=2  \\ -x+y-2z=3 \end{array}$$

 
       

Métodos iterativos

Jacobi

Vamos a aplicar el método de Jacobi para resolver el sistema

$$ \begin{array}{r} 2 x - y = 9  \\ x + 6y - 2z = 15  \\ 4x - 3y  + 8z=1 \end{array}$$

Despejamos la x de la primera ecuación, la y de la segunda y la z de la tercera y obtenemos la función:

J(x,y,z) = ( (9 + y)/2 , (15 - x + 2*z)/6 , (1 - 4*x + 3*y)/8 ) 
       

Si partimos de $x_0=0$, $y_0=0$, $z_0=0$, el vector $(x_1,y_1,z_1)$ dado por el método de Jacobi es:

x0=0.0 y0=0.0 z0=0.0 x1,y1,z1 = J(x0,y0,z0) (x1,y1,z1) 
       

Para obtener el siguiente paso, basta volver a aplicar la función a los valores obtenidos:

x2,y2,z2 = J(x1,y1,z1) (x2,y2,z2) 
       
x3,y3,z3 = J(x2,y2,z2) (x3,y3,z3) 
       

Podemos compararlo con la solución "exacta".

A=matrix(RDF,[[2,-1,0],[1,6,-2],[4,-3,8]]) b=vector(RDF,[9,15,1]) A.solve_right(b) 
       

Practica:

1. Aplica tres pasos del método de Jacobi al sistema. 

$\left.\begin{array}{r} 8x+z=1  \\ 2x-7y=2  \\ -x+y-6 z=3 \end{array}\right\}$

2.Automatización del método de Jacobi. Vamos a crear una función que reciba A, b y un vector x0 y devuelva x1 por el método de Jacobi.

  a) Partimos de la función sustitucion_regresiva. Cambiamos las variables de entrada a A, b y x0.

  b) Definimos como n el número de filas de A.

  c) El bucle for de la i se debería cambiar a [0..n-1] (no es necesario, pero ayuda para Gauss-Seidel). Es decir, recorremos todas las filas de la matriz

  d) Dentro del bucle, debemos aplicar el método, es decir, calcular $(b_i - \sum_{i\neq j} A_{ij} x_{0j})/A_{ii}$. Nótese que eso es sólo cambiar la matriz U_aumentada[i,n] por b[i], U_aumentada[i,j] por A[i,j], x por x0 y recorrer la lista tomando los números tales que i!=j.

3. (Opcional) Aplica tres pasos del método de Gauss-Seidel al sistema anterior (Pista: sustituir x por x0 en la función anterior)

 
       

Resolución de sistemas no lineales

Consideremos el sistema de ecuaciones

$\left.\begin{array}{r} x^3-y^2=0  \\ x^2+2y^2=1   \end{array}\right\}$

Vamos a encontrar una solución con el método de Newton.

En primer lugar, definimos la función $F(x,y)$ de modo que las soluciones al sistema sean un cero de $F$.

F(x,y)=(x^3-y^2,x^2+2*y^2-1);F 
       

Una solución al sistema es una solución a las dos ecuaciones. Podemos dibujar las ecuaciones usando implicit_plot. Donde se corten las dos curvas será la solución buscada.

curvas=implicit_plot(F(x,y)[0]==0,(x,0,1),(y,0,1))+implicit_plot(F(x,y)[1]==0,(x,0,1),(y,0,1),color='red') curvas 
       

Partimos del punto inicial $(1,1)$.

x0,y0=1.0,1.0 
       

Calculamos la matriz jacobiana en el punto $x_0$. Resolvemos el sistema $J x=-F(x_0,y_0)$ que nos da la diferencia entre $x_{n+1}$ y $x_n$. Para obtener $x_{n+1}$, sumamos la solución del sistema anterior a $x_n$.

J0=matrix(RDF,F.diff()(x0,y0)) J0 
       
-F(x0,y0) 
       
J0=matrix(RDF,F.diff()(x0,y0)) xd1,yd1=J0.solve_right(-F(x0,y0)); xd1,yd1 
       
x1,y1=x0+xd1,y0+yd1 x1,y1 
       

Repetimos el proceso para obtener una segunda aproximación.

J1=matrix(RDF,F.diff()(x1,y1)) xd2,yd2=J1.solve_right(-F(x1,y1)); x2,y2=x1+xd2,y1+yd2 x2,y2 
       

Y una tercera

J2=matrix(RDF,F.diff()(x2,y2)) xd3,yd3=J2.solve_right(-F(x2,y2)); x3,y3=x2+xd3,y2+yd3 x3,y3 
       

Pintamos los puntos obtenidos

curvas+points([ (x0,y0), (x1,y1), (x2,y2), (x3,y3) ],color='black') 
       

Practica:

1. Calcular mínimos, máximos y puntos de inflexión de $x^4+y^4-3 x^2 y+2 y^3+3 x y^2-4 x^2+x y-2 x-1$.

var('x,y') F(x,y)=x^4+y^4-3*x^2*y+2*y^3+3*x*y^2-4*x^2+x*y-2*x-1 contour_plot(F(x,y),(x,-3,3),(y,-3,3))+implicit_plot(F.diff()(x,y)[0],(x,-3,3),(y,-3,3),color=(0.3,1,0.3),linestyle='--')+implicit_plot(F.diff()(x,y)[1],(x,-3,3),(y,-3,3),color=(1,0.3,1),linestyle='--') 
       

Número de condición

El número de condición de una matriz nos mide cómo se propagan los errores al resolver el sistema. Para utilizarlo tenemos que ejecutar la siguiente línea.

from numpy.linalg.linalg import cond 
       

Partimos de una matriz "mal condicionada", es decir, con un número de condición alto.

B=Matrix(CDF,[[ 400,-201],[-800,401]]);print cond(B) 
       

Resolvemos el sistema para $b=(100,-200)$.

B.solve_right([100,-200]) 
       

Si perturbamos ligeramente el vector de términos independientes (del orden del 1%) obtenemos que las soluciones cambian notablemente.

B.solve_right([101,-198]) 
       

Si perturbamos ligeramente el vector de términos independientes y la matriz, las soluciones difieren mucho de las del sistema sin perturbar.

B2=Matrix(CDF,[[ 400,-201],[-800,401.999]]);B2.I*Matrix([[101],[-200]]) 
       

Difuminando y perfilando una imagen

Vamos a ver un ejemplo de cómo se puede aplicar el método iterativo de Gauss-Seidel. Concretamente, veremos cómo funciona (de modo básico) un modelo de difuminado y perfilado. Estos modelos se aplican por ejemplo a imágenes de satélite para eliminar los efectos de difusión de la atmósfera, al análisis de imágenes por ordenador, para retocar fotografías, etc.

Cada imagen viene cargada en una matriz (por simplicidad asumiremos que está en blanco y negro). Vamos a necesitar algunas funciones que nos pasen los datos de la matriz a un vector y viceversa. Para ello, basta ejecutar la siguiente línea.

def matrix_to_vector(A): n,m=A.dimensions(); v=vector(RDF,n*m) for i in range(0,n): for j in range(0,m): v[m*i+j]=A[i,j]; return v; def vector_to_matrix(v): n=sqrt(v.degree()); A=matrix(RDF,n); for i in range(0,n): for j in range(0,n): A[i,j]=v[n*i+j]; return A; 
       

Para pasar la matriz a un vector, numeramos las celdas comenzando por la primera fila:

$\left(\begin{array}{rrr}1 & 2 & 3 \\4 & 5 & 6 \\7 & 8 & 9\end{array}\right)\to(1,2,3,4,5,6,7,8,9)$

Definimos la matriz en la que tenemos la imagen. Cuando creamos la matriz estará a cero (color negro). Vamos a poner algunas posiciones a uno (blanco).

n=5; A=matrix(RDF,n) A[1,1:4]=Matrix([1,1,1]) A[3,1:4]=Matrix([1,1,1]) A[2,1]=1 show(matrix_plot(A), figsize=5) 
       

Creamos la matriz de difuminado (está libremente basado en el difuminado gaussiano). Para ello cada posición pasa a tomar "un poco del color de los vecinos". Estamos considerando cuatro vecinos (arriba, abajo, derecha e izquierda). Si A(i,j) es el color del punto (i,j), después de la transformación será: $$A(i,j)\to (1-4d) A(i,j)+d\,A(i-1,j)+d\,A(i+1,j)+d\,A(i,j-1)+d\,A(i,j+1)$$

Esto puede hacerse escribiendo la matriz $A$ en forma de vector y multiplicando el vector por cierta matriz (dada por la expresión anterior). Guardaremos en $N$ la matriz por la que hay que multiplicar (matriz de difuminado).

N=matrix(RDF,n*n); d=1/10; for i in range(0,n-1): for j in range(0,n): N[i+j*n,i+1+j*n]=d; # vecino a la derecha (i+1,j) N[i+1+j*n,i+j*n]=d; # vecino a la izquierda (i-1,j) N[i*n+j,(i+1)*n+j]=d; # vecino superior (i,j+1) N[(i+1)*n+j,i*n+j]=d; # vecino inferior (i,j-1) for i in range(0,n*n): N[i,i]=1-4*d; # posición original (i,j) 
       

Para difuminar la imagen, escribimos la matriz como vector y multiplicamos ese vector por la matriz de difuminado.

Av=matrix_to_vector(A);Av # Necesitamos considerar la imagen como un vector b=N*Av # Multiplicamos la imagen Av por la matriz N show(matrix_plot(vector_to_matrix(b)), figsize=5) # Para dibujar necesitamos volver a pasar a la forma de matriz 
       

Y podemos hacer el proceso inverso resolviendo el sistema $Nx=b$, para lo cual Gauss-Seidel es un método apropiado (a cada paso se perfila un poco más la imagen). Véamoslo en el ejemplo.

Construimos las matrices de Gauss-Seidel

U=matrix(RDF,n*n); for i in range(0,n*n): for j in range(i+1,n*n): U[i,j]=N[i,j]; DL=N-U; 
       

Y aplicamos el método

x0=copy(b) m=3 # Número de iteraciones del método M0=Matrix(RDF,m+1,n*n) # Aquí guardamos los pasos para pintarlo después M0[0,:]=x0; # guardamos la condición inicial for i in range(1,m+1): x0=(DL).solve_right(-U*x0+b) M0[i,:]=x0; # guardamos el vector en la fila i 
       

A continuación mostramos los pasos de Gauss-Seidel. Se puede observar que con 3 iteraciones prácticamente estamos en la matriz original. Esto no depende del número de puntos de la imagen, por lo que para matrices muy grandes este método es mucho más rápido que la eliminación gaussiana o que los métodos LU.

c0=matrix_plot(vector_to_matrix(vector(M0.row(0))), figsize=3) c1=matrix_plot(vector_to_matrix(vector(M0.row(1))), figsize=3) c2=matrix_plot(vector_to_matrix(vector(M0.row(2))), figsize=3) c3=matrix_plot(vector_to_matrix(vector(M0.row(3))), figsize=3) table([["Iteración 0", "Iteración 1","Iteración 2", "Iteración 3"],[c0,c1,c2,c3]], header_row=true)