s1Eva_IT2015_T4 Lingotes metales

Ejercicio: 1Eva_IT2015_T4 Lingotes metales

Se plantea que cada lingote debe aportar una proporción xi al lingote nuevo a ser fundido.

Se dispone de los compuestos de cada lingote por filas:

compuesto = np.array([[ 20, 50, 20, 10],
                      [ 30, 40, 10, 20],
                      [ 20, 40, 10, 30],
                      [ 50, 20, 20, 10]])
proporcion = np.array([ 27, 39.5, 14, 19.5])

El contenido final de cada componente basado en los aportes xi de cada lingote para cada componente.

Ejemplo para los 27 gramos de oro

20x_1 + 30x_2+ 20x_3 + 50x_4 = 27

se realiza el mismo procedimiento para los otros tipos de metal.

50x_1 + 40x_2+ 40x_3 + 20x_4 = 39.5 20x_1 + 10x_2+ 10x_3 + 20x_4 = 14 10x_1 + 20x_2+ 30x_3 + 10x_4 = 19.5

Las ecuaciones se escriben en la forma matricial Ax=B

\begin{bmatrix} 20 && 30&& 20 &&50 \\ 50 && 40 && 40 && 20 \\ 20 && 10 && 10 && 20 \\ 10 && 20 && 30 && 10 \end{bmatrix} = \begin{bmatrix} x_1 \\ x_2 \\x_3 \\ x_4 \end{bmatrix} = \begin{bmatrix} 27 \\ 39.5 \\ 14 \\ 19.5 \end{bmatrix}

Para resolver se plantea la matriz aumentada

\begin{bmatrix} 20 && 30&& 20 &&50 && 27\\ 50 && 40 && 40 && 20 && 39.5\\ 20 && 10 && 10 && 20 && 14 \\ 10 && 20 && 30 && 10 && 19.5 \end{bmatrix}

se pivotea por filas la matriz:

\begin{bmatrix} 50 && 40 && 40 && 20 && 39.5\\ 20 && 30&& 20 &&50 && 27\\ 10 && 20 && 30 && 10 && 19.5 \\ 20 && 10 && 10 && 20 && 14 \end{bmatrix}

Para eliminar hacia adelante:

\begin{bmatrix} 50 && 40 && 40 && 20 && 39.5 \\ 20 - 50\frac{20}{50} && 30-40\frac{20}{50} && 20-40\frac{20}{50} && 50-20\frac{20}{50} && 27-39.5\frac{20}{50}\\ 10 - 50\frac{10}{50} && 20-40\frac{10}{50} && 30-40\frac{10}{50} && 10-20\frac{10}{50} && 19.5-39.5\frac{10}{50} \\ 20 - 50\frac{20}{50} && 10-40\frac{20}{50} && 10-40\frac{20}{50} && 20-20\frac{20}{50} && 14-39.5\frac{20}{50} \end{bmatrix}

continuando con el desarrollo:

Elimina hacia adelante
[[50.  40.  40.  20.  39.5]
 [ 0.  14.   4.  42.  11.2]
 [ 0.  12.  22.   6.  11.6]
 [ 0.  -6.  -6.  12.  -1.8]]
Elimina hacia adelante
[[ 50.  40.  40.      20.  39.5]
 [  0.  14.   4.      42.  11.2]
 [  0.   0.  18.5714 -30.   2. ]
 [  0.   0.  -4.2857  30.   3. ]]
Elimina hacia adelante
[[ 50.  40.  40.      20.     39.5   ]
 [  0.  14.   4.      42.     11.2   ]
 [  0.   0.  18.5714 -30.      2.    ]
 [  0.   0.   0.      23.0769  3.4615]]
Elimina hacia adelante
[[ 50.  40.   40.      20.     39.5  ]
 [  0.  14.    4.      42.     11.2  ]
 [  0.   0.   18.5714 -30.      2.   ]
 [  0.   0.    0.      23.0769  3.4615]]
Elimina hacia atras
[[ 1.    0.    0.    0.    0.25]
 [ 0.    1.    0.    0.    0.25]
 [ 0.    0.    1.   -0.    0.35]
 [ 0.    0.    0.    1.    0.15]]
el vector solución X es:
[[0.25]
 [0.25]
 [0.35]
 [0.15]]
verificar que A.X = B
[[39.5]
 [27. ]
 [19.5]
 [14. ]]

Las proporciones de cada lingote a usar para el nuevo lingote que cumple con lo solicitado son:

[0.25, 0.25, 0.35, 0.15]


Algoritmo en python

usado para la solución es:

# Método de Gauss-Jordan ,
# Recibe las matrices A,B
# presenta solucion X que cumple: A.X=B
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
A = np.array([[ 50., 40, 40, 20],
              [ 20., 30, 20, 50],
              [ 10., 20, 30, 10],
              [ 20., 10, 10, 20]])
B1 = np.array([ 39.5, 27, 19.5, 14])

B = np.transpose([B1])

# PROCEDIMIENTO
casicero = 1e-15 # 0
AB = np.concatenate((A,B),axis=1)
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]

print('matriz aumentada: ')
print(AB)
# Gauss elimina hacia adelante
# tarea: verificar términos cero
for i in range(0,n,1):
    pivote = AB[i,i]
    adelante = i+1 
    for k in range(adelante,n,1):
        if (np.abs(pivote)>=casicero):
            factor = AB[k,i]/pivote
            AB[k,:] = AB[k,:] - factor*AB[i,:]
    print('Elimina hacia adelante')
    print(AB)

