Ejercicio: 1Eva_IIT2011_T2 Sistema de Ecuaciones, diagonal dominante
1. Desarrollo analítico
1.1 Solución iterativa usando el médodo de Jacobi
El sistema de ecuaciones
\begin{cases} -2x+5y+9z=1\\7x+y+z=6\\-3x+7y-z=-26\end{cases}
cambia a su forma matricial AX=B
\begin{bmatrix} -2 & 5 & 9 \\ 7 & 1 & 1 \\ -3 & 7 &-1 \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 1 \\ 6 \\ -26 \end{bmatrix}
para usarla en el algoritmo se intercambian filas, pivoteo, buscando hacerla diagonal dominante:
\begin{bmatrix} 7 & 1 & 1 \\ -3 & 7 &-1\\ -2 & 5 & 9 \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 6 \\ -26 \\1 \end{bmatrix}
las ecuaciones para el algoritmo se obtienen despejando una variable diferente en cada ecuación.
\begin{cases}x=(6-y-z)/(7)\\y=(-26+3x+z)/(7)\\z=(1+2x-5y)/(9)\end{cases}
Dado el vector inicial X(0) = [0, 0, 0], y una tolerancia de 0.0001, e desarrollan al menos 3 iteraciones:
Nota. El super índice entre paréntesis denota el número de la iteración. El subíndice denota el elemento del vector. El error se verifica para éste ejercicio como el mayor valor de la diferencia entre iteraciones consecutivas. Analogía al video del acople entre las aeronaves. Si una coordenada aún no es la correcta …. menor a lo tolerado, pues NO hay acople…
Iteración 1
\begin{cases}x=(6-(0)-(0))/(7)\\y=(-26+3(0)+(0))/(7)\\z=(1+2(0)-5(0))/(9)\end{cases}
X(1) = [6/7, -26/7, 1/9]
= [0.8571, -3,7142, 0.1111]
diferencia(1) = |X(1) - X(0)|
= |[0.8571,-3,7142,0.1111] - [0,0,0]|
= |[0.8571,-3,7142,0.1111]|
errado(1) = maximo|[0.8571,-3,7142,0.1111]|
= 26/7 ≅ 3.71,
error es más alto que lo tolerado
Iteración 2
\begin{cases}x=(6-(-26/7)-(1/9))/(7)\\y=(-26+3(6/7)+(1/9))/(7)\\z=(1+2(6/7)-5(-26/7))/(9)\end{cases}
X(2) = [1.3718, -3.3310 , 2,3650]
diferencia(2) = |X(2) - X(1)|
= |[1.3718, -3.3310 , 2,3650] - [0.8571, -3,7142, 0.1111]|
= |[0.5146, 0.3832, 2.2539]|
errado(2) ≅ 2.2538,
el error disminuye, pero es más alto que lo tolerado
Iteración 3
\begin{cases}x=(6-(-3.3310)-( 2,3650))/(7)\\y=(-26+3x+z)/(7)\\z=(1+2x-5y)/(9)\end{cases}
X(3) = [0.9951, -2.7884, 2.2665]
diferencia(3) = |[0.9951, -2.7884, 2.2665] - [1.3718, -3.3310 , 2,3650]|
= |[-0.3767, 0.5425, -0.0985]|
errado(3) ≅ 0.5425,
el error disminuye, pero es más alto que lo tolerado
Observación: Si el error disminuye en cada iteración, se puede intuir que se converge a la solución. Se puede continuar con la 4ta iteración…
2. Solución numérica usando Python
Usando el algoritmo desarrollado en clase se obtiene la respuesta con 13 iteraciones:
Xi, errado
[ 0.85714286 -3.71428571 0.11111111] 3.71428571429
[ 1.37188209 -3.33106576 2.36507937] 2.25396825397
[ 0.99514091 -2.78846777 2.26656589] 0.542597991578
[ 0.93170027 -2.96400162 1.8814023 ] 0.385163589245
[ 1.0117999 -3.04621384 1.96482318] 0.0834208883016
[ 1.01162724 -2.99996816 2.02829656] 0.0634733731283
[ 0.99595309 -2.99097453 2.00256614] 0.0257304176216
[ 0.99834406 -3.0013678 1.99408654] 0.0103932672897
[ 1.00104018 -3.00155447 2.0003919 ] 0.00630536414945
[ 1.00016608 -2.99949822 2.00109475] 0.00205624814582
[ 0.99977193 -2.99977243 1.99975814] 0.00133660433338
[ 1.00000204 -3.0001323 1.99982289] 0.000359867509479
[ 1.0000442 -3.00002443 2.00007395] 0.000251063280064
[ 0.99999292 -2.99997049 2.00002339] 5.39347670632e-05
iteraciones: 13
La gráfica muestra las coordenadas de las aproximaciones, observe el espiral que se va cerrando desde el punto inicial en verde al punto final en rojo.
Tarea: convierta el algoritmo a una función, así es más sencillo usar en otros ejercicios.
Se debe controlar el número de iteraciones para verificar la convergencia con iteramax.
NO es necesario almacenar los valores de los puntos, solo el último valor y el contador itera permite determinar si el sistema converge.
# 1ra Evaluación II Término 2011
# Tema 2. Sistema ecuaciones Jacobi
import numpy as np
# INGRESO
A = np.array([[ 7.0, 1, 1],
[-3.0, 7,-1],
[-2.0, 5, 9]])
B = np.array([6.0, -26.0, 1.0])
X = np.array([0.0, 0.0, 0.0],dtype=float)
tolera = 0.0001
iteramax = 100
# PROCEDIMIENTO
tamano = np.shape(A)
n = tamano[0]
m = tamano[1]
Xi1 = np.zeros(n, dtype=float)
errado = 2*tolera
itera = 0
while not(errado<=tolera or itera>=iteramax):
i = 0
while not(i>=n):
j = 0
nuevo = B [i]
while not(j>=m):
if (i!=j):
nuevo = nuevo - A[i,j]*X[j]
j = j+1
Xi1[i] = nuevo/A[i,i]
i = i+1
diferencia = np.abs(X - Xi1)
errado = np.max(diferencia)
X = np.copy(Xi1)
itera = itera +1
print(Xi1, errado)
# SALIDA
print('iteraciones: ', itera-1)
Tarea: Realice las modificaciones necesarias para usar el algoritmo de Gauss-Seidel. Luego compare resultados.
Revisión de resultados
Si se usan las ecuaciones sin cambiar a diagonal dominante, como fueron presentadas en el enunciado, el algoritmo Jacobi NO converge, el error aumenta, y muy rápido como para observar la espiral hacia afuera en una gráfica.
Xi, errado
[ -0.5 6. 26. ] 26.0
[ 131.5 -16.5 69.5] 132.0
[ 271. -984. -484.] 967.5
[-4638.5 -1407. -7675. ] 7191.0
[-38055.5 40150.5 4092.5] 41557.5
[ 118792. 262302. 395246.] 391153.5
[ 2434361.5 -1226784. 1479764. ] 2315569.5
[ 3591977.5 -18520288.5 -15890546.5] 17370310.5
....
