s2Eva_2023PAOI_T3 EDP elíptica, placa rectangular con frontera variable

Ejercicio: 2Eva_2023PAOI_T3 EDP elíptica, placa rectangular con frontera variable

Dada la EDP elíptica,

\frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{\partial y^2} = \Big( x^2 + y^2 \Big) e^{xy} 0 \lt x \lt 1 0 \lt y \lt 0.5

Se convierte a la versión discreta usando diferencias divididas centradas, según se puede mostrar con la gráfica de malla.

literal b

EDP eliptica rectangular frontera variable

literal a


\frac{u[i-1,j]-2u[i,j]+u[i+1,j]}{\Delta x^2} + + \frac{u[i,j-1]-2u[i,j]+u[i,j+1]}{\Delta y^2} = \Big( x^2 + y^2 \Big) e^{xy}

literal c

Se agrupan los términos Δx, Δy semejante a formar un λ al multiplicar todo por Δy2

\frac{\Delta y^2}{\Delta x^2}\Big(u[i-1,j]-2u[i,j]+u[i+1,j] \Big) + + \frac{\Delta y^2}{\Delta y^2}\Big(u[i,j-1]-2u[i,j]+u[i,j+1]\Big) = \Big( x^2 + y^2 \Big) e^{xy}\frac{\Delta y^2}{\Delta x^2}

los tamaños de paso en ambos ejes son de igual valor, se simplifica la ecuación

\lambda= \frac{\Delta y^2}{\Delta x^2} = 1

se simplifica el coeficiente en λ =1

u[i-1,j]-2u[i,j]+u[i+1,j] + + u[i,j-1]-2u[i,j]+u[i,j+1] = \Big( x^2 + y^2 \Big) e^{xy}

Se agrupan los términosiguales

u[i-1,j]-4u[i,j]+ u[i+1,j] + u[i,j-1] +u[i,j+1] = \Big( x^2 + y^2 \Big) e^{xy}

Se desarrollan las iteraciones para tres rombos y se genera el sistema de ecuacioens a resolver.
para i=1,j=1

u[0,1]-4u[1,1]+ u[2,1] + u[1,0] +u[1,2] = \Big( x[1]^2 + y[1]^2 \Big) e^{x[1]y[1]} 1-4u[1,1]+ u[2,1] + 1 + \frac{1}{8} = \Big( 0.25^2 + 0.25^2 \Big) e^{(0.25) (0.25)} -4u[1,1]+ u[2,1] = \Big( 0.25^2 + 0.25^2 \Big) e^{(0.25) (0.25)} - \frac{1}{8}

para i=2, j=1

u[1,1]-4u[2,1]+ u[3,1] + u[2,0] +u[2,2] = \Big( x[2]^2 + y[1]^2 \Big) e^{x[2]y[1]} u[1,1]-4u[2,1]+ u[3,1] + 1 + 0.25 = \Big( 0.5^2 + 0.25^2 \Big) e^{(0.5)(0.25)} u[1,1]-4u[2,1]+ u[3,1] = \Big( 0.5^2 + 0.25^2 \Big) e^{(0.5)(0.25)} -1.25

para i=3, j=1

u[2,1]-4u[3,1]+ u[4,1] + u[3,0] +u[3,2] = \Big( x[3]^2 + y[1]^2 \Big) e^{x[3]y[1]} u[2,1]-4u[3,1]+ 0.25 + 0 + \frac{3}{8} = \Big( 0.75^2 + 0.25^2 \Big) e^{(0.75)(0.25)} u[2,1]-4u[3,1] = \Big( 0.75^2 + 0.25^2 \Big) e^{(0.75)(0.25)} - 0.25 - \frac{3}{8}

con lo que se puede crear un sistema de ecuaciones y resolver el sistema para cada punto desconocido

\begin{pmatrix} -4 & 1 & 0 & \Big| & 0.008061807364732415 \\ 1 & -4 & 1 & \Big| & -0.8958911084166168 \\0 & 1 & -4 &\Big| & 0.1288939058881129 \end{pmatrix}

se obtiene los resultados para:

u[1,1] = 0.05953113
u[2,1] = 0.24618634
u[3,1] = 0.02932311