# Gauss-Jordan elimina hacia atras
ultfila = n-1
ultcolumna = m-1
for i in range(ultfila,0-1,-1):
    # Normaliza a 1 elemento diagonal
    AB[i,:] = AB[i,:]/AB[i,i]
    pivote = AB[i,i] # uno
    # arriba de la fila i
    atras = i-1 
    for k in range(atras,0-1,-1):
        if (np.abs(AB[k,i])>=casicero):
            factor = pivote/AB[k,i]
            AB[k,:] = AB[k,:]*factor - AB[i,:]
        else:
            factor= 'division para cero'
print('Elimina hacia atras')
print(AB)

X = AB[:,ultcolumna]

# Verifica resultado
verifica = np.dot(A,X)

# SALIDA
print('el vector solución X es:')
print(np.transpose([X]))

print('verificar que A.X = B')
print(np.transpose([verifica]))

Tarea: Revisar sobre la última pregunta.

s1Eva_IT2015_T2 Salida cardiaca

Ejercicio: 1Eva_IT2015_T2 Salida cardiaca

Solución presentada como introducción al tema de interpolación y solución de sistemas de ecuaciones.
No realiza el literal c, no se desarrolla el tema de integrales.

Note que el desarrollo del tema permite aumentar el grado del polinomio de interpolación, lo que se afecta al tamaño del sistema de ecuaciones (matriz).

Los valores obtenidos con la solución propuesta son:

solución para X por Gauss-Seidel
[[-0.42175867]
 [ 0.15610893]
 [ 0.02736763]]
verifica que A.X = B:
[[ -7.02831455e-05]
 [  1.50012970e+00]
 [  3.20000000e+00]]
polinomio interpolación, con puntos:  3
0.0273676337623498*x**2 + 0.156108926206362*x - 0.421758670607596
>>>

La gráfica para observar los datos experimentales es:

La gráfica con polinomio de interpolación de grado 2, con tres puntos:

instrucciones del problema para la solución por partes en python:

# 1ra Evaluación I Término 2015
# Tema 2. Flujo de sangre en corazón
# Tarea: parte c), no se ha realizado el áre bajo curva
#        falta calcular salida cardiaca.
import numpy as np
import matplotlib.pyplot as plt

# Gráfica de datos experimentales:
t = np.array([2,6,9,12,15,18])
y = np.array([0,1.5,3.2,4.1,3.4,2.0])

# SALIDA
plt.plot(t,y)
plt.title('datos del experimento: t vs concentración ')
plt.show()

# Sistema de ecuaciones para aproximar a polinomio grado 2
# para grado dos usa tres puntos,
# por ejemplo usando el punto [ 2, 0] del experimento
# a + b*(2) + c*(2**2) =  0
A = np.array([[1,2,2**2],
              [1,6,6**2],
              [1,9,9**2]])
B = np.array([[0],
              [1.5],
              [3.2]])

tolera = 0.0001
X = np.zeros(len(B), dtype = float)

# usando numpy para solucion de matrices
# Xnp = np.linalg.solve(A,B)
# print('solución para A.X=B con numpy')
# print(Xnp)

# algoritmo Gauss-Seidel
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 (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)
    itera = itera + 1
# Vector en columna
X =  np.transpose([X])
# No converge
if (itera>iteramax):
    X=0

Xgs = X

# Metodo numérico Gauss_Seidel
verifica = np.dot(A,Xgs)
print('solución para X por Gauss-Seidel')
print(Xgs)

# verificar resultado
print('verifica que A.X = B: ')
print(verifica)

# Observar interpolacion con polinomio creado
pt = lambda t: Xgs[0,0]+ Xgs[1,0]*t + Xgs[2,0]*t**2

ti = np.linspace(2,18,501)
pti = pt(ti)

plt.plot(ti,pti, label = 'interpolacion')
plt.plot(t,y,'*', label = 'datos experimento')
plt.title('interpolación con polinomio')
plt.legend()
plt.show()

# polinomio en sympy
import sympy as sp
x = sp.Symbol('x')
polinomio = 0
k = len(Xgs)
for i in range(0,k,1):
    polinomio = polinomio + Xgs[i,0]*x**i
print('polinomio interpolación, con puntos: ', k) 
print(polinomio)

s2Eva_IT2012_T3_MN EDO Taylor 2 Contaminación de estanque

Ejercicio: 2Eva_IT2012_T3_MN EDO Taylor 2 Contaminación de estanque

La ecuación a resolver con Taylor es:

s'- \frac{26s}{200-t} - \frac{5}{2} = 0

Para lo que se plantea usar la primera derivada:

s'= \frac{26s}{200-t}+\frac{5}{2}

con valores iniciales de s(0) = 0, h=0.1

La fórmula de Taylor para tres términos es:

s_{i+1}= s_{i}+s'_{i}h + \frac{s''_{i}}{2}h^2 + error

Para el desarrollo se compara la solución con dos términos, tres términos y Runge Kutta.

1. Solución con dos términos de Taylor

Iteraciones

i = 0, t0 = 0, s(0)=0

s'_{0}= \frac{26s_{0}}{200-t_{0}}+\frac{5}{2} = \frac{26(0)}{200-0}+\frac{5}{2} = \frac{5}{2} s_{1}= s_{0}+s'_{0}h = 0+ \frac{5}{2}*0.1= 0.25

