MUI_IIMA_1415_Numerico_2

1069 days ago by evalsanjuan

Resolución de sistemas lineales

Métodos iterativos

1. 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 
       
\newcommand{\Bold}[1]{\mathbf{#1}}\left(\left(\begin{array}{rrr}
2.0 & -1.0 & 0.0 \\
1.0 & 6.0 & -2.0 \\
4.0 & -3.0 & 8.0
\end{array}\right), \left(9.0,\,15.0,\,1.0\right)\right)
\newcommand{\Bold}[1]{\mathbf{#1}}\left(\left(\begin{array}{rrr}
2.0 & -1.0 & 0.0 \\
1.0 & 6.0 & -2.0 \\
4.0 & -3.0 & 8.0
\end{array}\right), \left(9.0,\,15.0,\,1.0\right)\right)

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

x0=vector(RDF,[0,0,0]);x0 
       
\newcommand{\Bold}[1]{\mathbf{#1}}\left(0.0,\,0.0,\,0.0\right)
\newcommand{\Bold}[1]{\mathbf{#1}}\left(0.0,\,0.0,\,0.0\right)

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 
       
\newcommand{\Bold}[1]{\mathbf{#1}}\left(\left(\begin{array}{rrr}
2.0 & 0.0 & 0.0 \\
0.0 & 6.0 & 0.0 \\
0.0 & 0.0 & 8.0
\end{array}\right), \left(\begin{array}{rrr}
0.0 & 0.0 & 0.0 \\
1.0 & 0.0 & 0.0 \\
4.0 & -3.0 & 0.0
\end{array}\right), \left(\begin{array}{rrr}
0.0 & -1.0 & 0.0 \\
0.0 & 0.0 & -2.0 \\
0.0 & 0.0 & 0.0
\end{array}\right)\right)
\newcommand{\Bold}[1]{\mathbf{#1}}\left(\left(\begin{array}{rrr}
2.0 & 0.0 & 0.0 \\
0.0 & 6.0 & 0.0 \\
0.0 & 0.0 & 8.0
\end{array}\right), \left(\begin{array}{rrr}
0.0 & 0.0 & 0.0 \\
1.0 & 0.0 & 0.0 \\
4.0 & -3.0 & 0.0
\end{array}\right), \left(\begin{array}{rrr}
0.0 & -1.0 & 0.0 \\
0.0 & 0.0 & -2.0 \\
0.0 & 0.0 & 0.0
\end{array}\right)\right)

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)) 
       
\newcommand{\Bold}[1]{\mathbf{#1}}\left(4.0,\,2.5,\,0.12\right)
\newcommand{\Bold}[1]{\mathbf{#1}}\left(4.0,\,2.5,\,0.12\right)

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)) 
       
\newcommand{\Bold}[1]{\mathbf{#1}}\left(6.0,\,1.8,\,-1.2\right)
\newcommand{\Bold}[1]{\mathbf{#1}}\left(6.0,\,1.8,\,-1.2\right)

Y una tercera.

x3=D.solve_right(-(L+U)*x2+b); show(x3.n(3)) 
       
\newcommand{\Bold}[1]{\mathbf{#1}}\left(5.0,\,1.2,\,-2.0\right)
\newcommand{\Bold}[1]{\mathbf{#1}}\left(5.0,\,1.2,\,-2.0\right)

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

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(20)) for k in [0..5]]], header = True ) 
       
Paso 0 Paso 1 Paso 2 Paso 3 Paso 4 Paso 5
\left(\begin{array}{r} 0.00000 \\ 0.00000 \\ 0.00000 \end{array}\right) \left(\begin{array}{r} 4.5000 \\ 2.5000 \\ 0.12500 \end{array}\right) \left(\begin{array}{r} 5.7500 \\ 1.7917 \\ -1.1875 \end{array}\right) \left(\begin{array}{r} 5.3958 \\ 1.1458 \\ -2.0781 \end{array}\right) \left(\begin{array}{r} 5.0729 \\ 0.90799 \\ -2.1432 \end{array}\right) \left(\begin{array}{r} 4.9540 \\ 0.94010 \\ -2.0710 \end{array}\right)
Paso 0 Paso 1 Paso 2 Paso 3 Paso 4 Paso 5
\left(\begin{array}{r} 0.00000 \\ 0.00000 \\ 0.00000 \end{array}\right) \left(\begin{array}{r} 4.5000 \\ 2.5000 \\ 0.12500 \end{array}\right) \left(\begin{array}{r} 5.7500 \\ 1.7917 \\ -1.1875 \end{array}\right) \left(\begin{array}{r} 5.3958 \\ 1.1458 \\ -2.0781 \end{array}\right) \left(\begin{array}{r} 5.0729 \\ 0.90799 \\ -2.1432 \end{array}\right) \left(\begin{array}{r} 4.9540 \\ 0.94010 \\ -2.0710 \end{array}\right)

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])) 
       
\newcommand{\Bold}[1]{\mathbf{#1}}\left(-0.0460069444444,\,-0.0598958333333,\,-0.0709635416667\right)
\newcommand{\Bold}[1]{\mathbf{#1}}\left(-0.0460069444444,\,-0.0598958333333,\,-0.0709635416667\right)

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(20)) for k in [0..5]]], header = True ) 
       
Paso 0 Paso 1 Paso 2 Paso 3 Paso 4 Paso 5
\left(\begin{array}{r} 0.00000 \\ 0.00000 \\ 0.00000 \end{array}\right) \left(\begin{array}{r} 4.5000 \\ 1.7500 \\ -1.4688 \end{array}\right) \left(\begin{array}{r} 5.3750 \\ 1.1146 \\ -2.1445 \end{array}\right) \left(\begin{array}{r} 5.0573 \\ 0.94227 \\ -2.0503 \end{array}\right) \left(\begin{array}{r} 4.9711 \\ 0.98805 \\ -1.9901 \end{array}\right) \left(\begin{array}{r} 4.9940 \\ 1.0043 \\ -1.9954 \end{array}\right)
Paso 0 Paso 1 Paso 2 Paso 3 Paso 4 Paso 5
\left(\begin{array}{r} 0.00000 \\ 0.00000 \\ 0.00000 \end{array}\right) \left(\begin{array}{r} 4.5000 \\ 1.7500 \\ -1.4688 \end{array}\right) \left(\begin{array}{r} 5.3750 \\ 1.1146 \\ -2.1445 \end{array}\right) \left(\begin{array}{r} 5.0573 \\ 0.94227 \\ -2.0503 \end{array}\right) \left(\begin{array}{r} 4.9711 \\ 0.98805 \\ -1.9901 \end{array}\right) \left(\begin{array}{r} 4.9940 \\ 1.0043 \\ -1.9954 \end{array}\right)

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])) 
       
\newcommand{\Bold}[1]{\mathbf{#1}}\left(-0.00597692418981,\,0.00431239752122,\,0.00460561116536\right)
\newcommand{\Bold}[1]{\mathbf{#1}}\left(-0.00597692418981,\,0.00431239752122,\,0.00460561116536\right)
 
       
Paso 0 Paso 1 Paso 2 Paso 3 Paso 4 Paso 5
\left(\begin{array}{r} 0.00 \\ 0.00 \\ 0.00 \end{array}\right) \left(\begin{array}{r} 2.2 \\ 1.0 \\ -0.31 \end{array}\right) \left(\begin{array}{r} 3.8 \\ 1.4 \\ -0.75 \end{array}\right) \left(\begin{array}{r} 4.5 \\ 1.5 \\ -1.1 \end{array}\right) \left(\begin{array}{r} 5.0 \\ 1.4 \\ -1.5 \end{array}\right) \left(\begin{array}{r} 5.0 \\ 1.2 \\ -1.6 \end{array}\right)
Paso 0 Paso 1 Paso 2 Paso 3 Paso 4 Paso 5
\left(\begin{array}{r} 0.00 \\ 0.00 \\ 0.00 \end{array}\right) \left(\begin{array}{r} 2.2 \\ 1.0 \\ -0.31 \end{array}\right) \left(\begin{array}{r} 3.8 \\ 1.4 \\ -0.75 \end{array}\right) \left(\begin{array}{r} 4.5 \\ 1.5 \\ -1.1 \end{array}\right) \left(\begin{array}{r} 5.0 \\ 1.4 \\ -1.5 \end{array}\right) \left(\begin{array}{r} 5.0 \\ 1.2 \\ -1.6 \end{array}\right)

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 
       
(x, y) |--> (x^3 - y^2, x^2 + 2*y^2 - 1)
(x, y) |--> (x^3 - y^2, x^2 + 2*y^2 - 1)

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 
       
[ 3.0 -2.0]
[ 2.0  4.0]
[ 3.0 -2.0]
[ 2.0  4.0]
-F(x0,y0) 
       
(0, -2)
(0, -2)
J0=matrix(RDF,F.diff()(x0,y0)) xd1,yd1=J0.solve_right(-F(x0,y0)); xd1,yd1 
       
(-0.25, -0.375)
(-0.25, -0.375)
x1,y1=x0+xd1,y0+yd1 x1,y1 
       
(0.750000000000000, 0.625000000000000)
(0.750000000000000, 0.625000000000000)

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 
       
(0.666666666666667, 0.537500000000000)
(0.666666666666667, 0.537500000000000)

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 
       
(0.657407407407407, 0.532890109101349)
(0.657407407407407, 0.532890109101349)

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) 
       
Iteración 0 Iteración 1 Iteración 2 Iteración 3
Iteración 0 Iteración 1 Iteración 2 Iteración 3
 
       

Ejercicios:

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

 

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

2. 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  \\ 5x+10exp(y)-2z=2  \\ -2x+y-6 sin(z)=3 \end{array}\right\}$$

 
       

Problema. Se tiene la siguiente matriz $M$, resultado de aplicar la matriz $N$ sobre cierta imagen.


M=matrix(RDF,[[0.0, 0.0, 0.0, 0.1, 0.0, 0.0, 0.0],[0.0, 0.0, 0.1, 0.7, 0.1, 0.0, 0.0],[0.0, 0.1, 0.2, 0.8, 0.2, 0.1, 0.0],[0.1, 0.7, 0.8, 1.0, 0.8, 0.7, 0.1],[0.0, 0.1, 0.2, 0.8, 0.2, 0.1, 0.0],[0.0, 0.0, 0.1, 0.7, 0.1, 0.0, 0.0],[0.0, 0.0, 0.0, 0.1, 0.0, 0.0, 0.0]]) n=7 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) 
       

Comparar las imágenes resultantes tras aplicar 3 pasos del método de Jacobi y tres pasos del método de Gauss-Seidel.