s1Eva_IIT2011_T2 Sistema de Ecuaciones, diagonal dominante

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.

Sistema Ecuaciones Diagonal Dominante 02

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
....