[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..
1. El método iterativo de Jacobi
Referencia: Burden 7.3 p334, Rodríguez 5.1 p154, Chapra 11.2.1p312
La analogía presentadas entre la «norma como distancia 3D» y el «error de acoplamiento de aeronaves», 
permite considerar desde un punto de partida o inicial las aproximaciones o iteraciones sucesivas hacia una solución del sistema de ecuaciones.
Las iteraciones pueden ser convergentes o no.
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.
Para describir el método iterativo de Gauss-Seidel, se usa un sistema de 3 incógnitas y 3 ecuaciones, diagonalmente dominante.
\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}Para facilitar la escritura del algoritmo, note el uso de índices ajustado a la descripción de arreglos en Python para la primera fila i=0 y primera columna j=0.
Semejante a despejar una variable de la ecuación para representar un plano, se plantea despejar una variable de cada ecuación. Se obtiene así los valores de cada xi, por cada por cada fila i:
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}}Observe que el patrón de cada fórmula para determinar x[i], tiene la forma:
x_i = \bigg(b_i - \sum_{j=0, j\neq i}^n A_{i,j}X_j\bigg) \frac{1}{A_{ii}}La parte de la sumatoria se realiza para cada término de A[i,j] en la fila i, excepto para el término de la diagonal A[i,i].
Si se tiene conocimiento del problema planteado y se puede «intuir» o «suponer» una solución para el vector X. Por ejemplo, iniciando con el vector cero, es posible calcular un nuevo vector X usando las ecuaciones para cada X[i] encontradas.
X = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}Con cada nuevo valor actualiza el el vector X, en Xnuevo, se calcula el vector diferencias entre X y el vector Xnuevo .
El error a prestar la atención es al mayor valor de las diferencias; se toma como condición para repetir la evaluación de cada vector.
nuevo = [ 0, 0, 0]
X = [x0, x1, x2]
diferencia = [|0 - x0|, |0 - x1|, |0 - x2|]
errado = max(diferencia)
Se observa los resultados de errado para cada iteración, relacionados con la convergencia. Si luego de «muchas» iteraciones se encuentra que (errado>tolera), se detiene el proceso.
[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..
2. Ejercicio
Referencia: Steven C. Chapra. Numerical Methods 7th Edition, Ejercicio 12.39 p339; 1Eva_IT2018_T3 Temperatura en nodos de placa.
La temperatura en los nodos de la malla de una placa se puede calcular con el promedio de las temperaturas de los 4 nodos vecinos de la izquierda, derecha, arriba y abajo. 
Una placa cuadrada de 3 m de lado tiene la temperatura en los nodos de los bordes como se indica en la figura,
a) Plantee el sistema de ecuaciones y resuelva con el método de Jacobi para encontrar los valores de a, b, c, d.
[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..
3. Desarrollo Analítico
Solución Propuesta: s1Eva_IT2018_T3 Temperatura en nodos de placa
a) Plantear el sistema de ecuaciones. Usando el promedio para cada nodo interior:
que reordenando se convierte en:
4a=150+c+b 4b=a+80+d 4c=a+160+d 4d=b+c+90simplificando:
4a-b-c= 150 a-4b+d = -80 a-4c+d = -160 b+c-4d = -90que a forma matricial se convierte en:
A = [[ 4, -1, -1, 0],
[ 1, -4, 0, 1],
[ 1, 0, -4, 1],
[ 0, 1, 1,-4]]
B = [150, -80,-160,-90]
Observación: la matriz A ya es diagonal dominante, no requiere pivotear por filas. Se aumentó el punto decimal a los valores de la matriz A y el vector B para que sean considerados como números reales.
El número de condición es: np.linalg.cond(A) = 3.0
que es cercano a 1 en un orden de magnitud, por lo que la solución matricial es «estable» y los cambios en los coeficientes afectan proporcionalmente a los resultados. Se puede aplicar métodos iterativos sin mayores inconvenientes.
b y c) método de Jacobi para sistemas de ecuaciones, con vector inicial
X_0 = [60,40,70,50]reemplazando los valores iniciales en cada ecuación sin cambios.
iteración 0
a=\frac{50+70+100+40}{4} = 65
iteración 1
a=\frac{50+67.5+100+47.5}{4} = 66.25
iteración 2
a=\frac{50+68.75+100+48.75}{4} = 66.875
siguiendo las iteraciones se debería llegar a:
>>> np.linalg.solve(A,B) array([67.5, 50. , 70. , 52.5])
[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]
..
4. Algoritmo en Python
Tarea: Realice el algoritmo de Jacobi en Python, incluyendo el pivoteo por filas de la matriz.
Resultados con el algoritmo
Matriz aumentada [[ 4. -1. -1. 0. 150.] [ 1. -4. 0. 1. -80.] [ 1. 0. -4. 1. -160.] [ 0. 1. 1. -4. -90.]] Pivoteo parcial: Pivoteo por filas NO requerido Iteraciones Jacobi itera,[X],errado 0 [60. 40. 70. 50.] 1.0 1 [65. 47.5 67.5 50. ] 7.5 2 [66.25 48.75 68.75 51.25] 1.25 3 [66.875 49.375 69.375 51.875] 0.625 4 [67.1875 49.6875 69.6875 52.1875] 0.3125 5 [67.34375 49.84375 69.84375 52.34375] 0.15625 6 [67.42188 49.92188 69.92188 52.42188] 0.078125 7 [67.46094 49.96094 69.96094 52.46094] 0.0390625 8 [67.48047 49.98047 69.98047 52.48047] 0.01953125 9 [67.49023 49.99023 69.99023 52.49023] 0.009765625 10 [67.49512 49.99512 69.99512 52.49512] 0.0048828125 11 [67.49756 49.99756 69.99756 52.49756] 0.00244140625 12 [67.49878 49.99878 69.99878 52.49878] 0.001220703125 13 [67.49939 49.99939 69.99939 52.49939] 0.0006103515625 14 [67.49969 49.99969 69.99969 52.49969] 0.00030517578125 15 [67.49985 49.99985 69.99985 52.49985] 0.000152587890625 16 [67.49992 49.99992 69.99992 52.49992] 7.62939453125e-05 numero de condición: 3.000000000000001 respuesta X: [67.49992 49.99992 69.99992 52.49992] iterado, incluyendo X0: 17
Instrucciones en Python
# 1Eva_2024PAOI_T2 Temperatura en nodos de placa cuadrada # Método de Jacobi 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: vertabla=True ''' A = np.array(A,dtype=float) B = np.array(B,dtype=float) X0 = np.array(X0,dtype=float) tamano = np.shape(A) n = tamano[0] m = tamano[1] diferencia = np.ones(n, dtype=float) errado = np.max(diferencia) X = np.copy(X0) xnuevo = np.copy(X0) tabla = [np.copy(X0)] itera = 0 if vertabla==True: print('Iteraciones Jacobi') print('itera,[X],errado') print(itera, xnuevo, errado) np.set_printoptions(precision) while not(errado<=tolera or itera>iteramax): for i in range(0,n,1): nuevo = B[i] for j in range(0,m,1): if (i!=j): # excepto diagonal de A nuevo = nuevo-A[i,j]*X[j] nuevo = nuevo/A[i,i] xnuevo[i] = nuevo diferencia = np.abs(xnuevo-X) errado = np.max(diferencia) X = np.copy(xnuevo) tabla = np.concatenate((tabla,[X]),axis = 0) itera = itera + 1 if vertabla==True: print(itera, xnuevo, errado) # No converge if (itera>iteramax): X=itera print('iteramax superado, No converge') return(X,tabla) 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 # INGRESO A = [[ 4, -1, -1, 0], [ 1, -4, 0, 1], [ 1, 0, -4, 1], [ 0, 1, 1, -4]] B = [150, -80,-160,-90] X0 = [60,40,70,50] tolera = 0.0001 iteramax = 100 verdecimal = 5 # PROCEDIMIENTO AB = pivoteafila(A,B,vertabla=True) # separa matriz aumentada en A y B [n,m] = np.shape(AB) A = AB[:,:m-1] B = AB[:,m-1] [X, puntos] = jacobi(A,B,X0,tolera,vertabla=True,precision=verdecimal) iterado = len(puntos) # numero de condicion ncond = np.linalg.cond(A) # SALIDA print('numero de condición:', ncond) print('respuesta X: ') print(X) print('iterado, incluyendo X0: ', iterado)
La gráfica de puntos por iteraciones en 3D mostrada al inicio se desarrolla en la propuesta de solución al ejercicio:
s1Eva_IIT2011_T2 Sistema de Ecuaciones, diagonal dominante
[ Método de Jacobi ] [ Ejercicio ] [ Analítico ] [ Algoritmo ]