Con respecto al separador decimal en este documento, ya que los resultados son automáticos y se convierten al formato del notebook sin manipulación, aparecerá "." en lugar de ",", ya que este corresponde al separador decimal por defecto del lenguaje de programación (python).
Igualmente el formateo de los vectores del resultado, corresponden a la forma string por defecto, por lo cual al iterar se indicara como Matriz un arreglo horizontal de datos, que corresponden a vectores columna de la solución.
- Edward Villegas (Curso de Métodos Numéricos, UdeM).El uso de métodos iterativos para sistemas matriciales $A\vec{x}=\vec{b}$ permite solventar la acumulación de error numérico que se presenta por la solución de métodos directos usando la aritmética de la máquina (se recordará de los métodos directos la abundancia de operaciones, principalmente de restas y divisiones). Por ello, se buscará la solución a traves de un esquema de punto fijo, $\vec{x}^{(k+1)}=T\vec{x}^{(k)}+\vec{c}$.
Tras dicha forma iterativa, para una iteración $k+1$, $\vec{x}^{(k+1)}=T^{k+1}\vec{x}^{(0)}+\sum\limits_{i=0}^kT^{i}\vec{c}$. Lo cual, si $k\rightarrow\infty$, solo es convergente si las potencias de $T$ tienden a anularse ($T^{k+1}\rightarrow 0_{K_{n,n}}$), lo que corresponde a una matriz convergente ($\rho(T)<1$).
Se propone considerar la matriz $A$ como la suma de tres matrices, de la forma que se expone a continuación. \begin{equation} A = D - L - U \end{equation}
\begin{equation} A = \begin{pmatrix} a_{0,0} & a_{0,1} & \ldots & a_{0,n-1}\\ a_{1,0} & a_{1,1} & \ldots & a_{1,n-1}\\ \vdots & \vdots & \ddots & \vdots\\ a_{n-1,0} & a_{n-1,1} & \ldots & a_{n-1,n-1}\\ \end{pmatrix} \end{equation}\begin{equation} D = \begin{pmatrix} a_{0,0} & 0 & \ldots & 0\\ 0 & a_{1,1} & \ldots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \ldots & a_{n-1,n-1}\\ \end{pmatrix} \end{equation}\begin{equation} U = \begin{pmatrix} 0 & -a_{0,1} & -a_{0, i} & -a_{0,n-1}\\ 0 & 0 & -a_{i,j} & -a_{1,n-1}\\ \vdots & \vdots & \ddots & -a_{i,n-1}\\ 0 & 0 & \ldots & 0\\ \end{pmatrix} \end{equation}\begin{equation} L = \begin{pmatrix} 0 & 0 & \ldots & 0\\ -a_{1,0} & 0 & \ldots & 0\\ -a_{i,0} & -a_{i,j} & \ddots & \vdots\\ -a_{n-1,0} & -a_{n-1,1} & -a_{n-1,j} & 0\\ \end{pmatrix} \end{equation}Debido a la gran cantidad de elementos ceros, una forma indicial de los métodos será mucho más eficiente que la forma matricial en el momento de la implementación (reducción de numero de operaciones y de almacenamiento de memoria).
El método de Jacobi propone el siguiente desarrollo matricial.
\begin{eqnarray} A\vec{x} & = & \vec{b} \\ (D - L - U)\vec{x} & = & \vec{b} \\ D\vec{x} & = & (L+U)\vec{x} + \vec{b} \\ \vec{x} & = & D^{-1}(L+U)\vec{x} + D^{-1}\vec{b} \\ \vec{x}^{(k+1)} & = & T_J \vec{x}^{(k)} + \vec{c}_J \quad \text{(Forma matricial)} \end{eqnarray}Esto puede ser llevado a una forma indicial, transformando los productos matriz-vector a su forma de producto punto fila-vector, y las sumas elemento a elemento. Igualmente, teniendo en cuenta que la inversa de una matriz diagonal es otra matriz diagonal en la que sus elementos son los inversos multiplicativos de los elementos originales.
\begin{equation} x_i^{(k+1)} = \frac{-\sum_{j=0}^{i-1}a_{ij}x_j^{(k)} -\sum_{j=i+1}^{n-1}a_{ij}x_j^{(k)} + b_i}{a_{ii}} \end{equation}Para este método es necesario que los $a_{ii} \neq 0$ para que $D$ (matriz diagonal) sea invertible. Si se realiza manualmente, tenga en cuenta que la inversa de una matriz diagonal, es una matriz diagonal con los inversos multiplicativos de los elementos originales. Y su multiplicación por matrices y vectores se obtiene rapidamente al multiplicar el elemento diagonal por la fila correspondiente de la otra matriz o vector.
Si $A$ es estrictamente diagonal dominante, este método converge a la solución única de $A\vec{x} = \vec{b}$.
Nótese que en la forma indicial del método de Jacobi, se observa que en la medida que se calcula una determinada componente $x_i$, las componentes anteriores ya son conocidas sus actualizaciones. De manera, que una mejor aproximación de esta componente es usar la información actualizada de las componentes anteriores a esta.
\begin{equation} x_i^{(k+1)} = \frac{-\sum_{j=0}^{i-1}a_{ij}x_j^{(k+1)} -\sum_{j=i+1}^{n-1}a_{ij}x_j^{(k)} + b_i}{a_{ii}} \end{equation}Esta forma indicial, puede ser llevada a la forma matricial si consideramos que
\begin{eqnarray} A\vec{x} & = & \vec{b} \\ (D-L-U)\vec{x} & = & \vec{b} \\ (D-L)\vec{x} & = & U\vec{x} + \vec{b} \\ \vec{x} & = & (D-L)^{-1}U\vec{x} + (D-L)^{-1}\vec{b} \\ \vec{x}^{(k+1)} & = & T_G\vec{x}^{(k)} + \vec{c}_G \end{eqnarray}Para este método es necesario que los $a_{ii} \neq 0$ para que $(D-L)$ (matriz triangular inferior) sea invertible.
Si $A$ es estrictamente diagonal dominante, este método converge a la solución única de $A\vec{x} = \vec{b}$.
El método de SOR (Successive Over-Relaxation) parte de la unión del método de Gauss-Seidel con el principio de punto fijo bajo la condición de convergencia, y toma una proporción de cada igualdad para formar una nueva forma iterativa.
\begin{equation} \begin{matrix} (1-\omega)\vec{x}^{(k+1)} & = & (1-\omega)\vec{x}^{(k)} \\ \omega \vec{x}^{(k+1)} & = & \omega (D-L)^{-1}(U\vec{x}^{(k)} + \vec{b}) \\ \hline \vec{x}^{(k+1)} & = & (1-\omega)\vec{x}^{(k)} + \omega (D-L)^{-1}(U\vec{x}^{(k)} + \vec{b}) \end{matrix} \end{equation}Luego, tras manipulación algebraica y aplicaciones de la condición de punto fijo en la convergencia, se sigue que
\begin{eqnarray} (D-L)\vec{x}^{(k+1)} - (1-\omega)(D-L)\vec{x}^{(k)} &=& \omega (U\vec{x}^{(k)} + \vec{b}) \\ D\vec{x}^{(k+1)} - L\vec{x}^{(k+1)} + L\vec{x}^{(k)} - \omega L\vec{x}^{(k)} &=& (D(1-\omega)+\omega U)\vec{x}^{(k)} + \omega\vec{b} \\ D\vec{x}^{(k+1)} - \omega L\vec{x}^{(k)} &=& (D(1-\omega)+\omega U)\vec{x}^{(k)} + \omega\vec{b} \\ (D - \omega L)\vec{x}^{(k+1)} & = & (D(1-\omega)+\omega U)\vec{x}^{(k)} + \omega\vec{b} \\ \vec{x}^{(k+1)} &=& (D - \omega L)^{-1}(D(1-\omega)+\omega U)\vec{x}^{(k)}\\ && + (D - \omega L)^{-1}\omega\vec{b} \\ \vec{x}^{(k+1)} &=& T_{SOR}\vec{x}^{(k)} + \vec{c}_{SOR} \end{eqnarray}De la misma forma, la forma indicial se obtiene como
\begin{equation} \begin{matrix} (1-\omega)\vec{x}_i^{(k+1)} & = & (1-\omega)\vec{x}_i^{(k)} \\ \omega \vec{x}_i^{(k+1)} & = & \omega \frac{-\sum_{j=0}^{i-1}a_{ij}x_j^{(k+1)} -\sum_{j=i+1}^{n-1}a_{ij}x_j^{(k)} + b_i}{a_{ii}} \\ \hline \vec{x}_i^{(k+1)} & = & (1-\omega)\vec{x}_i^{(k)} + \frac{-\sum_{j=0}^{i-1}a_{ij}x_j^{(k+1)} -\sum_{j=i+1}^{n-1}a_{ij}x_j^{(k)} + b_i}{a_{ii}} \end{matrix} \end{equation}Donde $0 \leq \omega \leq 2$, para que el método converja. Esto se puede asegurar si $A$ es matriz definida positiva. Sin embargo, la convergencia a la solución única de $A\vec{x}=\vec{b}$ no se asegura, y es fuertemente dependiente de la selección de $\omega$.
Para seleccionar el método más adecuado, se usa como criterio el radio espectral. A menor radio espectral, mayor velocidad de convergencia.
Para la definición del error, criterio de parada, se puede definir el error a tráves de las normas p, $||\vec{x}^{(k+1)}-\vec{x}^{(k)}||$ o $||A\vec{x}^{(k+1)}-\vec{b}||$. Estas definiciones son las versiones absolutas del error, las cuales también pueden definirse como relativas al dividir por la norma del vector verdadero. Entre las dos definiciones (una solo dependiente de las iteraciones, y la segunda dependiente de reconstruir el valor verdadero), la primera solo permite definir si el método esta convergiendo, mientras la segunda define si converge a la solución única.
In [1]:
from sympy import * # Paquete de calculo simbolico/exacto
init_printing() # Habilita la impresión de simbolos en formato LaTeX desde SymPy
from IPython.display import Latex # Habilita el uso LaTeX desde IPython
w = Symbol('w', real=True) # Se declara la variable simbolica para el parametro de relajación
In [2]:
A = Matrix([[4, 3, 0], [3, 4, -1], [0, -1, 4]])
b = Matrix([[24],[30],[-24]])
D = Matrix([[4 , 0 , 0],[0 , 4 , 0],[0, 0, 4]])
L = Matrix([[0, 0, 0], [-3, 0, 0], [0, 1, 0]])
U = Matrix([[0, -3, 0], [0, 0, 1], [0, 0, 0]])
T_J = D**(-1)*(U+L)
T_G = (D-L)**(-1)*U
T_SOR = (D - w*L)**(-1)*(D*(1-w)+w*U)
c_J = D**(-1)*b
c_G = (D-L)**(-1)*b
c_SOR = w*(D - w*L)**(-1)*b
In [3]:
Latex('\\begin{eqnarray}' + latex(A)+'\\vec{x}&=&' + latex(b) +
'\\\\ A&=&' + latex(A)+ '\\\\ D &=&' + latex(D)+'\\\\ L&=&' + latex(L)+
'\\\\ U&=&' + latex(U)+ '\\\\ T_J&=&' + latex(T_J) + '\\\\ T_G&=&' + latex(T_G) +
'\\\\ T_{SOR}&=&' + latex(T_SOR) +'\\\\ \\vec{c}_J&=&'+latex(c_J)+'\\\\ \\vec{c}_G&=&'+
latex(c_G)+'\\\\ \\vec{c}_{SOR}&=&'+latex(c_SOR)+'\\end{eqnarray}')
Out[3]:
A continuación obtenemos los autovalores de las matrices $T$, con el fin de calcular posteriormente el radio espectral (el cual permite determinar la convergencia del método).
In [4]:
eigen_J = T_J.eigenvals() # Autovalores de T_J
Latex('$'+latex(eigen_J)+'$')
Out[4]:
In [5]:
eigen_G = T_G.eigenvals()
Latex('$'+latex(eigen_G)+'$')
Out[5]:
Para el caso del método de SOR, es necesario asignar un valor para el parámetro de relajación, en este caso $\omega =$ 1,25.
In [6]:
T_SOR = (D - 1.25*L)**(-1)*(D*(1-1.25)+1.25*U)
c_SOR = 1.25*(D - 1.25*L)**(-1)*b
Latex('\\begin{eqnarray} T_{SOR} &=&' + latex(T_SOR) + '\\\\ \\vec{c}_{SOR}&=&' +
latex(c_SOR) + '\\end{eqnarray}')
Out[6]:
In [7]:
eigen_SOR = T_SOR.eigenvals()
Latex('$'+latex(eigen_SOR)+'$')
Out[7]:
Se define el radio espectral $\rho(T) = \max \lbrace |\lambda_i| \rbrace$.
In [8]:
rho_J = max([abs(i) for i in list(eigen_J.keys())])
rho_G = max([abs(i) for i in list(eigen_G.keys())])
rho_SOR = max([abs(i) for i in list(eigen_SOR.keys())])
In [9]:
Latex('\\begin{eqnarray} \\rho (T_J) & = &' + latex(rho_J) + '=' + latex(rho_J.evalf()) +
'\\\\ \\rho (T_G) & = &' + latex(rho_G) + '=' + latex(rho_G.evalf()) +
'\\\\ \\rho (T_{SOR}) & = &' + latex(rho_SOR) + '=' + latex(rho_SOR.evalf()) +
'\\end{eqnarray}')
Out[9]:
Observamos que la matriz $T$ con menor radio espectral es $T_{SOR}$, por lo que esperamos que de las iteraciones de los 3 métodos, este sea el que logre una convergencia más rápida, tal como se observa a continuación.
Solución exacta $\vec{x} = (3, 4, -5) $.
In [10]:
# Aplicación de los métodos iterativos
x_J = Matrix([[0], [0], [0]])
x_J_list = [x_J]
x_G = Matrix([[0], [0], [0]])
x_G_list = [x_G]
x_SOR = Matrix([[0], [0], [0]])
x_SOR_list = [x_SOR]
for i in range(11):
x_J = T_J*x_J + c_J
x_J_list.append(x_J)
x_G = T_G*x_G + c_G
x_G_list.append(x_G)
x_SOR = T_SOR*x_SOR + c_SOR
x_SOR_list.append(x_SOR)
In [11]:
# Método de Jacobi
M = x_J
for i in range(1, len(x_J_list)):
print(x_J_list[i].evalf())
In [12]:
# Método de Gauss Seidel
M = x_G
for i in range(1, len(x_G_list)):
print(x_G_list[i].evalf())
In [13]:
# Método de SOR
M = x_SOR
for i in range(1, len(x_SOR_list)):
print(x_SOR_list[i].evalf())