t1 =  t0+h = 0+0.1 = 0.1

i=1


s'_{1}= \frac{26s_{1}}{200-t_{1}}+\frac{5}{2} = \frac{26(0.25)}{200-0.1}+\frac{5}{2} = 2.5325 s_{2}= s_{1}+s'_{1}h = 0.25 + (2.5325)*0.1 = 0.5032

t2 =  t1+h = 0.1+0.1 = 0.2

i=2,

resolver como tarea


2. Resolviendo con Python

estimado
 [xi,yi Taylor,yi Runge-Kutta, diferencias]
[[ 0.0  0.0000e+00  0.0000e+00  0.0000e+00]
 [ 0.1  2.5000e-01  2.5163e-01 -1.6258e-03]
 [ 0.2  5.0325e-01  5.0655e-01 -3.2957e-03]
 [ 0.3  7.5980e-01  7.6481e-01 -5.0106e-03]
 [ 0.4  1.0197e+00  1.0265e+00 -6.7714e-03]
 [ 0.5  1.2830e+00  1.2916e+00 -8.5792e-03]
 [ 0.6  1.5497e+00  1.5601e+00 -1.0435e-02]
 [ 0.7  1.8199e+00  1.8322e+00 -1.2339e-02]
 [ 0.8  2.0936e+00  2.1079e+00 -1.4294e-02]
 [ 0.9  2.3710e+00  2.3873e+00 -1.6299e-02]
 [ 1.0  2.6519e+00  2.6703e+00 -1.8357e-02]
 [ 1.1  2.9366e+00  2.9570e+00 -2.0467e-02]
 [ 1.2  3.2250e+00  3.2476e+00 -2.2632e-02]
 [ 1.3  3.5171e+00  3.5420e+00 -2.4853e-02]
 [ 1.4  3.8132e+00  3.8403e+00 -2.7129e-02]
 [ 1.5  4.1131e+00  4.1426e+00 -2.9464e-02]
 [ 1.6  4.4170e+00  4.4488e+00 -3.1857e-02]
 [ 1.7  4.7248e+00  4.7592e+00 -3.4310e-02]
 [ 1.8  5.0368e+00  5.0736e+00 -3.6825e-02]
 [ 1.9  5.3529e+00  5.3923e+00 -3.9402e-02]
 [ 2.0  5.6731e+00  5.7152e+00 -4.2043e-02]]
error en rango:  0.04204310894163932


2. Algoritmo en Python

# EDO. Método de Taylor 3 términos 
# estima la solucion para muestras espaciadas h en eje x
# valores iniciales x0,y0
# entrega arreglo [[x,y]]
import numpy as np

def edo_taylor2t(d1y,x0,y0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,2),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0]
    x = x0
    y = y0
    for i in range(1,tamano,1):
        y = y + h*d1y(x,y) # + ((h**2)/2)*d2y(x,y)
        x = x+h
        estimado[i] = [x,y]
    return(estimado)

def rungekutta2(d1y,x0,y0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,2),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0]
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1 = h * d1y(xi,yi)
        K2 = h * d1y(xi+h, yi + K1)

        yi = yi + (K1+K2)/2
        xi = xi + h
        
        estimado[i] = [xi,yi]
    return(estimado)

# PROGRAMA PRUEBA
# 2Eva_IIT2016_T3_MN EDO Taylor 2, Tanque de agua

# INGRESO.
# d1y = y' = f, d2y = y'' = f'
d1y = lambda x,y: 26*y/(200-x)+5/2
x0 = 0
y0 = 0
h = 0.1
muestras = 20

# PROCEDIMIENTO
puntos = edo_taylor2t(d1y,x0,y0,h,muestras)
xi = puntos[:,0]
yi = puntos[:,1]

# Con Runge Kutta
puntosRK2 = rungekutta2(d1y,x0,y0,h,muestras)
xiRK2 = puntosRK2[:,0]
yiRK2 = puntosRK2[:,1]

# diferencias
diferencias = yi-yiRK2
error = np.max(np.abs(diferencias))
tabla = np.copy(puntos)
tabla = np.concatenate((puntos,np.transpose([yiRK2]),
                        np.transpose([diferencias])),
                       axis = 1)

# SALIDA
np.set_printoptions(precision=4)
print('estimado[xi,yi Taylor,yi Runge-Kutta, diferencias]')
print(tabla)
print('error en rango: ', error)

# Gráfica
import matplotlib.pyplot as plt
plt.plot(xi[0],yi[0],'o',
         color='r', label ='[x0,y0]')
plt.plot(xi[0:],yi[0:],'-',
         color='g',
         label ='y Taylor 2 términos')
plt.plot(xiRK2[0:],yiRK2[0:],'-',
         color='blue',
         label ='y Runge-Kutta 2Orden')
plt.axhline(y0/2)
plt.title('EDO: Taylor 2T vs Runge=Kutta 2Orden')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()

Usando Taylor con 3 términos

estimado
 [xi,        yi,        d1yi,      d2yi      ]
