Página web del curso

Contenidos

Método de Gauss-Seidel

Sea el sistema $A\mathbf{x}=\mathbf{b}$ donde

$$ A=\left(\begin{array}{cccc} 20 & 1 & 0 & 1\\ 1 & 20 & 3 & 1\\ 0 & 3 & 20 & 1\\ 1 & 0 & 1 & 20 \end{array}\right)\quad\mathbf{b}=\left(\begin{array}{c} 10\\ 7\\ 4\\ 6 \end{array}\right) $$

Resolver utilizando el método de Gauss-Seidel.

Utilizar como condición de parada $\left\Vert \mathbf{x}^{(k)}-\mathbf{x}^{(k-1)}\right\Vert _{\infty}\lt0.01$

Vamos a resolverlo usando python para hacer las cuentas

In [1]:
import numpy as np
np.set_printoptions(precision=3)

Escribimos los datos, que son la matriz de coeficientes y la de términos independientes

In [2]:
A = np.array([[20., 1, 0, 1],[1, 20, 3, 1],[0, 3, 20, 1],[1, 0, 1, 20]])
b = np.array([10., 7, 4, 6])

E inicializamos los dos vectores donde almacenaremos las dos últimas soluciones aproximadas

In [3]:
x0 = np.zeros(4)
x1 = np.zeros(4)

Si tenemos un sistema $4\times 4$ $$ \begin{array}{cll} a_{11}x_{1}+a_{12}x_{2}+a_{13}x_{3}+a_{14}x_{4} & = & b_{1}\\ a_{21}x_{1}+a_{22}x_{2}+a_{23}x_{3}+a_{24}x_{4} & = & b_{2}\\ a_{31}x_{1}+a_{32}x_{2}+a_{33}x_{3}+a_{34}x_{4} & = & b_{3}\\ a_{41}x_{1}+a_{42}x_{2}+a_{43}x_{3}+a_{44}x_{4} & = & b_{4} \end{array} $$

y despejamos la primera incógnita en la primera ecuación, la segunda incógnita en la segunda y así sucesivamente

$$ \begin{array}{cll} x_{1}& = &\dfrac{b_{1}-a_{12}x_{2}-a_{13}x_{3}-a_{14}x_{4}}{a_{11}} \\ x_{2} & = & \dfrac{b_{2}-a_{21}x_{1}-a_{23}x_{3}-a_{24}x_{4}}{a_{22}}\\ x_{3} & = & \dfrac{b_{3}-a_{31}x_{1}-a_{32}x_{2}-a_{34}x_{4}}{a_{33}}\\ x_{4} & = & \dfrac{b_{4}-a_{41}x_{1}-a_{42}x_{2}-a_{43}x_{3}}{a_{44}} \end{array} $$

Y obtenemos los valores de la primera iteración a partir de los valores iniciales, pero conforme vamos teniendo resultados nuevos de una variable, los actualizamos

$$ \begin{array}{cll} x_{1}^{\color{red}{(1)}} & = & \dfrac{b_{1}-a_{12}x_{2}^{(0)}-a_{13}x_{3}^{(0)}-a_{14}x_{4}^{(0)}}{a_{11}} \\ x_{2}^{\color{red}{(1)}} & = & \dfrac{b_{2}-a_{21}x_{1}^{\color{red}{(1)}}-a_{23}x_{3}^{(0)}-a_{24}x_{4}^{(0)}}{a_{22}}\\ x_{3}^{\color{red}{(1)}} & = & \dfrac{b_{3}-a_{31}x_{1}^{\color{red}{(1)}}-a_{32}x_{2}^{\color{red}{(1)}}-a_{34}x_{4}^{(0)}}{a_{33}}\\ x_{4}^{\color{red}{(1)}}& = & \dfrac{b_{4}-a_{41}x_{1}^{\color{red}{(1)}}-a_{42}x_{2}^{\color{red}{(1)}}-a_{43}x_{3}^{\color{red}{(1)}}}{a_{44}} \end{array} $$

Algoritmo

El algoritmo que hemos aplicado se expresa formalmente

  • Elegir una aproximación inicial $\mathbf{x}^{(0)}$
  • Para $k=1,2,\ldots,$MaxIter
    • Para $i=1,2,\ldots,n$, calcular $$x_{i}^{\color{red}{(k)}}=\dfrac{1}{a_{ii}}\left(b_{i}-\sum_{j=1}^{i-1}a_{ij}x_{j}^{\color{red}{(k)}}-\sum_{j=i+1}^{n}a_{ij}x_{j}^{(k-1)}\right)$$
    • Si se cumple el criterio de parada, tomar $\mathbf{x}^{(k)}$ como aproximación a la solución.

