Ejercicio: 1Eva_IIT2011_T2 Sistema de Ecuaciones, diagonal dominante
1. Desarrollo analítico
Solución iterativa usando el método 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 \\ y \\ z \end{bmatrix} = \begin{bmatrix} 1 \\ 6 \\ -26 \end{bmatrix}
y luego como matriz aumentada y realizar el pivoteo parcial por filas, buscando hacerla diagonal dominante:
\begin{pmatrix} 7 & 1 & 1 &6\\ -3 & 7 &-1 & -26\\ -2 & 5 & 9 &1 \end{pmatrix}
las ecuaciones para el algoritmo se obtienen despejando una variable diferente en cada ecuación.
x_{i+1}=\frac{6-y_i-z_i}{7}
y_{i+1}=\frac{-26+3x_i+z_i}{7}
z_{i+1}=\frac{1+2x_i-5y_i}{9}
Dado el vector inicial X0= [0, 0, 0], y una tolerancia de 0.0001, e desarrollan al menos 3 iteraciones:
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…
itera = 0
x_1=\frac{6-(0)-(0)}{7} = 0.8571
y_1=\frac{-26+3(0)+(0)}{7} = -3,7142
z_1=\frac{1+2(0)-5(0)}{9} =0.111
diferencia = [ 0.8571 -0 , -3,7142 -0 ,0.1111- 0 ]
= [ 0.8571 , -3,7142,0.1111 ]
errado = max |diferencia| = 3,7142
error es más alto que lo tolerado
X_1= [ 0.8571 , -3.7142 ,0.1111 ]
Itera = 1
x_2=\frac{6-(-3,7142) -0.1111}{7} = 1.3718
y_2=\frac{-26+3(0.8571)+0.1111}{7} = -3.3310
z_2=\frac{1+2(0.8571)-5(-3,7142)}{9} = 2.3650
diferencia = [ 1.3718-0.8571 , -3.3310 -(-3.7142) , 2.3650 - 0.1111 ]
= [0.5146 , 0.3832, 2.2539 ]
errado = max |diferencia| = 2.2539
el error disminuye, pero es más alto que lo tolerado
X_2= [ 1.3718 , -3.3310, 2.3650 ]
Itera = 2
x_3 = \frac{6-(-3.3310)-2.3650}{7} = 0.9951
y_3 = \frac{-26+3(1.3718)+2.3650}{7} = -2.7884
z_3 = \frac{1+2(1.3718)-5(-3.3310)}{9} = 2.2665
diferencia = [ 0.9951 - 1.3718 , -2.7884 -(-3.3310) , 2.2665- 2.3650]
diferencia = [ -0.3766 , 0.5425 , -0.0985]
errado = max |diferencia| = 0.5425
el error disminuye, pero es más alto que lo tolerado
X_3= [ 0.9951 , -2.7884 , 2.2665]
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:
Matriz aumentada
[[ -2. 5. 9. 1.]
[ 7. 1. 1. 6.]
[ -3. 7. -1. -26.]]
Pivoteo parcial:
1 intercambiar filas: 0 y 1
2 intercambiar filas: 1 y 2
AB
Iteraciones Jacobi
itera,[X],errado
0 [0. 0. 0.] 1.0
1 [ 0.85714286 -3.71428571 0.11111111] 3.7142857142857144
2 [ 1.37188209 -3.33106576 2.36507937] 2.2539682539682544
3 [ 0.99514091 -2.78846777 2.26656589] 0.5425979915775843
4 [ 0.93170027 -2.96400162 1.8814023 ] 0.3851635892452223
5 [ 1.0117999 -3.04621384 1.96482318] 0.0834208883015708
6 [ 1.01162724 -2.99996816 2.02829656] 0.06347337312830126
7 [ 0.99595309 -2.99097453 2.00256614] 0.025730417621558477
8 [ 0.99834406 -3.0013678 1.99408654] 0.010393267289710018
9 [ 1.00104018 -3.00155447 2.0003919 ] 0.006305364149447268
10 [ 1.00016608 -2.99949822 2.00109475] 0.002056248145824391
11 [ 0.99977193 -2.99977243 1.99975814] 0.0013366043333828959
12 [ 1.00000204 -3.0001323 1.99982289] 0.0003598675094789172
13 [ 1.0000442 -3.00002443 2.00007395] 0.0002510632800638568
14 [ 0.99999292 -2.99997049 2.00002339] 5.3934767063168465e-05
respuesta de A.X=B :
[ 0.99999292 -2.99997049 2.00002339]
iterado, incluyendo X0: 15
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.
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.
# 1Eva_IIT2011_T2 Sistema de Ecuaciones, diagonal dominante
# MATG1052 Métodos Numéricos
# propuesta: edelros@espol.edu.ec,
import numpy as np
def jacobi(A,B,X0, tolera, iteramax=100, vertabla=False):
''' Método de Jacobi, tolerancia, vector inicial X0
para mostrar iteraciones: vertabla=True
'''
A = np.array(A,dtype=float)
B = np.array(B,dtype=float)
X0 = np.array(X0,dtype=float)
tamano = np.shape(A)
n = tamano[0]
m = tamano[1]
diferencia = np.ones(n, dtype=float)
errado = np.max(diferencia)
X = np.copy(X0)
xnuevo = np.copy(X0)
tabla = [np.copy(X0)]
itera = 0
if vertabla==True:
print('Iteraciones Jacobi')
print('itera,[X],errado')
print(itera, xnuevo, errado)
while not(errado<=tolera or itera>iteramax):
for i in range(0,n,1):
nuevo = B[i]
for j in range(0,m,1):
if (i!=j): # excepto diagonal de A
nuevo = nuevo-A[i,j]*X[j]
nuevo = nuevo/A[i,i]
xnuevo[i] = nuevo
diferencia = np.abs(xnuevo-X)
errado = np.max(diferencia)
X = np.copy(xnuevo)
tabla = np.concatenate((tabla,[X]),axis = 0)
itera = itera + 1
if vertabla==True:
print(itera, xnuevo, errado)
# No converge
if (itera>iteramax):
X=itera
print('iteramax superado, No converge')
return(X,tabla)
def pivoteafila(A,B,vertabla=False):
'''
Pivotea parcial por filas, entrega matriz aumentada AB
Si hay ceros en diagonal es matriz singular,
Tarea: Revisar si diagonal tiene ceros
'''
A = np.array(A,dtype=float)
B = np.array(B,dtype=float)
# Matriz aumentada
nB = len(np.shape(B))
if nB == 1:
B = np.transpose([B])
AB = np.concatenate((A,B),axis=1)
if vertabla==True:
print('Matriz aumentada')
print(AB)
print('Pivoteo parcial:')
# Pivoteo por filas AB
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
# Para cada fila en AB
pivoteado = 0
for i in range(0,n-1,1):
# columna desde diagonal i en adelante
columna = np.abs(AB[i:,i])
dondemax = np.argmax(columna)
# dondemax no es en diagonal
if (dondemax != 0):
# intercambia filas
temporal = np.copy(AB[i,:])
AB[i,:] = AB[dondemax+i,:]
AB[dondemax+i,:] = temporal
pivoteado = pivoteado + 1
if vertabla==True:
print(' ',pivoteado, 'intercambiar filas: ',i,'y', dondemax+i)
if vertabla==True:
if pivoteado==0:
print(' Pivoteo por filas NO requerido')
else:
print('AB')
return(AB)
# PROGRAMA Búsqueda de solucion
# INGRESO --------
# INGRESO
A = [[-2, 5, 9],
[ 7, 1, 1],
[-3, 7,-1]]
B = [1, 6,-26]
X0 = [0, 0, 0]
tolera = 0.0001
iteramax = 100
# PROCEDIMIENTO --------
AB = pivoteafila(A,B,vertabla=True)
# separa matriz aumentada en A y B
[n,m] = np.shape(AB)
A = AB[:,:m-1]
B = AB[:,m-1]
[X, puntos] = jacobi(A,B,X0,tolera,vertabla=True)
iterado = len(puntos)
# SALIDA --------
print('respuesta de A.X=B : ')
print(X)
print('iterado, incluyendo X0: ', iterado)
# GRAFICA DE ITERACIONES
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
figura = plt.figure() # Una hoja para dibujar
grafica = figura.add_subplot(111,projection = '3d')
# Puntos de la iteración
pxi = puntos[:,0]
pyi = puntos[:,1]
pzi = puntos[:,2]
grafica.plot(pxi,pyi,pzi, color = 'magenta',
label = 'puntos[i]')
grafica.scatter(pxi[0],pyi[0],pzi[0], color = 'green', marker='o')
grafica.scatter(pxi[iterado-1],pyi[iterado-1],pzi[iterado-1],
color = 'red', marker='o')
grafica.set_title('Puntos de iteración')
grafica.set_xlabel('eje x')
grafica.set_ylabel('eje y')
grafica.set_zlabel('eje z')
grafica.legend()
plt.show()
Gráfica de planos
Dado que el sistema de ecuaciones es de tres incógnitas, la solución se puede interpretar como la intersección de los planos formados por cada una de las ecuaciones en el espacio X,Y,Z.
se comenta la última instrucción del algoritmo anterior, #plt.show(), y se añaden las siguientes líneas al algoritmo anterior para obtener la gráfica:
# Tema 2. Sistema de ecuaciones 3x3
# Concepto como interseccion de Planos
##import matplotlib.pyplot as plt
##from mpl_toolkits.mplot3d import Axes3D
ax = -5 # Intervalo X
bx = 5
ay = ax-2 # Intervalo Y
by = bx+2
muestras = 11
# PROCEDIMIENTO --------
# Ecuaciones de planos
z0 = lambda x,y: (-A[0,0]*x - A[0,1]*y + B[0])/A[0,2]
z1 = lambda x,y: (-A[1,0]*x - A[1,1]*y + B[1])/A[1,2]
z2 = lambda x,y: (-A[2,0]*x - A[2,1]*y + B[2])/A[2,2]
xi = np.linspace(ax,bx, muestras)
yi = np.linspace(ay,by, muestras)
Xi, Yi = np.meshgrid(xi,yi)
Z0 = z0(Xi,Yi)
Z1 = z1(Xi,Yi)
Z2 = z2(Xi,Yi)
# solución al sistema
punto = np.linalg.solve(A,B)
# SALIDA
print('respuesta de A.X=B : ')
print(punto)
# Interseccion entre ecuación 1 y 2
# PlanoXZ, extremo inferior de y
Aa = np.copy(A[0:2,[0,2]])
Ba = np.copy(B[0:2])
Ba = Ba-ay*A[0:2,1]
pta = np.linalg.solve(Aa,Ba)
pa = np.array([ay])
pxza = np.array([pta[0],ay,pta[1]])
# PlanoXZ, extremo superior de y
Bb = np.copy(B[0:2])
Bb = Bb-by*A[0:2,1]
ptb = np.linalg.solve(Aa,Bb)
pb = np.array([by])
pxzb = np.array([ptb[0],by,ptb[1]])
# GRAFICA de planos
fig3D = plt.figure()
graf3D = fig3D.add_subplot(111, projection='3d')
graf3D.plot_wireframe(Xi,Yi,Z0,
color ='blue',
label='Z1')
graf3D.plot_wireframe(Xi,Yi,Z1,
color ='green',
label='Z2')
graf3D.plot_wireframe(Xi,Yi,Z2,
color ='orange',
label='Z3')
# recta intersección planos 1 y 2
graf3D.plot([pxza[0],pxzb[0]],
[pxza[1],pxzb[1]],
[pxza[2],pxzb[2]],
label='Sol 1y2',
color = 'violet',
linewidth = 4)
# Punto solución del sistema 3x3
graf3D.scatter(punto[0],punto[1],punto[2],
color = 'red',
marker='o',
label ='punto',
linewidth = 6)
graf3D.set_title('Sistema de ecuaciones 3x3')
graf3D.set_xlabel('x')
graf3D.set_ylabel('y')
graf3D.set_zlabel('z')
graf3D.set_xlim(ax, bx)
graf3D.set_xlim(ay, by)
graf3D.legend()
graf3D.view_init(45, 45)
# rotacion de ejes
for angulo in range(45, 360+45, 5 ):
graf3D.view_init(45, angulo)
plt.draw()
plt.pause(.001)
plt.show()
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 pivorear y 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
....
