s1Eva_IIT2007_T2 Aplicar Gauss-Seidel 6×6

Ejercicio: 1Eva_IIT2007_T2 Aplicar Gauss-Seidel 6×6

1. Desarrollo Analítico

– Verificar que la matriz es diagonal dominante

No es necesario realizar el pivoteo por filas, ya la matriz tiene la diagonal dominante.

A = [[7.63, 0.30, 0.15,  0.50, 0.34, 0.84],
     [0.38, 6.40, 0.70,  0.90, 0.29, 0.57],
     [0.83, 0.19, 8.33,  0.82, 0.34, 0.37],
     [0.50, 0.68, 0.86, 10.21, 0.53, 0.70],
     [0.71, 0.30, 0.85,  0.82, 5.95, 0.55],
     [0.43, 0.54, 0.59,  0.66, 0.31, 9.25]]

B = [ -9.44, 25.27, -48.01, 19.76, -23.63, 62.59]

– revisar el número de condición

cond(A) = || A||.||A^-1||

El número de condición no es «muy alto»,  los valores de la diagonal son los mayores en toda la fila, por lo que el sistema converge.

>>> np.linalg.cond(A)
2.04518539662910119

– realizar iteraciones con las expresiones:

$x_0 = \Big(-9.44 -0.3 x_1 -0.15x_2 -0.50 x_3 -0.34 x_4 -0.84 x_5\Big)\frac{1}{7.63}$ $x_1 = \Big( 25.27 -0.38 x_0 -0.70 x_2 -0.90 x_3 -0.29 x_4 -0.57 x_5 \Big) \frac{1}{6.40}$ $x_2 = \Big(-48.01 -0.83x_0 -0.19 x_1 -0.82x_3 -0.34x_4 -0.37x_5 \Big)\frac{1}{8.33}$ $x_3 = \Big(19.76 -0.50 x_0 -0.68 x_1 -0.86 x_2 -0.53 x_4 -0.70x_5 \Big)\frac{1}{10.21}$ $x_4 = \Big( -23.63 - 0.71 x_0 -0.30 x_1 -0.85 x_2 -0.82 x_3 -0.55x_5 \Big)\frac{1}{5.95}$ $x_5 = \Big( 62.59 - 0.43x_0 -0.54x_1 - 0.59 x_2 -0.66 x_3 -0.31 x_4 \Big)\frac{1}{9.25}$

Dado que no se establece en el enunciado el vector inicial, se usará el vector cero. La tolerancia requerida es 10^-5
X0 = [ 0. 0. 0. 0. 0. 0.] = [x₀,x₁,x₂,x₃,x₄,x₅]

iteración 0

$x_0 = \Big(-9.44 -0.3 (0) -0.15(0) -0.50 (0)$ $-0.34(0) -0.84 (0)\Big)\frac{1}{7.63} = -1.23722$ $x_1 = \Big( 25.27 -0.38 (-1.23722) -0.70 (0) -0.90 (0)$ $-0.29 (0) -0.57 (0) \Big) \frac{1}{6.40} = 4.0219$ $x_2 = \Big(-48.01 -0.83 (-1.23722) -0.19 (4.0219) -0.82(0)$ $-0.34(0) -0.37(0) \Big)\frac{1}{8.33} = -5.73196$ $x_3 = \Big( 19.76 -0.50 (-1.23722) -0.68 (4.0219) -0.86 (-5.73196)$ $-0.53 (0) -0.70(0) \Big)\frac{1}{10.21} = 2.21089$ $x_4 = \Big( -23.63 - 0.71 (-1.23722) -0.30 (4.0219) -0.85 (-5.73196)$ $-0.82 (2.21089) -0.55(0) \Big)\frac{1}{5.95} = -3.51242$ $x_5 = \Big( 62.59 - 0.43(-1.23722) -0.54(4.0219) - 0.59 (-5.73196)$ $-0.66 (2.21089) -0.31 (-3.51242) \Big)\frac{1}{9.25} = 6.91478$ $X_1 = [-1.23722, 4.0219, -5.73196, 2.21089, -3.51242, 6.91478 ]$ $diferencia = X1 - [0,0,0,0,0,0]$ $errado = max(|diferencia|) = 6.91478$

iteración 1

$x_0 = \Big(-9.44 -0.3 (4.0219) -0.15(4.0219) -0.50 (2.21089)$ $-0.34 (-3.51242) -0.84(6.91478) \Big)\frac{1}{7.63} = -2.0323$ $x_1 = \Big( 25.27 -0.38 (-2.0323) -0.70 ( -5.73196) -0.90 (2.21089)$ $-0.29 (-3.51242) -0.57 (6.91478) \Big) \frac{1}{6.40} = 3.92844$ $x_2 = \Big(-48.01 -0.83(-2.0323) -0.19 (3.92844) -0.82(2.21089)$ $-0.34(-3.51242) -0.37(6.91478) \Big)\frac{1}{8.33} = -6.03203$ $x_3 = \Big(19.76 -0.50 (-2.0323) -0.68 (3.92844) -0.86 (-6.03203)$ $-0.53 (-3.51242) -0.70(6.91478) \Big)\frac{1}{10.21} = 1.98958$ $x_4 = \Big( -23.63 - 0.71 (-2.0323) -0.30 (3.92844) -0.85 (-6.03203)$ $-0.82 (1.98958) -0.55(6.91478) \Big)\frac{1}{5.95} = -3.97865$ $x_5 = \Big( 62.59 - 0.43(-2.0323) -0.54(3.92844) - 0.59 (-6.03203)$ $-0.66 (1.98958) -0.31 (-3.97865) \Big)\frac{1}{9.25} = 7.00775$ $X_2 = [-2.0323, 3.92844, -6.03203, 1.98958, -3.97865, 7.00775]$ $diferencia = X2 - [-1.23722, 4.0219, -5.73196, 2.21089, -3.51242, 6.91478 ]$ $errado = max(|diferencia|) = 0.79507$

