GDIDP_AMAT_17_Practica2

614 days ago by etlopez17

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

A=matrix(RDF,[[3,1,2],[6,3,2],[-3,0,-8]]);A 
       
[ 3.0  1.0  2.0]
[ 6.0  3.0  2.0]
[-3.0  0.0 -8.0]
[ 3.0  1.0  2.0]
[ 6.0  3.0  2.0]
[-3.0  0.0 -8.0]

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

B=3*A;B 
       
[  9.0   3.0   6.0]
[ 18.0   9.0   6.0]
[ -9.0   0.0 -24.0]
[  9.0   3.0   6.0]
[ 18.0   9.0   6.0]
[ -9.0   0.0 -24.0]
A*B+A^2 
       
[ 36.0  24.0 -32.0]
[120.0  60.0   8.0]
[ 60.0 -12.0 232.0]
[ 36.0  24.0 -32.0]
[120.0  60.0   8.0]
[ 60.0 -12.0 232.0]

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 
       
[-1.0  1.0  2.0]
[ 6.0  3.0  2.0]
[-3.0  0.0 -8.0]
[-1.0  1.0  2.0]
[ 6.0  3.0  2.0]
[-3.0  0.0 -8.0]

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] 
       
-1.0
-1.0
A[0,2] 
       
2.0
2.0
A[2,2] 
       
-8.0
-8.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;A 
       
[ 1.0  1.0  2.0]
[ 6.0  3.0  2.0]
[-3.0  0.0 -8.0]
[ 1.0  1.0  2.0]
[ 6.0  3.0  2.0]
[-3.0  0.0 -8.0]
A[1,0]=A[1,0]/3;A 
       
[ 1.0  1.0  2.0]
[ 2.0  3.0  2.0]
[-3.0  0.0 -8.0]
[ 1.0  1.0  2.0]
[ 2.0  3.0  2.0]
[-3.0  0.0 -8.0]

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,:] 
       
[2.0 3.0 2.0]
[2.0 3.0 2.0]
A.row(1) 
       
(2.0, 3.0, 2.0)
(2.0, 3.0, 2.0)
A[:,2] 
       
[ 2.0]
[ 2.0]
[-8.0]
[ 2.0]
[ 2.0]
[-8.0]
A.column(2) 
       
(2.0, 2.0, -8.0)
(2.0, 2.0, -8.0)

Podemos modificar una columna entera, bien accediendo a las posiciones, bien usando set_column

A[:,2]=[[1],[3],[4]];A 
       
[ 1.0  1.0  1.0]
[ 2.0  3.0  3.0]
[-3.0  0.0  4.0]
[ 1.0  1.0  1.0]
[ 2.0  3.0  3.0]
[-3.0  0.0  4.0]
A.set_column(2,(1,2,3));A 
       
[ 1.0  1.0  1.0]
[ 2.0  3.0  2.0]
[-3.0  0.0  3.0]
[ 1.0  1.0  1.0]
[ 2.0  3.0  2.0]
[-3.0  0.0  3.0]

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

A[[0,2],[0,1]] 
       
[ 1.0  1.0]
[-3.0  0.0]
[ 1.0  1.0]
[-3.0  0.0]

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} $$

 
       

 

Resolución de sistemas lineales

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 
       
[  3.0   1.0   2.0]
[  6.0   4.0   2.0]
[ -3.0   6.0 -17.0]
[  3.0   1.0   2.0]
[  6.0   4.0   2.0]
[ -3.0   6.0 -17.0]
b=vector([1,2,3]);b 
       
(1, 2, 3)
(1, 2, 3)
Aa=A.augment(b);Aa 
       
[  3.0   1.0   2.0   1.0]
[  6.0   4.0   2.0   2.0]
[ -3.0   6.0 -17.0   3.0]
[  3.0   1.0   2.0   1.0]
[  6.0   4.0   2.0   2.0]
[ -3.0   6.0 -17.0   3.0]

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 
       
[  6.0   4.0   2.0   2.0]
[  3.0   1.0   2.0   1.0]
[ -3.0   6.0 -17.0   3.0]
[  6.0   4.0   2.0   2.0]
[  3.0   1.0   2.0   1.0]
[ -3.0   6.0 -17.0   3.0]

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 
       
[  6.0   4.0   2.0   2.0]
[  0.0  -1.0   1.0   0.0]
[ -3.0   6.0 -17.0   3.0]
[  6.0   4.0   2.0   2.0]
[  0.0  -1.0   1.0   0.0]
[ -3.0   6.0 -17.0   3.0]

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 
       
