s1Eva_2024PAOI_T3 Tasas de natalidad de reemplazo en Ecuador

Ejercicio: 1Eva_2024PAOI_T3 Tasas de natalidad en países desarrollados

Los datos de la tabla que se requieren usar para la tabla son los años 1965, 1980, 1995 y 2010.

Promedio Hijos por mujer en Ecuador. INEC 29 Agosto 2019 [4]
tasa 7.32 4.89 3.50 2.79
año 1965 1980 1995 2010
k 0 3 6 9

literal a

Para cuatro datos o muestras el grado del polinomio que se puede plantear es 3. La ecuación se ordena siguiendo el modelo para la matriz de Vandermonde, primer coeficiente es del termino de mayor grado.

p(x) = ax^3+bx^2+cx+d

Por el orden de magnitud del eje x comparado con los valores de las constantes, se usan los valores del vector k que corresponde a los índices de la tabla. El valor del año se obtiene al hacer un cambio de escala que inicia en el año t0 = 1965 y tamaño de paso Δx = 5 años.

el sistema de ecuaciones usando cada punto será:

p(0) = a(0)^3+b(0)^2+c(0)+d = 7.32 p(3) = a(3)^3+b(3)^2+c(3)+d = 4.89 p(6) = a(6)^3+b(6)^2+c(6)+d = 3.50 p(9) = a(9)^3+b(9)^2+c(9)+d = 2.79

literal b

El sistema de ecuaciones para la Matriz de Vandermonde D y el vector de constantes B en la forma de matriz aumentada D|B:

\left[ \begin{array}{cccc | c} (0)^3 & (0)^2 & (0) & 1 & 7.32 \\ (3)^3 & (3)^2 & (3) &1 &4.89\\ (6)^3 & (6)^2 & (6) &1 & 3.50 \\ (9)^3 & (9)^2 & (9) & 1 & 2.79 \end{array} \right]

literal c

Resolviendo el sistema de ecuaciones con el algoritmo y usando la instrucción np.linalg.solve(A,B)

>>> np.linalg.solve(A,B)
array([-2.22222222e-03,  7.77777778e-02,
       -1.02333333e+00,  7.32000000e+00])

literal d

p(x) = -0.00222 x^3 + 0.07777 x^2 - 1.0233 x + 7.32

Para el ajuste de escala de años, se usa:

año = x*5 +1965

la gráfica usando el algoritmo, muestra que el polinomio pasa por los puntos seleccionados y existe un error entre los puntos que no participan en la construcción del polinomio.

tasa reemplazo colapso Ec v1

literal e

El año 2020 en x se escala como

2020 = x*5 +1965

x = 11

o revisando la tabla proporcionada en el ejercicio

p(11) = -0.00222 (11)^3 + 0.07777(11)^2 - 1.0233 (11)x + 7.32 = 2.5166

el error se calcula usando el valor medido para el año 2020

error = | 2.51- 2.35| = 0.16

literal f

p(x) = -0.00222 x^3 + 0.07777 x^2 - 1.0233 x + 7.32 = 2.51

para resolverlo se requiere usar cualquiera de los algoritmos de la unidad 2.

>>> import scipy.optimize as opt
>>> opt.bisect(px,11,30,xtol=0.001)
20.312713623046875
>>>

que corresponde al año = 20.31*5+1965 = 2066.55

Sin embargo la interpolación se usa para estimar valores dentro del intervalo de muestras. El concepto de extrapolación corresponde a otro capítulo. Sin embargo la pregunta se la considera para evaluar la comprensión de la unidad 2.


Algoritmo con Python

El resultado con el algoritmo se observa como

xi [0 3 6 9]
fi [7.32 4.89 3.5  2.79]
Matriz Vandermonde: 
[[  0.   0.   0.   1.]
 [ 27.   9.   3.   1.]
 [216.  36.   6.   1.]
 [729.  81.   9.   1.]]
los coeficientes del polinomio: 
[-2.22222222e-03  7.77777778e-02 
 -1.02333333e+00  7.32000000e+00]
Polinomio de interpolación: 
-0.00222222224*x**3 + 0.077777778*x**2 - 1.02333333*x + 7.32

 formato pprint
                       3                      2                            
- 0.002222222224*x  + 0.0777777778*x  - 1.023333333*x + 7.32
>>> px(11)
2.5166666666667066

Instrucciones en Python

# 1Eva_2024PAOI_T3 Tasas de natalidad en países desarrollados
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

# INGRESO
tasa = [7.32, 6.39, 5.58, 4.89,4.31,3.85,3.50,3.26,3.03,2.79,2.54,2.35]
anio = [1965, 1970, 1975, 1980,1985,1990,1995,2000,2005,2010,2015,2020]
k    = [   0,     1,   2,    3,   4,   5,   6,   7,   8,   9,  10,  11]
cual = [0,3,6,9]
tasamin = 2.1

# PROCEDIMIENTO
fi = np.take(tasa,cual)
xi = np.take(k,cual)

# El polinomio de interpolación
# muestras = tramos+1
muestras = 25

# PROCEDIMIENTO
# Convierte a arreglos numpy 
xi = np.array(xi)
B = np.array(fi)
n = len(xi)

# Matriz Vandermonde D
D = np.zeros(shape=(n,n),dtype =float)
for i in range(0,n,1):
    for j in range(0,n,1):
        potencia = (n-1)-j # Derecha a izquierda
        D[i,j] = xi[i]**potencia

# Aplicar métodos Unidad03. Tarea
# Resuelve sistema de ecuaciones A.X=B
coeficiente = np.linalg.solve(D,B)

# Polinomio en forma simbólica
x = sym.Symbol('x')
polinomio = 0
for i in range(0,n,1):
    potencia = (n-1)-i   # Derecha a izquierda
    termino = coeficiente[i]*(x**potencia)
    polinomio = polinomio + termino

# Polinomio a forma Lambda
# para evaluación con vectores de datos xin
px = sym.lambdify(x,polinomio)

# Para graficar el polinomio en [a,b]
a = np.min(xi)
b = np.max(xi)
xin = np.linspace(a,b,muestras)
yin = px(xin)
   
# SALIDA
print('xi', xi)
print('fi',fi)
print('Matriz Vandermonde: ')
print(D)
print('los coeficientes del polinomio: ')
print(coeficiente)
print('Polinomio de interpolación: ')
print(polinomio)
print('\n formato pprint')
sym.pprint(polinomio)

# Grafica
plt.plot(anio,tasa,'o')
plt.plot(xi*5+anio[0],fi,'o',color='red', label='[xi,fi]')
plt.plot(xin*5+anio[0],yin,color='red', label='p(x)')
plt.xlabel('xi')
plt.ylabel('fi')
plt.legend()
plt.axhline(0)
plt.axhline(2.1,linestyle='-.')
plt.xlabel('anio')
plt.ylabel('tasa')
plt.title(polinomio)
plt.show()

 

s1Eva_2024PAOI_T2 Temperatura en nodos de placa cuadrada

Placa Cuadrada Calentada Nodos01Ejercicio: 1Eva_2024PAOI_T2 Temperatura en nodos de placa cuadrada

literal a

El planteamiento de las ecuaciones según se cita en la solución iterativa es el promedio de los nodos adjacentes:

a = \frac{100+40+b+c}{4} b = \frac{40+20+a+d}{4} c = \frac{100+80+a+d}{4} d = \frac{80+20+c+b}{4}

reordenando las ecuaciones para su forma maticial:

4a -b -c = 100+40 4b -a -d= 40+20 4c -a -d= 100+80 4d -c -b= 80+20

literal b

la forma matricial del sistema de ecuaciones es:

\begin{bmatrix} 4 & -1 &-1 & 0 \\ -1 & 4 & 0 &-1 \\ -1 & 0 & 4 &-1 \\ 0 & -1 &-1 & 4 \end{bmatrix} \begin{bmatrix} a \\ b \\ c \\ d \end{bmatrix} = \begin{bmatrix} 140 \\ 60 \\ 180 \\ 100 \end{bmatrix}

matriz aumentada:

\left[ \begin{array} {cccc|c} 4 & -1 &-1 & 0 & 140\\ -1 & 4 & 0 &-1 & 60 \\ -1 & 0 & 4 &-1 & 180 \\ 0 & -1 &-1 & 4 & 100 \end{array}\right]

al revisar el procedimiento de pivoteo parcial por filas se muestra que ya se encuentra pivoteada la matriz. Por lo que no es necesario realizar cambios en las filas.

