Resolución de sistemas lineales: métodos directos¶
import numpy as np
Precisamos el formato en que queremos imprimir
np.set_printoptions(precision = 2) # solo dos decimales
np.set_printoptions(suppress = True) # no usar notación exponencial
Vamos a resolver sistemas de ecuaciones lineales con el mismo número de ecuaciones que de incógnitas con matriz de coeficientes no singular que, por lo tanto, tiene solución única.
Dados los números $a_{ij}$ y $b_{j}$ para $i,j=1,2,\ldots,n$ se trata de hallar los números $x_{1},x_{2},\ldots,x_{n}$ que verifican las $n$ ecuaciones lineales
$$ \begin{array}{cll} a_{11}x_{1}+a_{12}x_{2}+\ldots+a_{1n}x_{n} & = & b_{1},\\ a_{21}x_{1}+a_{22}x_{2}+\ldots+a_{2n}x_{n} & = & b_{2},\\ \vdots\\ a_{n1}x_{1}+a_{n2}x_{2}+\ldots+a_{nn}x_{n} & = & b_{n}. \end{array} $$Sustitución regresiva¶
Análogamente, para el sistema
$$ \mathbf{U}\mathbf{x}=\mathbf{b} $$con $\mathbf{U}$ de la forma
$$ \mathbf{U}=\left(\begin{array}{ccccc} u_{11} & u_{12} & u_{13} & \ldots & u_{1n}\\ 0 & u_{22} & u_{23} & \ldots & u_{2n}\\ 0 & 0 & u_{33} & \ldots & u_{3n}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & 0 & u_{nn} \end{array}\right) $$la $i$-ésima ecuación es, para $i=1,2,\ldots,n$,
$$ u_{ii}x_{i}+u_{ii+1}x_{i+1}+\cdots+u_{in}x_{n}=b_{i} $$Por lo tanto, para $i=1,2,\ldots,n$
$$ x_{i}=\frac{b_{i}-u_{ii+1}x_{i+1}-\cdots-u_{in}x_{n}}{u_{ii}}= \frac{1}{u_{ii}}\left(b_{i}-\sum_{j=i+1}^{n}u_{ij}x_{j}\right), $$Algoritmo¶
- $x_{n}=\dfrac{b_{n}}{u_{nn}}$
- Para $i=n-1,\ldots,1$ $$ x_{i}=\frac{1}{u_{ii}}\left(b_{i}-\sum_{j=i+1}^{n}u_{ij}x_{j}\right). $$
Escribir una función sust_regre(U,b) cuyo argumento de salida $x$ es la solución del sistema $Ux=b$, con U una matriz triangular superior y b una matriz columna.
Usarla para resolver el sistema
$$\left(\begin{array}{rrr} 2 & 1 & 1\\ 0 & 2 & 1\\ 0 & 0 & 2 \end{array}\right) \left(\begin{array}{c} x_1\\ x_2\\ x_3 \end{array}\right)= \left(\begin{array}{r} 9\\ 4\\ 4 \end{array}\right) $$U = np.array([[2., 1, 1], [0, 2, 1], [0, 0, 2]])
b = np.array([9., 4, 4])
- Comprobar la solución calculando $Ux$ y comparándolo con $b$.
- Comprobar que la función también funciona bien con el sistema $Ux=b$ de matrices $L$ y $b$ generadas como sigue
n = 5 np.random.seed(2) U = np.random.random((n,n)) U = np.triu(U) # Haz cero los elementos bajo la diagonal b = np.random.random(n)
Nota
- Las dimensiones de una matriz se pueden obtener con
n_filas,n_columnas = A.shapeque devuelve una tupla con el número de filas y columnas. - El número de elementos de un vector se pueden obtener con
n = len(b) - Si queremos que
xtenga las mismas dimensiones queb, podemos inicializarlox = np.zeros(len(b))ox = np.zeros_like(b)
%run Ejercicio1.py
Gauss¶
La idea del método de Gauss es transformar el sistema original en un sistema con matriz de coeficientes triangular superior que tenga la misma solución. Para ello, usaremos las siguientes operaciones:
- Sumamos una fila multiplicada por una constante a otra $f_i\rightarrow f_i +\lambda f_j,\quad i\ne j$
- Intercambiamos filas.
Una vez que la matriz del sistema es triangular, podemos resolver el sistema por sustitución regresiva.
En general, el método debe incluir el intercambio de filas, pero empezaremos escribiendo el código sin intercambio de filas.
Algoritmo¶
- Para $k=1,\ldots,n-1$
- Si $a_{kk}=0$, acabar.
- Caso contrario, para $i=k+1,\ldots,n,$
- $f = \dfrac{a_{ik}}{a_{kk}}$
- $b_{i}=b_{i}-f\,b_{k}$
- $a\,\mathrm{fila}_i \leftarrow a\,\mathrm{fila}_i - f\times a\,\mathrm{fila}_k$
Escribir una función triang(A,b) cuyos argumentos de salida son la matriz triangular superior U y el vector c, obtenidos a partir de la matriz del sistema A y su término independiente b.
Usarla con las matrices del sistem $Ax=b$
$$\left(\begin{array}{rrr} 2 & 1 & 1\\ 1 & 2 & 1\\ 1 & 1 & 2 \end{array}\right) \left(\begin{array}{c} x_1\\ x_2\\ x_3 \end{array}\right)= \left(\begin{array}{r} 2\\ 4\\ 6 \end{array}\right) $$- Escribir los argumentos de salida
Uyc. - Resolver el sistema equivalente $Ux=c$ y comprobar que la solución del sistema es correcta, calculando $Ax$ y comparándolo con $b$.
- Comprobar que el código también funciona correctamente para resolver el sistema $Ax=b$ cuyas matrices $A$ y $b$ han sido generadas de la manera siguiente
n = 5 np.random.seed(3) A = np.random.random((n,n)) b = np.random.random(n)
Nota
- Podemos extraer y operar sobre la fila
ide una matrizAconA[i,:]o simplementeA[i].
%run Ejercicio2.py
Matriz Inversa por Gauss-Jordan¶
Para calcular la inversa de una matriz $A_{n\times n}$ por Gauss-Jordan, seguiremos los pasos siguientes:
Escribimos una matriz $n\times 2n$ que consiste en una matriz $M$ formada por la concatenación de las matrices $A$ a la izquierda y la matriz identidad $I$ de dimensión $n\times n$ a la derecha, es decir $M = [A\,|\,I]$
Mediante operaciones por filas, transformamos la matriz $A$ en una matriz $I$. Entonces la matriz de la derecha se transformará en la matriz inversa de $A$. Es decir, a partir de $M = [A\,|\,I]$ llegaremos a $M = [I\,|\,A^{-1}]$. Las operaciones por filas serán:
Para $k = 1,\ldots,n$
- Si $m_{kk}=0$, acabar.
- Caso contrario,
- $m\,\mathrm{fila}_k \leftarrow \dfrac{m\,\mathrm{fila}_k}{m_{kk}}$
- Para $i = 1,\ldots,n \quad i \ne k$
$m\,\mathrm{fila}_i \leftarrow m\,\mathrm{fila}_i - m_{ik}\times m\,\mathrm{fila}_k$
Escribir una función gaussjordan(A) cuyo argumento de entrada es la matriz cuadrada A y cuyo argumento de salida es la matriz inversa de A y que calcule la inversa utilizando el método de Gauss-Jordan sin pivote.
Calcular la inversa de las matrices $A$ del ejercicio anterior.
Notas
- La matriz identidad $n \times n$ se puede obtener con
I = np.eye(n). - Dos matrices
AyBse pueden concatenar horizontalmente (si tienen el mismo número de filas) conM = np.concatenate((A,B),axis=1). - Podemos extraer y operar sobre la fila
ide una matrizMconM[i,:]o simplemente,M[i]. - Podemos extraer las últimas columnas de una matriz desde la columna
nconM[:,n:].
%run Ejercicio3.py
Sustitución progresiva¶
Para el sistema
$$ \mathbf{L}\mathbf{x}=\mathbf{b}, $$con $\mathbf{L}$ de la forma
$$ \mathbf{L}=\left(\begin{array}{ccccc} l_{11} & 0 & 0 & \ldots & 0\\ l_{21} & l_{22} & 0 & \ldots & 0\\ l_{31} & l_{32} & l_{33} & \ldots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ l_{n1} & l_{n2} & l_{n3} & \ldots & l_{nn} \end{array}\right) $$la ecuación $i$-ésima es, para $i=1,2,\ldots,n$,
$$ l_{i1}x_{1}+l_{i2}x_{2}+\cdots+l_{ii-1}x_{i-1}+l_{ii}x_{i}=b_{i}. $$Por lo tanto, para $i=1,2,\ldots,n$
$$ x_{i}=\frac{b_{i}-l_{i1}x_{1}-\cdots-l_{ii-1}x_{i-1}}{l_{ii}}= \frac{1}{l_{ii}}\left(b_{i}-\sum_{j=1}^{i-1}l_{ij}x_{j}\right). $$Algoritmo¶
- $x_{1}=\dfrac{b_{1}}{l_{11}}$
- Para $i=2,\ldots,n$ $$ x_{i}=\frac{1}{l_{ii}}\left(b_{i}-\sum_{j=1}^{i-1}l_{ij}x_{j}\right). $$
Escribir una función sust_prog(L,b) cuyo argumento de salida x es la solución del sistema $Lx=b$, con L una matriz triangular inferior y b una matriz columna.
Usarla para resolver el sistema
$$\left(\begin{array}{rrr} 2 & 0 & 0\\ 1 & 2 & 0\\ 1 & 1 & 2 \end{array}\right) \left(\begin{array}{c} x_1\\ x_2\\ x_3 \end{array}\right)= \left(\begin{array}{r} 2\\ -1\\ 0 \end{array}\right) $$- Comprobar la solución calculando $Lx$ y comparándolo con $b$.
- Comprobar que la función también funciona bien con el sistema $Lx=b$ de matrices $L$ y $b$ generadas como sigue
n = 5 np.random.seed(1) L = np.random.random((n,n)) L = np.tril(L) # hace cero los elementos por encima de la diagonal b = np.random.random(n)
%run Ejercicio4.py
Factorización LU¶
Expresaremos $\mathbf{A}$ como el producto de dos matrices: una matriz triangular inferior, $\mathbf{L}$, y una matriz triangular superior, $\mathbf{U}$,
$$ \mathbf{A}=\mathbf{L}\mathbf{U} $$Si tenemos el sistema
$$ \mathbf{A}\mathbf{x}=\mathbf{b}, $$donde $\mathbf{A}$ admite la factorización $\mathbf{L}\mathbf{u}$
$$ \mathbf{A}\mathbf{x}=\mathbf{b}\Longleftrightarrow \mathbf{L}\mathbf{U}\mathbf{x}=\mathbf{b} $$Una vez factorizada la matriz $\mathbf{A}$ solucionaremos el sistem $\mathbf{A}\mathbf{x}=\mathbf{b}$ en dos pasos
- Resolver $\mathbf{L}\mathbf{y}=\mathbf{b}$, con sustitución progresiva para obtener $\mathbf{y}$.
- Resolver $\mathbf{U}\mathbf{x}=\mathbf{y}$, con sustitución regresiva para obtener $\mathbf{x}$.
De momento, no estamos considerando el intercambio de filas.
Usar la factorización $LU$ de $A$ para resolver el sistem $Ax=b$, donde
$$ A=\left(\begin{array}{rrr} 2 & 1 & 1\\ 1 & 2 & 1\\ 1 & 1 & 2 \end{array}\right)\quad\quad b=\left(\begin{array}{c} 2\\ 4\\ 6 \end{array}\right) $$- Crear la función
LU(A)que da la factorización $LU$ de $A$. - Obtener y escribir las matrices $L$ y $U$.
- Usar las funciones
sust_regreysust_progpara resolver los sistemas resultantes. - Comprobar la solución calculando $Ax$ y comparándolo con $b$.
- Las matrices se pueden inicializar:
U = np.zeros((n,n))(matriz de ceros $n\times n$, como $A$) oU = np.copy(A)(matriz copia de $A$) ¡Ojo! no valeU = A.L = np.eye(n)(matriz identidad $n\times n$)
Comprobar que la función funciona bien con el sistema $Ax=b$ cuyas matrices $A$ y $b$ se generan de la forma siguiente
n = 5
np.random.seed(4)
A = np.random.random((n,n))
b = np.random.random(n)
%run Ejercicio5.py