[[0.         0.         2.5        0.325     ]
 [0.1        0.251625   2.53272761 0.32958302]
 [0.2        0.50654568 2.56591685 0.33423301]
 [0.3        0.76480853 2.59957447 0.33895098]
 [0.4        1.02646073 2.63370731 0.34373796]
 [0.5        1.29155015 2.66832233 0.348595  ]
 [0.6        1.56012536 2.70342658 0.35352316]
 [0.7        1.83223563 2.73902723 0.35852351]
 [0.8        2.10793097 2.77513155 0.36359715]
 [0.9        2.38726211 2.81174694 0.36874519]
 [1.         2.67028053 2.84888087 0.37396876]
 [1.1        2.95703846 2.88654098 0.37926901]
 [1.2        3.24758891 2.92473497 0.3846471 ]
 [1.3        3.54198564 2.96347069 0.39010422]
 [1.4        3.84028323 3.00275611 0.39564157]
 [1.5        4.14253705 3.04259931 0.40126036]
 [1.6        4.44880328 3.08300849 0.40696184]
 [1.7        4.75913894 3.12399199 0.41274727]
 [1.8        5.07360187 3.16555827 0.41861793]
 [1.9        5.39225079 3.2077159  0.42457511]
 [2.         5.71514526 0.         0.        ]]

 

s1Eva_IT2012_T2_MN Modelo Leontief

Ejercicio: 1Eva_IT2012_T2_MN Modelo Leontief

Planteamiento

X – TX = D

A(I-T) = D

(I-T)X = D

para el algoritmo:

A = I – T

B = D


Algoritmo en Python

Resultados del algoritmo

respuesta X: 
[[158.56573701]
 [288.73225044]
 [323.87373581]]
verificar A.X=B: 
[[ 79.99999997]
 [139.99999998]
 [200.        ]]
>>> itera
8
>>>

Instrucciones en Python

# Método de Gauss-Seidel
# solución de sistemas de ecuaciones
# por métodos iterativos

import numpy as np

# INGRESO
T = np.array([[0.40, 0.03, 0.02],
              [0.06, 0.37, 0.10],
              [0.12, 0.15, 0.19]])
A = np.identity(3) - T

B = np.array([80.0, 140.0, 200.0],dtype=float)

X0 = np.array([200.0,200.0,200.0])

tolera = 0.00001
iteramax = 100

# PROCEDIMIENTO
# Gauss-Seidel
tamano = np.shape(A)
n = tamano[0]
m = tamano[1]
#  valores iniciales
X = np.copy(X0)
diferencia = np.ones(n, dtype=float)
errado = 2*tolera

itera = 0
while not(errado<=tolera or itera>iteramax):
    # por fila
    for i in range(0,n,1):
        # por columna
        suma = 0 
        for j in range(0,m,1):
            # excepto diagonal de A
            if (i!=j): 
                suma = suma-A[i,j]*X[j]
        
        nuevo = (B[i]+suma)/A[i,i]
        diferencia[i] = np.abs(nuevo-X[i])
        
        X[i] = nuevo
    errado = np.max(diferencia)
    itera = itera + 1

# Respuesta X en columna
X = np.transpose([X])

# revisa si NO converge
if (itera>iteramax):
    X=0
# revisa respuesta
verifica = np.dot(A,X)

# SALIDA
print('respuesta X: ')
print(X)
print('verificar A.X=B: ')
print(verifica)

s1Eva_IIT2011_T2 Sistema de Ecuaciones, diagonal dominante

Ejercicio: 1Eva_IIT2011_T2 Sistema de Ecuaciones, diagonal dominante

1. Desarrollo analítico

1.1 Solución iterativa usando el médodo 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_0 \\ x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 1 \\ 6 \\ -26 \end{bmatrix}

para usarla en el algoritmo se intercambian filas, pivoteo, buscando hacerla diagonal dominante:

\begin{bmatrix} 7 & 1 & 1 \\ -3 & 7 &-1\\ -2 & 5 & 9 \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 6 \\ -26 \\1 \end{bmatrix}

las ecuaciones para el algoritmo se obtienen despejando una variable diferente en cada ecuación.

\begin{cases}x=(6-y-z)/(7)\\y=(-26+3x+z)/(7)\\z=(1+2x-5y)/(9)\end{cases}

Dado el vector inicial X(0) = [0, 0, 0], y una tolerancia de 0.0001, e desarrollan al menos 3 iteraciones:

Nota. El super índice entre paréntesis denota el número de la iteración. El subíndice denota el elemento del vector. 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…

Iteración 1

\begin{cases}x=(6-(0)-(0))/(7)\\y=(-26+3(0)+(0))/(7)\\z=(1+2(0)-5(0))/(9)\end{cases}
X(1) = [6/7, -26/7, 1/9] 
    = [0.8571, -3,7142, 0.1111]
diferencia(1) = |X(1) - X(0)|
             = |[0.8571,-3,7142,0.1111] - [0,0,0]|
             = |[0.8571,-3,7142,0.1111]|
errado(1) = maximo|[0.8571,-3,7142,0.1111]|
         = 26/7 ≅ 3.71, 
error es más alto que lo tolerado

Iteración 2

\begin{cases}x=(6-(-26/7)-(1/9))/(7)\\y=(-26+3(6/7)+(1/9))/(7)\\z=(1+2(6/7)-5(-26/7))/(9)\end{cases}
X(2) = [1.3718, -3.3310 ,  2,3650]
diferencia(2) = |X(2) - X(1)|
             = |[1.3718, -3.3310 , 2,3650] - [0.8571, -3,7142, 0.1111]| 
             = |[0.5146, 0.3832, 2.2539]|
errado(2) ≅ 2.2538, 
el error disminuye, pero es más alto que lo tolerado

Iteración 3