literal c

Las expresiones a partir de la matriz aumentada y pivoteada para usar en un método iterativo como Gauss-Seidel son:

a = 0.25(b+c+140) b = 0.25*(a+d+60) c = 0.25(a+d+180) d = 0.25(c+b+100)

el vector inicial debe considerar que las temperaturas serán medias entre los  bordes de la placa

X0 = [ (100+40)/2, (40+20)/2, (100+80)/2, (80+20)/2 ]

X0 = [ 70, 30, 90, 50 ]

literal d iteraciones

itera = 0

a = 0.25(30+90+140)= 65

b = 0.25*(65+50+60) = 43.75

c = 0.25(65+50+180) = 73.75

d = 0.25(73.75+43.75+100) = 54.375

error = max|[65-70, 43.75-30, 73.75-90, 54.375-50]|

error = max|[ -5, 13.75, -16.25, 4.375| = 16.25

X1 = [ 65, 43.75, 73.75, 54.375 ]

itera = 1

a = 0.25(43.75+73.75+140)= 64.375

b = 0.25*(64.375+54.375+60) = 44.687

c = 0.25(64.375+54.375+180) = 74.687

d = 0.25(74.687+44.687+100) = 54.843

error = max|[64.375-65, 44.687-43.75, 74.687-73.75, 54.843-54.375]|

error = max|[-0.625, 0.937, 0.937, 0.468| = 0.937

X2 = [ 64.375, 44.687, 74.687, 54.843 ]

itera = 2

a = 0.25(44.687+74.687+140)= 64.843

b = 0.25*(64.843 +54.843+60) = 44.921

c = 0.25(64.843+54.843+180) = 74.921

d = 0.25(74.921+44.921+100) = 54.960

error = max|[64.843-64.375, 44.921-44.687, 74.921-74.687, 54.960-54.843]|

error = max|[0.468, 0.234, 0.234, 0.117| = 0.468

X3 = [ 64.843, 44.921, 74.921, 54.960 ]

literal e

El método converge pues los errores en cada iteración disminuyen. Los valores obtenidos en el resultado son acorde a las temperaturas esperadas en los nodos.

Algoritmos con Python

Los resultados obtenidos son:

X 0 = [65.    43.75  73.75  54.375]
X 1 = [64.375   44.6875  74.6875  54.84375]
X 2 = [64.84375   44.921875  74.921875  54.9609375]
X 3 = [64.9609375  44.98046875 74.98046875 54.99023438]
X 4 = [64.99023438 44.99511719 74.99511719 54.99755859]
X 5 = [64.99755859 44.9987793  74.9987793  54.99938965]
X 6 = [64.99938965 44.99969482 74.99969482 54.99984741]
X 7 = [64.99984741 44.99992371 74.99992371 54.99996185]
X 8 = [64.99996185 44.99998093 74.99998093 54.99999046]
X 9 = [64.99999046 44.99999523 74.99999523 54.99999762]
respuesta X: 
[[64.99999046]
 [44.99999523]
 [74.99999523]
 [54.99999762]]
verificar A.X=B: 
[[139.99997139]
 [ 59.99999285]
 [179.99999285]
 [100.        ]]
>>>

el algoritmo usado es:

# 1Eva_2024PAOI_T2 Temperatura en nodos de placa cuadrada
# Método de Gauss-Seidel
import numpy as np
# INGRESO
A = np.array([[4,-1,-1,0],
              [-1,4,0,-1],
              [-1,0,4,-1],
              [0,-1,-1,4]],dtype=float)
B = np.array([140,60,180,100],dtype=float)
X0  = np.array([70,30,90,50],dtype=float)

tolera = 0.0001
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
    print('X',itera,'=',X)
    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_2024PAOI_T1 Vaciado de reservorio semiesférico

Ejercicio: 1Eva_2024PAOI_T1 Vaciado de reservorio semiesférico

literal a. Planteamiento

A partir de la solución expresión propuesta y simplificada del ejercicio, se reemplaza los valores de R, K, a y g que son constantes. Como se da el valor de t=(15 min)(60 s/min), se unifica la unidad de medida de tiempo a segundos pues g=9.6 m/s2.


-\frac{4}{3}R h^{3/2} + \frac{2}{5} h^{5/2} = -\frac{Ka}{\pi} t \sqrt{2g}
-\frac{4}{3}(3) h^{3/2} + \frac{2}{5} h^{5/2} = -\frac{(0.85)(0.01)}{\pi} (15)(60) \sqrt{2(9.8)}

la función para encontrar la raíz se reordena como f(h)=0

f(h) = -\frac{4}{3}(3) h^{3/2} + \frac{2}{5} h^{5/2} +\frac{(0.85)(0.01)}{\pi} (15)(60) \sqrt{2(9.8)} f(h) = -4h^{3/2} + 0.4 h^{5/2} +10.7805

tanque semiesferatanque semiesfera v2Usando la gráfica proporcionada del reservorio semiesférico, el valor de h no puede superar la altura R. Al vaciar el reservorio h=0.

Por lo que el intervalo a usar es [0,3]. Se verifica la existencia de una raíz al con la gráfica y Python.

b. Método de Newton Raphson

El método requiere la derivada f'(h)

f(h) = -4h^{3/2} + 0.4 h^{5/2} +10.7805 f'(h) = -6\frac{ h^{3/2}}{h} + \frac{h^{5/2}}{h} f'(h) = -6 h^{1/2}+ h^{3/2}

El valor inicial puede ser por ejemplo, el centro del intervalo:

x_0= \frac{a+b}{2} = \frac{0+3}{2} = 1.5

itera=0

f(1.5) = -4(1.5)^{3/2} + 0.4 (1.5)^{5/2} +10.7805= 4.5343 f'(1.5) = -6\frac{ 1.5^{3/2}}{1.5} + \frac{1.5^{5/2}}{1.5} = -5.5113 x_1= x_0 - \frac{f(1.5)}{f'(1.5)}= 1.5-\frac{4.5343}{-5.5113} =2.3227 error_1 =|x1-x0| = |2.3227-1.5| =0.8227

itera=1

f(2.3227) = -4(2.3227)^{3/2} + 0.4 (2.3227)^{5/2} +10.7805= -0.09019 f'(2.3227) = -6\frac{ 2.3227^{3/2}}{2.3227} + \frac{2.3227^{5/2}}{2.3227} = -5.6043 x_2= 2.3227-\frac{-0.09019}{-5.6043} = 2.3066 error_2 = |2.3066-2.3227| =0.01611

itera=2

f(2.3227) = -4(2.3066)^{3/2} + 0.4 (2.3066)^{5/2} +10.7805= -0.0007104 f'(2.3227) = -6\frac{ 2.3227^{3/2}}{2.3227} + \frac{2.3227^{5/2}}{2.3227} = -5.6093 x_3= 2.3227-\frac{-0.0007104}{-5.6093} = 2.3066 error_3 = |2.3066-2.3066| =0.000007241

literal c, convergencia

considerando tolerancia de 10-3 (milimétrica), solo serán necesarias 3 iteraciones.

El método converge al mostrarse que los errores disminuyen en cada iteración.

El resultado se encuentra dentro del intervalo y es x = 2.3066

Algoritmo con Python

resultados:

fh:
         ____          ____                   
        /  3          /  5                    
- 4.0*\/  h   + 0.4*\/  h   + 10.7805172327811
dfh:
         ____          ____
        /  3          /  5 
  6.0*\/  h     1.0*\/  h  
- ----------- + -----------
       h             h     
['xi', 'xnuevo', 'tramo']
[[1.5000e+00 2.3227e+00 8.2272e-01]
 [2.3227e+00 2.3066e+00 1.6118e-02]
 [2.3066e+00 2.3066e+00 7.2412e-06]]
raiz en:  2.306612665792577
con error de:  7.241218404896443e-06

>>> f(1.5)
4.534318388683996
>>> df(1.5)
-5.5113519212621505
>>> f(2.3227)
-0.09019959114247378
>>> df(2.3227)
-5.604354799446854
>>> f(2.3066)
7.104683901282272e-05
>>> df(2.3066)
-5.609349350102558

instrucciones:

# 1Eva_2024PAOI_T1 Vaciado de reservorio semiesférico
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

# INGRESO
h = sym.Symbol('h')
t = 15*60
R = 3 
g = 9.8 
area = 0.01
K = 0.85 
fh = -(4/3)*R*sym.sqrt(h**3)+(2/5)*sym.sqrt(h**5)+(K*area/np.pi)*np.sqrt(2*g)*(t)

# Grafica de f(h)
a = 0
b = 3
muestras = 15

# PROCEDIMIENTO
hi = np.linspace(a,b,muestras)
# para evaluación numérica con numpy
f = sym.lambdify(h,fh)
fi = f(hi)

# derivada con sympy
dfh = sym.diff(fh,h,1)
df = sym.lambdify(h,dfh)

# SALIDA
print('fh:')
sym.pprint(fh)
print('dfh:')
sym.pprint(dfh)

plt.plot(hi,fi)
plt.axhline(0, color='black')
plt.xlabel('h')
plt.ylabel('f')
plt.title(fh)
plt.show()

## Método de Newton-Raphson

# INGRESO
fx  = f
dfx = df
x0 = (a+b)/2 # mitad del intervalo (ejemplo)
tolera = 0.001

# PROCEDIMIENTO
tabla = []
tramo = abs(2*tolera)
xi = x0
while (tramo>=tolera):
    xnuevo = xi - fx(xi)/dfx(xi)
    tramo  = abs(xnuevo-xi)
    tabla.append([xi,xnuevo,tramo])
    xi = xnuevo

# convierte la lista a un arreglo.
tabla = np.array(tabla)
n = len(tabla)

# SALIDA
print(['xi', 'xnuevo', 'tramo'])
np.set_printoptions(precision = 4)
print(tabla)
print('raiz en: ', xi)
print('con error de: ',tramo)

 

s1Eva_2023PAOI_T3 Recoger los restos de sumergible

ejercicio: 1Eva_2023PAOI_T3 Recoger los restos de sumergible

literal a

Para realizar el planteamiento del ejercicio, se realiza la gráfica para observar mejor los puntos.

recorrido de polinomio

Se podría observar los puntos y plantear un primer recorrido como el siguiente

Rov Sumergible Polinomio01a

los puntos corresponden a los índices j = [0 ,2, 5, 6, 11]:

# INGRESO
xj = np.array([0.2478, 0.6316, 1.3802, 2.4744, 2.7351, 2.8008,
               3.5830, 3.6627, 3.7213, 4.2796, 4.3757, 4.6794])
fj = np.array([1.8108, 0.1993, 4.7199, 2.4529, 2.3644, 3.5955,
               2.6558, 4.3657, 4.1932, 4.6998, 4.7536, 1.5673])

# selecionando puntos
j = [0, 2, 5, 6, 11]
xi = xj[j]
fi = fj[j]

literal b

Partiendo del algoritmo de Diferencias divididas de Newton, y añadiendo vectores xj y fj para todos los puntos dados en el ejercicio.

Tabla Diferencia Dividida
[['i   ', 'xi  ', 'fi  ', 'F[1]', 'F[2]', 'F[3]', 'F[4]', 'F[5]']]
[[ 0.      0.2478  1.8108  2.569  -1.3163  0.3389 -0.0561  0.    ]
 [ 1.      1.3802  4.7199 -0.7915 -0.1861  0.09    0.      0.    ]
 [ 2.      2.8008  3.5955 -1.2014  0.111   0.      0.      0.    ]
 [ 3.      3.583   2.6558 -0.9928  0.      0.      0.      0.    ]
 [ 4.      4.6794  1.5673  0.      0.      0.      0.      0.    ]]
dDividida: 
[ 2.569  -1.3163  0.3389 -0.0561  0.    ]
polinomio: 
2.56896856234546*x - 0.0561488275125463*(x - 3.583)*(x - 2.8008)*(x - 1.3802)*(x - 0.2478) 
+ 0.338875729488637*(x - 2.8008)*(x - 1.3802)*(x - 0.2478) 
- 1.31628089036375*(x - 1.3802)*(x - 0.2478) + 1.17420959025079
polinomio simplificado: 
-0.0561488275125463*x**4 + 0.788728905753656*x**3 
- 3.98331084054791*x**2 + 7.41286537452727*x 
+ 0.206696844067185
>>>

usando la tabla de diferencias divididas se plantea la expresión para el polinomio:

p(x) = 1.8108 + 2.569 (x-0.2478) + + (-1.3163)(x-0.2478)(x-1.3802) + + 0.3389 (x-0.2478)(x-1.3802)(x-2.8008) + + (-0.0561)(x-0.2478)(x-1.3802)(x-2.8008)(x-3.583)

el algoritmo simplifica la expresión al siguiente polinomio,

p(x) = -0.0561 x^4 + 0.7887 x^3 - 3.9833 x^2 + 7.4128 x + 0.2066

literal c

Se usa el algoritmo Interpolación polinómica de Lagrange con Python para desarrollar el polinomio a partir de los puntos no usados en el literal anterior.

j = [3,4,7,8,9,10]

Para evitar oscilaciones grandes, se excluye el punto con índice 1, y se obtiene una opción mostrada:

Rov Sumergible Polinomio 02a

con el siguiente resultado:

    valores de fi:  [2.4529 2.3644 4.3657 4.1932 4.6998 4.7536]
divisores en L(i):  [-1.32578996  0.60430667 -0.02841115  0.02632723 -0.09228243  0.13986517]

Polinomio de Lagrange, expresiones
-1.85014223337671*(x - 4.3757)*(x - 4.2796)*(x - 3.7213)*(x - 3.6627)*(x - 2.7351) 
+ 3.9125829809921*(x - 4.3757)*(x - 4.2796)*(x - 3.7213)*(x - 3.6627)*(x - 2.4744) 
- 153.661523787291*(x - 4.3757)*(x - 4.2796)*(x - 3.7213)*(x - 2.7351)*(x - 2.4744) 
+ 159.272361555669*(x - 4.3757)*(x - 4.2796)*(x - 3.6627)*(x - 2.7351)*(x - 2.4744) 
- 50.928437685792*(x - 4.3757)*(x - 3.7213)*(x - 3.6627)*(x - 2.7351)*(x - 2.4744) 
+ 33.9870187307222*(x - 4.2796)*(x - 3.7213)*(x - 3.6627)*(x - 2.7351)*(x - 2.4744)

Polinomio de Lagrange: 
-9.26814043907703*x**5 + 163.708008152209*x**4 
- 1143.16428460497*x**3 + 3941.13583467614*x**2 
- 6701.74628786718*x + 4496.64364271972
>>> 

Para presentar la parte analítica puede seleccionar menos puntos, se busca mostrar la aplicación del algoritmo, no necesariamente cuan largo puede ser.

p(x) = + 2.4529\frac{(x - 4.3757)(x - 4.2796)(x - 3.7213)(x - 3.6627)(x - 2.7351) }{(2.4744 - 4.3757)(2.4744 - 4.2796)(2.4744 - 3.7213)(2.4744 - 3.6627)(2.4744 - 2.7351)} + 2.3644\frac{ (x - 4.3757)(x - 4.2796)(x - 3.7213)(x - 3.6627)(x - 2.4744) }{(2.7351 - 4.3757)(2.7351 - 4.2796)(2.7351 - 3.7213)(2.7351 - 3.6627)(2.7351 - 2.4744)} +4.3657\frac{(x - 4.3757)(x - 4.2796)(x - 3.7213)(x - 2.7351)(x - 2.4744) }{(3.6627 - 4.3757)(3.6627 - 4.2796)(3.6627 - 3.7213)(3.6627 - 2.7351)(3.6627 - 2.4744)} + 4.1932 \frac{(x - 4.3757)(x - 4.2796)(x - 3.6627)(x - 2.7351)(x - 2.4744) }{(3.7213 - 4.3757)(3.7213 - 4.2796)(3.7213 - 3.6627)(3.7213 - 2.7351)(3.7213 - 2.4744)} +4.6998\frac{(x - 4.3757)(x - 3.7213)(x - 3.6627)(x - 2.7351)(x - 2.4744) }{(4.2796 - 4.3757)(4.2796 - 3.7213)(4.2796 - 3.6627)(4.2796 - 2.7351)(4.2796 - 2.4744)} + 4.7536 \frac{(x - 4.2796)(x - 3.7213)(x - 3.6627)(x - 2.7351)(x - 2.4744)}{(4.3757 - 4.2796)(4.3757 - 3.7213)(4.3757 - 3.6627)(4.3757 - 2.7351)(4.3757 - 2.4744)}

literal d

las gráficas se han presentado en el planteamiento para justificar los criterios usados.

literal e

Los errores de encuentran como la diferencia de los puntos no usados en el polinomio, y evaluando el polinomio en cada coordenada x.

Como las propuestas de polinomio pueden ser variadas se obtendrán diferentes respuestas para cada literal.

Se revisa la aplicación de conceptos en las propuestas.

 

Algoritmos usados

literal b. Diferencias divididas de Newton

# 1Eva_2023PAOI_T3 Recoger los restos de sumergible
# Polinomio interpolación
# Diferencias Divididas de Newton
# Tarea: Verificar tamaño de vectores,
#        verificar puntos equidistantes en x
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO
xj = np.array([0.2478, 0.6316, 1.3802, 2.4744, 2.7351, 2.8008,
               3.5830, 3.6627, 3.7213, 4.2796, 4.3757, 4.6794])
fj = np.array([1.8108, 0.1993, 4.7199, 2.4529, 2.3644, 3.5955,
               2.6558, 4.3657, 4.1932, 4.6998, 4.7536, 1.5673])

# selecionando puntos
j = [0, 2, 5, 6, 11]
xi = xj[j]
fi = fj[j]


# PROCEDIMIENTO

# Tabla de Diferencias Divididas Avanzadas
titulo = ['i   ','xi  ','fi  ']
n = len(xi)
ki = np.arange(0,n,1)
tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
tabla = np.transpose(tabla)

# diferencias divididas vacia
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((tabla,dfinita), axis=1)

# Calcula tabla, inicia en columna 3
[n,m] = np.shape(tabla)
diagonal = n-1
j = 3
while (j < m):
    # Añade título para cada columna
    titulo.append('F['+str(j-2)+']')

    # cada fila de columna
    i = 0
    paso = j-2 # inicia en 1
    while (i < diagonal):
        denominador = (xi[i+paso]-xi[i])
        numerador = tabla[i+1,j-1]-tabla[i,j-1]
        tabla[i,j] = numerador/denominador
        i = i+1
    diagonal = diagonal - 1
    j = j+1

# POLINOMIO con diferencias Divididas
# caso: puntos equidistantes en eje x
dDividida = tabla[0,3:]
n = len(dfinita)

# expresión del polinomio con Sympy
x = sym.Symbol('x')
polinomio = fi[0]
for j in range(1,n,1):
    factor = dDividida[j-1]
    termino = 1
    for k in range(0,j,1):
        termino = termino*(x-xi[k])
    polinomio = polinomio + termino*factor

# simplifica multiplicando entre (x-xi)
polisimple = polinomio.expand()

# polinomio para evaluacion numérica
px = sym.lambdify(x,polisimple)

# Puntos para la gráfica
muestras = 101
a = np.min(xi)
b = np.max(xi)
pxi = np.linspace(a,b,muestras)
pfi = px(pxi)

# SALIDA
np.set_printoptions(precision = 4)
print('Tabla Diferencia Dividida')
print([titulo])
print(tabla)
print('dDividida: ')
print(dDividida)
print('polinomio: ')
print(polinomio)
print('polinomio simplificado: ' )
print(polisimple)

# Gráfica
plt.plot(xj,fj,'o',color='red')
plt.plot(xi,fi,'o', label = 'Puntos')
##for i in range(0,n,1):
##    plt.axvline(xi[i],ls='--', color='yellow')

plt.plot(pxi,pfi, label = 'Polinomio')
plt.legend()
plt.grid()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Diferencias Divididas - Newton')
plt.show()

literal c. Polinomio de Lagrange

# 1Eva_2023PAOI_T3 Recoger los restos de sumergible
# Interpolacion de Lagrange
# divisoresL solo para mostrar valores
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO , Datos de prueba
xj = np.array([0.2478, 0.6316, 1.3802, 2.4744, 2.7351, 2.8008,
               3.5830, 3.6627, 3.7213, 4.2796, 4.3757, 4.6794])
fj = np.array([1.8108, 0.1993, 4.7199, 2.4529, 2.3644, 3.5955,
               2.6558, 4.3657, 4.1932, 4.6998, 4.7536, 1.5673])

# selecionando puntos
#j = [0, 2, 5, 6, 11]
j = [3,4,7,8,9,10]
xi = xj[j]
fi = fj[j]

# PROCEDIMIENTO
# Polinomio de Lagrange
n = len(xi)
x = sym.Symbol('x')
polinomio = 0
divisorL = np.zeros(n, dtype = float)
for i in range(0,n,1):
    
    # Termino de Lagrange
    numerador = 1
    denominador = 1
    for j  in range(0,n,1):
        if (j!=i):
            numerador = numerador*(x-xi[j])
            denominador = denominador*(xi[i]-xi[j])
    terminoLi = numerador/denominador

    polinomio = polinomio + terminoLi*fi[i]
    divisorL[i] = denominador

# simplifica el polinomio
polisimple = polinomio.expand()

# para evaluación numérica
px = sym.lambdify(x,polisimple)

# Puntos para la gráfica
muestras = 101
a = np.min(xi)
b = np.max(xi)
pxi = np.linspace(a,b,muestras)
pfi = px(pxi)

# SALIDA
print('    valores de fi: ',fi)
print('divisores en L(i): ',divisorL)
print()
print('Polinomio de Lagrange, expresiones')
print(polinomio)
print()
print('Polinomio de Lagrange: ')
print(polisimple)

# Gráfica
plt.plot(xj,fj,'o',color='red')
plt.plot(xi,fi,'o', label = 'Puntos')
plt.plot(pxi,pfi, label = 'Polinomio')
plt.legend()
plt.grid()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Interpolación Lagrange')
plt.show()

 

s1Eva_2023PAOI_T2 Productos en combo por subida de precio

Ejercicio: 1Eva_2023PAOI_T2 Productos en combo por subida de precio

literal a

Para el ejercicio solo es posible plantear tres ecuaciones, se muestran datos solo para tres semanas. Se tienen 4 incógnitas, por lo que tendría infinitas soluciones y no se podría resolver.

500 a + 600 b + 400 c + 90 d = 1660 800 a + 450 b + 300 c + 100 d = 1825 400 a + 300 b + 600 c + 80 d = 1430

Sin embargo desde el punto de vista práctico se puede usar una variable libre, considerando algunos criterios como:

– que en la nota se da el precio para la docena de huevos a 2.5 USD
– que la variable para huevos es la de menor cantidad se puede incurrir en un menor error si el precio ha variado en los últimos días respecto a lo que se podría haber tenido como referencia.

500 a + 600 b + 400 c = 1660 - 90 d 800 a + 450 b + 300 c = 1825 - 100 d 400 a + 300 b + 600 c = 1430 - 80 d

que al sustituir con d = 2.5

500 a + 600 b + 400 c = 1660 - 90 (2.5) =1435 800 a + 450 b + 300 c = 1825 - 100 (2.5) = 1575 400 a + 300 b + 600 c = 1430 - 80 (2.5) = 1230

Nota: el estudiante puede realizar una justificación diferente para el uso de otra variable libre.

literal b

La forma matricial del sistema de ecuaciones se convierte a:

\begin{pmatrix} 500 & 600 & 400 & \Big| & 1435 \\ 800 & 450 & 300 & \Big| & 1575 \\ 400 & 300 & 600 &\Big| & 1230 \end{pmatrix}

literal c

pivoteando la matriz aumentada, se encuentra que para la primera columna se pivotean la fila 1 y 2, por ser 800 el valor de mayor magnitud:

\begin{pmatrix} 800 & 450 & 300 & \Big| & 1575 \\ 500 & 600 & 400 & \Big| & 1435 \\ 400 & 300 & 600 &\Big| & 1230 \end{pmatrix}

luego, observando la segunda columna desde el elemento de la diagonal hacia abajo, se observa que el mayor valor en magnitud ya se encuentra en la diagonal. No se requieren intercambios de fila para la tercera columna.

literal d

Para el método iterativo de Gauss-Seidel se requiere el vector inicial, tomando las observaciones de la nota. Se da un precio de 50 USD por un quintal métrico de 100Kg, por lo que se estima que el precio por kilogramo ronda 50/100= 0.5 USD. Se indica además que los demás valores se encuentran en el mismo orden de magnitud, por lo que se tiene como vector inicial

X0 = [0.5, 0.5, 0.5]

Para un método iterativo se despejan las ecuaciones como:

a = \frac {1575 - 450 b - 300 c}{800} b = \frac {1435 -500 a - 400 c}{600} c = \frac {1230 - 400 a - 300 b}{600}

Para el cálculo del error se considera que se busca un precio, lo regular es considerar el centavo, es decir una tolerancia 10-2. Para el desarrollo se considera usar cuatro decimales por si se considera truncar o redondear.

itera = 0

a = \frac {1575 - 450 (0.5) - 300 (0.5)}{800} = 1.5 b = \frac {1435 -500(1.5) - 400 (0.5)}{600} = 0.8083 c = \frac {1230 - 400(1.5) - 300 (0.8083)}{600} =0.6458

errado = max|[1.5-0.5, 0.8083-0.5, 0.6458-0.5]|
errado = max|[1, 0.3083, 0.1458]| = 1

itera = 1

X1 = [1.5, 0.8083, 0.6458]

a = \frac {1575 - 450 (0.8083) - 300 (0.6458)}{800} = 1.2719 b = \frac {1435 -500(1.2719) - 400 (0.6458)}{600} = 0.9012 c = \frac {1230 - 400 (1.2719) - 300 (0.9012)}{600} = 0.7514

errado = max|[1.2719-1.5, 0.9012-0.8083, 0.7514-0.6458]|
errado = max|[-0.2281, 0.0929, 0.1056]| = 0.2281

el error disminuye entre iteraciones.

itera = 2

X2 = [1.2719 , 0.9012, 0.7514]

a = \frac {1575 - 450 (0.9012) - 300 (0.7514)}{800} = 1.1805 b = \frac {1435 -500 (1.1805) - 400 (0.7514)}{600} = 0.9069 c = \frac {1230 - 400 (1.1805) - 300 (0.9069)}{600} = 0.8095

errado = max|[ 1.1805-1.2719 , 0.9069-0.9012, 0.8095-0.7514]|
errado = max|[-0.0914, 0.0057, 0.1056]| = 0.1056

el error disminuye entre iteraciones.

literal e

Si el error entre iteraciones disminuye, se considera que el método converge.

usando el algoritmo Método de Gauss-Seidel con Python, se tiene como resultado en 13 iteraciones:

[1.5        0.80833333 0.64583333] 1.0
[1.271875   0.90121528 0.75147569] 0.2281249999999999
[1.18001302 0.90733869 0.80965531] 0.09186197916666661
[1.15475125 0.88960375 0.83536396] 0.025708648304880177
[1.1550864  0.87218536 0.84384972] 0.01741839593330019
[1.16170209 0.86101511 0.84502438] 0.011170246538965367
[1.16754486 0.85536303 0.84395525] 0.005842764350612484
[1.17112508 0.85309227 0.84270381] 0.0035802211655899807
[1.17287167 0.85247107 0.84185002] 0.0017465903731879173
[1.17354127 0.85248226 0.84139802] 0.0006695985845832642
[1.17370447 0.85264759 0.84120656] 0.00019146597344232852
[1.17368327 0.8527929  0.84114804] 0.00014530950399282982
[1.17362348 0.85288174 0.84114348] 8.884049018109685e-05
respuesta X: 
[[1.17362348]
 [0.85288174]
 [0.84114348]]
verificar A.X=B: 
[[1575.03861029]
 [1434.99817609]
 [1230.        ]]
iteraciones 13
>>> 

instrucciones en Python

Se debe recordar que las matrices y vectores deben ser tipo real (float) .

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

import numpy as np

# INGRESO
A = np.array([[800, 450, 300],
              [500, 600, 400],
              [400, 300, 600]], dtype=float)

B = np.array([1575, 1435,1230], dtype=float)

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

tolera = 0.0001
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)
    print(X,errado)
    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)
