En esta ocasión exploraremos una implementación del método de Jacobi para la aproximación de solución a sistemas de ecuaciones lineales.
Empecemos declarando un sistema de ecuaciones lineales en su forma matricial:
In [ ]:
from numpy import matrix
In [ ]:
A = matrix([[72, 0, 0, 9, 0, 0],
[ 0, 2.88, 0, 0, 0, -4.5],
[ 0, 0, 18, 9, 0, 0],
[ 9, 0, 9, 12, 0, 0],
[ 0, 0, 0, 0, 33, 0],
[ 0, -4.5, 0, 0, 0, 33]])
b = matrix([[2],
[0.5],
[1],
[0],
[1.2],
[5]])
Con lo que, trivialmente, podemos ver su solución por medio de la implementación mas simple:
In [ ]:
A.I*b
Sin embargo, deseamos implementar un método diferente, por lo que empezaremos obteniendo información acerca del sistema, por ejemplo, si queremos saber la dimensión de la matriz $A$, y por consecuente la dimensión del sistema, tendremos:
In [ ]:
A.shape
Ya sabemos que esta matriz es cuadrada, por lo que podemos obtener el primer elemento solamente con:
In [ ]:
A.shape[0]
Recordando el método de Jacobi, tenemos que empezar con una aproximación $0$ en todas las variables, por lo que definimos un arreglo de ceros, con la funcion zeros, y lo convertimos en un objeto de tipo matrix:
In [ ]:
from numpy import zeros, copy
In [ ]:
x0 = matrix(zeros(A.shape[0])).T
x0
Lo siguiente que voy a necesitar es un arreglo con todos los elementos de la diagonal principal, para dividir cada una de las formulas de la aproximación con estos elementos:
In [ ]:
A.diagonal()
In [ ]:
divisores = copy(A.diagonal().T)
divisores
Una vez que obtengo la diagonal, voy a sustituir todos los elementos de esta por $0$, de manera que cuando multiplique cada una de las filas de esta matriz, por el valor de la aproximación de $x$, $x_0$, obtendré la mayoría de los elementos de cada una de las formulas:
$$ 0 x_1 + a_{12}x_2 + a_{13}x_3+\ldots+a_{1n}x_n $$
In [ ]:
from numpy import fill_diagonal
In [ ]:
fill_diagonal(A, 0)
In [ ]:
A
In [ ]:
-A*x0
Tomando en cuenta que nuestra primera aproximación es toda $0$, este vector es completamente correcto, ahora tenemos que dividir cada uno de estos elementos por el vector diagonal:
$$ \frac{0 x_1 + a_{12}x_2 + a_{13}x_3+\ldots+a_{1n}x_n}{a_{11}} $$
In [ ]:
from numpy import divide
In [ ]:
divide(-A*x0, divisores)
Por ultimo tenemos que restar cada uno de estos elementos a la division del vector $b$ con la diagonal, de tal manera que tengamos:
$$ \frac{b}{a_{11}} -\frac{0 x_1 + a_{12}x_2 + a_{13}x_3+\ldots+a_{1n}x_n}{a_{11}} $$
In [ ]:
x1 = divide(b, divisores) - divide(A*x0, divisores)
x1
lo cual corresponde a nuestra primer aproximación, por lo que lo hemos guardado en x1.
Si ahora obtenemos la diferencia entre las primeras dos aproximaciones, y lo guardamos en la variable dif:
In [ ]:
dif = x1 - x0
dif
podemos obtener el cuadrado de cada una de estas diferencias con dif.T*dif, es decir:
y a este valor, tan solo tenemos que sacarle la raiz cuadrada:
In [ ]:
from numpy import sqrt
In [ ]:
e1 = sqrt((dif.T*dif)[0,0])
e1
por lo que obtenemos nuestro error inicial. Ahora tan solo tenemos que seguir haciendo iteraciones hasta obtener una buena aproximación:
In [ ]:
x2 = divide(b, divisores) - divide(A*x1, divisores)
x2
In [ ]:
dif = x2 - x1
e2 = sqrt((dif.T*dif)[0,0])
e2
In [ ]:
x3 = divide(b, divisores) - divide(A*x2, divisores)
x3
In [ ]:
dif = x3 - x2
e3 = sqrt((dif.T*dif)[0,0])
e3
Este valor ya es menor que $0.1$, por lo que terminaremos el ejemplo aqui...
Implemente la función que calcule el método de Jacobi, utilizando las partes necesarias de la explicación anterior, considere el siguiente codigo como un inicio:
def jacobi(A, b, error):
from numpy import matrix, zeros, copy, fill_diagonal, divide, sqrt
xs = []
# Inicialize valores a utilizar y arreglos a llenar, como el anterior
while True:
# Escribe tu codigo aqui
if sqrt(xs[-1].T*xs[-1] - xs[-2].T*xs[-2]) < error:
return xs[-1]