[  6.0   4.0   2.0   2.0]
[  0.0  -1.0   1.0   0.0]
[  0.0   8.0 -16.0   4.0]
[  6.0   4.0   2.0   2.0]
[  0.0  -1.0   1.0   0.0]
[  0.0   8.0 -16.0   4.0]

Tenemos que intercambiar de nuevo la segunda y tercera fila.

Aa.swap_rows(1,2);Aa 
       
[  6.0   4.0   2.0   2.0]
[  0.0   8.0 -16.0   4.0]
[  0.0  -1.0   1.0   0.0]
[  6.0   4.0   2.0   2.0]
[  0.0   8.0 -16.0   4.0]
[  0.0  -1.0   1.0   0.0]

Y hacemos "ceros"

Aa[2,:]=Aa[2,:]-Aa[1,:]*Aa[2,1]/Aa[1,1];Aa 
       
[  6.0   4.0   2.0   2.0]
[  0.0   8.0 -16.0   4.0]
[  0.0   0.0  -1.0   0.5]
[  6.0   4.0   2.0   2.0]
[  0.0   8.0 -16.0   4.0]
[  0.0   0.0  -1.0   0.5]

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

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

Y finalmente la primera.

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

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]; for j in [i+1..n-1]: x[i]=x[i]-U_aumentada[i,j]*x[j] x[i]=x[i]/U_aumentada[i,i] return x 
       
sustitucion_regresiva(Aa) 
       
(0.833333333333, -0.5, -0.5)
(0.833333333333, -0.5, -0.5)

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) 
       
(0.833333333333, -0.5, -0.5)
(0.833333333333, -0.5, -0.5)

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}$$

 
       

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 
       
(
[ 2.0 -1.0  0.0]                  
[ 1.0  6.0 -2.0]                  
[ 4.0 -3.0  8.0], (9.0, 15.0, 1.0)
)
(
[ 2.0 -1.0  0.0]                  
[ 1.0  6.0 -2.0]                  
[ 4.0 -3.0  8.0], (9.0, 15.0, 1.0)
)

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

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

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 
       
(
[2.0 0.0 0.0]  [ 0.0  0.0  0.0]  [ 0.0 -1.0  0.0]
[0.0 6.0 0.0]  [ 1.0  0.0  0.0]  [ 0.0  0.0 -2.0]
[0.0 0.0 8.0], [ 4.0 -3.0  0.0], [ 0.0  0.0  0.0]
)
(
[2.0 0.0 0.0]  [ 0.0  0.0  0.0]  [ 0.0 -1.0  0.0]
[0.0 6.0 0.0]  [ 1.0  0.0  0.0]  [ 0.0  0.0 -2.0]
[0.0 0.0 8.0], [ 4.0 -3.0  0.0], [ 0.0  0.0  0.0]
)

Para hacer sustitución progresiva:

%auto def sustitucion_progresiva(L_aumentada): n=L_aumentada.nrows() x=vector(RDF,n) for i in [0,..,(n-1)]: x[i]=L_aumentada[i,n]; for j in [0..i-1]: x[i]=x[i]-L_aumentada[i,j]*x[j] x[i]=x[i]/L_aumentada[i,i] return x 
       

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

x1=sustitucion_progresiva(D.augment(-(L+U)*x0+b)); x1 
       
(4.5, 2.5, 0.125)
(4.5, 2.5, 0.125)

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

x2=sustitucion_progresiva(D.augment(-(L+U)*x1+b)); x2 
       
(5.75, 1.79166666667, -1.1875)
(5.75, 1.79166666667, -1.1875)

Y una tercera.

x3=sustitucion_progresiva(D.augment(-(L+U)*x2+b)); x3 
       
(5.39583333333, 1.14583333333, -2.078125)
(5.39583333333, 1.14583333333, -2.078125)

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=sustitucion_progresiva(D.augment(-(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=sustitucion_progresiva((D+L).augment(-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: 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\}$

 
       

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, -2.00000000000000)
(0, -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)
J1 
       
[1.6875  -1.25]
[   1.5    2.5]
[1.6875  -1.25]
[   1.5    2.5]
F(x1,y1) 
       
(0.0312500000000000, 0.343750000000000)
(0.0312500000000000, 0.343750000000000)

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