print('iteraciones',itera)

s1Eva_2023PAOI_T1 Desacople de cohete de dos etapas

Ejercicio: 1Eva_2023PAOI_T1 Desacople de cohete de dos etapas

literal a. Planteamiento

La ecuación a usar según el enunciado es y usando los valores dados es:

v = u \ln\Big(\frac{m_0}{m_0-qt}\Big) - gt 800 = 1870 \ln\Big(\frac{195000}{195000-2500 t}\Big) - 9.8 t

con lo que la función para buscar la raíz es:

f(t) = 1870 \ln\Big(\frac{195000}{195000-2500 t}\Big) - 9.8 t -800

Literal b. Intervalo de búsqueda

Para el intervalo de búsqueda se puede usar una gráfica e interpretar el punto a buscar alrededor de 800 m/s. Que de la gráfica se observa que un intervalo alrededor de 35 sería válido para el método de la Bisección. Para otros métodos abiertos, también es posible deducir un punto t0.

La validación se muestra con la primera iteración al evaluar f(30) que es negativo y f(40) que es de signo positivo.

Desacople de cohete de dos etapasInstrucciones en Python

# 1Eva_2023PAOI_T1 Desacople de cohete de dos etapas
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
v = lambda t: 1870*np.log(195000/(195000-2500*t))-9.8*t
a = 0
b = 40
tramos = 51

