Introducción a python para computación numérica III¶
import numpy as np
Si tenemos el sistema lineal $Ax=b$ donde la matriz de coeficientes $A$ y la de términos independientes $b$ son
n = 7
A1 = np.diag(np.ones(n))*3
A2 = np.diag(np.ones(n-1),1)
A = A1 + A2 + A2.T
b = np.arange(n,2*n)*1.
print('A')
print(A)
print('b')
print(b)
Podemos resolver el sistema con el comando genérico
x = np.linalg.solve(A,b)
print('x')
print(x)
Veamos un método específico para resolver esta clase de sistemas: los sistemas lineales tridiagonales
El método de Gauss para sistemas lineales tridiagonales¶
Un sistema tridiagonal tiene por matriz de coeficientes una matriz tridiagonal, que es una matriz donde todos los elementos son cero excepto los de la diagonal principal y las adyacentes por encima y por debajo. Por ejemplo, el sistema lineal siguiente
tiene por matriz de coeficientes y términos independientes
$$ A = \left( \begin{matrix} a_{11} & a_{12} & 0 & 0 & 0 & 0& 0\\ a_{21} & a_{22} & a_{23} & 0 & 0 & 0& 0\\ 0 & a_{32} & a_{33} & a_{34} & 0 & 0& 0\\ 0 & 0 & a_{43} & a_{44} & a_{45} & 0& 0\\ 0 & 0 & 0 & a_{54} & a_{55} & a_{56}& 0\\ 0 & 0 & 0 & 0 & a_{65} & a_{66}& a_{67}\\ 0 & 0 & 0 & 0 & 0 & a_{76}& a_{77} \end{matrix} \right) \quad \quad b = \left( \begin{matrix} b_1\\ b_2\\ b_3\\ b_4\\ b_5\\ b_6\\ b_7 \end{matrix} \right) $$Por Gauss, resolvemos el sistema en dos pasos:
- Triangularización de la matriz aumentada.
- Sustitución regresiva.
Triangularización¶
Tenemos que hacer ceros por debajo de la diagonal principal utilizando operaciones por filas. Para el ejemplo que estamos viendo
Dado un sistema tridiagonal, usando un bucle for, escribir una función At, bt = triangulariza(A,b) que devuelva las matrices At y bt que son las matrices transformadas de A y b tras el proceso de triangularización.
- Usar, para que la impresión de las matrices sea más legible
np.set_printoptions(precision = 2) # solo dos decimales
np.set_printoptions(suppress = True) # no usar notación exponencial
Extraer las dimensiones de
Aconshape. Por ejemplom, n = A.shape. También se puede obtener el número de elementos de un vector conlen. Por ejemplon = len(b).Inicializar la matriz triangularizada
At = np.copy(A)si no queremos modificar la matrizA. Analogamente conb.Triangularizar el sistema $Ax=b$ donde las matrices
Aybse generan con el código
n = 7
A1 = np.diag(np.ones(n))*3
A2 = np.diag(np.ones(n-1),1)
A = A1 + A2 + A2.T
b = np.arange(n,2*n)*1.
- Sin modificar la función, triangularizar también el sistema de matrices
Aybgeneradas con el código
n = 8
np.random.seed(3)
A1 = np.diag(np.random.rand(n))
A2 = np.diag(np.random.rand(n-1),1)
A = A1 + A2 + A2.T
b = np.random.rand(n)
%run Ejercicio1.py
Despejamos las incógnitas de abajo a arriba, empezando por la última
Usando un bucle for, escribir una función x = sust_reg(At,bt) que resuelva un sistema triangular como el anterior por sustitución regresiva. Utilizar como entrada las matrices transformadas obtenidas con la función de triangularización.
Nota:
- Se puede usar un bucle con paso negativo. Por ejemplo donde
for k in range(n1,n2,-1):
n1es el valor inicialn2el valor siguiente al final y-1el paso.
%run Ejercicio2.py
Almacenamiento de los coeficientes en una matriz rectangular de tres columnas¶
Almacenemos ahora los coeficientes en una matriz rectangular. El sistema era
Si guardamos los coeficientes en una matriz rectangular
$$ A_r = \left( \begin{matrix} 0 & a_{11} & a_{12}\\ a_{21} & a_{22} & a_{23}\\ a_{32} & a_{33} & a_{34}\\ a_{43} & a_{44} & a_{45}\\ a_{54} & a_{55} & a_{56}\\ a_{65} & a_{66} & a_{67}\\ a_{76} & a_{77} & 0 \end{matrix} \right)\quad = \quad \left( \begin{matrix} 0 & r_{12} & r_{13}\\ r_{21} & r_{22} & r_{23}\\ r_{31} & r_{32} & r_{33}\\ r_{41} & r_{42} & r_{43}\\ r_{51} & r_{52} & r_{53}\\ r_{61} & r_{62} & r_{63}\\ r_{71} & r_{72} & 0 \end{matrix} \right) $$Ahora, teniendo en cuenta la nueva forma de almacenar los coeficientes, será
Y la matriz triangulada queda
Y ahora la sustitución regresiva será
En el segundo sistema, la matriz de coeficientes es tridiagonal. Este tipo de matrices dispersas (con muchos ceros) son habituales en problemas de la física (ecuación del calor, ecuación de onda) que se resuelven con diferencias finitas.
Por todo esto, se adaptan algoritmos ya existentes a estos casos particulares. Para ello, por ejemplo, en lugar de almacenar la matriz $A$ del segundo sistema como una matriz cuadrada, podemos almacenar sólo las diagonales en una matriz rectangular $A_r$ de la siguiente forma
n = 7
Ar = np.zeros((n,3))
Ar[:,0] = np.concatenate((np.array([0]),np.ones((n-1),)))
Ar[:,1] = np.ones((n),)*3
Ar[:,2] = np.concatenate((np.ones((n-1),),np.array([0])))
b = np.arange(n,2*n)*1.
Modificar el código de las funciones de los dos ejercicios anteriores de forma que admitan la matriz $A$ en la forma $A_r$.
- Sin modificar las funciones, aplicarlas a los datos
n = 8
np.random.seed(3)
Ar = np.zeros((n,3))
Ar[:,1] = np.random.rand(n)
Ar[:,0] = np.concatenate((np.array([0]),np.random.rand(n-1)))
Ar[0:n-1,2] = Ar[1:n,0]
b = np.random.rand(n)
%run Ejercicio3.py
Producto de matrices¶
Utilizaremos tres matrices numpy
- $A$, matriz $m\times n$
- $B$, matriz $n\times p$
- $C$, matriz $m\times p$
Como vamos a multiplicar $A$ por $B$, el número de columnas de $A$, $n$, ha de ser igual al número de filas de $B$. Y la matriz resultante $C$, ha de tener las mismas filas que $A$ y las mismas columnas que $B$.
Un elemento de la matriz producto $C$ viene dado por
$$ c_{ij} = a_{i1}\,b_{1j}+a_{i2}\,b_{2j}+\dots+a_{in}\,b_{nj} = \sum_{k=1}^{n}a_{ik}\,b_{kj}\quad \mathrm{para} \; i = 1,\ldots,m \quad j = 1,\ldots,p \quad (1) $$
Utilizando tres bucles for, uno con el índice i de la fila, otro con el índice j de la columna y otro índice k para la suma acumulada, construir una función C = multiplicar_3bucles(A,B) que multiplique las matrices A y B y devuelva la matriz producto C.
- Extraer las dimensiones de
AyBconshape. Por ejemplom, n = A.shape - Inicializar la matriz
C = np.zeros((m,p))siendomel número de filas deAypel número de columnas deB - Cada elemento de la matriz
Cha sido inicializado a cero y se construye como una suma acumulada. Es decir, empezamos con $c_{ij}=0$. Luego $$k= 1\qquad c_{ij}\leftarrow c_{ij}+a_{i1}b_{1j}$$ $$k= 2\qquad c_{ij}\leftarrow c_{ij}+a_{i2}b_{2j}$$ $$\ldots$$ $$k= n\qquad c_{ij}\leftarrow c_{ij}+a_{in}b_{nj}$$ - Utilizar la función para multiplicar las matrices
A = np.array([[-3.,2],[-2,0],[-4,4],[4,-4]]) B = np.array([[4.,-3,1],[-2,1,1]])
- Sin cambiar la función, utilizarla para multiplicar las matrices
A1 = np.array([[-3.,2],[-2,0],[-4,4],[4,-4],[1,1]]) B1 = np.array([[4.,-3,1,1],[-2,1,1,1]])
%run Ejercicio4.py
Podemos obtener el producto de forma más eficaz quitando el bucle interior y sustituyéndolo por una sola línea de código.
- La fila
ideAesA[i,:], que es un vector de n elementos. - La columna
jdeBesB[:,j], que también es un vector de n elementos. - Los multiplicamos elemento a elemento
A[i,:]*B[:,j]y obtenemos un nuevo vector de n elementos. - Sumamos los elementos de este vector, producto de dos vectores,
np.sum(A[i,:]*B[:,j])
Construir una función multiplicar_2bucles(A,B) que multiplique las matrices A y B y devuelva la matriz producto C usando solo dos bucles. Multiplicar otra vez las matrices del Ejercicio 4.
%run Ejercicio5.py
Podemos aumentar el grado de vectorización obteniendo una columna (o fila) de la matriz C por línea de código. Así quitaríamos otro bucle.
B[:,j]es un vector de n componentes que contiene la columnajdeB.- En
D = A*B[:,j]cada fila deDes el producto elemento a elemento de cada fila deAcon el vectorB[:,j]. Por lo tanto,Dtiene el mismo tamaño queA. C[:,j] = np.sum(D,axis=1)es una matriz column donde cada elemento es la suma de los elementos de una fila deDy el el producto matricial de la matrizApor la matriz columnaB[:,j].
A = np.array([[-3.,2],[-2,0],[-4,4],[4,-4]])
B = np.array([[4.,-3,1],[-2,1,1]])
m = A.shape[0]
p = B.shape[1]
C = np.zeros((m,p))
for j in range(p):
C[:,j] = np.sum(A*B[:,j],axis=1)
print(C)
%run Ejercicio6.py
Comparemos el tiempo de ejecución de las tres funciones de multiplicar matrices, el comando numpy np.dot(A,B) y A @ B. Usar las matrices $A_{m\times n}$ y $B_{n\times p}$ con elementos obtenidos alatoriamente. Usar el comando time.time() visto en la práctica anterior para calcular el tiempo de ejecución. np.random.rand(m,n) crea una matriz de elementos aleatorios distribuidos uniformemente en $[0,1]$.
m = 300; n = 200; p = 400
A = np.random.rand(m,n)
B = np.random.rand(n,p)
%run Ejercicio7.py