s1Eva_IIT2019_T3 Circuito eléctrico

Ejercicio: 1Eva_IIT2019_T3 Circuito eléctrico

Las ecuaciones del problema son:

55 I_1 - 25 I_4 =-200 -37 I_3 - 4 I_4 =-250 -25 I_1 - 4 I_3 +29 I_4 =100 I_2 =-10

Planteo del problema en la forma A.X=B

A = [[ 55.0, 0,  0, -25],
     [  0  , 0,-37,  -4],
     [-25  , 0, -4,  29],
     [  0  ,  1, 0,   0]]

B = [-200,-250,100,-10]

El ejercicio se puede simplificar con una matriz de 3×3 dado que una de las corrientes I2 es conocida con valor -10, queda resolver el problema para
[I1 ,I3 ,I4 ]

A = [[ 55.0,   0, -25],
     [  0  , -37,  -4],
     [-25  ,  -4,  29]]

B = [-200,-250,100]

conformando la matriz aumentada

[[  55.    0.  -25. -200.]
 [   0.  -37.   -4. -250.]
 [ -25.   -4.   29.  100.]]

que se pivotea por filas para acercar a matriz diagonal dominante:

[[  55.    0.  -25. -200.]
 [   0.  -37.   -4. -250.]
 [ -25.   -4.   29.  100.]]

Literal a

Para métodos directos se aplica el método de eliminación hacia adelante.

Usando el primer elemento en la diagonal se convierten en ceros los números debajo de la posición primera de la diagonal

[[  55.    0.  -25.         -200.      ]
 [   0.  -37.   -4.         -250.      ]
 [   0.   -4.   17.636363      9.090909]]

luego se continúa con la segunda columna:

[[  55.    0.  -25.         -200.      ]
 [   0.  -37.   -4.         -250.      ]
 [   0.    0.   18.068796     36.117936]]

y para el método de Gauss se emplea sustitución hacia atrás
se determina el valor de I4

18.068796 I_4 = 36.11793612 I_4 =\frac{36.11793612}{18.068796}= 1.99891216 -37 I_3 -4 I_4 = -250 -37 I_3= -250 + 4 I_4 I_3=\frac{-250 + 4 I_4}{-37} I_3=\frac{-250 + 4 (1.99891216)}{-37} = 6.54065815

y planteando se obtiene el último valor

55 I_1 +25 I_4 = -200 55 I_1 = -200 -25 I_4 I_1 = \frac{-200 -25 I_4}{55} I_1 = \frac{-200 -25(1.99891216)}{55} = -2.7277672

con lo que el vector solución es:

[-2.7277672   6.54065815  1.99891216]

sin embargo, para verificar la respuesta se aplica A.X=B

verificar que A.X = B, obteniendo nuevamente el vector B.
[-200.  -250.  100.]]

literal b

La norma de la matriz infinito se determina como:

||x|| = max\Big[ |x_i| \Big]

considere que en el problema el término en A de magnitud mayor es 55.
El vector suma de filas es:

[[| 55|+|  0|+|-25|],    [[80],
 [|  0|+|-37|+| -4|],  =  [41],
 [[-25|+| -4|+| 29|]]     [58]]

por lo que la norma ∞ ejemplo ||A||∞ 
es el maximo de suma de filas: 80

para revisar la estabilidad de la solución, se observa el número de condición

>>> np.linalg.cond(A)
4.997509004325602

En éste caso no está muy alejado de 1. De resultar alejado del valor ideal de uno,  la solución se considera poco estable. Pequeños cambios en la entrada del sistema generan grandes cambios en la salida.

Tarea: Matriz de transición de Jacobi


Literal c

En el método de Gauss-Seidel acorde a lo indicado, se inicia con el vector cero. Como no se indica el valor de tolerancia para el error, se considera tolera = 0.0001

las ecuaciones para el método son:

I_1 =\frac{-200 + 25 I_4}{55} I_3 = \frac{-250+ 4 I_4}{-37} I_4 =\frac{100 +25 I_1 + 4 I_3}{29}

Como I2 es constante, no se usa en las iteraciones

I_2 =-10

teniendo como resultados de las iteraciones:

Matriz aumentada
[[  55.    0.  -25. -200.]
 [   0.  -37.   -4. -250.]
 [ -25.   -4.   29.  100.]]
Pivoteo parcial:
  Pivoteo por filas NO requerido
Iteraciones Gauss-Seidel
itera,[X]
      diferencia,errado
0 [0. 0. 0.] 2e-05
1 [-3.6363636  6.7567568  1.2454461]
   [3.6363636 6.7567568 1.2454461] 6.756756756756757
2 [-3.0702518  6.6221139  1.7149021]
   [0.5661119 0.1346428 0.469456 ] 0.5661118513783094
3 [-2.8568627  6.5713619  1.891858 ]
   [0.2133891 0.050752  0.1769559] 0.2133891067340583
4 [-2.7764282  6.5522316  1.9585594]
   [0.0804345 0.0191304 0.0667014] 0.08043447732439457
5 [-2.7461094  6.5450206  1.9837016]
   [0.0303188 0.007211  0.0251423] 0.030318816370094925
6 [-2.7346811  6.5423025  1.9931787]
   [0.0114283 0.0027181 0.0094771] 0.011428316023939011
7 [-2.7303733  6.541278   1.996751 ]
   [0.0043078 0.0010246 0.0035723] 0.004307767346479974
8 [-2.7287495  6.5408918  1.9980975]
   [0.0016238 0.0003862 0.0013465] 0.001623761494915943
9 [-2.7281375  6.5407462  1.9986051]
   [0.0006121 0.0001456 0.0005076] 0.0006120575185017962
10 [-2.7279068  6.5406913  1.9987964]
   [2.3070778e-04 5.4871039e-05 1.9131760e-04] 0.00023070777766820427
11 [-2.7278198  6.5406707  1.9988685]
   [8.6962544e-05 2.0682983e-05 7.2114885e-05] 8.696254366213907e-05
12 [-2.727787   6.5406629  1.9988957]
   [3.2779493e-05 7.7962038e-06 2.7182845e-05] 3.277949307367578e-05
13 [-2.7277747  6.5406599  1.998906 ]
   [1.2355839e-05 2.9386860e-06 1.0246249e-05] 1.235583874370505e-05
14 [-2.72777    6.5406588  1.9989098]
   [4.6573860e-06 1.1077026e-06 3.8622013e-06] 4.6573859666665385e-06
numero de condición: 4.997509004325604
respuesta con Gauss-Seidel
[-2.72777    6.5406588  1.9989098]
>>>

con lo que el vector resultante es:

respuesta con Gauss-Seidel
[-2.72777 6.5406588 1.9989098]

que para verificar, se realiza la operación A.X
observando que el resultado es igual a B

[[-200.00002751]
 [-249.9999956 ]
 [ 100.0000125 ]]


Solución alterna


Usando la matriz de 4×4, los resultados son iguales para las corrientes
[I1 ,I2 , I3 ,I4 ]. Realizando la matriz aumentada,

[[  55.    0.    0.  -25. -200.]
 [   0.    0.  -37.   -4. -250.]
 [ -25.    0.   -4.   29.  100.]
 [   0.    1.    0.    0.  -10.]]

que se pivotea por filas para acercar a matriz diagonal dominante:

[[  55.    0.    0.  -25. -200.]
 [   0.    1.    0.    0.  -10.]
 [   0.    0.  -37.   -4. -250.]
 [ -25.    0.   -4.   29.  100.]]

Literal a

Para métodos directos se aplica el método de eliminación hacia adelante.

Usando el primer elemento  en la diagonal.

[[  55.     0.     0.   -25.         -200.        ]
 [   0.     1.     0.     0.          -10.        ]
 [   0.     0.   -37.    -4.         -250.        ]
 [   0.     0.    -4.    17.63636364    9.09090909]]

para el segundo no es necesario, por debajo se encuentran valores cero.
Por lo que se pasa al tercer elemento de la diagonal

[[  55.     0.     0.     -25.         -200.        ]
 [   0.     1.     0.      0.          -10.        ]
 [   0.     0.   -37.     -4.         -250.        ]
 [   0.     0.     0.     18.06879607    36.11793612]]

y para el método de Gauss se emplea sustitución hacia atras.
para x4:

18.06879607 x_4 = 36.11793612 x_4 = 1.99891216

para x3:

-37 x_3 -4 x_3 = -250 37 x_3 = 250-4 x_4 = 250-4(1.99891216) x_3 = 6.54065815

como ejercicio, continuar con x1, dado que x2=-10

55 x_1 + 25 x_4 = -200

El vector solución obtenido es:

el vector solución X es:
[[ -2.7277672 ]
 [-10.        ]
 [  6.54065815]
 [  1.99891216]]

sin embargo, para verificar la respuesta se aplica A.X=B.

[[-200.]
 [-250.]
 [ 100.]
 [ -10.]]

Se revisa el número de condición de la matriz:

>>> np.linalg.cond(A)
70.21827416891405

Y para éste caso, el número de condición se encuentra alejado del valor 1, contrario a la respuesta del la primera forma de solución con la matriz 3×3. De resultar alejado del valor ideal de uno, la solución se considera poco estable. Pequeños cambios en la entrada del sistema generan grandes cambios en la salida.


Algoritmo en Python

Presentado por partes para revisión:

Para el método de Gauss, los resultados del algoritmo se muestran como:

Matriz aumentada
[[  55.    0.  -25. -200.]
 [   0.  -37.   -4. -250.]
 [ -25.   -4.   29.  100.]]
Pivoteo parcial:
  Pivoteo por filas NO requerido
Elimina hacia adelante:
 fila 0 pivote:  55.0
   factor:  0.0  para fila:  1
   factor:  -0.45454545454545453  para fila:  2
 fila 1 pivote:  -37.0
   factor:  0.10810810810810811  para fila:  2
 fila 2 pivote:  18.06879606879607
[[  55.            0.          -25.         -200.        ]
 [   0.          -37.           -4.         -250.        ]
 [   0.            0.           18.06879607   36.11793612]]
solución: 
[-2.7277672   6.54065815  1.99891216]

Instrucciones en Python usando las funciones creadas en la unidad:

# 1Eva_IIT2019_T3 Circuito eléctrico
# Método de Gauss
# Solución a Sistemas de Ecuaciones
# de la forma A.X=B
import numpy as np

def pivoteafila(A,B,vertabla=False):
    '''
    Pivotea parcial por filas
    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)

def gauss_eliminaAdelante(AB,vertabla=False, casicero = 1e-15):
    ''' Gauss elimina hacia adelante
    tarea: verificar términos cero
    '''
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    if vertabla==True:
        print('Elimina hacia adelante:')
    for i in range(0,n,1):
        pivote = AB[i,i]
        adelante = i+1
        if vertabla==True:
            print(' fila',i,'pivote: ', pivote)
        for k in range(adelante,n,1):
            if (np.abs(pivote)>=casicero):
                factor = AB[k,i]/pivote
                AB[k,:] = AB[k,:] - factor*AB[i,:]
                if vertabla==True:
                    print('   factor: ',factor,' para fila: ',k)
            else:
                print('  pivote:', pivote,'en fila:',i,
                      'genera division para cero')
    if vertabla==True:
        print(AB)
    return(AB)

def gauss_sustituyeAtras(AB,vertabla=False, casicero = 1e-15):
    ''' Gauss sustituye hacia atras
    '''
    tamano = np.shape(AB)
    n = tamano[0]
    m = tamano[1]
    # Sustitución hacia atras
    X = np.zeros(n,dtype=float) 
    ultfila = n-1
    ultcolumna = m-1
    for i in range(ultfila,0-1,-1):
        suma = 0
        for j in range(i+1,ultcolumna,1):
            suma = suma + AB[i,j]*X[j]
        X[i] = (AB[i,ultcolumna]-suma)/AB[i,i]
    return(X)

# INGRESO
A = [[ 55.0,   0, -25],
     [  0  , -37,  -4],
     [-25  ,  -4,  29]]

B = [-200,-250,100]

# PROCEDIMIENTO
AB = pivoteafila(A,B,vertabla=True)

AB = gauss_eliminaAdelante(AB,vertabla=True)

X = gauss_sustituyeAtras(AB,vertabla=True)

# SALIDA
print('solución: ')
print(X)

literal c

Resultados con el algoritmo de Gauss Seidel

Matriz aumentada
[[  55.    0.  -25. -200.]
 [   0.  -37.   -4. -250.]
 [ -25.   -4.   29.  100.]]
Pivoteo parcial:
  Pivoteo por filas NO requerido
Iteraciones Gauss-Seidel
itera,[X]
      diferencia,errado
0 [0. 0. 0.] 2e-05
1 [-3.6363636  6.7567568  1.2454461]
   [3.6363636 6.7567568 1.2454461] 6.756756756756757
2 [-3.0702518  6.6221139  1.7149021]
   [0.5661119 0.1346428 0.469456 ] 0.5661118513783094
3 [-2.8568627  6.5713619  1.891858 ]
   [0.2133891 0.050752  0.1769559] 0.2133891067340583
4 [-2.7764282  6.5522316  1.9585594]
   [0.0804345 0.0191304 0.0667014] 0.08043447732439457
5 [-2.7461094  6.5450206  1.9837016]
   [0.0303188 0.007211  0.0251423] 0.030318816370094925
6 [-2.7346811  6.5423025  1.9931787]
   [0.0114283 0.0027181 0.0094771] 0.011428316023939011
7 [-2.7303733  6.541278   1.996751 ]
   [0.0043078 0.0010246 0.0035723] 0.004307767346479974
8 [-2.7287495  6.5408918  1.9980975]
   [0.0016238 0.0003862 0.0013465] 0.001623761494915943
9 [-2.7281375  6.5407462  1.9986051]
   [0.0006121 0.0001456 0.0005076] 0.0006120575185017962
10 [-2.7279068  6.5406913  1.9987964]
   [2.3070778e-04 5.4871039e-05 1.9131760e-04] 0.00023070777766820427
11 [-2.7278198  6.5406707  1.9988685]
   [8.6962544e-05 2.0682983e-05 7.2114885e-05] 8.696254366213907e-05
12 [-2.727787   6.5406629  1.9988957]
   [3.2779493e-05 7.7962038e-06 2.7182845e-05] 3.277949307367578e-05
13 [-2.7277747  6.5406599  1.998906 ]
   [1.2355839e-05 2.9386860e-06 1.0246249e-05] 1.235583874370505e-05
14 [-2.72777    6.5406588  1.9989098]
   [4.6573860e-06 1.1077026e-06 3.8622013e-06] 4.6573859666665385e-06
numero de condición: 4.997509004325604
respuesta con Gauss-Seidel
[-2.72777    6.5406588  1.9989098]
>>> 

Instrucciones en Python

# 1Eva_IIT2019_T3 Circuito eléctrico
# Algoritmo Gauss-Seidel
# solución de matrices
# métodos iterativos
# Referencia: Chapra 11.2, p.310,
#      Rodriguez 5.2 p.162
import numpy as np

def gauss_seidel(A,B,X0,tolera, iteramax=100,
                 vertabla=False, precision=4):
    ''' Método de Gauss Seidel, 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 = 2*tolera*np.ones(n, dtype=float)
    errado = np.max(diferencia)
    X = np.copy(X0)

    itera = 0
    if vertabla==True:
        print('Iteraciones Gauss-Seidel')
        print('itera,[X]')
        print('      diferencia,errado')
        print(itera, X, errado)
        np.set_printoptions(precision)
    while (errado>tolera and itera<iteramax):
        for i in range(0,n,1):
            xi = B[i]
            for j in range(0,m,1):
                if (i!=j):
                    xi = xi-A[i,j]*X[j]
            xi = xi/A[i,i]
            diferencia[i] = np.abs(xi-X[i])
            X[i] = xi
        errado = np.max(diferencia)
        itera = itera + 1
        if vertabla==True:
            print(itera, X)
            print('  ', diferencia, errado)
    # No converge
    if (itera>iteramax):
        X=itera
        print('iteramax superado, No converge')
    return(X)

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)

# INGRESO
A = [[ 55.0,   0, -25],
     [  0  , -37,  -4],
     [-25  ,  -4,  29]]

B = [-200,-250,100]

X0  = [0.,0.,0.]

tolera = 0.00001
iteramax = 100
verdecimal = 7

# PROCEDIMIENTO
# numero de condicion
ncond = np.linalg.cond(A)

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]

respuesta = gauss_seidel(A,B,X0, tolera,
                         vertabla=True, precision=verdecimal)

# SALIDA
print('numero de condición:', ncond)
print('respuesta con Gauss-Seidel')
print(respuesta)

.