# PROCEDIMIENTO
ti = np.linspace(a,b,tramos)
vi = v(ti)

# SALIDA
plt.plot(ti,vi)
plt.xlabel('ti')
plt.ylabel('vi')
plt.title('Velocidad vertical vs tiempo')
plt.grid()
plt.show()

literal c. Desarrollo con algoritmo de Bisección

itera =0, intervalo [30,40]

f(30) = 1870 \ln\Big(\frac{195000}{195000-2500 (30)}\Big) - 9.8 (30) -800 = -186.100 f(40) = 1870 \ln\Big(\frac{195000}{195000-2500 (40)}\Big) - 9.8 (40) -800 = 152.759 c= \frac{a+b}{2}= \frac{30+40}{2} = 35 f(35) = 1870 \ln\Big(\frac{195000}{195000-2500 (35)}\Big) - 9.8 (35) -800 = -29.399

error = tramo = |40-30| = 10

como f(a) y f(c) son del mismo signo, el intervalo nuevo será [35,40]

itera =1 , intervalo [35,40]

c= \frac{35+40}{2} = 37.5 f(37.5) = 1870 \ln\Big(\frac{195000}{195000-2500 (37.5)}\Big) - 9.8 (37.5) -800 = 58.111

error = tramo = |40-35| = 5

