GIIT_AMA_18_Numérico_Tutorial_2

111 days ago by etlopez18

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 
       
(1.0, 2.0, 3.0, 4.0)
(1.0, 2.0, 3.0, 4.0)

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 
       
(49.0, -6.0, -41.0)
(49.0, -6.0, -41.0)

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] 
       
4.0
4.0
A[0,3] 
       
6.0
6.0

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,:] 
       
[ 6.0 -1.0  2.0 -4.0]
[ 6.0 -1.0  2.0 -4.0]
A[:,2] 
       
[ 5.0]
[ 2.0]
[-2.0]
[ 5.0]
[ 2.0]
[-2.0]

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 $(2,1)$ (la (1,0) en Sage).
  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 
       
(0.0, 0.0, -0.2201834862385321)
(0.0, 0.0, -0.2201834862385321)

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 
       
(0.0, 0.018348623853211062, -0.2201834862385321)
(0.0, 0.018348623853211062, -0.2201834862385321)

Y finalmente la primera.

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

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) 
       
(0.2844036697247706, 0.018348623853211062, -0.2201834862385321)
(0.2844036697247706, 0.018348623853211062, -0.2201834862385321)

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) 
       
(0.2844036697247706, 0.018348623853211062, -0.2201834862385321)
(0.2844036697247706, 0.018348623853211062, -0.2201834862385321)

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) 
       
(4.50000000000000, 2.50000000000000, 0.125000000000000)
(4.50000000000000, 2.50000000000000, 0.125000000000000)

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) 
       
(5.75000000000000, 1.79166666666667, -1.18750000000000)
(5.75000000000000, 1.79166666666667, -1.18750000000000)
x3,y3,z3 = J(x2,y2,z2) (x3,y3,z3) 
       
(5.39583333333333, 1.14583333333333, -2.07812500000000)
(5.39583333333333, 1.14583333333333, -2.07812500000000)

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(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]]) 
       

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.0],[0.0]]) # Para guardar las soluciones for i in [1..5]: 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]) for k in [0..5]]], header = True ) 
       
Paso 0 Paso 1 Paso 2 Paso 3 Paso 4 Paso 5
Paso 0 Paso 1 Paso 2 Paso 3 Paso 4 Paso 5

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(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]]) 
       

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 [1..5]: x0=(D+L).solve_right(-U*x0+b) M0=M0.augment(x0) # Guardamos las soluciones html.table( [ [("Paso {pi}".format(pi=k)) for k in [0..5]],[(M0[:,k]) for k in [0..5]]], header = True ) 
       
Paso 0 Paso 1 Paso 2 Paso 3 Paso 4 Paso 5
Paso 0 Paso 1 Paso 2 Paso 3 Paso 4 Paso 5

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:

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 
       
(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.000000000000000, -2.00000000000000)
(-0.000000000000000, -2.00000000000000)
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)

Pintamos los puntos obtenidos

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