iteración 2

Desarrollar como tarea.

---

2. Algoritmo en Python

La tabla de aproximaciones sucesivas para el vector X es:

Pivoteo por filas NO requerido
Iteraciones Gauss-Seidel
itera,[X],[errado]
0 [0. 0. 0. 0. 0. 0.] [1.]
1 [-1.23722  4.0219  -5.73196  2.21089 -3.51242  6.91478] [6.91478]
2 [-2.0323   3.92844 -6.03203  1.98958 -3.97865  7.00775] [0.79507]
3 [-1.99768  4.00317 -6.00049  1.99808 -4.00082  6.9999 ] [0.07473]
4 [-1.99994  4.00037 -5.99979  2.      -4.00005  6.99996] [0.00281]
5 [-2.00001  3.99998 -6.       2.00001 -4.       7.     ] [0.00038]
6 [-2.  4. -6.  2. -4.  7.] [1.60155e-05]
7 [-2.  4. -6.  2. -4.  7.] [1.74554e-06]
Matriz aumentada y pivoteada:
[[  7.63   0.3    0.15   0.5    0.34   0.84  -9.44]
 [  0.38   6.4    0.7    0.9    0.29   0.57  25.27]
 [  0.83   0.19   8.33   0.82   0.34   0.37 -48.01]
 [  0.5    0.68   0.86  10.21   0.53   0.7   19.76]
 [  0.71   0.3    0.85   0.82   5.95   0.55 -23.63]
 [  0.43   0.54   0.59   0.66   0.31   9.25  62.59]]
Vector Xi: 
[-2.  4. -6.  2. -4.  7.]
>>> 

el gráfico de los iteraciones vs errores es:

La tabla obtiene aplicando la función de Gauss-Seidel, tomando como vector inicial el vector de ceros.

Tarea: X=TX+C

Instrucciones en Python

# 1Eva_IIT2007_T2 Aplicar Gauss-Seidel
# Algoritmo Gauss-Seidel
import numpy as np

def pivoteafila(A,B,vertabla=False):
    '''
    Pivotea parcial por filas
    Si hay ceros en diagonal es matriz singular,
    Tarea: Revisar si diagonal tiene ceros
    '''
    # Matriz aumentada
    nB = len(np.shape(B))
    if nB == 1:
        B = np.transpose([B])
    AB  = np.concatenate((A,B),axis=1)
    M = np.copy(AB)
    
    # Pivoteo por filas AB
    tamano = np.shape(M)
    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(M[i:,i])
        dondemax = np.argmax(columna)
        
        # dondemax no es en diagonal
        if (dondemax != 0):
            # intercambia filas
            temporal = np.copy(M[i,:])
            M[i,:] = M[dondemax+i,:]
            M[dondemax+i,:] = temporal

            pivoteado = pivoteado + 1
            if vertabla==True:
                print(pivoteado, 'intercambiar: ',i,dondemax)
    if vertabla==True and pivoteado==0:
        print('Pivoteo por filas NO requerido')
    return(M)

def gauss_seidel(A,B,tolera,X0, iteramax=100, vertabla=False, precision=5):
    ''' Método de Gauss Seidel, tolerancia, vector inicial X0
        para mostrar iteraciones: vertabla=True
    '''
    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],[errado]')
        np.set_printoptions(precision)
        print(itera, X, np.array([errado]))
    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, np.array([errado]))        
    
    if (itera>iteramax): # No converge
        X = itera
        print('iteramax superado, No converge')
    return(X)

# Programa de prueba #######
# INGRESO
A = np.array([[7.63, 0.30, 0.15,  0.50, 0.34, 0.84],
              [0.38, 6.40, 0.70,  0.90, 0.29, 0.57],
              [0.83, 0.19, 8.33,  0.82, 0.34, 0.37],
              [0.50, 0.68, 0.86, 10.21, 0.53, 0.70],
              [0.71, 0.30, 0.85,  0.82, 5.95, 0.55],
              [0.43, 0.54, 0.59,  0.66, 0.31, 9.25]])

B = np.array([-9.44,25.27,-48.01,19.76,-23.63,62.59])

tolera = 0.00001
X = np.zeros(len(A), dtype=float)

# PROCEDIMIENTO
n = len(A)
A = np.array(A,dtype=float)
B = np.array(B,dtype=float)


AB = pivoteafila(A,B, vertabla=True)
A = AB[:,:n]
B = AB[:,n]
respuesta = gauss_seidel(A,B, tolera, X, vertabla=True)

# SALIDA
print('Matriz aumentada y pivoteada:')
print(AB)
print('Vector Xi: ')
print(respuesta)

En el caso de la norma infinito, para la matriz A, se puede usar el algoritmo desarrollado en clase.
Como valor para verificar su algoritmo, se obtuvo:

>>> np.linalg.norm(A, np.inf)
13.479999999999999

Tarea: incluir la norma infinito para T