>>> import numpy as np
>>> (0.25**2+0.25**2)*np.exp(0.25*0.25) - 1/8
0.008061807364732415
>>> (0.5**2+0.25**2)*np.exp(0.5*0.25) - 1.25
-0.8958911084166168
>>> (0.75**2+0.25**2)*np.exp(0.75*0.25) - 0.25 -3/8
0.1288939058881129
>>> A=[[-4,1,0],[1,-4,1],[0.0,1.0,-4.0]]
>>> B = [0.008061807364732415, -0.8958911084166168, 0.1288939058881129]
>>> np.linalg.solve(A,B)
array([0.05953113, 0.24618634, 0.02932311])

2Eva_2023PAOI_T3 EDP elíptica, placa rectangular con frontera variable

2da Evaluación 2023-2024 PAO I. 29/Agosto/2023

Tema 3 (35 puntos) Aproxime la solución de la Ecuación Diferencial Parcial

\frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{\partial y^2} = \Big( x^2 + y^2 \Big) e^{xy} 0 \lt x \lt 1 0 \lt y \lt 0.5

Con las condiciones de frontera:

u(0,y)=1, u(1,y)= y, 0≤y≤0.5
u(x,0)=1, u(x,0.5)=x/2, 0≤x≤1

Aproxime la solución con tamaños de paso Δx = 0.25, Δy = 0.25
Utilice diferencias finitas centradas para las variables independientes x,y

a. Plantee las ecuaciones para usar un método numérico en un nodo i,j

b. Realice la gráfica de malla,

c. desarrolle y obtenga el modelo discreto para u(xi,tj)

d. Realice al menos tres iteraciones en el eje tiempo.

e. Estime el error de u(xi,tj) y adjunte los archivos del algoritmo y resultados.

Rúbrica: Aproximación de las derivadas parciales (5 puntos), construcción de la malla (10), construcción del sistema lineal (15), resolución del sistema (5 puntos).

Referencia: 2Eva_IT2012_T3 EDP elíptica, placa rectangular

s2Eva_2023PAOI_T2 Péndulo vertical amortiguado

Ejercicio: 2Eva_2023PAOI_T2 Péndulo vertical amortiguado

literal a

\frac{d^2 \theta}{dt^2} = -\mu \frac{d\theta}{ dt}-\frac{g}{L}\sin (\theta)

Se simplifica su forma a:

\frac{d\theta}{dt}= z = f_t(t,\theta,z) \frac{d^2\theta }{dt^2}= z' = -\mu z -\frac{g}{L}\sin (\theta) = g_t(t,\theta,z)

se usan los valores dados: g = 9.81 m/s2, L = 2 m

f_t(t,\theta,z) = z g_t(t,\theta,z) = - 0.5 z -\frac{9.81}{2}\sin (\theta)

y los valores iniciales para la tabla: θ(0) = π/4 rad, θ’ (0) = 0 rad/s, se complementan los valores en la medida que se aplica el desarrollo.

ti θ(ti) θ'(ti)=z
0 π/4 0
0.2 0.7161 -0.6583
0.4 0.5267 -0.1156
0.6 0.2579 -0.1410

literal b

Iteración 1:  ti = 0 ; yi = π/4 ; zi = 0

K1y = h * ft(ti,yi,zi) 
    = 0.2*(0) = 0
