GT_MM_1314_Sistemas_Lineales

1425 days ago by trinidad

Resolución de sistemas lineales

0. Manejo de matrices con Sage

En el tutorial vimos que Sage puede trabajar con matrices y vectores. Vamos a profundizar un poco sobre las operaciones que podemos utilizar.

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).

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

Podemos multiplicar matrices por un escalar y sumar o multiplicar matrices, siempre que las dimensiones cuadren.

B=3*A;B 
       
A*B+A^2 
       

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

B=A

En ese caso, Sage no crea una matriz B que contenga los mismos datos que A, sino que B pasa a ser otro nombre para la matriz A (Sage es orientado a objetos)

B=A B[0,0]=-1 A 
       

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,2)$ el de la esquina superior derecha, etc.

A[0,0] 
       
A[0,2] 
       
A[2,2] 
       

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;A 
       
A[1,0]=A[1,0]/3;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$. También se pueden utilizar las funciones row y column, aunque a diferencia de la orden : no devuelven una matriz, sino un vector.

A[1,:] 
       
A.row(1) 
       
A[:,2] 
       
A.column(2) 
       
A[:,2]=[[1],[3],[4]];A 
       

Por último, podemos quedarnos con un menor de la matriz indicando las filas y columnas que queremos:

A[[0,2],[0,1]] 
       

Practica: Crea la matriz

 

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

  1. Modifica la matriz para que en la posición superior izquierda haya un 3.
  2. Guarda en una variable la segunda fila.
  3. Obtén el resultado de sumar la segunda y tercera columnas
  4. Modifica la tercera fila para que pase a valer la suma de la primera y la tercera fila de la matriz, es decir:

 

$$ A=\begin{pmatrix} 3.0 & 1.0 & 2.0 \\ 6.0 & 4.0 & 2.0 \\ 0.0 & 7.0 & -15.0 \end{pmatrix} $$

 
       

1. Eliminación gaussiana con pivote parcial

Vamos a aplicar la eliminación gaussiana con pivote parcial a la siguiente matriz:

A=matrix(RDF,[[3,1,2],[6,4,2],[-3,6,-17]]);A 
       
b=vector([1,2,3]);b 
       
Aa=A.augment(b);Aa 
       

Intercambiamos las dos primeras filas para que el mayor valor esté en la posición (1,1) de la matriz.

Hay que tener en cuenta que Sage comienza numerando por (0,0).

Para intercambiar filas utilizaremos swap_rows, que recibe las posiciones de dos filas y las intercambia.

Aa.swap_rows(0,1);Aa 
       

Ahora "hacemos cero" debajo del elemento de la diagonal.

La segunda fila pasa a ser la segunda fila menos la primera dividida entre dos.

Aa[1,:]=Aa[1,:]-Aa[0,:]*(1/2);Aa 
       

Hacemos lo mismo con la tercera fila. Ahora el multiplicador será $-1/2$. Nótese que se puede hacer que lo calcule automaticamente:

Aa[2,:]=Aa[2,:]-Aa[0,:]*Aa[2,0]/Aa[0,0];Aa 
       

Tenemos que intercambiar de nuevo la segunda y tercera fila.

Aa.swap_rows(1,2);Aa 
       

Y hacemos "ceros"

Aa[2,:]=Aa[2,:]-Aa[1,:]*Aa[2,1]/Aa[1,1];Aa 
       

Sustitución regresiva

Para obtener la solución, despejamos "de abajo a arriba".

x=vector(RDF,[0,0,0]) #Creamos el vector donde guardaremos las soluciones 
       

Despejamos la tercera incógnita (x[2]) de la tercera línea de la matriz ampliada.

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

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

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

Y finalmente la primera.

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

Automatizando el proceso

En Sage podemos resolver diréctamente el sistema $Ax=b$ mediante eliminación gaussiana con pivoteo utilizando el comando solve_right.

A.solve_right(b) 
       

Practica:

Resuelve el siguiente sistema mediante eliminación gaussiana con pivote. Hazlo primero "a mano" y después utilizando el método solve_right.

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

 
       

