Contenidos
Método de Jacobi¶
Sea el sistema $A\mathbf{x}=\mathbf{b}$ donde
Resolver utilizando el método de Jacobi.
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
import numpy as np
np.set_printoptions(precision=3)
Escribimos los datos, que son la matriz de coeficientes y la de términos independientes
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
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
$$ \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}^{(0)}-a_{23}x_{3}^{(0)}-a_{24}x_{4}^{(0)}}{a_{22}}\\ x_{3}^{{\color{red}{(1)}}} & = & \dfrac{b_{3}-a_{31}x_{1}^{(0)}-a_{32}x_{2}^{(0)}-a_{34}x_{4}^{(0)}}{a_{33}}\\ x_{4}^{{\color{red}{(1)}}}& = & \dfrac{b_{4}-a_{41}x_{1}^{(0)}-a_{42}x_{2}^{(0)}-a_{43}x_{3}^{(0)}}{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}^{(k-1)}-\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)}}$
k = 1
x1[0] = (b[0] - A[0,1] * x0[1] - A[0,2] * x0[2] - A[0,3] * x0[3]) / A[0,0]
x1[1] = (b[1] - A[1,0] * x0[0] - A[1,2] * x0[2] - A[1,3] * x0[3]) / A[1,1]
x1[2] = (b[2] - A[2,0] * x0[0] - A[2,1] * x0[1] - A[2,3] * x0[3]) / A[2,2]
x1[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
print ('x1 = ', x1)
print ('x1-x0 = ', x1-x0)
print ('max(abs(x1-x0)) = ', max(abs(x1-x0)))
Como la norma es mayor que la $0.01$, seguimos iterando
Iteración 2
A partir de $\mathbf{x^{(1)}}$ calculamos $\mathbf{x^{(2)}}$ (aunque los guardamos respectivamente en x0 y x1)
k = 2
x0 = np.copy(x1)
x1[0] = (b[0] - A[0,1] * x0[1] - A[0,2] * x0[2] - A[0,3] * x0[3]) / A[0,0]
x1[1] = (b[1] - A[1,0] * x0[0] - A[1,2] * x0[2] - A[1,3] * x0[3]) / A[1,1]
x1[2] = (b[2] - A[2,0] * x0[0] - A[2,1] * x0[1] - A[2,3] * x0[3]) / A[2,2]
x1[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
print ('x1 = ', x1)
print ('x1-x0 = ', x1-x0)
print ('max(abs(x1-x0)) = ', max(abs(x1-x0)))
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)
k = 3
x0 = np.copy(x1)
x1[0] = (b[0] - A[0,1] * x0[1] - A[0,2] * x0[2] - A[0,3] * x0[3]) / A[0,0]
x1[1] = (b[1] - A[1,0] * x0[0] - A[1,2] * x0[2] - A[1,3] * x0[3]) / A[1,1]
x1[2] = (b[2] - A[2,0] * x0[0] - A[2,1] * x0[1] - A[2,3] * x0[3]) / A[2,2]
x1[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
print ('x1 = ', x1)
print ('x1-x0 = ', x1-x0)
print ('max(abs(x1-x0)) = ', max(abs(x1-x0)))
Como la norma es mayor que la $0.01$, seguimos iterando
Iteración 4
A partir de $\mathbf{x^{(3)}}$ calculamos $\mathbf{x^{(4)}}$ (aunque los guardamos respectivamente en x0 y x1)
k = 4
x0 = np.copy(x1)
x1[0] = (b[0] - A[0,1] * x0[1] - A[0,2] * x0[2] - A[0,3] * x0[3]) / A[0,0]
x1[1] = (b[1] - A[1,0] * x0[0] - A[1,2] * x0[2] - A[1,3] * x0[3]) / A[1,1]
x1[2] = (b[2] - A[2,0] * x0[0] - A[2,1] * x0[1] - A[2,3] * x0[3]) / A[2,2]
x1[3] = (b[3] - A[3,0] * x0[0] - A[3,1] * x0[1] - A[3,2] * x0[2]) / A[3,3]
Imprimimos el resultado
print ('x1 = ', x1)
print ('x1-x0 = ', x1-x0)
print ('max(abs(x1-x0)) = ', max(abs(x1-x0)))
Ahora ya se verifica que la norma del vector diferencia es menor que $0.01$ y paramos. Nuestra solución aproximada es
print ('x1 = ', x1)
La solución exacta es
x = np.linalg.solve(A,b);
print(x)
Y el error
print ('error = ', max(abs(x1-x)))