Publicado en por Edison Del Rosario

3.6.1 Método de Gauss-Seidel con Python

Referencia: Ejemplo 01 Chapra Ejemplo 11.3  p311 pdf335

\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 fué implementado el procedimiento.

Planteamiento

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}

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

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 = [2.61, -2.79, 7.00]

diferencias = [2.61, -2.79, 7.00] – [0,0,0]
diferencias = [2.61, -2.79, 7.00]
error = ||diferencias|| , se usa la norma de max(diferencias)
error = 7.00

itera = 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 = [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 = ||diferencias||, se usa la norma de max(diferencias)

error = 0.38

itera = 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.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 = ||diferencias||, se usa la norma de max(diferencias)

error = 0.01

El 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 raiz.

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.

Ejemplo 01 Chapra Ejemplo 11.3 p311 pdf335

\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 fué 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.