[ Gauss-Seidel ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]
..
Referencia: Chapra 11.2 p310, Burden 7.3 p337, Rodríguez 5.2 p162
El método de Gauss-Seidel realiza operaciones semejantes al método de Jacobi. 
El método de Gauss-Sidel también usa el vector inicial X0, la diferencia consiste en que la actualización del vector X en cada iteración se realiza por cada nuevo valor del vector que se ha calculado. Por lo que las iteraciones llegan más rápido al punto cuando el método converge.
[ Gauss-Seidel ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]
..
2. Ejercicio
Referencia: Chapra ejemplo 11.3 p311
El ejemplo de referencia, ya presenta una matriz pivoteada por filas, por lo que no fue necesario desarrollar esa parte en el ejercicio.
\begin{cases} 3 x_0 -0.1 x_1 -0.2 x_2 = 7.85 \\ 0.1 x_0 +7x_1 -0.3 x_2 = -19.3 \\ 0.3 x_0 -0.2 x_1 +10 x_2 = 71.4 \end{cases}Encuentre la solución usando el método de Gauss-Seidel y tolerancia de 0.00001
[ Gauss-Seidel ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]
..
3. Desarrollo analítico
Para el desarrollo con lápiz y papel, se despeja una variable de cada ecuación, se empieza con la primera expresión para obtener x0
3 x_0 -0.1 x_1 -0.2 x_2 = 7.85 3 x_0 = 7.85 + 0.1 x_1 + 0.2 x_2 x_0 = \frac{7.85 + 0.1 x_1 + 0.2 x_2}{3}con la segunda ecuación se obtiene x1
0.1 x_0 +7x_1 -0.3 x_2 = -19.3 7x_1 = -19.3 - 0.1 x_0 +0.3 x_2 x_1 = \frac{ -19.3 - 0.1 x_0 +0.3 x_2}{7}usando la tercera ecuación se obtiene x2
0.3 x_0 - 0.2 x_1 + 10 x_2 = 71.4 10 x_2 = 71.4 - 0.3 x_0 + 0.2 x_1 x_2 = \frac{71.4 - 0.3 x_0 + 0.2 x_1}{10}Las ecuaciones listas para usar en el algoritmo son:
x_0 = \frac{7.85 + 0.1 x_1 + 0.2 x_2}{3} x_1 = \frac{ -19.3 - 0.1 x_0 +0.3 x_2}{7} x_2 = \frac{71.4 - 0.3 x_0 + 0.2 x_1}{10}Vector inicial
Como el vector inicial no se indica en el enunciado, se considera usar el vector de ceros para iniciar las iteraciones.
X = [0,0,0]
con tolerancia de 0.00001
Iteraciones
Con las ecuaciones obtenidas en el planteamiento, se desarrolla usando el vector inicial presentado, hasta que el |error|<tolera.
El método de Gauss -Seidel actualiza el vector de respuestas Xi+1 luego de realizar cada cálculo. es decir, aprovechas las aproximaciones de cada iteración en el momento que se realizan. Observe la diferencia con el método de Jacobi que espera a terminar con la iteración para actualizar Xi+1
itera = 0
x_0 = \frac{7.85 + 0.1 (0) + 0.2 (0)}{3} = 2.61 x_1 = \frac{ -19.3 - 0.1 (2.61) +0.3 (0)}{7} = -2.79 x_2 = \frac{71.4 - 0.3 (2.61) + 0.2 (-2.79)}{10} = 7.00 X_1 = [2.61, -2.79, 7.00] diferencias = [2.61, -2.79, 7.00] - [0,0,0] diferencias = [2.61, -2.79, 7.00] error = max |diferencias|= 7.00itera = 1
X = [2.61, -2.79, 7.00]
x_0 = \frac{7.85 + 0.1 (-2.79) + 0.2 (7.00)}{3} = 2.99 x_1 = \frac{ -19.3 - 0.1 (2.99) +0.3 (7.00)}{7} = -2.49 x_2 = \frac{71.4 - 0.3 (2.99) + 0.2 (-2.49)}{10} = 7.00 X_1 = [2.99, -2.49, 7.00] diferencias = [2.99, -2.49, 7.00] - [2.61, -2.79, 7.00] diferencias = [0.38, 0.3 , 0. ] error = max |diferencias| = 0.38itera = 2
X = [2.99, -2.49, 7.00]
x_0 = \frac{7.85 + 0.1 (-2.49) + 0.2 (7.00)}{3} = 3.00 x_1 = \frac{ -19.3 - 0.1 (3.00) +0.3 (7.00)}{7} = -2.5 x_2 = \frac{71.4 - 0.3 (3.00) + 0.2 (-2.5)}{10} = 7.00 X_3 = [3.00, -2.50, 7.00] diferencias = [3.00, -2.50, 7.00] - [2.99, -2.49, 7.00] diferencias = [ 0.01, -0.01, 0. ] error = max |diferencias|= 0.01El error aún es mayor que tolera, por lo que es necesario continuar con las iteraciones.
Observaciones
El error disminuye en cada iteración, por lo que el método converge hacia la raíz.
El ejercicio por tener solo 3 incógnitas, la solución se puede interpretar como la intersección de planos en el espacio.
[ Gauss-Seidel ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]
..
4. Algoritmo con Python
Con la descripción dada para el método de Gauss-Seidel, se muestra una forma básica de implementar el algoritmo.
Referencia: Chapra ejemplo 11.3 p311
\begin{cases} 3 x_0 -0.1 x_1 -0.2 x_2 = 7.85 \\ 0.1 x_0 +7x_1 -0.3 x_2 = -19.3 \\ 0.3 x_0 -0.2 x_1 +10 x_2 = 71.4 \end{cases}El ejemplo de referencia, ya presenta una matriz pivoteada por filas, por lo que no fue implementado el procedimiento. Para generalizar el algoritmo, incluya como tarea aumentar el procedimiento de pivoteo por filas.
# Método de Gauss-Seidel # solución de sistemas de ecuaciones # por métodos iterativos import numpy as np # INGRESO A = np.array([[3. , -0.1, -0.2], [0.1, 7 , -0.3], [0.3, -0.2, 10 ]]) B = np.array([7.85,-19.3,71.4]) X0 = np.array([0.,0.,0.]) tolera = 0.00001 iteramax = 100 # PROCEDIMIENTO # Gauss-Seidel tamano = np.shape(A) n = tamano[0] m = tamano[1] # valores iniciales X = np.copy(X0) diferencia = np.ones(n, dtype=float) errado = 2*tolera itera = 0 while not(errado<=tolera or itera>iteramax): # por fila for i in range(0,n,1): # por columna suma = 0 for j in range(0,m,1): # excepto diagonal de A if (i!=j): suma = suma-A[i,j]*X[j] nuevo = (B[i]+suma)/A[i,i] diferencia[i] = np.abs(nuevo-X[i]) X[i] = nuevo errado = np.max(diferencia) itera = itera + 1 # Respuesta X en columna X = np.transpose([X]) # revisa si NO converge if (itera>iteramax): X=0 # revisa respuesta verifica = np.dot(A,X) # SALIDA print('respuesta X: ') print(X) print('verificar A.X=B: ') print(verifica)
que da como resultado:
respuesta X: [[ 3. ] [-2.5] [ 7. ]] verificar A.X=B: [[ 7.84999999] [-19.3 ] [ 71.4 ]] >>>
que es la respuesta del problema obtenida durante el desarrollo del ejemplo con otros métodos.
Tarea
Completar el algoritmo para usar una matriz diagonal dominante, usando al menos el pivoteo parcial por filas.
[ Gauss-Seidel ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]
..
5. Algoritmo como función
Se agrupa las instrucciones como función. Recuerde que la matriz AB tiene que pivotearse por filas antes de usar en el algoritmo. Se obtienen los siguientes resultados:
Matriz aumentada
[[ 3. -0.1 -0.2 7.85]
[ 0.1 7. -0.3 -19.3 ]
[ 0.3 -0.2 10. 71.4 ]]
Pivoteo parcial:
Pivoteo por filas NO requerido
Iteraciones Gauss-Seidel
itera,[X]
diferencia,errado
0 [0. 0. 0.] 2e-05
1 [ 2.6166667 -2.7945238 7.0056095]
[2.6166667 2.7945238 7.0056095] 7.005609523809525
2 [ 2.9905565 -2.4996247 7.0002908]
[0.3738898 0.2948991 0.0053187] 0.3738898412698415
3 [ 3.0000319 -2.499988 6.9999993]
[0.0094754 0.0003633 0.0002915] 0.00947538997430053
4 [ 3.0000004 -2.5 7. ]
[3.1545442e-05 1.2043402e-05 7.0549522e-07] 3.154544153582961e-05
5 [ 3. -2.5 7. ]
[3.5441370e-07 3.5298562e-08 1.1338384e-08] 3.544137046063156e-07
numero de condición: 3.335707415629964
respuesta con Gauss-Seidel
[ 3. -2.5 7. ]
>>>
Instrucciones en Python como función
Se ha incorporado la función de pivoteo por filas.
# Algoritmo Gauss-Seidel # solución de matrices # métodos iterativos # Referencia: Chapra 11.2, p.310, # Rodriguez 5.2 p.162 import numpy as np def gauss_seidel(A,B,X0,tolera, iteramax=100, vertabla=False, precision=4): ''' Método de Gauss Seidel, 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 = 2*tolera*np.ones(n, dtype=float) errado = np.max(diferencia) X = np.copy(X0) itera = 0 if vertabla==True: print('Iteraciones Gauss-Seidel') print('itera,[X]') print(' diferencia,errado') print(itera, X, errado) np.set_printoptions(precision) while (errado>tolera and itera<iteramax): for i in range(0,n,1): xi = B[i] for j in range(0,m,1): if (i!=j): xi = xi-A[i,j]*X[j] xi = xi/A[i,i] diferencia[i] = np.abs(xi-X[i]) X[i] = xi errado = np.max(diferencia) itera = itera + 1 if vertabla==True: print(itera, X) print(' ', diferencia, errado) # No converge if (itera>iteramax): X=itera print('iteramax superado, No converge') 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 de prueba ------------ # INGRESO A = [[3. , -0.1, -0.2], [0.1, 7 , -0.3], [0.3, -0.2, 10 ]] B = [7.85,-19.3,71.4] X0 = [0.,0.,0.] tolera = 0.00001 iteramax = 100 verdecimal = 7 # PROCEDIMIENTO # numero de condicion ncond = np.linalg.cond(A) 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] respuesta = gauss_seidel(A,B,X0, tolera, vertabla=True, precision=verdecimal) # SALIDA print('numero de condición:', ncond) print('respuesta con Gauss-Seidel') print(respuesta)
[ Gauss-Seidel ] [ Ejercicio ] [ Analítico ] [ Algoritmo ] [ Función ]