como f(c) y f(b) son del mismo signo, el intervalo nuevo será [35,37.5]

itera =2 , intervalo [35,37.5]

c= \frac{37.5+35}{2} = 36.25 f(36.25) = 1870 \ln\Big(\frac{195000}{195000-2500 (36.25)}\Big) - 9.8 (36.25) -800 = 13.518

error = tramo = |37.5-35| = 2.5

como f(c) y f(b) son del mismo signo, el intervalo nuevo será [35,36.25]

[ i, a,    c,     b,     f(a),    f(c),   f(b),   tramo]
  1 30.000 35.000 40.000 -186.100 -29.399 152.759 10.000 
  2 35.000 37.500 40.000 -29.399   58.111 152.759  5.000 
  3 35.000 36.250 37.500 -29.399   13.518 58.111   2.500 
  4 35.000 35.625 36.250 -29.399   -8.144 13.518   1.250 
  5 35.625 35.938 36.250 -8.144     2.635 13.518   0.625 
  6 35.625 35.781 35.938 -8.144    -2.767 2.635    0.312 
  7 35.781 35.859 35.938 -2.767    -0.069 2.635    0.156 
  8 35.859 35.898 35.938 -0.069     1.282 2.635    0.078 
raiz:  35.8984375
>>>

literal d. tolerancia y errores

La tolerancia depende de la escala a la que se mide y el instrumento de medición. Si consideramos décimas de segundo la tolerancia será de 10-1. ó 0.1

Los errores entre iteraciones se muestran en el literal anterior.

literal e. convergencia

Los errores en cada iteración disminuye, lo que muestra que el método converge. Luego de 8 iteraciones se encuentra el tiempo ti a usar como 35.8 s.

Se adjunta el algoritmo de la bisección ajustado para el ejercicio. La gráfica es complementaria a la presentada en el literal b, que puede ser incorporada al mismo algoritmo para la presentación.

Instrucciones en Python

# Algoritmo de Bisección
# 1Eva_2023PAOI_T1 Desacople de cohete de dos etapas
import numpy as np

# INGRESO
fx = lambda t: 1870*np.log(195000/(195000-2500*t))-9.8*t -800
a = 30
b = 40
tolera = 0.1