2. Factorización LU

Vamos a utilizar la factorización LU que tiene Sage. Esta vez no vamos a definir la matriz del tipo RDF porque en dicho caso sólo permite la factorización LU con pivote y queremos obtener también la factorización sin pivote.

A=matrix([[ 6, 6, 6 ],[2, 8, 11 ],[1, 13, 7]]);A 
       

Para factorizar sin pivote, utilizamos el método LU, con el argumento "pivot='nonzero'".

P, L, U = A.LU(pivot='nonzero');P,L,U 
       

Y si queremos utilizar el método con pivote, no ponemos ningún argumento o ponemos "pivot='partial'".

P, L, U = A.LU(pivot='partial');P,L,U 
       

Resolución de sistemas mediante la factorización LU

Una vez tenemos factorizada la matriz $A$, vamos a resolver el sistema $Ax=b$.

A=matrix(RDF,[[ 6, 6, 6 ],[2, 8, 11 ],[1, 13, 7]]);A b=vector([1,2,3]) P, L, U = A.LU() A,b 
       

En primer lugar, el sistema es no singular pues su determinante es no nulo. Lo podemos calcular -en valor absoluto- como el producto de los determinantes de  $U$ y $L$. Para obtener el signo exacto, tendríamos que calcular el determinante de $P$ que es $+1$ o $-1$.

abs(det(A)),prod(L.diagonal())*prod(U.diagonal()) 
       

Resolvemos $Ly=Pb$ (lo que se debería hacer por sustitución regresiva)

y=var('y') y=L.solve_right(P*b);y 
       

y para terminar, resolvemos $Ux=y$ (de nuevo por sustitución regresiva).

x=U.solve_right(y);x 
       

Comprobamos la solución:

A*x 
       

Cholesky

Cuando la matriz es simétrica y definida positiva, la factorización LU más rápida es la de Cholesky. Es decir, obtener una matriz triangular inferior $L$ tal que $A=L L^\top$.

Tomamos una matriz de ejemplo:

A = Matrix(CDF,[[4,4,6],[4,5,7],[6,7,74]]);A 
       

Sage nos permite calcular diréctamente su descomposición de Cholesky.

L=A.cholesky_decomposition();L 
       

Comprobemos que es correcta:

A,L*L.transpose() 
       

3. Métodos iterativos

Jacobi

Vamos a aplicar el método de Jacobi para resolver el sistema $Ax=b$ dado por la matriz y el vector siguientes:

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

Partimos de una solución inicial aproximada. Usualmente se toma el vector con todas sus posiciones ceros.

x0=vector(RDF,[0,0,0]);x0 
       

Construimos las matrices que necesitamos para aplicar el método.

D=matrix(RDF,[[2,0,0],[0,6,0],[0,0,8]]) L=matrix(RDF,[[0,0,0],[1,0,0],[4,-3,0]]) U=matrix(RDF,[[0,-1,0],[0,0,-2],[0,0,0]]) D,L,U 
       

Para obtener una aproximación de la solución, debemos resolver $D x=-(L+U)x0+b$.

x1=D.solve_right(-(L+U)*x0+b); show(x1.n(3)) 
       

Repitendo el proceso, ahora con $x1$, obtenemos una nueva aproximación de la solución.

x2=D.solve_right(-(L+U)*x1+b); show(x2.n(3)) 
       

Y una tercera.

x3=D.solve_right(-(L+U)*x2+b); show(x3.n(3)) 
       

Automatizando el método

Podemos automatizar el proceso, usando un blucle for.

Volvemos a definir las matrices:

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

Y utilizamos un bucle "for" para repetir el proceso. Usamos una matriz M0 para ir guardando los vectores obtenidos y mostrarlos luego en forma de tabla (no es necesario para utilizar el método). Nótese que se ha redondeado a sólo dos cifras decimales.

