Ejercicio: 1Eva_IIT2008_T1 Distribuidores de productos
literal a
Siguiendo las instrucciones del enunciado, el promedio de precios del nodo A, se conforma de los precios en los nodos aledaños menos el costo de transporte.
precio en X1 para A = precio en nodoA – costo de transporteA
siguiendo el mismo procedimiento,
precio en X1 para A: (3.1-0.2)
precio de X1 para B: (2.8-0.3)
precio de X1 para C: (2.7-0.4)
precio de X1 para X2: (X2-0.1)
precio de X1 para X3: (X3-0.5)
Se continua con el mismo proceso para los siguientes nodos:
x_2 = \frac{1}{4} \Big[ (3.2-0.5)+(3.4-0.3) +(x_1-0.1)+(x_3-0.2)\Big] 4x_2 = (3.2-0.5)+(3.4-0.3) +(x_1-0.1)+(x_3-0.2) 4x_2 = 2.7+3.1 +x_1+x_3-0.3 -x_1+4x_2-x_3 = 5.5Para X3
x_3 = \frac{1}{4} \Big[ (3.3-0.3)+(2.9-0.2) +(x_1-0.5)+(x_2-0.2)\Big] 4x_3 = 3.0+2.7+x_1+x_2-0.7 4x_3 = 5+x_1+x_2 -x_1-x_2+4x_3= 5El sistema de ecuaciones se convierte en:
\begin{cases} 5x_1-x_2-x_3 = 7.1 \\ -x_1+4x_2-x_3 = 5.5 \\-x_1-x_2+4x_3= 5 \end{cases}Para resolver se cambia a la forma Ax=B
\begin{bmatrix} 5 && -1 && -1 \\ -1 && 4 && -1 \\ -1 && -1 && 4 \end{bmatrix} . \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 7.1 \\ 5.5 \\5 \end{bmatrix}En los métodos directos se usa la forma de matriz aumentada
\begin{bmatrix} 5 && -1 && -1 && 7.1 \\ -1 && 4 && -1 && 5.5 \\ -1 && -1 && 4 && 5 \end{bmatrix}pivoteo: no es necesario, pues la matriz ya está ordenada de forma diagonalmente dominante.
Eliminación hacia adelante
i=0, f=0
\begin{bmatrix} 5 && -1 && -1 && 7.1 \\ -1-5\big(\frac{-1}{5}\big) && 4 -(-1)\big(\frac{-1}{5}\big) && -1 -(-1)\big(\frac{-1}{5}\big) && 5.5-7.1\big(\frac{-1}{5}\big) \\ -1-5\big(\frac{-1}{5}\big) && -1-(-1)1\big(\frac{-1}{5}\big) && 4-(-1)\big(\frac{-1}{5}\big) && 5-7.1\big(\frac{-1}{5}\big) \end{bmatrix} \begin{bmatrix} 5 && -1 && -1 && 7.1 \\ 0 && 3.8 && -1.2 && 6.92 \\ 0 && -1.2 && 3.8 && 6.42 \end{bmatrix}Eliminación hacia adelante para i = 1, j=1
Elimina hacia adelante [[ 5. -1. -1. 7.1 ] [ 0. 3.8 -1.2 6.92 ] [ 0. 0. 3.42105263 8.60526316]]
Eliminación hacia atrás, continuando el desarrollo de forma semejante a los pasos anteriores se obtiene:
Elimina hacia atrás [[ 1. -0. -0. 2.44615385] [ 0. 1. -0. 2.61538462] [ 0. 0. 1. 2.51538462]] el vector solución X es: [[2.44615385] [2.61538462] [2.51538462]] verificar que A.X = B [[7.1] [5.5] [5. ]]
literal b
Para el literal b se usa como referencia el número de condición:
>>> np.linalg.cond(A) 2.5274158815808474
El número de condición es cercano a 1, dado que la matriz A es diagonalmente dominante pues los valores mayores de la fila se encuentran en la diagonal. Como el número de condición es cercano a 1 el sistema converge usando métodos iterativos.
La selección del vector inicial para las iteraciones siguiendo el enunciado del problema, se evita el vector cero dado que el precio de un producto para una fabrica no puede ser cero. Se observa los valores de los precios, y se encuentra que el rango de existencia en los nodos es [ 2.7, 3.4] que restando el valor de transporte podrían ser un valor menor a 2.7. Por lo que un buen vector inicial será [2,2,2]
literal c
Se plantean las ecuaciones para el método de Jacobi a partir del sistema de ecuaciones, a partir del pivoteo por filas:
\begin{cases} 5x_1-x_2-x_3 = 7.1 \\ -x_1+4x_2-x_3 = 5.5 \\-x_1-x_2+4x_3= 5 \end{cases} x_1 = \frac{7.1 +x_2 + x_3}{5} x_2 = \frac{ 5.5 +x_1 + x_3}{4} x_3 = \frac{5 +x_1 + x_2}{4}Si consideramos que el costo mínimo podría ser 2, el precio debería ser mayor x0 = [2,2,2]
itera = 0
x_1 = \frac{7.1 +(2) + (2)}{5} = 2.22 x_2 = \frac{ 5.5 +(2) + (2)}{4} = 2.375 x_3 = \frac{5 +(2) + (2)}{4} =2.25 errado = max|[2.22-2, 2.375-2, 2.25-2]| = max |[0.22, 0.375, 0.25]| = 0.375 X_1 = [2.22, 2.375, 2.25]|itera = 1
x_1 = \frac{7.1 +(2.375) + (2.25)}{5} = 2.345 x_2 = \frac{ 5.5 +(2.22) + (2.25)}{4} = 2.4925 x_3 = \frac{5 +(2.22) + (2.375)}{4} =2.39875 errado = max|[2.345-2.22, 2.4925-2.375, 2.39875 - 2.25]| = 0.14874999999999972 X_2 = [2.345, 2.4925, 2.39875]itera = 2
x_1 = \frac{7.1 +(2.4925) + (2.39875)}{5} = 2.39825 x_2 = \frac{ 5.5 +(2.345) + (2.39875)}{4} = 2.5609375 x_3 = \frac{5 +(2.345) + (2.4925)}{4} = 2.459375 errado = max|[2.39825 - 2.345, 2.5609375- 2.4925, 2.459375 - 2.39875]| = 0.06843749999999993 X_3 = [2.39825, 2.5609375, 2.459375 ]El error disminuye entre iteraciones, por lo que el método converge.
Los datos de las iteraciones usando el algoritmo son:
Matriz aumentada [[ 5. -1. -1. 7.1] [-1. 4. -1. 5.5] [-1. -1. 4. 5. ]] Pivoteo parcial: Pivoteo por filas NO requerido itera,[X],errado 0 [2. 2. 2.] 1.0 1 [2.22 2.375 2.25 ] 0.375 2 [2.345 2.4925 2.39875] 0.14874999999999972 3 [2.39825 2.5609375 2.459375 ] 0.06843749999999993 4 [2.4240625 2.58940625 2.48979687] 0.030421875000000043 5 [2.43584062 2.60346484 2.50336719] 0.014058593750000181 6 [2.44136641 2.60980195 2.50982637] 0.006459179687499983 7 [2.44392566 2.61279819 2.51279209] 0.0029962402343750583 8 [2.44511806 2.61417944 2.51418096] 0.001388874511718985 9 [2.44567208 2.61482476 2.51482437] 0.0006453167724607134 10 [2.44592983 2.61512411 2.51512421] 0.00029983517456066977 11 [2.44604966 2.61526351 2.51526348] 0.0001393951034542873 12 [2.4461054 2.61532829 2.51532829] 6.480845146183967e-05 numero de condición: 2.5274158815808474 respuesta con Jacobi [2.4461054 2.61532829 2.51532829]
Algoritmo en Python
# 1Eva_IIT2008_T1 Distribuidores de productos # Método de Jacobi import numpy as np def jacobi(A,B,tolera,X0,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) itera = 0 if vertabla==True: 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) itera = itera + 1 if vertabla==True: print(itera, xnuevo, 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) if vertabla==True: if pivoteado==0: print(' Pivoteo por filas NO requerido') else: print('AB') return(AB) # INGRESO A = [[ 5, -1, -1], [ -1, 4, -1], [ -1, -1, 4]] B = [7.1, 5.5,5] X0 = [2,2,2] tolera = 0.0001 iteramax = 100 # PROCEDIMIENTO # numero de condicion ncond = np.linalg.cond(A) # Pivoteo parcial por filas AB = pivoteafila(A,B,vertabla=True) n = len(A) A = AB[:,:n] B = AB[:,n] respuesta = jacobi(A,B,tolera,X0,vertabla=True) # SALIDA print('numero de condición:', ncond) print('respuesta con Jacobi') print(respuesta)