\begin{cases}x=(6-(-3.3310)-( 2,3650))/(7)\\y=(-26+3x+z)/(7)\\z=(1+2x-5y)/(9)\end{cases}
X(3) = [0.9951, -2.7884, 2.2665]
diferencia(3) = |[0.9951, -2.7884, 2.2665] - [1.3718, -3.3310 ,  2,3650]|
             = |[-0.3767, 0.5425, -0.0985]|
errado(3) ≅ 0.5425, 
el error disminuye, pero es más alto que lo tolerado

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:

Xi, errado
[ 0.85714286 -3.71428571  0.11111111] 3.71428571429
[ 1.37188209 -3.33106576  2.36507937] 2.25396825397
[ 0.99514091 -2.78846777  2.26656589] 0.542597991578
[ 0.93170027 -2.96400162  1.8814023 ] 0.385163589245
[ 1.0117999  -3.04621384  1.96482318] 0.0834208883016
[ 1.01162724 -2.99996816  2.02829656] 0.0634733731283
[ 0.99595309 -2.99097453  2.00256614] 0.0257304176216
[ 0.99834406 -3.0013678   1.99408654] 0.0103932672897
[ 1.00104018 -3.00155447  2.0003919 ] 0.00630536414945
[ 1.00016608 -2.99949822  2.00109475] 0.00205624814582
[ 0.99977193 -2.99977243  1.99975814] 0.00133660433338
[ 1.00000204 -3.0001323   1.99982289] 0.000359867509479
[ 1.0000442  -3.00002443  2.00007395] 0.000251063280064
[ 0.99999292 -2.99997049  2.00002339] 5.39347670632e-05
iteraciones:  13

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.

Tarea: convierta el algoritmo a una función, así es más sencillo usar en otros ejercicios.
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.

# 1ra Evaluación II Término 2011
# Tema 2. Sistema ecuaciones Jacobi
import numpy as np

# INGRESO
A = np.array([[ 7.0, 1, 1],
              [-3.0, 7,-1],
              [-2.0, 5, 9]])

B = np.array([6.0, -26.0, 1.0])

X = np.array([0.0, 0.0,  0.0],dtype=float)
tolera = 0.0001
iteramax = 100

# PROCEDIMIENTO
tamano = np.shape(A)
n = tamano[0]
m = tamano[1]
Xi1 = np.zeros(n, dtype=float)
errado = 2*tolera
itera = 0
while not(errado<=tolera or itera>=iteramax):
    i = 0
    while not(i>=n):
        j = 0
        nuevo = B [i]
        while not(j>=m):
            if (i!=j):
                nuevo = nuevo - A[i,j]*X[j]
            j = j+1
        Xi1[i] = nuevo/A[i,i]
        i = i+1
    diferencia = np.abs(X - Xi1)
    errado = np.max(diferencia)
    X = np.copy(Xi1)
    itera = itera +1
    print(Xi1, errado)

# SALIDA
print('iteraciones: ', itera-1)

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

s1Eva_IT2011_T3_MN Precios unitarios en factura, k

Ejercicio: 1Eva_IT2011_T3_MN Precios unitarios en factura, k

Las ecuaciones basadas en las sumas de cantidad.preciounitario representan el valor pagado en cada factura.

Siendo Xi el precio unitario de cada material:

2x_1 + 5x_2 + 4x_3 = 35 3x_1 + 9x_2 + 8x_3 = k 5x_1 + 3x_2 + x_3 = 17

se escriben en la forma matricial Ax=B

\begin{bmatrix} 2 && 5 && 4 \\ 3 && 9 && 8 \\ 5 && 3 && 1 \end{bmatrix}.\begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix}=\begin{bmatrix} 35 \\ k \\ 17 \end{bmatrix}

luego se escribe la matriz aumentada:

\begin{bmatrix} 2 && 5 && 4 && 35\\ 3 && 9 && 8 && k\\ 5 && 3 && 1 && 17\end{bmatrix}

se pivotea por filas buscando tener una matriz diagonal dominante,

\begin{bmatrix} 5 && 3 && 1 && 17 \\ 3 && 9 && 8 && k\\ 2 && 5 && 4 && 35\\\end{bmatrix}

Luego se usa el procedimiento de eliminación hacia adelante,

\begin{bmatrix} 5 && 3 && 1 && 17 \\ 3-5\frac{3}{5} && 9-3\frac{3}{5} && 8-1\frac{3}{5} && k-17\frac{3}{5} \\ 2-5\frac{2}{5} && 5-3\frac{2}{5} && 4-1\frac{2}{5} && 35-17\frac{2}{5} \end{bmatrix} \begin{bmatrix} 5 && 3 && 1 && 17 \\ 0 && \frac{36}{5} && \frac{37}{5} && k-\frac{51}{5} \\ 0 && \frac{19}{5} && \frac{18}{5} && \frac{141}{5} \end{bmatrix} \begin{bmatrix} 5 && 3 && 1 && 17 \\ 0 && 36 && 37 && 5k-51 \\ 0 && 19 && 18 && 141 \end{bmatrix} \begin{bmatrix} 5 && 3 && 1 && 17 \\ 0 && 36 && 37 && 5k-51 \\ 0 && 19-36\frac{19}{36} && 18-37\frac{19}{36} && 141-(5k-51)\frac{19}{36} \end{bmatrix} \begin{bmatrix} 5 && 3 && 1 && 17 \\ 0 && 36 && 37 && 5k-51 \\ 0 && 0 && \frac{-55}{36} && \frac{6045-5k}{36} \end{bmatrix}