x0=vector(RDF,[0,0,0]) M0=matrix(RDF,[[0],[0],[0]]) # Para guardar las soluciones for i in range(1,6): x0=D.solve_right(-(L+U)*x0+b) M0=M0.augment(x0) # Guardamos x0 html.table( [ [("Paso {pi}".format(pi=k)) for k in [0..5]],[(M0[:,k].n(14)) for k in [0..5]]], header = True ) 
       

La verdadera solución del sistema es (5,1,-2), a la cual se aproxima la sucesión que define el método. Los errores que hemos cometido en este caso son:

show(x0-vector(RDF,[5,1,-2])) 
       

Gauss-Seidel

El método de Gauss-Seidel converge más rapidamente:

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

Utilizamos un bucle "for" para automatizar el procesor. 

x0=vector(RDF,[0,0,0]) M0=matrix(RDF,[[0],[0],[0]]) # Aquí vamos guardando los pasos para mostrarlos al final. for i in range(1,6): x0=(D+L).solve_right(-U*x0+b) M0=M0.augment(x0) # Guardamos x0 html.table( [ [("Paso {pi}".format(pi=k)) for k in [0..5]],[(M0[:,k].n(14)) for k in [0..5]]], header = True ) 
       

Y con el método de Gauss-Seidel los errores son menores, para el mismo número de iteraciones.

show(x0-vector(RDF,[5,1,-2])) 
       

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

 

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

 
       

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 
       

Finalmente, dibujamos los puntos obtenidos para comprobar que convergen a la solución buscada.

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

Practica:

Aproximar un cero del sistema

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

mediante tres pasos del método de Newton, tomando como valor inicial $(1,1)$.

 
       

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=vector(RDF,n*n) 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) html.table([["Iteración 0", "Iteración 1","Iteración 2", "Iteración 3"],[c0,c1,c2,c3]], header=True) 
       

Ejercicios:

Crea una variable dni=0."número de DNI". Si sois un grupo, usad la media de todos los DNI.

1.Resolver mediante el método de Gauss el sistema. Comprobar la solución obtenida.

$$\left.\begin{array}{r} 3x+z=1  \\ 2x+(4+dni) y=2  \\ -x+y-6 z=3 \end{array}\right\}$$

2. Resolver mediante LU el sistema:

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

3. Aplicar tres pasos del método de Gauss-Seidel al sistema:

$$\left.\begin{array}{r} (12+dni)x-y+z=1  \\ 5x+4y-2z=2  \\ -2x+y-6 z=3 \end{array}\right\}$$

Discutir si podemos asegurar la convergencia.

4. Aplicar tres pasos del método de Newton para, partiendo de $(0,0,0)$, aproximar un cero del sistema:

$$\left.\begin{array}{r} (12+dni)x-y^3-exp(z)=1  \\ -sin(x)+10 y-2z=2\\ -2x+-\cos(y)+16\sin(z)=3  \end{array}\right\}$$

Si quieres representar las ecuaciones y los puntos, deber usar las versiones 3d de parametric_plot y points.

 
       

Problema. Se mide la altura de un objeto en cada instante de tiempo desde el minuto $0$ al minuto $19$, obteniendo los siguientes datos:

altura=[4.75351711919536, 7.56752950275510, 10.9960497733263, 14.0971301046088, 17.4130376577309, 20.9649523428482, 24.3568624902812, 27.3541737270660, 30.8184613374497, 33.8515985570556, 37.8654187115242, 41.3977351632576, 44.1241410048809, 47.7356466165504, 52.0740780452740, 55.3145040396807, 58.4873686846329, 62.2134285713830, 66.0439512440366, 70.6808944275451] 
       

Vamos a obtener una recta $y=at+b$ tal que la suma de los errores en las alturas sea mínima.

1. Dibujar en una gráfica los puntos $(t,altura(t))$, usando el comando points.

2. Definir la función "suma de las distancias vertical entre la recta y la curva" (también valdría la diferencia de alturas al cuadrado).

3. Minimizar la función del apartado anterior utilizando el Método de Newton (hasta que no varíen los 4 primeros decimales) sobre la derivada (será una función de $a$ y $b$).

4. Dibujar la recta obtenida y los puntos.