# PROCEDIMIENTO
tabla = []
tramo = b-a

fa = fx(a)
fb = fx(b)
i = 1
while (tramo>tolera):
    c = (a+b)/2
    fc = fx(c)
    tabla.append([i,a,c,b,fa,fc,fb,tramo])
    i = i + 1
                 
    cambia = np.sign(fa)*np.sign(fc)
    if (cambia<0):
        b = c
        fb = fc
    else:
        a=c
        fa = fc
    tramo = b-a
c = (a+b)/2
fc = fx(c)
tabla.append([i,a,c,b,fa,fc,fb,tramo])
tabla = np.array(tabla)

raiz = c

# SALIDA
np.set_printoptions(precision = 4)
print('[ i, a, c, b, f(a), f(c), f(b), tramo]')
# print(tabla)

# Tabla con formato
n=len(tabla)
for i in range(0,n,1):
    unafila = tabla[i]
    formato = '{:.0f}'+' '+(len(unafila)-1)*'{:.3f} '
    unafila = formato.format(*unafila)
    print(unafila)
    
print('raiz: ',raiz)

s1Eva_2022PAOII_T3 Trayectoria de dron con polinomios

Ejercicio: 1Eva_2022PAOII_T3 Trayectoria de dron con polinomios

La variable independiente para la trayectoria es tiempo, con datos en el vector de ti.

ti  = [0, 1, 2, 3, 4]
xti = [2, 1, 3, 4, 2]
yti = [0, 1, 5, 1, 0]

Considerando que los puntos marcan posiciones por donde debe pasar el dron y se define la trayectoria, se usarán todos los puntos. Cada polinomio será de grado 4 al incluir los 5 puntos disponibles para cada eje.

Nota: podría usar polinomios de menor grado, siempre que considere que se debe completar la trayectoria y regresar al punto de salida.

px(t) = 2\frac{(t-1)(t-2)(t-3)(t-4)}{(0-1)(0-2)(0-3)(0-4)} + 1 \frac{(t-0)(t-2)(t-3)(t-4)}{(1-0)(1-2)(1-3)(1-4)} + 3 \frac{(t-0)(t-1)(t-3)(t-4)}{(2-0)(2-1)(2-3)(2-4)} + 4 \frac{(t-0)(t-1)(t-2)(t-4)}{(3-0)(3-1)(3-2)(3-4)} + 2 \frac{(t-0)(t-1)(t-2)(t-3)}{(4-0)(4-1)(4-2)(4-3)}

simplificando con el algoritmo:

px(t) = \frac{1}{12}t^4 - \frac{7}{6}t^3 + \frac{53}{12}t^2 - \frac{13}{3}t + 2

Realizando lo mismo con el algoritmo para polinomio de Lagrange se obtiene:

py(t) = \frac{11}{12}t^4 - \frac{22}{3}t^3 + \frac{205}{12}t^2 - \frac{29}{3}t

se muestra la gráfica de trayectorias por cada eje vs tiempo

Observaciones: La trayectoria usada tiene el mismo punto de salida como de retorno. La trayectoria presenta lóbulos que podrían ser reducidos y minimizar uso de recursos como bateria. Considere usar trazadores cúbicos y observe la misma gráfica de trayectorias x(t) vs y(t).

Resultado con el algoritmo:

Polinomio de Lagrange x: 
x**4/12 - 7*x**3/6 + 53*x**2/12 - 13*x/3 + 2
Polinomio de Lagrange y: 
11*x**4/12 - 22*x**3/3 + 205*x**2/12 - 29*x/3

Algoritmo en Python

# Interpolacion de Lagrange
# divisoresL solo para mostrar valores
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO , Datos de prueba
ti  = [0,1,2,3,4]
xti = [2,1,3,4,2]
yti = [0,1,5,1,0]

# PROCEDIMIENTO
x = sym.Symbol('x')

def interpola_lagrange(xi,yi):
    '''
    Interpolación con método de Lagrange
    resultado: polinomio en forma simbólica
    '''
    # PROCEDIMIENTO
    n = len(xi)
    x = sym.Symbol('x')
    # Polinomio
    polinomio = 0
    for i in range(0,n,1):
        # Termino de Lagrange
        termino = 1
        for j  in range(0,n,1):
            if (j!=i):
                termino = termino*(x-xi[j])/(xi[i]-xi[j])
        polinomio = polinomio + termino*yi[i]
    # Expande el polinomio
    polinomio = polinomio.expand()
    return(polinomio)

# para ejex
polinomiox = interpola_lagrange(ti,xti)
polisimplex = polinomiox.expand()
px = sym.lambdify(x,polisimplex)

# para ejey
polinomioy = interpola_lagrange(ti,yti)
polisimpley = polinomioy.expand()
py = sym.lambdify(x,polisimpley)

# Puntos para la gráfica
muestras = 101
a = np.min(ti)
b = np.max(ti)
ti = np.linspace(a,b,muestras)
pxi = px(ti)
pyi = py(ti)

# SALIDA
print('Polinomio de Lagrange x: ')
print(polisimplex)
print('Polinomio de Lagrange y: ')
print(polisimpley)

# Gráfica
figura, enplano = plt.subplots()
plt.scatter(xti,yti, color='red')
plt.plot(pxi,pyi)
plt.ylabel('y(t)')
plt.xlabel('x(t)')
plt.title('trayectoria 2D')
plt.grid()

figura, entiempo = plt.subplots()
plt.plot(ti,pxi, label = 'px')
plt.plot(ti,pyi, label = 'py')
plt.legend()
plt.title('posicion en tiempo')
plt.xlabel('t')
plt.ylabel('p(t)')
plt.grid()

plt.show()

s1Eva_2022PAOII_T2 Admisión universitaria – cupos por recursos

Ejercicios: 1Eva_2022PAOII_T2 Admisión universitaria – cupos por recursos

Se requiere determinar la distribución de cupos en base a los costos relativos al promedio por estudiante para docencia, infraestructura y servicios mostrados en la tabla.

Costo referencial /carrera Mecatrónica Computación Civil Matemáticas
Docencia 1.5 0.9 0.6 0.7
Infraestructura 0.8 1.4 0.4 0.5
Servicios 0.45 0.55 1.1 0.5

Con los datos del total de recursos relativos al promedio por estudiante disponibles son docencia 271, infraestructura 250 y servicios 230.

1.5 a + 0.9 b + 0.6 c + 0.7 d = 271 0.8 a + 1.4 b + 0.4 c + 0.5 d = 250 0.45 a + 0.55 b + 1.1 c + 0.5 d = 230

se indica que en carreras como matemáticas de baja demanda, se establece el cupo de 10,

1.5 a + 0.9 b + 0.6 c + 0.7 (10) = 271 0.8 a + 1.4 b + 0.4 c + 0.5(10) = 250 0.45 a + 0.55 b + 1.1 c + 0.5(10) = 230

el sistema se convierte en:

1.5 a + 0.9 b + 0.6 c = 271 - 0.7 (10) 0.8 a + 1.4 b + 0.4 c = 250 - 0.5(10) 0.45 a + 0.55 b + 1.1 c = 230 - 0.5(10)

Para usar un método iterativo se convierte a matriz aumentada:

\begin{pmatrix} 1.5 & 0.9 & 0.6 & \Big| & 264 \\ 0.8 & 1.4 & 0.4 & \Big| & 245 \\ 0.45 & 0.55 & 1.1 &\Big| & 225 \end{pmatrix}

con pivoteo parcial por filas, la matriz aumentada se mantiene igual, pues los valores de la diagonal ya son los mayores posibles según el algoritmo.

Para un método iterativo se despeja una ecuación por cada incógnita.

a = \frac{1}{1.5}(264 - 0.9 b - 0.6 c) b = \frac{1}{1.4}(245 - 0.8 a - 0.4 c) c = \frac{}{1.1}(225 -0.45 a - 0.55 b)

Para los valores iniciales se consideran números mayores que cero, pues existen recursos para los cupos. No se admiten cupos negativos.

X_0 = [50,50,50]

Las iteraciones para el método iterativo de Gauss-Seidel

itera = 0

a = \frac{1}{1.5}(264 - 0.9 (50) - 0.6 (50)) = 126 b = \frac{1}{1.4}(245 - 0.8 (126) - 0.4 (50)) = 88.714 c = \frac{}{1.1}(225 -0.45 (126) - 0.55 (88.714)) =108.642 diferencia = [126-50, 88.714-50, 108.42-50] diferencia = [76, 38.714, 58.642] errado = max|[76, 38.714, 58.642]| =76 X = [126, 88.714, 108.42]