multiplicando la última fila por 36,

\begin{bmatrix} 5 && 3 && 1 && 17 \\ 0 && 36 && 37 && 5k-51 \\ 0 && 0 && -55 && 6045-5k \end{bmatrix}

con lo que se pueden obtener cada precio unitario en función de k.
Como variante, se continua siguiendo el procedimieno de Gauss, dejando como tarea el uso de Gauss-Jordan

-55x_3 = 6045-5k x_3 = -\frac{6045-5k}{55} 36 x_2 + 37 x_3 = 5k-51 x_2 = \frac{1}{36}(5k-51 - 37 x_3) x_2 = \frac{1}{36} \Big( 5k-51 - 37 \big(-\frac{6045-5k}{55}\big) \Big) 5x_1 + 3 x_2 +x_3 = 17 x_1 = \frac{1}{5} \Big[ 17 - 3 x_2 - x_3 \Big] x_1 = \frac{1}{5} \Big[17-3\frac{1}{36} \Big( 5k-21 - 37 \big(-\frac{6045-5k}{55}\big) \Big) - \Big( -\frac{6045-5k}{55} \Big) \Big]

para luego simplificar las expresiones (tarea).

En el literal c se indica que el valor de k es 65, con lo que se requiere sustituir en la solución el valor de K para encontrar los precios unitarios.

\begin{bmatrix} 5 && 3 && 1 && 17 \\ 3 && 9 && 8 && 65\\ 2 && 5 && 4 && 35\\\end{bmatrix}

Se encuentra que:

el vector solución X es:

[[-0.18181818]
 [ 5.18181818]
 [ 2.36363636]]

Lo que muestra que debe existir un error en el planteamiento del enunciado, considerando que los precios NO deberían ser negativos como sucede con el primer precio unitario de la respuesta.

que es lo que suponemos ser trata de corregir en el literal d, al indicar que se cambie en la matriz el valor de 5 por 5.1. Los resultados en éste caso son más coherentes con el enunciado. Todas las soluciones son positivas.

A = np.array([[ 5.1, 3  , 1],
              [ 3. , 9  , 8],
              [ 2. , 5.1, 4]])
B1 = np.array([ 17, 65, 35])

el vector solución X es:
[[0.33596838]
 [3.88669302]
 [3.62648221]]

El error relativo de los precios encontrados entre las ecuaciones planteadas es:

diferencia = [0.335968-0.181818,
              3.886693-5.181818, 
              3.626482-2.363636]
           = [0.154150, -1.295125, 1.262845]
error(dolares) = max|diferencia| = 1.295125
Por las magnitudes de los precios, el error se aprecia
usando el error relativo referenciado 
al mayor valor de la nueva solución
error relativo = 1.295125/3.886693 = 0.333220
es decir de aproximadamente 33%

Para revisar otra causa del error se analiza el número de condición de la matriz:

>>> A
array([[5.1, 3. , 1. ],
       [3. , 9. , 8. ],
       [2. , 5.1, 4. ]])
>>> np.linalg.cond(A)
60.28297696795716

El número de condición resulta lejano a 1, por lo que para éste problema:
pequeños cambios en la matriz de entrada producen grandes cambios en los resultados.

por ejemplo: un 0.1/5= 0.02 que es un 2% de variación en la entrada produce un cambio del 33% en el resultado.

s1Eva_IT2011_T1_MN Fondo de inversión

Ejercicio: 1Eva_IT2011_T1_MN Fondo de inversión

Se desarrolla para la función:

C(t)=Ate^{-t/3}

siguiento las instrucciones por partes, se obtienen los siguientes resultados:
los valores resultantes:

derivada de la función: 
-x*exp(-x/3)/3 + exp(-x/3)
valor maximo en : 
3.0000001192092896
A para un millon: 
906093.94282
el valor de t para la meta es: 
11.0779035867

las instrucciones en python para observar la función son:

# 1ra Evaluación I Término 2011
# Tema 1. Fondo de Inversion
import numpy as np
import matplotlib.pyplot as plt

ft = lambda t: t*np.exp(-t/3)

a=0
b=20
tolera = 0.000001
muestras = 101
meta = 0.25
# PROCEDIMIENTO
# Observar la función entre [a,b]
ti = np.linspace(a,b,muestras)
fti = ft(ti)

# Salida
# Gráfica
plt.plot(ti,fti, label='f(t)')
plt.axhline(meta, color = 'r')
plt.axhline(0, color = 'k')
plt.legend()
plt.show()


Para desarrollar el literal a), donde se busca el valor máximo, usando la derivada de la función cuando existe el cruce por cero.
Para la derivada se usa la forma simbólica de la función, que se convierte a forma numérica lambda para evaluarla de forma más fácil.

# Literal a) usando derivada simbólica
import sympy as sp
x = sp.Symbol('x')
fxs = x*sp.exp(-x/3)
dfxs = fxs.diff(x,1)

# convierte la expresión a lambda
dfxn = sp.utilities.lambdify(x,dfxs,'numpy')
dfxni = dfxn(ti)
print('derivada de la función: ')
print(dfxs)
# Gráfica de la derivada.
plt.plot(ti,dfxni, label='df(t)/dt')
plt.axhline(0, color = 'k')
plt.legend()
plt.show()

derivada de la función: 
-x*exp(-x/3)/3 + exp(-x/3)

