[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfica ]
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,
[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfica ]
2. Ejercicio
Referencia: 1Eva_IIT2011_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.
[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfica ]
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.]\begin{cases} 7x+y+z=6\\-3x+7y-z=-26\\-2x+5y+9z=1\end{cases}
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]\text{diferencias}=[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]\text{diferencias}= |[0.51474, -0.38322, 2.25397]| \text{errado}=max |X_2-X_1|= 2.2539
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]\text{diferencias}=|[-0.3767, 0.5426, 0.0985]| \text{errado}=max |X_2-X_1|= 0.5425
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: 14
[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfica ]
..
4. 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)
[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfica ]
..
5. Algoritmo en Python como función
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 3×3 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') 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)
[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfica ]
..
5.1 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()
[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfica ]
5.2 Gráfica de iteraciones para sistema de 3×3
Si el sistema es de 3×3, 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()
[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Gráfica ]