itera = 1

a = \frac{1}{1.5}(264 - 0.9 (88.714) - 0.6 (108.42)) = 79.314 b = \frac{1}{1.4}(245 - 0.8 (79.314) - 0.4 (108.42)) = 98.637 c = \frac{}{1.1}(225 -0.45 (79.314) - 0.55 (98.637)) =122.780 diferencia = [79.314-126, 88, 98.637-88.714, 122.780-108.42] diferencia = [46.685,9.922, 14.137] errado = max| [46.685,9.922, 14.137] | = 46.685

el error disminuye en la iteración

X = [79.314 , 98.637, 122.780]

itera = 2

a = \frac{1}{1.5}(264 - 0.9 (79.314) - 0.6 (122.780)) = 67.705 b = \frac{1}{1.4}(245 - 0.8 (67.705) - 0.4 (122.780)) = 101.230 c = \frac{}{1.1}(225 -0.45 (67.705) - 0.55 (101.230)) =126.232 diferencia = [67.705-79.314, 101.230-98.637, 126.232-122.780] diferencia = [-11.608, 2.594, 3.451] errado = max| [-11.608, 2.594, 3.451] | = 11.608

el error disminuye en la iteración, se considera que el método converge

X = [67.705 , 101.230, 126.232]

con el algoritmo se tiene como resultado:

[126.          88.71428571 108.64285714]
[76.         38.71428571 58.64285714]

[ 79.31428571  98.63673469 122.78033395]
[46.68571429  9.92244898 14.13747681]

[ 67.7058256  101.23086138 126.23218611]
[11.60846011  2.59412669  3.45185216]

[ 64.76860873 101.92302755 127.08769174]
[2.93721688 0.69216617 0.85550564]

[ 64.01110677 102.11145563 127.30336487]
[0.75750196 0.18842808 0.21567312]

[ 63.81178067 102.16373537 127.3587675 ]
[0.1993261  0.05227973 0.05540263]

[ 63.75825178 102.17849398 127.37328637]
[0.05352889 0.01475862 0.01451887]

[ 63.74358906 102.18272443 127.37716953]
[0.01466272 0.00423045 0.00388316]

[ 63.73949753 102.18395297 127.37822907]
[0.00409153 0.00122854 0.00105954]

[ 63.73833659 102.18431364 127.37852366]
[0.00116094 0.00036067 0.0002946 ]

[ 63.73800235 102.18442047 127.37860699]
[3.34240232e-04 1.06824300e-04 8.33224905e-05]

respuesta X: 
[[ 63.73800235]
 [102.18442047]
 [127.37860699]]
verificar A.X=B: 
[[264.00014614]
 [245.00003333]
 [225.        ]]
>>>

se interpreta la respuesta como la parte entera de la solución:

cupos = [ 63, 102 , 127]

s1Eva_2022PAOII_T1 Esfera flotando en agua

Ejercicio: 1Eva_2022PAOII_T1 Esfera flotando en agua

Según el principio de Arquímedes, la fuerza de flotación o empuje es igual al peso de el fluido desplazado por la porción sumergida de un objeto.

F_{empuje} = F_{peso} \rho_{agua} V_{sumergido} \text{ } g = \rho_{esfera}V_{esfera} \text{ } g V_{sumergido} = \frac{\rho_{esfera}}{\rho_{agua}}V_{esfera} V_{esfera} - V_{sobreagua} = \frac{\rho_{esfera}}{\rho_{agua}}V_{esfera} V_{sobreagua} = \Big( 1- \frac{\rho_{esfera}}{\rho_{agua}}\Big) V_{esfera} V_{esfera} = \frac{4}{3}\pi r^3 \frac{\pi h^2}{3}(3r-h) = \Big( 1- \frac{\rho_{esfera}}{\rho_{agua}}\Big) \frac{4}{3}\pi r^3 h^2(3r-h) = \Big( 1- \frac{\rho_{esfera}}{\rho_{agua}}\Big) 4 r^3

El planteamiento para la búsqueda de raíces es f(x) = 0, que para este caso será:

f(h) = h^2(3r-h) - \Big( 1 - \frac{\rho_{esfera}}{\rho_{agua}}\Big) 4 r^3 = 0

usando los valores dados para el ejercicio, r=1 y ρesfera = 200 Kg/m3 y ρagua    = 1000 kg/m3 se tiene que:

f(h) = h^2(3-h) - \Big( 1 - \frac{200}{1000}\Big) 4 f(h) = h^2(3-h) - \frac{16}{5}

Se observa la gráfica de f(h) en el intervalo de h entre[0,2] interpretado como totalmente sumergida y totalmente flotando sobre el agua, confirmando que existe una raíz

Para el caso de aplicar el método del punto fijo se plantea que x=g(x),

h = g(h) h^2(3-h) = \frac{16}{5}

con lo que se puede plantear dos ecuaciones al despejar h

h = \sqrt{ \frac{16}{5(3-h)}} h = 3-\frac{16}{5 h^2}

Iteraciones de la primera ecuación

itera = 0 ; h = h0 = 0.5 ;

g(h) = \sqrt{ \frac{16}{5(3-0.5)}} = 1.1313 tramo = |1.1313-0.5|=0.6313

itera = 1 ; h = 1.1313 ;

g(h) = \sqrt{ \frac{16}{5(3-1.1313)}} = 1.3086 tramo = |1.3086-1.1313| = 0.1772

itera = 2 ; h = 1.3086 ;

g(h) = \sqrt{ \frac{16}{5(3-1.3086)}} = 1.3754 tramo = |1.3754-1.3086| = 0.0668

Observando los errores o tramos en cada iteración se tiene que se reduce, el método converge.


resultados.txt

x,g(x),tramo
0.5 1.131370849898476 0.631370849898476
1.131370849898476 1.308619626317284 0.17724877641880799
1.308619626317284 1.3754802083033437 0.06686058198605971
1.3754802083033437 1.4035002223557855 0.02802001405244181
1.4035002223557855 1.4157629993958152 0.012262777040029649
1.4157629993958152 1.4212317895316 0.005468790135784829
1.4212317895316 1.4236912066694054 0.0024594171378053975
1.4236912066694054 1.424801422465215 0.0011102157958096104
1.424801422465215 1.4253034412081806 0.0005020187429656264
raiz: 1.4253034412081806

Algoritmo en Python

# Algoritmo de punto fijo
# [a,b] intervalo de búsqueda
# error = tolera

import numpy as np
import matplotlib.pyplot as plt

def puntofijo(gx,a,tolera, iteramax = 15):
    i = 1 # iteración
    b = gx(a)
    tramo = abs(b-a)
    print('x,g(x),tramo')
    print(a,b,tramo)
    while(tramo>=tolera and i<=iteramax ):
        a = b
        b = gx(a)
        tramo = abs(b-a)
        print(a,b,tramo)
        i = i + 1
    respuesta = b
    
    # Validar respuesta
    if (i>=iteramax ):
        respuesta = np.nan
    return(respuesta)

# PROGRAMA ---------
# INGRESO
fx = lambda h: h**2*(3-h)-16/5
gx = lambda h: np.sqrt(16/(5*(3-h)))

#fx = lambda h: h**2*(3-h)-16/5
#gx = lambda h: 3-16/(5*(h**2))

x0 = 0.5
tolera = 0.001
iteramax = 50  # itera máximo
a = 0     # intervalo
b = 2
muestras = 51  # gráfico

# PROCEDIMIENTO
respuesta = puntofijo(gx,x0,tolera)

# SALIDA
print('raiz:',respuesta)

hi = np.linspace(a,b,muestras)
fi = fx(hi)
gi = gx(hi)
plt.plot(hi,fi,label='f(h)')
plt.plot(hi,gi,label='g(h)')
plt.plot(hi,hi,label='Identidad')
plt.axhline(0,color='grey')
plt.grid()
plt.xlabel('h')
plt.ylabel('f(h)')
plt.title('esfera sumergida')
plt.legend()
plt.show()

s1Eva_2022PAOI_T3 Interpolar crecimiento de contagios

Ejercicio: 1Eva_2022PAOI_T3 Interpolar crecimiento de contagios

Día del mes 1 8 15 22
Contagios 1 5.6 27 43.5

