1. El método iterativo de Jacobi
Referencia: Burden 7.3 p334, Rodríguez 5.1 p154, Chapra 11.2.1p312
Los métodos iterativos para sistemas de ecuaciones, son semejantes al método de punto fijo para búsqueda de raíces, requieren un punto inicial para la búsqueda "de la raíz o solución" que satisface el sistema. Se plantea considerar aproximaciones o iteraciones de forma sucesiva desde un punto de partida o inicial X0, hacia una solución del sistema de ecuaciones.
Las iteraciones pueden ser o no, convergentes.
Para describir el método iterativo de Jacobi, se usa como ejemplo: un sistema de 3 incógnitas y 3 ecuaciones, diagonalmente dominante, planteado en su forma matricial.
\begin{cases} a_{0,0} x_0+a_{0,1}x_1+a_{0,2} x_2 = b_{0} \\ a_{1,0} x_0+a_{1,1}x_1+a_{1,2} x_2 = b_{1} \\ a_{2,0} x_0+a_{2,1}x_1+a_{2,2} x_2 = b_{2} \end{cases}La expresión para encontrar nuevos valores usando la matriz de coeficientes A, el vector de constantes B y como incógnitas X es:
x_i = \bigg(b_i - \sum_{j=0, j\neq i}^n a_{i,j}x_j\bigg) \frac{1}{a_{ii}}El método de Jacobi, reescribe las expresiones despejando la incógnita de la diagonal en cada ecuación. El patrón que se observa en las ecuaciones despejadas, se describe en la fórmula para xi de la fila i, que se usará en el algoritmo en Python.
x_0 = \frac{b_{0} -a_{0,1}x_1 -a_{0,2} x_2 }{a_{0,0}} x_1 = \frac{b_{1} -a_{1,0} x_0 -a_{1,2} x_2}{a_{1,1}} x_2 = \frac{b_{2} -a_{2,0} x_0 - a_{2,1} x_1}{a_{2,2}}La convergencia del método iterativo se puede revisar usando el número de condición de la matriz de coeficientes A. Un resultado cercano a 1, implica que el sistema será convergente.
cond(A) = ||A|| ||A-1||
np.linalg.cond(A)
Usando de forma experimental los ejercicios de las evaluaciones de la sección matriciales iterativos, el número de condición no debería superar el valor de 7 para que el sistema sea convergente. Puede revisar que si es mayor, el método tiende a no ser convergente. Aunque los textos del libro no lo indican directamente,
2. Ejercicio
Referencia: 1Eva2011TII_T2 Sistema de Ecuaciones, diagonal dominante
Considere el siguiente sistema de ecuaciones A.X=B dado por:
\begin{cases} -2x+5y+9z=1\\7x+y+z=6\\-3x+7y-z=-26\end{cases}Luego del pivoteo parcial por filas, el sistema de ecuaciones a usar para el método de Jacobi es:
\begin{cases} 7x+y+z=6\\-3x+7y-z=-26\\-2x+5y+9z=1\end{cases}Desarrolle el ejercicio usando el método de Jacobi, tomando como vector inicial X0=[0,0,0] y tolerancia de 0.0001.
3. Desarrollo Analítico
Para el método de Jacobi, se usan las ecuaciones reordenas con el pivoteo parcial por filas, para que en lo posible sea diagonal dominante que mejora la convergencia. El resultado del algoritmo a partir de las ecuaciones presentadas al inicio es:
Pivoteo parcial por filas AB:
[[ 7. 1. 1. 6.]
[ -3. 7. -1. -26.]
[ -2. 5. 9. 1.]]
A:
[[ 7. 1. 1.]
[-3. 7. -1.]
[-2. 5. 9.]]
B:
[ 6. -26. 1.]Se despejan la incógnita de la diagonal en cada ecuación.
x=\frac{6-y-z}{7} y= \frac{-26+3x+z}{7} z=\frac{1+2x-5y}{9}El vector inicial X0=[0,0,0] representa un punto de partida, valores estimados o semilla para buscar la solución a las variables x0,y0,z0.
Como las ecuaciones son el comportamiento del 'sistema', se usan con X0 para obtener un nuevo y mejor valor estimado para la solución.
itera = 0
X_0=[x_0,y_0,z_0] X_0=[0,0,0] x_1=\frac{6-(0)-(0)}{7} =\frac{6}{7}=0.8571 y_1= \frac{-26+3(0)+(0)}{7} =\frac{-26}{7}=-3.7142 z_1=\frac{1+2(0)-5(0)}{9} =\frac{1}{9}=0.1111 X_1=[0.8571, -3.7142, 0.1111]Para revisar si se está más cerca de la solución, se calculan las diferencias, tramos o desplazamientos entre los valores nuevos y anteriores:
\text{diferencias}=X_1-X_0 =[0.8571-0, -3.7142-0, 0.1111-0] [0.8571, -3.7142, 0.1111]
-[0. , 0. , 0. ]
__________________________
[0.8571, -3.7142, 0.1111]Para estimar el valor del error, se presenta un valor simplificado o escalar de las diferencias (norma del error). Por simplicidad se usa el valor de magnitud máxima de los errores. Se usa como variable 'errado', para evitar el conflicto de nombre de variable error usada como palabra reservada en Python:
\text{errado}= max \Big|\text{diferencias}\Big| \text{errado}= max \Big|[0.8571, -3.7142, 0.1111]\Big| = 3.7142Al comparar el valor errado con la tolerancia, se observa que el error aún es mayor que lo tolerado. Por lo que se continúa con otra iteración.
itera = 1
X_1=[x_1,y_1,z_1] X_1=[0.8571, -3.7142, 0.1111] x_2=\frac{6-(-3.7142)-(0.1111)}{7} = 1.3718 y_2= \frac{-26+3(0.8571)+(0.1111)}{7} =-3.3310 z_2=\frac{1+2(0.8571)-5(-3.7142)}{9} =2.3650 X_2=[1.3718, -3.3310, 2.3650] \text{diferencias}=|X_2-X_1| =|[0.8571-1.3718, -3.3310-(-3.7142), 2.3650-0.1111]| [1.3718, -3.3310, 2.3650]
-[0.8571, -3.7142, 0.1111]
__________________________
[0.51474, -0.38322, 2.25397]itera = 2
X_2=[1.3718, -3.3310, 2.3650] x_3=\frac{6-(-3.3310)-(2.3650)}{7} = 0.9951 y_3= \frac{-26+3(1.3718)+(2.3650)}{7} = -2.7884 z_3=\frac{1+2(1.3718)-5(-3.3310)}{9} = 2.2665 X_3=[0.99514, -2.7884, 2.2665] \text{diferencias}=|X_3-X_2| =|[0.99514-1.3718, -2.7884-(-3.3310), 2.2665-2.3650]| [ 0.9951, -2.7884, 2.2665]
-[ 1.3718, -3.3310, 2.3650]
__________________________
[-0.3767, 0.5426, 0.0985]
En cada iteración el valor de error disminuye, por lo que se estima que el método converge.
Al usar el algoritmo, de los resultados se observa que se detiene luego de 14 iteraciones, los errores son menores a la tolerancia:
Iteraciones Jacobi
itera,[X]
,errado,|diferencia|
0 [0. 0. 0.]
nan
1 [ 0.85714 -3.71429 0.11111]
3.7142857142857144 [0.85714 3.71429 0.11111]
2 [ 1.37188 -3.33107 2.36508]
2.2539682539682544 [0.51474 0.38322 2.25397]
3 [ 0.99514 -2.78847 2.26657]
0.5425979915775843 [0.37674 0.5426 0.09851]
4 [ 0.9317 -2.964 1.8814]
0.3851635892452223 [0.06344 0.17553 0.38516]
...
0.0003598675094789 [2.30116e-04 3.59868e-04 6.47473e-05]
13 [ 1.00004 -3.00002 2.00007]
0.0002510632800638 [4.21600e-05 1.07871e-04 2.51063e-04]
14 [ 0.99999 -2.99997 2.00002]
5.393476706316e-05 [5.12763e-05 5.39348e-05 5.05593e-05]
Metodo de Jacobi
numero de condición: 1.9508402675105447
X: [ 0.99999 -2.99997 2.00002]
errado: 5.3934767063168465e-05
iteraciones: 144. Algoritmo en Python
Para el algoritmo se usa el sistema de ecuaciones en su forma matricial A.X=B. Se asume que se ha aplicado el pivoteo parcial por filas y matriz aumentada (luego se añade como una función). Se asegura que A, X y B sean arreglos de números reales.
Las instrucciones inician seleccionando la primera ecuación, primera fila, i=0. El numerador de la expresión se acumula en suma, que inicia con B[i]. Luego se restan todas las combinaciones de A[i,j]*X[j], menos en la diagonal. El nuevo valor de X[i], o Xnuevo, es suma dividido para el coeficiente de la diagonal A[i,i]
i=0 # una ecuacion
suma = B[i] # numerador
for j in range(0,m,1): # un coeficiente
if (i!=j): # excepto diagonal de A
suma = suma-A[i,j]*X[j]
Xnuevo[i] = suma/A[i,i]
Se repite el proceso para las siguientes ecuaciones, cambiando el valor de 'i'
En cada iteración se determinan las diferencias |X[i+1]-X[i]| como los errores de iteración para cada variable.
diferencia = abs(Xnuevo-X)
errado = np.max(diferencia)
Como las diferencias es un vector y la tolerancia solo un número (escalar), se simplifica usando la norma np.max(diferencia).
De ser necesario, se repite la iteración, actualizando el vector X.
# Método de Jacobi
import numpy as np
# INGRESO
# NOTA: Para AB, se ha usado pivoteo parcial por filas
A = [[ 7, 1, 1],
[-3, 7, -1],
[-2, 5, 9]]
B = [6, -26, 1]
X0 = [0,0,0]
tolera = 0.0001
iteramax = 100
# PROCEDIMIENTO
# Matrices como arreglo, numeros reales
A = np.array(A,dtype=float)
B = np.array(B,dtype=float)
X0 = np.array(X0,dtype=float)
tamano = np.shape(A) # tamaño A
n = tamano[0]
m = tamano[1]
# valores iniciales
diferencia = np.ones(n, dtype=float)
errado = 2*tolera # np.max(diferencia)
itera = 0
X = np.copy(X0)
Xnuevo = np.copy(X0)
while errado>tolera and itera<=iteramax:
for i in range(0,n,1): # una ecuacion
suma = B[i] # numerador
for j in range(0,m,1): # un coeficiente
if (i!=j): # excepto diagonal de A
suma = suma-A[i,j]*X[j]
Xnuevo[i] = suma/A[i,i]
diferencia = abs(Xnuevo-X)
errado = np.max(diferencia)
print(itera, X)
print(' ',errado,diferencia)
X = np.copy(Xnuevo) # actualiza X
itera = itera + 1
if (itera>iteramax): # No converge
X = np.nan
print('No converge,iteramax superado')
# numero de condicion
ncond = np.linalg.cond(A)
# SALIDA
print('Metodo de Jacobi')
print('numero de condición:', ncond)
print('X: ',X)
print('errado:',errado)
print('iteraciones:', itera)
5. Algoritmo como función en Python
Las instrucciones del algoritmo, empaquetadas como función, añaden algunas facilidades para el desarrollo del ejercicio. Antes de usar el método de Jacobi, se aplica el pivoteo parcial por filas para mejorar de ser posible la convergencia.
En la función jacobi() se añaden las instrucciones para mostrar los resultados parciales en cada iteración. Se añade una tabla para revisión de valores al final del algoritmo, o para realizar la gráfica de errores por iteración.
Si el sistema es de 3x3 se puede añadir una gráfica de los resultados por iteración, mostrando la trayectoria descrita por los resultados parciales en 3D.
# 1Eva_IIT2011_T2 Sistema de Ecuaciones, diagonal dominante
# Metodos Numericos
import numpy as np
def jacobi(A,B,X0, tolera, iteramax=100,
vertabla=False, precision=4):
''' Método de Jacobi, tolerancia, vector inicial X0
para mostrar iteraciones y tabla: vertabla=True
'''
# Matrices como arreglo, numeros reales
A = np.array(A,dtype=float)
B = np.array(B,dtype=float)
X0 = np.array(X0,dtype=float)
tamano = np.shape(A) # tamaño A
n = tamano[0]
m = tamano[1]
# valores iniciales
diferencia = 2*tolera*np.ones(n, dtype=float)
errado = 2*tolera # np.max(diferencia)
tabla = [np.copy(X0)] # tabla de iteraciones
tabla = np.concatenate((tabla,[[np.nan]]),
axis=1) # errado
if vertabla==True:
print('Iteraciones Jacobi')
print('itera,[X]')
print(' ,errado,|diferencia|')
print(0,X0)
print(' ',np.nan)
np.set_printoptions(precision)
itera = 0 # Jacobi
X = np.copy(X0)
Xnuevo = np.copy(X0)
while errado>tolera and itera<=iteramax:
for i in range(0,n,1): # una ecuacion
suma = B[i]
for j in range(0,m,1): # un coeficiente
if (i!=j): # excepto diagonal de A
suma = suma-A[i,j]*X[j]
Xnuevo[i] = suma/A[i,i]
diferencia = abs(Xnuevo-X)
errado = np.max(diferencia)
X = np.copy(Xnuevo)
Xfila = np.concatenate((Xnuevo,[errado]),axis=0)
tabla = np.concatenate((tabla,[Xfila]),axis=0)
itera = itera + 1
if vertabla==True:
print(itera, Xnuevo)
print(' ',errado,diferencia)
# No converge
if (itera>iteramax):
X = np.nan
print('No converge,iteramax superado')
if vertabla==True:
X = [X,tabla]
return(X)
def pivoteafila(A,B,vertabla=False):
'''
Pivotea parcial por filas, entrega matriz aumentada AB
Si hay ceros en diagonal es matriz singular,
Tarea: Revisar si diagonal tiene ceros
'''
A = np.array(A,dtype=float)
B = np.array(B,dtype=float)
# Matriz aumentada
nB = len(np.shape(B))
if nB == 1:
B = np.transpose([B])
AB = np.concatenate((A,B),axis=1)
if vertabla==True:
print('Matriz aumentada')
print(AB)
print('Pivoteo parcial:')
# Pivoteo por filas AB
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
# Para cada fila en AB
pivoteado = 0
for i in range(0,n-1,1):
# columna desde diagonal i en adelante
columna = np.abs(AB[i:,i])
dondemax = np.argmax(columna)
# dondemax no es en diagonal
if (dondemax != 0):
# intercambia filas
temporal = np.copy(AB[i,:])
AB[i,:] = AB[dondemax+i,:]
AB[dondemax+i,:] = temporal
pivoteado = pivoteado + 1
if vertabla==True:
print(' ',pivoteado, 'intercambiar filas: ',i,'y', dondemax+i)
if vertabla==True:
if pivoteado==0:
print(' Pivoteo por filas NO requerido') else: print('AB') print(AB) return(AB)
# PROGRAMA Búsqueda de solucion --------
# INGRESO
A = [[ 7, 1, 1],
[-3, 7, -1],
[-2, 5, 9]]
B = [6, -26, 1]
X0 = [0,0,0]
tolera = 0.0001
iteramax = 100
verdecimal = 5
# PROCEDIMIENTO
AB = pivoteafila(A,B,vertabla=True)
n,m = np.shape(AB)
A = AB[:,:n] # separa en A y B
B = AB[:,n]
[X, tabla] = jacobi(A,B,X0,tolera,iteramax,
vertabla=True,
precision=verdecimal)
n_itera = len(tabla)-1
errado = tabla[-1,-1]
# numero de condicion
ncond = np.linalg.cond(A)
# SALIDA
print('Metodo de Jacobi')
print('numero de condición:', ncond)
print('X: ',X)
print('errado:',errado)
print('iteraciones:', n_itera)
6. Gráficas de iteraciones
Se muestra la gráfica de errores vs iteraciones.
# GRAFICA de iteraciones errados = tabla[:,n] import matplotlib.pyplot as plt # grafica de error por iteracion fig2D = plt.figure() graf = fig2D.add_subplot(111) graf.plot(errados) graf.set_xlabel('itera') graf.set_ylabel('|error|') graf.set_title('Método de Jacobi: error[itera]') graf.grid() plt.show()
6.1 Gráfica de iteraciones para sistema de 3x3
Si el sistema es de 3x3, se muestra una gráfica de la trayectoria de resultados parciales, semejante a la presentada para la descripción del concepto. Permite observar si la espiral es convergente o no.
# Grafica de iteraciones x,y,z # solo para 3x3 xp = tabla[:,0] yp = tabla[:,1] zp = tabla[:,2] fig3D = plt.figure() graf3D = fig3D.add_subplot(111,projection='3d') graf3D.plot(xp,yp,zp) graf3D.plot(xp[0],yp[0],zp[0],'o',label='X[0]') graf3D.plot(xp[n_itera],yp[n_itera],zp[n_itera], 'o',label='X['+str(n_itera)+']') graf3D.set_xlabel('x') graf3D.set_ylabel('y') graf3D.set_zlabel('z') graf3D.set_title('Método de Jacobi: iteraciones X') # graf3D.view_init(45,45) plt.legend() plt.show()