Iteración 1

A partir de $\mathbf{x^{(0)}}$ calculamos $\mathbf{x^{(1)}}$, pero para ir actualizando los valores de una coordenada en cuanto lo hemos actualizado, una estrategia sería usar un sólo vector durante la iteración. No obstante, necesitamos guardar la iteración anterior para comprobar el criterio de parada

In [4]:
k = 1
xant = np.copy(x0)
x0[0] = (b[0] - A[0,1] * x0[1] - A[0,2] * x0[2] - A[0,3] * x0[3]) / A[0,0]
x0[1] = (b[1] - A[1,0] * x0[0] - A[1,2] * x0[2] - A[1,3] * x0[3]) / A[1,1]
x0[2] = (b[2] - A[2,0] * x0[0] - A[2,1] * x0[1] - A[2,3] * x0[3]) / A[2,2]
x0[3] = (b[3] - A[3,0] * x0[0] - A[3,1] * x0[1] - A[3,2] * x0[2]) / A[3,3]

Imprimimos la nueva iteración, el vector diferencia y la norma de dicho vector diferencia, que nos va a valer para controlar el error

In [5]:
print ('x0                = ', x0)
print ('x0-xant           = ', x0-xant)
print ('max(abs(x0-xant)) = ', max(abs(x0-xant)))

x0                =  [0.5   0.325 0.151 0.267]
x0-xant           =  [0.5   0.325 0.151 0.267]
max(abs(x0-xant)) =  0.5

Como la norma es mayor que la $0.01$, seguimos iterando

Iteración 2

A partir de $\mathbf{x^{(1)}}$ calculamos $\mathbf{x^{(2)}}$

In [6]:
k = 2
xant = np.copy(x0)
x0[0] = (b[0] - A[0,1] * x0[1] - A[0,2] * x0[2] - A[0,3] * x0[3]) / A[0,0]
x0[1] = (b[1] - A[1,0] * x0[0] - A[1,2] * x0[2] - A[1,3] * x0[3]) / A[1,1]
x0[2] = (b[2] - A[2,0] * x0[0] - A[2,1] * x0[1] - A[2,3] * x0[3]) / A[2,2]
x0[3] = (b[3] - A[3,0] * x0[0] - A[3,1] * x0[1] - A[3,2] * x0[2]) / A[3,3]

Imprimimos la nueva iteración, el vector diferencia y la norma de dicho vector diferencia, que nos va a valer para controlar el error

In [7]:
print ('x0                = ', x0)
print ('x0-xant           = ', x0-xant)
print ('max(abs(x0-xant)) = ', max(abs(x0-xant)))

x0                =  [0.47  0.29  0.143 0.269]
x0-xant           =  [-0.03  -0.035 -0.008  0.002]
max(abs(x0-xant)) =  0.03457828125000001

Como la norma es mayor que la $0.01$, seguimos iterando

Iteración 3

A partir de $\mathbf{x^{(2)}}$ calculamos $\mathbf{x^{(3)}}$ (aunque los guardamos respectivamente en x0 y x1)

In [8]:
k = 3
xant = np.copy(x0)
x0[0] = (b[0] - A[0,1] * x0[1] - A[0,2] * x0[2] - A[0,3] * x0[3]) / A[0,0]
x0[1] = (b[1] - A[1,0] * x0[0] - A[1,2] * x0[2] - A[1,3] * x0[3]) / A[1,1]
x0[2] = (b[2] - A[2,0] * x0[0] - A[2,1] * x0[1] - A[2,3] * x0[3]) / A[2,2]
x0[3] = (b[3] - A[3,0] * x0[0] - A[3,1] * x0[1] - A[3,2] * x0[2]) / A[3,3]

Imprimimos el nuevo resultado y la norma del vector diferencia

In [9]:
print ('x0                = ', x0)
print ('x0-xant           = ', x0-xant)
print ('max(abs(x0-xant)) = ', max(abs(x0-xant)))

x0                =  [0.472 0.291 0.143 0.269]
x0-xant           =  [ 1.634e-03  1.052e-03 -2.522e-04 -6.911e-05]
max(abs(x0-xant)) =  0.0016343965429687635

Ahora ya se verifica que la norma del vector diferencia es menor que $0.01$ y paramos. Nuestra solución aproximada es

In [10]:
print ('x0 = ', x0)

x0 =  [0.472 0.291 0.143 0.269]

La solución exacta es

In [11]:
x = np.linalg.solve(A,b);
print(x)

[0.472 0.292 0.143 0.269]

Y el error

In [12]:
print ('error = ', max(abs(x0-x)))

error =  5.146953503498697e-05