K1z = h * gt(ti,yi,zi) 
    = 0.2*(-0.5(0) -(9.81/2)sin (π/4) = -0.6930
        
K2y = h * ft(ti+h, yi + K1y, zi + K1z)
    = 0.2*(0-0.6930)= -0.1386
K2z = h * gt(ti+h, yi + K1y, zi + K1z)
    = 0.2*(-0.5(0-0.6930) -(9.81/2)sin(π/4-0) 
    = -0.6237

yi = yi + (K1y+K2y)/2 
   = π/4+ (0+(-0.1386))/2 = 0.7161
zi = zi + (K1z+K2z)/2 
   = 0+(-0.6930-0.6237)/2 = -0.6583
ti = ti + h = 0 + 0.2 = 0.2

estimado[i] = [0.2,0.7161,-0.6583]

Iteración 2:  ti = 0.2 ; yi = 0.7161 ; zi = -0.6583

K1y = h * ft(ti,yi,zi) 
    = 0.2*(-0.6583) = -0.1317
K1z = h * gt(ti,yi,zi) 
    = 0.2*(- 0.5 ( -0.6583) -(9.81/2)sin (0.7161) 
    = -0.5775
        
K2y = h * ft(ti+h, yi + K1y, zi + K1z)
    = 0.2*(-0.6583 -0.5775)= -0.2472
K2z = h * gt(ti+h, yi + K1y, zi + K1z)
    = 0.2*(- 0.5 (-0.6583 -0.5775) -(9.81/2)sin(0.7161-0.1317) 
    = -0.4171

yi = yi + (K1y+K2y)/2 
   = 0.7161 + (-0.1317-0.2472)/2 = 0.5267
zi = zi + (K1z+K2z)/2 
   = -0.6583+(-0.5775-0.4171)/2 = -0.1156
ti = ti + h = 0.2 + 0.2 = 0.4

estimado[i] = [0.4,0.5267,-0.1156]

Iteración 3:  ti = 0.4 ; yi = 0.5267 ; zi = -1.156

K1y = h * ft(ti,yi,zi) 
    = 0.2*(-1.156) = -0.2311
K1z = h * gt(ti,yi,zi) 
    = 0.2*(- 0.5(-1.156) -(9.81/2)sin (0.5267) 
    = -0.3771
        
K2y = h * ft(ti+h, yi + K1y, zi + K1z)
    = 0.2*(-1.156 -0.3771)= -0.3065
K2z = h * gt(ti+h, yi + K1y, zi + K1z)
    = 0.2*(- 0.5 ( -1.156 -0.3771) -(9.81/2)sin(0.5267-0.2311) 
    = -0.1322

yi = yi + (K1y+K2y)/2 
   = 0.5267 + (-0.2311-0.3065)/2 = 0.2579
zi = zi + (K1z+K2z)/2 
   = -1.156+(-0.3771-0.1322)/2 = -1.410
ti = ti + h = 0.4 + 0.2 = 0.6

estimado[i] = [0.6,0.2579,-1.410]

literal c

resultados del algoritmo:

[ t, 		 y, 	 dyi/dti=z,  K1y,	 K1z,	    K2y,      K2z]
[[ 0.000e+00  7.854e-01  0.000e+00  0.000e+00  0.000e+00  0.000e+00   0.000e+00]
 [ 2.000e-01  7.161e-01 -6.583e-01  0.000e+00 -6.930e-01 -1.386e-01  -6.237e-01]
 [ 4.000e-01  5.267e-01 -1.156e+00 -1.317e-01 -5.775e-01 -2.472e-01  -4.171e-01]
 [ 6.000e-01  2.579e-01 -1.410e+00 -2.311e-01 -3.771e-01 -3.065e-01  -1.322e-01]
 [ 8.000e-01 -3.508e-02 -1.377e+00 -2.820e-01 -1.089e-01 -3.038e-01    1.756e-01]
...

péndulo amortiguado 03

con h=0.2 se tienen 1/0.2 = 5 tramos por segundo, por lo que para 10 segundo serán 50 tramos. La cantidad de muestras = tramos + 1(valor inicial) = 51

con lo que se puede usar el algoritmo en EDO Runge-Kutta d2y/dx2

literal d

Se observa que la respuesta es oscilante y amortiguada en magnitud como se esperaba según el planteamiento. Con el tiempo se estabilizará en cero.

Instrucciones en Python

# 2Eva_2023PAOI_T2 Péndulo vertical amortiguado
import numpy as np
import matplotlib.pyplot as plt

def rungekutta2_fg(f,g,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,7),dtype=float)
    # incluye el punto [x0,y0,z0]
    estimado[0] = [x0,y0,z0,0,0,0,0]
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1y = h * f(xi,yi,zi)
        K1z = h * g(xi,yi,zi)
        
        K2y = h * f(xi+h, yi + K1y, zi + K1z)
        K2z = h * g(xi+h, yi + K1y, zi + K1z)

        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        xi = xi + h
        
        estimado[i] = [xi,yi,zi,K1y,K1z,K2y,K2z]
    return(estimado)

# INGRESO theta = y
g = 9.8
L  = 2
ft = lambda t,y,z: z
gt = lambda t,y,z: -0.5*z +(-g/L)*np.sin(y)

t0 = 0
y0 = np.pi/4
z0 = 0
h = 0.2
muestras = 51

# PROCEDIMIENTO
tabla = rungekutta2_fg(ft,gt,t0,y0,z0,h,muestras)

# SALIDA
np.set_printoptions(precision=3)
print(' [ t, \t\t y, \t dyi/dti=z, K1y,\t K1z,\t K2y,\t K2z]')
print(tabla)

# Grafica
ti = np.copy(tabla[:,0])
yi = np.copy(tabla[:,1])
zi = np.copy(tabla[:,2])
plt.subplot(121)
plt.plot(ti,yi)
plt.grid()
plt.xlabel('ti')
plt.title('yi')
plt.subplot(122)
plt.plot(ti,zi, color='green')
plt.xlabel('ti')
plt.title('dyi/dti')
plt.grid()
plt.show()

2Eva_2023PAOI_T2 Péndulo vertical amortiguado

2da Evaluación 2023-2024 PAO I. 29/Agosto/2023

Tema 2 (35 puntos) Una mejor aproximación a un péndulo oscilante con un ángulo θ más amplio y con un coeficiente de amortiguamiento μ se expresa con una ecuación diferencial ordinaria de segundo orden.

\frac{d^2 \theta}{dt^2} = -\mu \frac{d\theta}{ dt}-\frac{g}{L}\sin (\theta)

g = 9.81 m/s2
L = 2 m
θ(0) = π/4 rad
θ’ (0) = 0 rad/s

El péndulo se suelta desde el reposo, desde un ángulo de π/4 respecto al eje vertical. El coeficiente de amortiguamiento μ=0.5 es proporcional a la velocidad angular.

a. Realice el planteamiento del ejercicio usando Runge-Kutta de 2do Orden

b. Desarrolle tres iteraciones para θ(t) con tamaño de paso h=0.2

c. Usando el algoritmo, aproxime la solución entre t=0 a t=10 s, adjunte sus resultados en la evaluación.

d. Realice una observación sobre el movimiento estimado del péndulo a lo largo del tiempo.

Rúbrica: literal a (5 puntos), literal b (15 puntos), literal c (10 puntos), literal d (5 puntos)

Referencia: 2Eva_IT2019_T2 Péndulo vertical

Vista general de ecuaciones diferenciales I Capítulo 1, 6min 54s. 3Blue1Brown 31-Marzo-2023.

s2Eva_2023PAOI_T1 Material para medalla de academia

Ejercicio: 2Eva_2023PAOI_T1 Material para medalla de academia

medalla area con integral numerico f(x) = 2-8\Big( \frac{1}{2} - x \Big)^2

0 \le x \lt 1 g(x) = -\Big( 1-x\Big)\ln \Big( 1- x \Big)

Para f(x) se usará Simpson de 1/3 que requiere al menos dos tramos para aplicar:

a. Realice el planteamiento de las ecuaciones para el ejercicio.

I\cong \frac{h}{3}[f(x_0)+4f(x_1) + f(x_2)]

b. Describa el criterio usado para determinar el número de tramos usado en cada caso.

h = \frac{b-a}{2} = \frac{1-0}{2} = 0.5

c. Desarrolle las expresiones completas del ejercicio para cada función.

I_{fx}\cong \frac{0.5}{3}[f(0)+4f(0.5) + f(1)] f(0) = 2-8\Big( \frac{1}{2} - (0) \Big)^2 = 0 f(0.5) = 2-8\Big( \frac{1}{2} - (0.5) \Big)^2 = 2 f(1) = 2-8\Big( \frac{1}{2} - (1) \Big)^2 = 0 I_{fx} = \frac{1}{6}[0+4(2) + 0] = \frac{8}{6} = \frac{4}{3} = 1.3333

cota de error O(h5) = O(0.55)= O(0.03125)

Para g(x) se usará Simpson de 3/8 que requiere al menos tres tramos para aplicar:

I\cong \frac{3h}{8}[f(x_0)+3f(x_1) +3 f(x_2)+f(x_3)] h = \frac{b-a}{3} = \frac{1-0}{3} = 0.3333 I_{gx}\cong \frac{3(0.3333)}{8}[f(0)+3f(0.3333) +3 f(0.6666)+f(1)] g(0) = -\Big( 1-0\Big)\ln \Big( 1- 0 \Big) = 0 g(0.3333) = -\Big( 1-0.3333\Big)\ln \Big( 1- 0.3333 \Big) = 0.2703 g(0.6666) = -\Big( 1-0.6666\Big)\ln \Big( 1- 0.6666 \Big) = 0.3662 g(0.9999) = -\Big( 1-0.9999\Big)\ln \Big( 1- 0.9999 \Big) = 0

para la evaluación numérica de 1 se usa un valor muy cercano desplazado con la tolerancia aplicada.

I_{gx}\cong \frac{3(0.3333)}{8}[0+3(0.2703) + 3(0.3662)+0] = 0.2387

d. Indique el resultado obtenido para el área requerida y la cota de error
Area = I_{fx} – I_{gx} = 1.3333 – 0.2387 = 1.0945

cota de error = O(0.03125) + O(0.00411) = 0.03536

e. Encuentre el valor del tamaño de paso si se requiere una cota de error de 0.00032

Si el factor de mayor error es de Simpson 1/3, se considera como primera aproximación que:

cota de error O(h5) = O(0.00032), h = (0.00032)(1/5) = 0.2
es decir el número de tramos es de al menos (b-a)/tramos = 0.2 , tramos = 5.
El número de tramos debe ser par en Simpson de 1/3, por lo que se toma el entero mayor tramos=6 y el tamaño de paso recomendado es al menos 1/6. EL error al aplicar 3 veces la formula es 3(O((1/6)5)) = 0.0003858.

Lo que podría indicar que es necesario al menos dos tramos adicionales con h=1/8 y error O(0,00012) que cumple con el requerimiento.

Se puede aplicar el mismo criterio para Simpson 3/8 y se combinan los errores para verificar que cumplen con el requerimiento.

Algoritmo con Python

Resultados

Ifx:  1.3333332933359998
Igx:  0.238779092876627
Area:  1.094554200459373

medalla area con integral numerico

Instrucciones en Python usando las funciones

# 2Eva_2023PAOI_T1 Material para medalla de academia
import numpy as np
import matplotlib.pyplot as plt

def integrasimpson13_fi(xi,fi,tolera = 1e-10):
    ''' sobre muestras de fi para cada xi
        integral con método de Simpson 1/3
        respuesta es np.nan para tramos desiguales,
        no hay suficientes puntos.
    '''
    n = len(xi)
    i = 0
    suma = 0
    while not(i>=(n-2)):
        h = xi[i+1]-xi[i]
        dh = abs(h - (xi[i+2]-xi[i+1]))
        if dh<tolera:# tramos iguales
            unS13 = (h/3)*(fi[i]+4*fi[i+1]+fi[i+2])
            suma = suma + unS13
        else:  # tramos desiguales
            suma = 'tramos desiguales'
        i = i + 2
    if i<(n-1): # incompleto, faltan tramos por calcular
        suma = 'tramos incompletos, faltan '
        suma = suma ++str((n-1)-i)+' tramos'
    return(suma)

def integrasimpson38_fi(xi,fi,tolera = 1e-10):
    ''' sobre muestras de fi para cada xi
        integral con método de Simpson 3/8
        respuesta es np.nan para tramos desiguales,
        no hay suficientes puntos.
    '''
    n = len(xi)
    i = 0
    suma = 0
    while not(i>=(n-3)):
        h  = xi[i+1]-xi[i]
        h1 = (xi[i+2]-xi[i+1])
        h2 = (xi[i+3]-xi[i+2])
        dh = abs(h-h1)+abs(h-h2)
        if dh<tolera:# tramos iguales
            unS38 = fi[i]+3*fi[i+1]+3*fi[i+2]+fi[i+3]
            unS38 = (3/8)*h*unS38
            suma = suma + unS38
        else:  # tramos desiguales
            suma = 'tramos desiguales'
        i = i + 3
    if (i+1)<n: # incompleto, tramos por calcular
        suma = 'tramos incompletos, faltan '
        suma = suma +str(n-(i+1))+' tramos'
    return(suma)

# INGRESO
fx = lambda x: 2-8*(0.5-x)**2
gx = lambda x: -(1-x)*np.log(1-x)
a = 0
b = 1-1e-4
muestras1 = 2+1
muestras2 = 3+1

# PROCEDIMIENTO
xi1 = np.linspace(a,b,muestras1)
xi2 = np.linspace(a,b,muestras2)
fi = fx(xi1)
gi = gx(xi2)

Ifx = integrasimpson13_fi(xi1,fi)
Igx = integrasimpson38_fi(xi2,gi)
Area = Ifx - Igx

# SALIDA
print('Ifx: ', Ifx)
print('Igx: ', Igx)
print('Area: ', Area)

plt.plot(xi1,fi,'ob',label='f(x)')
plt.plot(xi2,gi,'or', label='g(x)')
plt.grid()
plt.legend()
plt.xlabel('xi')

# curvas suave con mas muestras (no en evaluación)
xi = np.linspace(a,b,51)
fxi = fx(xi)
gxi = gx(xi)
plt.fill_between(xi,fxi,gxi,color='navajowhite')
plt.plot(xi,fxi,color='blue',linestyle='dotted')
plt.plot(xi,gxi,color='red',linestyle='dotted')

plt.show()

2Eva_2023PAOI_T1 Material para medalla de academia

2da Evaluación 2023-2024 PAO I. 29/Agosto/2023

Tema 1 (30 puntos) medalla area con integral numerico
Una academia encarga a un joyero un modelo de medalla cuyo costo unitario se determina por el área descrita entre las funciones f(x) y g(x) presentadas.

Se considera que el grosor de la medalla es único e independiente de la forma de la medalla.

f(x) = 2-8\Big( \frac{1}{2} - x \Big)^2 0 \le x \lt 1 g(x) = -\Big( 1-x\Big)\ln \Big( 1- x \Big)

Para el desarrollo numérico, use diferentes métodos de Simpson para cada función.

a. Realice el planteamiento de las ecuaciones para el ejercicio.

b. Describa el criterio usado para determinar el número de tramos usado en cada caso.

c. Desarrolle las expresiones completas del ejercicio para cada función.

d. Indique el resultado obtenido para el área requerida y la cota de error.

e. Encuentre el valor del tamaño de paso si se requiere una cota de error de 0.00032

Nota: en Python ln(x) se escribe np.log(x).

Rúbrica: literal a (5 puntos), literal b (5 puntos), literal c (10 puntos), literal d (5 puntos), literal e (5 puntos)

Referencia: Star Trek https://intl.startrek.com/
¿A quien se le ocurrió crear la moneda? | Discovery en Español Youtube.com 8 nov 2016.

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()

 

1Eva_2023PAOI_T3 Recoger los restos de sumergible

1ra Evaluación 2023-2024 PAO I. 4/Julio/2023

Tema 3 (35 puntos) Suponga que se requiere planificar la operación de recolección de restos del accidente del sumergible ocurrido en junio del 2023 usando un robot autónomo sumergible. rov sumergible

Existen limitaciones debidas a la profundidad del lecho submarino, las fuertes corrientes que fluyen desde el oeste, limitada energía en las baterías del robot, etc.

El recorrido del robot se puede programar usando un polinomio para el intervalo de las coordenadas presentadas. recorrido de polinomio

Se estima poder realizar tan solo dos recorridos desde el oeste hacia el este.

a. Plantee el ejercicio usando interpolación polinómica para realizar un recorrido para la recolección de restos, describiendo los criterios usados para la selección de los puntos.

b. Desarrolle un primer polinomio usando el método de diferencias divididas de Newton

c. Excluyendo los puntos usados en el literal anterior, desarrolle otro polinomio usando el método de Lagrange.

d. Grafique los resultados obtenidos para cada caso.

e. Determine el error de cada polinomio como la distancia entre cada punto por el que no pasa el polinomio y la coordenada equivalente en el polinomio. Realice observaciones a los resultados.

xi = 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])
fi = 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])

Nota: Los valores de la tabla se establecen para el ejercicio y no corresponden a una referencia publicada.

Rúbrica: literal a (5 puntos), literal b (10 puntos), literal c (10 puntos), literal d (5 puntos), literal e (5 puntos)

Referencia: [1] La Guardia Costera de EE.UU. anuncia una investigación oficial sobre la implosión del sumergible Titán. RTVE.es / EFE. 26.06.2023 https://www.rtve.es/noticias/20230626/eeuu-anuncia-investigacion-oficial-sumergible-titan/2450440.shtml

[2] Así será la recuperación de los restos del sumergible Titán, según un experto militar. CNN en Español. 23 jun 2023.

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)

1Eva_2023PAOI_T2 Productos en combo por subida de precio

1ra Evaluación 2023-2024 PAO I. 4/Julio/2023

Tema 2 (35 puntos) Un restaurante usa materia prima como: arroz, cebolla, papa, huevos, leche, etc.mercado alimentos

Sin embargo los precios han presentado variaciones en el mercado presentando como causas: las lluvias fuertes, el fenómeno del niño, daños en las cosechas y variaciones en insumos por la guerra Ucrania.

Un proveedor mayorista ofrece “combo” de productos por un valor total registrado en varias semanas como:

arroz [Kg] cebolla [Kg] papa [Kg] huevo [docena] Total Pagado [USD]
Semana 1 500 600 400 90 1660
Semana 2 800 450 300 100 1825
Semana 3 400 300 600 80 1430

Nota: Los valores de la tabla se establecen para el ejercicio y no corresponden a una referencia publicada.

Se requiere conocer el precio unitario de los insumos para estimar el presupuesto operativo semanal del restaurante.

a. Plantee el ejercicio usando los datos de la tabla. Considere justificar el uso de una variable libre.

b. Establezca la forma matricial del sistema de ecuaciones y como matriz aumentada

c. De ser necesario realice el pivoteo parcial por filas

d. Use el método iterativo de Gauss-Seidel, realizando al menos 3 iteraciones. Estime la tolerancia, y justifique.

e. Comente sobre la convergencia del método y justifique sus observaciones usando los errores entre iteraciones.

Nota: considere el precio de 50 USD para el quintal métrico (100Kg) de arroz, y precios menores del mismo orden de magnitud para los demás productos. La docena de huevo a 2,5 USD.

 Rúbrica: literal a (2 puntos), literal b (5 puntos), literal c (3puntos), literal d (15 puntos), literal e (5 puntos)

Referencia: [1] Papa, arroz, huevos, cebolla y gas escasean en las islas Galápagos por problemas de abastecimiento. Eluniverso.com. 9 Junio 2023. https://www.eluniverso.com/noticias/ecuador/escasez-de-productos-en-islas-galapagos-persiste-nota/

[2] El Gobierno importará arroz para controlar que el precio no suba por especulación o escasez. Ecuavisa. 16 junio 2023. https://www.ecuavisa.com/noticias/ecuador/el-gobierno-importara-arroz-para-controlar-que-el-precio-no-suba-por-especulacion-o-escasez-CK5396253

[3] Estos son los alimentos de la canasta básica que se han encarecido por las lluvias. Ecuavisa. 05 junio 2023. https://www.ecuavisa.com/la-noticia-a-fondo/estos-son-los-alimentos-de-la-canasta-basica-que-se-han-encarecido-por-las-lluvias-FX5303236

[4] Cebolla, arroz, piña y tomate son los productos que han subido su precio en mercados de Quito. La Hora. Junio 23, 2023. https://www.lahora.com.ec/pais/cebolla-arroz-pina-y-tomate-son-los-productos-que-han-subido-su-precio-en-mercados-de-quito/

[5] Escasez de arroz ya se refleja en los mercados de Ecuador | Televistazo | Ecuavisa. 14 jun 2023