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} $$Con from __future__ import division obligamos a que el resultado de la división con / sea siempre float. Si no lo incluímos, la división de dos enteros con / sería también un entero.
Importamos los módulos numpy y pprint. El módulo pprint proporciona una forma de imprimir datos más ordenada.
from __future__ import division
import numpy as np
import pprint
Para el sistema $$ L\mathbf{x}=\mathbf{b}, $$ con $L$ de la forma $$ 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}x_{i}=b_{i}. $$
Por lo tanto, para $i=1,2,\ldots,n$, si $l_{ii}$ es no nulo, entonces $$ 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). $$
Si $l_{ii}= 0$, acabar.
En caso contrario $$ 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$.
# define las matrices A y b usando arrays numpy
L = np.array([[2, 0, 0], [1, 2, 0], [1, 1, 2]])
b = np.array([2, -1, 0])
Comprobar que la función también funciona correctamente con el sistema $Lx=b$ de matrices $L$ y $b$ generadas de la siguiente forma
n = 5
np.random.seed(1) # para que los números aleatorios se repitan
L = np.random.random((n,n)) # genera una matriz aleatoria nxn
L = np.tril(L) # hacer cero los elementos por encima de la diagonal
b = np.random.random((n,)) # genera un vector aleatorio de dimensión n
print 'L'
pprint.pprint(L)
print 'b'
pprint.pprint(b)
%run Ejercicio1.py
Análogamente, para el sistema
$$ U\mathbf{x}=\mathbf{b}' $$con $U$ de la forma
$$ 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$, si $l_{ii}$ no se anula
$$ 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), $$
- Si $u_{ii}=0$, acabar.
- Caso contrario,
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) $$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,))
print 'U'
pprint.pprint(U)
print 'b'
pprint.pprint(b)
%run Ejercicio2.py
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:
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.
- Si $a_{kk}=0$, acabar.
- Caso contrario, $i=k+1,\ldots,n,\quad \quad b_{i}=b_{i}-\dfrac{a_{ik}}{a_{kk}}b_{k}. $
- Para $i=k+1,\ldots,n, \quad a_{ij}=a_{ij}-\dfrac{a_{ik}}{a_{kk}}a_{ik}. $
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) $$U y c.n = 5
np.random.seed(3)
A = np.random.random((n,n))
b = np.random.random((n,))
print '\n A \n'
pprint.pprint(A)
print '\n b \n'
pprint.pprint(b)
%run Ejercicio3.py
Expresaremos $A$ como el producto de dos matrices: una matriz triangular inferior, $L$, y una matriz triangular superior, $U$,
$$ A=LU. $$Si tenemos el sistema
$$ A\mathbf{x}=\mathbf{b}, $$donde $A$ admite la factorización $LU$
$$ A\mathbf{x}=\mathbf{b}\Longleftrightarrow LU\mathbf{x}=\mathbf{b} $$Una vez factorizada la matriz $A$ solucionaremos el sistem $Ax=b$ en dos pasos
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 & 2\\ 1 & 1 & 2 \end{array}\right)\quad\quad b=\left(\begin{array}{c} 2\\ 4\\ 6 \end{array}\right) $$LU(A) que da la factorización $LU$ de $A$.sust_regre y sust_prog para resolver los sistemas resultantes.n = 5
np.random.seed(4) # for the random numbers to be repeatable
A = np.random.random((n,n)) # generate a random nxn matrix
b = np.random.random((n,)) # generate a random vector of dimension n
print '\n A \n'
pprint.pprint(A)
print '\n b \n'
pprint.pprint(b)
%run Ejercicio4.py
Los métodos iterativos son, en general, más eficientes que los métodos directos para resolver sistemas de ecuaciones lineales poco densos (con muchos elementos de la matriz de coeficientes igual a cero).
Reescribimos el sistema
$$ A\mathbf{x}=\mathbf{b}\Longleftrightarrow\mathbf{x}=B\mathbf{x}+\mathbf{c} $$donde $B$ es una matriz de orden $n$ y $c$ es un vector columna de dimensión $n$.
Entonces, si $\mathbf{x}^{(0)}$ es el punto inicial, para $k=1,2,\ldots$
$$ \mathbf{x}^{(k)}=B\mathbf{x}^{(k-1)}+\mathbf{c} $$
- Para $i=1,2,\ldots,n$ calcular $$ x_{i}^{(k)}=\frac{1}{a_{ii}}(b_{i}-\sum_{j=1}^{i-1}a_{ij}x_{j}^{(k-1)}-\sum_{j=i+1}^{n}a_{ij}x_{j}^{(k-1)}) $$
- Si se cumple el criterio de parada, tomar $\mathbf{x}^{(k)}$ como solución aproximada.
Escribir una jacobi(A,b,tol) para resolver $Ax=b$
por el método de Jacobi.
tol = 1.e-6 con norma 2. Es decir, usaremos como criterio de parad np.linalg.norm(x-xant) < tol, con x la solución aproximada en esta iteración y xant la solución aproximada en la iteración anterior. Tomar como punto inicial $\mathbf{x^{(0)}=0}.$x y el número de iteraciones num_iter.xs = np.linalg.solve(A,b).n = 5
np.random.seed(5)
A = np.random.random((n,n)) + 2*np.eye(n)
b = np.random.random((n,))
print '\n A \n'
pprint.pprint(A)
print '\n b \n'
pprint.pprint(b)
%run Ejercicio5.py
- Para $i=1,2,\ldots,n$ calcular $$ x_{i}^{(k)}=\frac{1}{a_{ii}}(b_{i}-\sum_{j=1}^{i-1}a_{ij}x_{j}^{(k)}-\sum_{j=i+1}^{n}a_{ij}x_{j}^{(k-1)}) $$
- Si se satisface el criterio de parada, tomar $\mathbf{x}^{(k)}$ como solución aproximada
Escribir una función gauss_seidel(A,b,tol) para resolver $Ax=b$
por el método de Gauss-Seidel.
%run Ejercicio6.py