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 = np.array([[ 55, 0, 0,-25],
[ 0, 0, -37, -4],
[-25, 0, -4, 29],
[ 0, 1, 0, 0]],dtype=float)
B = np.array([[-200],
[-250],
[ 100],
[ -10]],dtype=float)
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 = np.array([[ 55, 0,-25],
[ 0, -37, -4],
[-25, -4, 29]],dtype=float)
B = np.array([[-200],
[-250],
[ 100]],dtype=float)
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 atras
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:
iteración: 1
Xnuevo: [-3.63636364 6.75675676 1.99891216]
diferencia: [3.63636364 6.75675676 1.99891216]
error: 6.756756756756757
iteración: 2
Xnuevo: [-2.7277672 6.54065815 1.99891216]
diferencia: [0.90859643 0.21609861 0. ]
error: 0.9085964348406557
iteración: 3
Xnuevo: [-2.7277672 6.54065815 1.99891216]
diferencia: [0. 0. 0.]
error: 0.0
con lo que el vector resultante es:
Método de Gauss-Seidel
respuesta de A.X=B :
[[-2.7277672 ]
[ 6.54065815]
[ 1.99891216]]
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:
# Método de Gauss,
# Recibe las matrices A,B
# presenta solucion X que cumple: A.X=B
import numpy as np
# INGRESO
A = np.array([[ 55, 0,-25],
[ 0, -37, -4],
[-25, -4, 29]],dtype=float)
B = np.array([[-200],
[-250],
[ 100]],dtype=float)
# PROCEDIMIENTO
# Matriz Aumentada
casicero = 1e-15 # 0
AB = np.concatenate((A,B),axis=1)
print('aumentada: ')
print(AB)
# pivotea por fila AB
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
for i in range(0,n-1,1):
# columna desde diagonal i en adelante
columna = np.abs(AB[i:,i])
dondemax = np.argmax(columna)
# revisa dondemax no está en la diagonal
if (dondemax != 0):
# intercambia fila
temporal = np.copy(AB[i,:])
AB[i,:] = AB[dondemax+i,:]
AB[dondemax+i,:] = temporal
print('pivoteada: ')
print(AB)
# Gauss elimina hacia adelante
# tarea: verificar términos cero
print('Elimina hacia adelante')
for i in range(0,n,1):
pivote = AB[i,i]
adelante = i+1
print('i: ', i)
for k in range(adelante,n,1):
print('k: ',k)
if (np.abs(pivote)>=casicero):
factor = AB[k,i]/pivote
AB[k,:] = AB[k,:] - factor*AB[i,:]
print(AB)
print('Elimina hacia adelante')
print(AB)
# 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]
X = np.transpose([X])
# Verifica resultado
verifica = np.dot(A,X)
# SALIDA
print('por sustitución hacia atras')
print('el vector solución X es:')
print(X)
print('verificar que A.X = B')
print(verifica)
# -------------------------------------------------
# # Algoritmo Gauss-Seidel,
# matrices, métodos iterativos
# Referencia: Chapra 11.2, p.310, pdf.334
# Rodriguez 5.2 p.162
# ingresar iteramax si requiere más iteraciones
import numpy as np
def gauss_seidel(A,B,tolera,X,iteramax=100):
tamano = np.shape(A)
n = tamano[0]
m = tamano[1]
diferencia = np.ones(n, dtype=float)
errado = np.max(diferencia)
itera = 0
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]
diferencia[i] = np.abs(nuevo-X[i])
X[i] = nuevo
errado = np.max(diferencia)
print(' iteración: ', itera)
print(' Xnuevo: ', X)
print(' diferencia: ', diferencia)
print(' error: ' ,errado)
itera = itera + 1
# Vector en columna
X = np.transpose([X])
# No converge
if (itera>iteramax):
X=0
return(X)
# Programa de prueba #######
tolera = 0.0001
# PROCEDIMIENTO
# usando la pivoteada por filas
n = len(B)
Xgs = np.zeros(n, dtype=float)
respuesta = gauss_seidel(AB[:,:n],AB[:,n],tolera,Xgs)
verificags = np.dot(A,respuesta)
# SALIDA
print(' Método de Gauss-Seidel')
print('respuesta de A.X=B : ')
print(respuesta)
print('verificar A.X: ')
print(verificags)
