Resolución de sistemas lineales: métodos iterativos¶
import numpy as np
Establecemos el formato
np.set_printoptions(precision = 8) # ocho 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} $$Métodos iterativos¶
Los métodos iterativos son, en general, más eficientes que los métodos directos para resolver sistemas de ecuaciones lineales poco densos o dispersos (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} $$y generamos una sucesión de vectores (matrices columna) $\left\{ \mathbf{x}^{(m)} \right\}$ que, caso que converja, converge a la solución del sistema $\mathbf{x}=\left(x_1,x_2,\ldots,x_n\right)^T.$
Método de Jacobi¶
Si tenemos un sistema $4\times 4$
$$ \begin{array}{cll} a_{11}x_{1}+a_{12}x_{2}+a_{13}x_{3}+a_{14}x_{4} & = & b_{1}\\ a_{21}x_{1}+a_{22}x_{2}+a_{23}x_{3}+a_{24}x_{4} & = & b_{2}\\ a_{31}x_{1}+a_{32}x_{2}+a_{33}x_{3}+a_{34}x_{4} & = & b_{3}\\ a_{41}x_{1}+a_{42}x_{2}+a_{43}x_{3}+a_{44}x_{4} & = & b_{4} \end{array} $$y despejamos la primera incógnita en la primera ecuación, la segunda incógnita en la segunda y así sucesivamente
$$ \begin{array}{cll} x_{1}& = &\dfrac{b_{1}-a_{12}x_{2}-a_{13}x_{3}-a_{14}x_{4}}{a_{11}} \\ x_{2} & = & \dfrac{b_{2}-a_{21}x_{1}-a_{23}x_{3}-a_{24}x_{4}}{a_{22}}\\ x_{3} & = & \dfrac{b_{3}-a_{31}x_{1}-a_{32}x_{2}-a_{34}x_{4}}{a_{33}}\\ x_{4} & = & \dfrac{b_{4}-a_{41}x_{1}-a_{42}x_{2}-a_{43}x_{3}}{a_{44}} \end{array} $$Y obtenemos los valores de la primera iteración a partir de los valores iniciales
$$ \begin{array}{cll} x_{1}^{{\color{red}{(1)}}} & = & \dfrac{b_{1}-a_{12}x_{2}^{(0)}-a_{13}x_{3}^{(0)}-a_{14}x_{4}^{(0)}}{a_{11}} \\ x_{2}^{{\color{red}{(1)}}} & = & \dfrac{b_{2}-a_{21}x_{1}^{(0)}-a_{23}x_{3}^{(0)}-a_{24}x_{4}^{(0)}}{a_{22}}\\ x_{3}^{{\color{red}{(1)}}} & = & \dfrac{b_{3}-a_{31}x_{1}^{(0)}-a_{32}x_{2}^{(0)}-a_{34}x_{4}^{(0)}}{a_{33}}\\ x_{4}^{{\color{red}{(1)}}}& = & \dfrac{b_{4}-a_{41}x_{1}^{(0)}-a_{42}x_{2}^{(0)}-a_{43}x_{3}^{(0)}}{a_{44}} \end{array} $$Si generalizamos tenemos el siguiente algoritmo
Algoritmo¶
- Elegir una aproximación inicial $\mathbf{x}^{(0)}$
- Para $k=1,2,\ldots,\mathrm{MaxIter}$
- Para $i=1,2,\ldots,n$, calcular $$x_{i}^{{\color{red}{(k)}}}=\dfrac{1}{a_{ii}}\left(b_{i}-\sum_{j=1}^{i-1}a_{ij}x_{j}^{(k-1)}-\sum_{j=i+1}^{n}a_{ij}x_{j}^{(k-1)}\right)$$
- Si se cumple el criterio de parada, tomar $\mathbf{x}^{(k)}$ como aproximación a la solución.
Escribir una jacobi(A,b,tol,maxiter=1000) para resolver $Ax=b$
por el método de Jacobi.
- Resolver
$$
\begin{cases}
5x_{1}+x_{2}-x_{3}-x_{4} & =1\\
x_{1}+4x_{2}-x_{3}+x_{4} & =1\\
x_{1}+x_{2}-5x_{3}-x_{4} & =1\\
x_{1}+x_{2}+x_{3}-4x_{4} & =1
\end{cases}
$$
por el método de Jacobi, usando una tolerancia absoluta
tol = 1.e-6con norma 2. Es decir, usaremos como criterio de paradanp.linalg.norm(x-xant) < tol, conxla solución aproximada en esta iteración yxantla solución aproximada en la iteración anterior. Tomar como punto inicial $\mathbf{x^{(0)}=0}.$A = np.array([[5.,1,-1,-1],[1,4,-1,1],[1,1,-5,-1],[1,1,1,-4]]) b = np.array([1.,1,1,1])
- Como argumento de salida, dar la solución
xy el número de iteracionesnum_iter. - Comprobar que la solución es correcta comparándola con
xs = np.linalg.solve(A,b). - Comprobar que también funciona correctamente para resolver el sistema $Ax=b$ de matrices $A$ y $b$ generadas de la forma siguiente
n = 9 A1 = np.diag(np.ones(n))*2 A2 = np.diag(np.ones(n-1),1) A = A1 + A2 + A2.T b = np.concatenate((np.arange(1,6),np.arange(4,0,-1)))*1.
Nota
Para asignar un vector o matriz se debe hacer una copia
xant = np.copy(x)
%run Ejercicio1.py
Método de Gauss-Seidel¶
Si tenemos un sistema $4\times 4$ $$ \begin{array}{cll} a_{11}x_{1}+a_{12}x_{2}+a_{13}x_{3}+a_{14}x_{4} & = & b_{1}\\ a_{21}x_{1}+a_{22}x_{2}+a_{23}x_{3}+a_{24}x_{4} & = & b_{2}\\ a_{31}x_{1}+a_{32}x_{2}+a_{33}x_{3}+a_{34}x_{4} & = & b_{3}\\ a_{41}x_{1}+a_{42}x_{2}+a_{43}x_{3}+a_{44}x_{4} & = & b_{4} \end{array} $$
y despejamos la primera incógnita en la primera ecuación, la segunda incógnita en la segunda y así sucesivamente
$$ \begin{array}{cll} x_{1}& = &\dfrac{b_{1}-a_{12}x_{2}-a_{13}x_{3}-a_{14}x_{4}}{a_{11}} \\ x_{2} & = & \dfrac{b_{2}-a_{21}x_{1}-a_{23}x_{3}-a_{24}x_{4}}{a_{22}}\\ x_{3} & = & \dfrac{b_{3}-a_{31}x_{1}-a_{32}x_{2}-a_{34}x_{4}}{a_{33}}\\ x_{4} & = & \dfrac{b_{4}-a_{41}x_{1}-a_{42}x_{2}-a_{43}x_{3}}{a_{44}} \end{array} $$Y obtenemos los valores de la primera iteración a partir de los valores iniciales, pero conforme vamos teniendo resultados nuevos de una variable, los actualizamos
$$ \begin{array}{cll} x_{1}^{\color{red}{(1)}} & = & \dfrac{b_{1}-a_{12}x_{2}^{(0)}-a_{13}x_{3}^{(0)}-a_{14}x_{4}^{(0)}}{a_{11}} \\ x_{2}^{\color{red}{(1)}} & = & \dfrac{b_{2}-a_{21}x_{1}^{\color{red}{(1)}}-a_{23}x_{3}^{(0)}-a_{24}x_{4}^{(0)}}{a_{22}}\\ x_{3}^{\color{red}{(1)}} & = & \dfrac{b_{3}-a_{31}x_{1}^{\color{red}{(1)}}-a_{32}x_{2}^{\color{red}{(1)}}-a_{34}x_{4}^{(0)}}{a_{33}}\\ x_{4}^{\color{red}{(1)}}& = & \dfrac{b_{4}-a_{41}x_{1}^{\color{red}{(1)}}-a_{42}x_{2}^{\color{red}{(1)}}-a_{43}x_{3}^{\color{red}{(1)}}}{a_{44}} \end{array} $$Si generalizamos tenemos el siguiente algoritmo
Algoritmo¶
El algoritmo que hemos aplicado se expresa formalmente
- Elegir una aproximación inicial $\mathbf{x}^{(0)}$
- Para $k=1,2,\ldots,$MaxIter
- Para $i=1,2,\ldots,n$, calcular $$x_{i}^{\color{red}{(k)}}=\dfrac{1}{a_{ii}}\left(b_{i}-\sum_{j=1}^{i-1}a_{ij}x_{j}^{\color{red}{(k)}}-\sum_{j=i+1}^{n}a_{ij}x_{j}^{(k-1)}\right)$$
- Si se cumple el criterio de parada, tomar $\mathbf{x}^{(k)}$ como aproximación a la solución.
Escribir una función gauss_seidel(A,b,tol,maxiter=1000) para resolver $Ax=b$
por el método de Gauss-Seidel.
- Resolver los mismos sistemas del ejercicio anterior usando Gauss-Seidel con la misma tolerancia y criterio de parada.
- Escribir las mismos argumentos de salida que en el ejercicio anterior y comparar el número de iteraciones.
%run Ejercicio2.py
Método de relajación¶
Algoritmo¶
- Elegir el punto inicial $\mathbf{x}^{(0)}$
- Para $k=1,2,\ldots,MaxIter$
- Para $i=1,2,\ldots,n$ calcular $$ x_{i}^{\color{red}{(k)}}=\frac{\omega}{a_{ii}}\left(b_{i}-\sum_{j=1}^{i-1}a_{ij}x_{j}^{\color{red}{(k)}}-\sum_{j=i+1}^{n}a_{ij}x_{j}^{(k-1)}\right)+(1-\omega) \, x_i^{(k-1)} $$
- Si se satisface el criterio de parada, tomar $\mathbf{x}^{(k)}$ como solución aproximada
Si converge, entonces $\omega \in (0,2).$
Escribir una función relajacion(A,b,w,tol,maxiter=1000) para resolver $Ax=b$
por el método de relajación.
- Resolver el segundo sistema del ejercicio anterior usando el método de relajación con $\omega = 1.5$ y la misma tolerancia y criterio de parada.
- Escribir las mismos argumentos de salida que en el ejercicio anterior y comparar el número de iteraciones.
%run Ejercicio3.py
Ejercicios propuestos¶
Relajación: estudio del parámetro $\omega$¶
Escribir un programa que calcule el valor de $\omega$ óptimo para encontrar la solución del sistema del ejercicio 3 con el método de sobrerrelajación. Usar maxiter = 1000 y tol = 1.e-12. Calcular y dibujar el número de iteraciones para los valores de $\omega$ w = np.arange(0.1,2.1,0.01).
Nota
- Usar
k = np.argmin(num_iter)para encontrar el índice del valor óptimo.
%matplotlib inline
%run Ejercicio4.py
Sistemas tridiagonales¶
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 = 9
Ar = np.zeros((n,3))
Ar[:,0] = np.concatenate((np.array([0]),np.ones((n-1),)))
Ar[:,1] = np.ones((n),)*2
Ar[:,2] = np.concatenate((np.ones((n-1),),np.array([0])))
b = np.concatenate((np.arange(1,6),np.arange(4,0,-1)))*1.
Modificar el código de las funciones de los tres ejercicios anteriores de forma que admitan la matriz $A$ en la forma $A_r$ (matriz rectangular).
- Resolver el sistema con la misma tolerancia y criterio de parada que en los ejercicios anteriores.
- Escribir los mismos argumentos de salida que en los ejercicios anteriores y comparar el número de iteraciones.
- Tomar, para el método de relajación, $\omega = 1.5$.
%run Ejercicio5.py