Se busca la raiz con algún método, por ejemplo bisección.

# Busca el máximo en dfxni
def biseccion(funcionx,a,b,tolera):
    fa = funcionx(a)
    fb = funcionx(b)
    tramo = np.abs(b-a)
    while (tramo>=tolera):
        c = (a+b)/2
        fc = funcionx(c)
        cambia = np.sign(fa)*np.sign(fc)
        if (cambia<0):
            b = c
            fb = fc
        else:
            a = c
            fa = fc
        tramo = np.abs(b-a)
    return(c)

# usa función para encontrar el máximo
raizmax = biseccion(dfxn, a, b, tolera)
verifica =  dfxn(raizmax)
print('valor maximo en : ')
print(raizmax)

# que el máximo sea un millon
tmax = raizmax
A = (1000000)/ft(tmax)
print('A para un millon: ')
print(A)
valor maximo en : 3.0000001192092896
A para un millon: 906093.94282

Para el literal b, se busca la laiz usando el metodo de Newton-Raphson como se indica en el enunciado.
En la función nueva se usa el valor de A encontrado, y la meta establecida.
Se obtiene la misa derivada del problema anterio multiplicada por A, por ser solo un factor que multiplica a la función original. El valor de meta es una constante, que se convierte en cero al derivar.

# literal b), buscar cumplir meta de 0.25 millones
def newtonraphson(funcionx, fxderiva, c, tolera):
    tramo = abs(2*tolera)
    while (tramo>=tolera):
        xnuevo = c - funcionx(c)/fxderiva(c)
        tramo = abs(xnuevo-c)
        c = xnuevo
    return(c)

ft1 = lambda t: A*t*np.exp(-t/3) - 250000
# usar método de newton,
# puede usar la misma derivada multiplicada por A
dft1s = A*(fxs.diff(x,1))
dft1 = sp.utilities.lambdify(x,dft1s,'numpy')
c = 10

raiz4 = newtonraphson(ft1, dft1, c, tolera)
ft1i = ft1(ti)

print('el valor de t para la meta es: ')
print(raiz4)

# Gráfica
plt.plot(ti,ft1i, label = 'f(t) con A>1')
plt.axhline(meta, color = 'r')
plt.axhline(0, color = 'k')
plt.axvline(raiz4, color = 'm')
plt.legend()
plt.show()

el valor de t para la meta es: 11.0779035867

s1Eva_IT2011_T1 Encontrar α en integral

Ejercicio: 1Eva_IT2011_T1 Encontrar α en integral

Desarrollo Analítico

Se iguala la ecuación al valor buscado = 10, y se resuelve

\int_{\alpha}^{2\alpha} x e^{x}dx = 10

siendo: μ = x , δv = ex, δu = δx , v = ex

\int u dv = uv - \int v \delta u xe^x \Big|_{\alpha}^{2 \alpha} - \int_{\alpha}^{2\alpha} e^{x}dx - 10 = 0 2\alpha e^{2 \alpha} -\alpha e^{\alpha} - (e^{2\alpha} - e^{\alpha}) - 10 = 0 (2\alpha-1)e^{2 \alpha}+ (1-\alpha) e^{\alpha} - 10 = 0

la función a usar en el método es

f(\alpha) = (2\alpha-1)e^{2 \alpha}+ (1-\alpha)e^{\alpha} -10

Se obtiene la derivada para el método de Newton Raphson

f'(\alpha) = 2e^{2 \alpha} + 2(2\alpha-1)e^{2 \alpha} - e^{\alpha} + (1-\alpha) e^{\alpha} f'(\alpha) = (2 + 2(2\alpha-1))e^{2 \alpha} +(-1 + (1-\alpha)) e^{\alpha} f'(\alpha) = 4\alpha e^{2 \alpha} -\alpha e^{\alpha}

la fórmula para el método de Newton-Raphson

\alpha_{i+1} = \alpha_i - \frac{f(\alpha)}{f'(\alpha)}

se prueba con α0 = 1, se evita el valor de cero por la indeterminación que se da por f'(0) = 0

iteración 1:

f(1) = (2(1)-1)e^{2(1)}+ (1-(1))e^{(1)} -10 f'(1) = 4(1) e^{2 (1)} -(1) e^{(1)} \alpha_{1} = 1- \frac{f(1)}{f'(1)}

\alpha_{1} = 1.0973
error = 0.0973

iteración 2:
f(2) = (2(1.0973)-1)e^{2(1.0973)}+ (1-(1.0973))e^{(1.0973)} -10

f'(2) = 4(1.0973) e^{2 (1.0973)} -(1.0973) e^{(1.0973)} \alpha_{2} = 1.0973 - \frac{f(1.0973)}{f'(1.0973)}

\alpha_{2} = 1.0853
error = 0.011941

iteración 3:
f(3) = (2(1.0853)-1)e^{2(1.0853)}+ (1-(1.0853))e^{(1.0853)} -10

f'(3) = 4(1.0853) e^{2 (1.0853)} -(1.0853) e^{(1.0853)} \alpha_{3} = 1.0853- \frac{f(1.0853)}{f'(1.0853)}

\alpha_{3} = 1.0851
error = 0.00021951