a) Realice el planteamiento del sistema de ecuaciones que se usaría usando el método de interpolación polinómica.

El modelo de polinomio de grado máximo que se puede obtener es grado 3:

p_3(t) = a_0 + a_1 t + a_2 t^2 + a_3 t^3

por lo que usando los valores de los puntos dados en la tabla:

p_3(1) = a_0 + a_1 (1) + a_2 (1)^2 + a_3 (1)^3 = 1 p_3(8) = a_0 + a_1 (8) + a_2 (8)^2 + a_3 (8)^3 = 5.6 p_3(15) = a_0 + a_1 (15) + a_2 (15)^2 + a_3 (15)^3 = 27 p_3(22) = a_0 + a_1 (22) + a_2 (22)^2 + a_3 (22)^3 = 43.5

b) Realice el planteamiento del sistema de ecuaciones en su forma matricial y muestre la matriz aumentada.

\begin{pmatrix} 1 & 1 & 1^2 & 1^2\\ 1 & 8 & 8^2 & 8^3 \\ 1 & 15 & 15^2 & 15^3 \\ 1 & 22 & 22^2 & 22^3 \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \end{pmatrix} \begin{pmatrix} 1 \\ 56 \\ 27 \\43.5 \end{pmatrix}

matriz aumentada,

\begin{pmatrix} 1 & 1 & 1^2 & 1^2 & 1 \\ 1 & 8 & 8^2 & 8^3 & 56 \\ 1 & 15 & 15^2 & 15^3 & 27\\ 1 & 22 & 22^2 & 22^3 & 43.5\end{pmatrix}

c) Desarrolle el pivoteo parcial por filas, indicando las operaciones realizadas en éste proceso

pivoteo parcial por filas

\begin{pmatrix} 1 & 1 & 1^2 & 1^2 & 1 \\ 1 & 22 & 22^2 & 22^3 & 43.5 \\ 1 & 15 & 15^2 & 15^3 & 27 \\ 1 & 8 & 8^2 & 8^3 & 56\end{pmatrix}

d) Usando el método directo de Gauss-Jordan, muestre las expresiones necesarias para el algoritmo.

eliminación hacia adelante

\begin{pmatrix} 1 & 1 & 1^2 & 1^2 & 1 \\ 1-1 & 22-1 & 22^2-1^2 & 22^3 -1^3& 43.5 - 1\\ 1-1 & 15-1 & 15^2 -1^2& 15^3 -1^3& 27 -1\\ 1-1 & 8-1 & 8^2 -1^2& 8^3 -1^3& 56-1\end{pmatrix} \begin{pmatrix} 1 & 1 & 1 & 1 & 1 \\ 0 & 21 & 483 & 10647& 42.5 \\ 0 & 14 & 224 & 3376& 26\\ 0 & 7 & 63& 511 & 55\end{pmatrix} \begin{pmatrix} 1 & 1 & 1 & 1 & 1 \\ 0 & 21 & 483 & 10647& 42.5 \\ 0 & 14-\frac{14}{21} 21 & 224-\frac{14}{21}483& 3376 -\frac{14}{21}10647& 26-\frac{14}{21}42.5\\ 0 & 7-\frac{7}{21}21 & 63-\frac{7}{21}483& 511-\frac{7}{21}10647 & 55-\frac{7}{21}42.5\end{pmatrix} \begin{pmatrix} 1 & 1 & 1 & 1 & 1 \\ 0 & 21 & 483 & 10647& 42.5 \\ 0 & 0 &-98 & -3722& -2.33 \\ 0 & 0 & -98 & -3038 & 40.83 \end{pmatrix} \begin{pmatrix} 1 & 1 & 1 & 1 & 1 \\ 0 & 21 & 483 & 10647& 42.5 \\ 0 & 0 &-98 & -3722& -2.33 \\ 0 & 0 & -98-\frac{98}{-98}98 & -3038 -\frac{98}{-98}3722& 40.83 -\frac{98}{-98}(-2.33)\end{pmatrix} \begin{pmatrix} 1 & 1 & 1 & 1 & 1 \\ 0 & 21 & 483 & 10647& 42.5 \\ 0 & 0 &-98 & -3722& -2.33 \\ 0 & 0 & 0 & 686 & -7.23\end{pmatrix}

realizando el proceso de eliminación hacia atrás, semejante al método anterior se obtiene

\begin{pmatrix} 1 & 0 & 0 & 0 & 2.98\\ 0 & 1 & 0 & 0& -2.39 \\ 0 & 0 & 1 & 0 & 0.424 \\ 0 & 0 & 0 &1 & -0.0105\end{pmatrix}

con lo que el vector resultado es:

X= [2.98, -2.39, 0.42, -0.0105]

El polinomio de interpolación resultante es:

p(t)= 2.98 -2.39 t + 0.42 t^2 -0.0105 t^3

e) Para el día 19 se encuentra que el valor correspondiente a contagios es de 37%. Estime el error presentado del modelo para ese día.

p(19)= 2.98 -2.39 (19) + 0.42 (19)^2 -0.0105 (19)^3 = 38.42 error = |38.42-37| = 1.42

f) Desarrolle el ejercicio usando otro método para encontrar el polinomio de interpolación.

usando diferencias finitas

Tabla Diferencia Finita
[['i', 'xi', 'fi', 'df1', 'df2', 'df3', 'df4']]
[[  0.    1.    1.    4.6  16.8 -21.7   0. ]
 [  1.    8.    5.6  21.4  -4.9   0.    0. ]
 [  2.   15.   27.   16.5   0.    0.    0. ]
 [  3.   22.   43.5   0.    0.    0.    0. ]]

polinomio:

p(t) = 1+\frac{4.6}{1! (7)}(t-1) + + \frac{16.8}{2!(7^2)}(t-1)(t-8) + +\frac{-21.7}{3!(7^3}(t-1)(t-8)(t-15)

Algoritmo en Python

Para literal f

# Polinomio interpolación
# Diferencias finitas avanzadas
# Tarea: Verificar tamaño de vectores,
#        verificar puntos equidistantes en x

import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

# INGRESO , Datos de prueba
xi = np.array([1,8,15,22],dtype=float)
fi = np.array([1,5.6,27,43.5],dtype=float)

# PROCEDIMIENTO

# Tabla de Diferencias Finitas
titulo = ['i','xi','fi']
n = len(xi)
ki = np.arange(0,n,1)
tabla = np.concatenate(([ki],[xi],[fi]),axis=0)
tabla = np.transpose(tabla)

# diferencias finitas vacia
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((tabla,dfinita), axis=1)

# Calcula tabla, inicia en columna 3
[n,m] = np.shape(tabla)
diagonal = n-1
j = 3
while (j < m):
    # Añade título para cada columna
    titulo.append('df'+str(j-2))
    # cada fila de columna
    i = 0
    while (i < diagonal):
        tabla[i,j] = tabla[i+1,j-1]-tabla[i,j-1]
        i = i+1
    diagonal = diagonal - 1
    j = j+1

# POLINOMIO con diferencias Finitas avanzadas
# caso: puntos equidistantes en eje x
h = xi[1] - xi[0]
dfinita = tabla[0,3:]
n = len(dfinita)

# expresión del polinomio con Sympy
x = sym.Symbol('x')
polinomio = fi[0]
for j in range(1,n,1):
    denominador = np.math.factorial(j)*(h**j)
    factor = dfinita[j-1]/denominador
    termino = 1
    for k in range(0,j,1):
        termino = termino*(x-xi[k])
    polinomio = polinomio + termino*factor

# simplifica multiplicando entre (x-xi)
polisimple = polinomio.expand()

# polinomio para evaluacion numérica
px = sym.lambdify(x,polisimple)

# Puntos para la gráfica
muestras = 101
a = np.min(xi)
b = np.max(xi)
pxi = np.linspace(a,b,muestras)
pfi = px(pxi)

# SALIDA
print('Tabla Diferencia Finita')
print([titulo])
print(tabla)
print('dfinita: ')
print(dfinita)
print('polinomio: ')
print(polinomio)
print('polinomio simplificado: ' )
print(polisimple)

# Gráfica
plt.plot(xi,fi,'o', label = 'Puntos')
##for i in range(0,n,1):
##    plt.axvline(xi[i],ls='--', color='yellow')
plt.plot(pxi,pfi, label = 'Polinomio')

plt.legend()
plt.xlabel('xi')
plt.ylabel('fi')
plt.title('Interpolación polinómica')
plt.show()