[  xi,          xnuevo,      f(xi),      f'(xi),       tramo ]
[  1.           1.0973      -2.6109      26.8379       0.0973]
[  1.0973e+00   1.0853e+00   4.3118e-01   3.6110e+01   1.1941e-02]
[  1.0853e+00   1.0851e+00   7.6468e-03   3.4836e+01   2.1951e-04]
[  1.0851e+00   1.0851e+00   2.5287e-06   3.4813e+01   7.2637e-08]
raiz:  1.08512526549

se obtiene el valor de la raíz con 4 iteraciones, con error de aproximación de 7.2637e-08

Desarrollo con Python

s1Eva_IT2010_T2_MN Uso de televisores

Ejercicio: 1Eva_IT2010_T2_MN Uso de televisores

El enunciado indica encontrar el máximo y luego el mínimo, por lo que la curva bajo análisis es la derivada de la función dp(x)/dx.

Adicionalmente, para encontrar los puntos se requiere usar el método de Newton-Raphson que corresponden a las raíces de dp(x)/dx. La función bajo análisis ahora es la derivada y para el método se su la derivada: d2p(x)/dx2.

Al usar el computador para las fórmulas, se usa la forma simbólica de la función p(x), para obtener dpx y d2px.

primera derivada: 
-3.13469387755102*x*exp(-8*x/7) + 2.74285714285714*exp(-8*x/7) + 13.7142857142857*exp(-24*x/7)*sin(12*x/7) - 6.85714285714286*exp(-24*x/7)*cos(12*x/7)
segunda derivada: 
(3.58250728862974*x - 6.26938775510204 - 35.265306122449*exp(-16*x/7)*sin(12*x/7) + 47.0204081632653*exp(-16*x/7)*cos(12*x/7))*exp(-8*x/7)

La gráfica requiere la evaluación las funciones, que por simplicidad de evaluación, su formas simbólicas se convierten a su forma ‘lambda’.

Con la gráfica se verifica que la raiz de dp(x)/dx (en naranja) pasa por el máximo y mínimo de p(x) (en azul).

que se obtienen con las siguientes instrucciones en python:

# 1ra Evaluación I Término 2010
# tema 2. encendido tv
# Tarea: aplicar el método de Newton-Raphson
# solo se muestra la función y sus derivadas 1 y 2
import numpy as np
import sympy as sp
import matplotlib.pyplot as plt

# función bajo análisis en forma simbólica
x = sp.Symbol('x')
pxs = (1/2.5)*(-10*sp.sin(12*x/7)*sp.exp(-24*x/7) + (48*x/7)*sp.exp(-8*x/7)+0.8)

# derivadas
dpxs = pxs.diff(x,1)
d2pxs = pxs.diff(x,2)

# SALIDA
print('primera derivada: ')
print(dpxs)
print('segunda derivada: ')
print(d2pxs)

# conversion a lambda
pxn = sp.utilities.lambdify(x,pxs, 'numpy')
dpxn = sp.utilities.lambdify(x,dpxs, 'numpy')
d2pxn = sp.utilities.lambdify(x,d2pxs, 'numpy')

# observar gráfica
a = 0
b = 4
muestras
tolera = 0.0001

xi = np.linspace(a,b, muestras)
pxi = pxn(xi)
dpxi = dpxn(xi)
d2pxi = d2pxn(xi)

# Gráfica
plt.plot(xi,pxi, label = 'pxi')
plt.plot(xi,dpxi, label = 'dpxi')
plt.plot(xi,d2pxi, label = 'd2pxi')
plt.axhline(0)
plt.legend()
plt.show()

# Tarea: encontrar la raiz de dpxn
# usando el método de Newton-Raphson

s1Eva_IT2010_T1_MN Demanda y producción sin,log

Ejercicio: 1Eva_IT2010_T1_MN Demanda y producción sin,log

Desarrollo Analítico

Para la demanda, el intervalo de existencia es [0,3]

demanda(t) = sin(t)

Para la oferta, el intervalo de existencia inicia en 1, limitado por la demanda [1,3]

oferta(t) = ln(t)

la oferta satisface la demanda cuando ambas son iguales

demanda(t) = oferta(t) sin(t) = ln(t)

por lo que el tiempo t se encuentra con algun método para determinar la raiz de:

sin(t) - ln(t) = 0 f(t) = sin(t) - ln(t)

Observe que las curvas de oferta y demanda se intersectan en el mismo punto en el eje x que la función f(t).

Use un método para encontrar el valor t que satisface la ecuación.


Algoritmo en Python

instrucciones para la gráfica, no para el algoritmo de búsqueda de raiz que es tarea.

# 1Eva_IT2010_T1_MN Demanda y producción sin,log
# Solo para analizar el problema
# Tarea: Añadir algoritmo de buscar raiz.
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# demanda
ad = 0
bd = 3
muestrasd = 31
# oferta
ao = 1
bo = 3
muestras0  = 21

demanda = lambda t: np.sin(t)
oferta = lambda t: np.log(t)
f = lambda t: demanda(t)-oferta(t)

# PROCEDIMIENTO
tid = np.linspace(ad,bd,muestrasd)
demandai = demanda(tid)

tio = np.linspace(ao,bo,muestras0)
ofertai = oferta(tio)

fi = f(tio)

# SALIDA
plt.plot(tid,demandai, label='demanda')
plt.plot(tio,ofertai, label ='oferta')
plt.plot(tio,fi,label='f(t)= demanda-oferta')
plt.axhline(0,color='black')
plt.axvline(2.2185, color = 'magenta')
plt.xlabel('tiempo')
plt.ylabel('unidades')
plt.legend()
plt.grid()
plt.show()