Edison Del Rosario – Página 4 – Métodos numéricos

2023-02-072023-02-07

3Eva_2022PAOII_T4 Recesión económica, PIB y diferenciación numérica

3ra Evaluación 2022-2023 PAO II. 7/febrero/2023

Tema 4. (20 puntos) Las recesiones económicas se caracterizan por presentar una disminución del consumo, de la inversión y de la producción de bienes y servicios. Lo cual provoca, a su vez, que se despidan trabajadores y por tanto, aumente el desempleo. La recesión es también conocida como el periodo de «vacas flacas».

La opinión emitida por Julius Shiskin en un artículo publicado el 28 de agosto de 1975 en el diario New York Times en torno a dos trimestres consecutivos de caída del PIB como plazo definitorio para el considerar una recesión.

a. Plantee y describa los métodos de diferenciación numérica que usen dos y tres puntos para primera derivada.

b. Realice tres iteraciones con los métodos numéricos seleccionados. Describa el tamaño de paso usado en cada método.

c. Compare los resultados del literal anterior y escriba sus observaciones respecto a las cotas de error.

d. Determine los periodos de “recesión económica” para los datos proporcionados entre el año 2012 y 20122. Liste acorde a lo definido por J. Shiskin. Use los algoritmos y adjunte los archivos usados (py,txt,png).

trimestres = ['2012.I', '2012.II', '2012.III', '2012.IV',
 '2013.I', '2013.II', '2013.III', '2013.IV', '2014.I',
 '2014.II', '2014.III', '2014.IV', '2015.I', '2015.II',
 '2015.III', '2015.IV', '2016.I', '2016.II', '2016.III',
 '2016.IV', '2017.I', '2017.II', '2017.III', '2017.IV',
 '2018.I', '2018.II', '2018.III', '2018.IV', '2019.I',
 '2019.II', '2019.III', '2019.IV', '2020.I', '2020.II',
 '2020.III', '2020.IV', '2021.I', '2021.II', '2021.III',
 '2021.IV', '2022.I', '2022.II', '2022.III']

PIB trimestral = [21.622937, 21.908844, 22.106937,
 22.285826, 23.019786, 23.441324, 24.238576, 24.429973,
 24.829431, 25.540887, 25.9404, 25.415613, 25.052739,
 25.086195, 24.779738, 24.371709, 24.913573, 24.926186,
 24.910741, 25.187196, 26.000261, 25.99355, 25.960907,
 26.341144, 26.510612, 26.761827, 27.078404, 27.211165,
 26.914897, 27.058331, 27.054758, 27.080023, 26.314576,
 23.110752, 24.64388, 25.221916, 25.412756, 26.20682,
 26.828611, 27.717679, 28.372038, 28.74092, 29.334581]

Rúbrica: literal a (3 puntos), literal b (9 puntos). literal c(3 puntos) literal d (5 puntos)

Referencia: Recesión económica. Wikipedia. https://es.wikipedia.org/wiki/Recesi%C3%B3n
Recesión económica. economipedia. https://economipedia.com/definiciones/recesion-economica.html

Boletín de Cuentas Nacionales Trimestrales No. 121, valores constantes USD 2007 y corrientes, período : 2000.I – 2022.IIIIT Banco Central del Ecuador (2022) https://contenido.bce.fin.ec/documentos/PublicacionesNotas/Catalogo/CuentasNacionales/Indices/c121122022.htm

2023-02-072023-02-11

3Eva_2022PAOII_T3 EDO cabezal lector en disco duro

3ra Evaluación 2022-2023 PAO II. 7/febrero/2023

Tema 3. (35 puntos) El objetivo de un sistema de Disco duro es posicionar con precisión el dispositivo de lectura en la pista buscada y moverse entre una pista y otra.

Se requiere identificar el plato de disco, el sensor y el controlador.

El disco duro usa un motor DC de imán permanente para posicionar el brazo lector con el sensor en un extremo. Un resorte metálico se usa para permitir que el cabezal flote sobe el disco a una distancia menor a 100nm.

El cabezal toma lectura del flujo magnético y da una señal al amplificador.

Suponiendo que dispone del dispositivo de lectura de precisión, una aproximación del modelo de control del motor con K_a=40, se supone que el brazo es rígido con parámetros como los mostrados, el sistema se puede aproximar como un sistema de orden 2 en el dominio s o en su forma de ecuación diferencial.

Y(s)(s^2+20s+5K_a )=X(s)5K_a

\frac{\delta^2}{\delta t^2 } y(t) + 20 \frac{\delta}{\delta t} y(t) + 5 K_a y(t) = x(t) 5 K_a

y(0) = 0 y’(0) = 0

x(t) = \begin{cases} 0 & t\lt 0 \\ 1 & t≥0 \end{cases}

Encuentre la respuesta del sistema y(t) ante una señal de entrada x(t), con las condiciones iniciales dadas.

a. Plantee la solución usando el método de Runge-Kutta de 2do orden.
b. Desarrolle tres iteraciones para Δt = 0.01
c. Estime el error del modelo usado
d. Realice la gráfica para y(t) para el intervalo de [0,1] segundos. Adjunte los archivos de los algoritmos.py usados para los cálculos, los resultados.txt y gráfica.png

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

Referencia: Bishop R. & Dorf R. (2017) 13th Edition. 2.10 sequential Design example: Disk Drive read system. p122.
How do Hard Disk Drives Work? Branch Education. 22 diciembre 2022.

https://youtu.be/wtdnatmVdIg

2023-02-072023-02-11

3Eva_2022PAOII_T2 Globo meteorológico espía distancia mas cercana

3ra Evaluación 2022-2023 PAO II. 7/febrero/2023

Tema 2. (20 puntos) Se requiere hacer el seguimiento a la trayectoria del globo aerostático del tema anterior, para descartar las sospechas de espionaje.

Dadas las coordenadas de un lugar considerado como de seguridad nacional p1(x,y)=[25,50] , se requiere revisar la distancia más cercana de la trayectoria y(x) del globo al punto de “interés”.

Se compararía la distancia mínima con el alcance las cámaras y sensores encontrados en los escombros del globo derribado.

Usando la trayectoria obtenida como resultado del tema anterior, se requiere:

a. Plantee el ejercicio describiendo los criterios usados, el método numérico y una tolerancia a usar.

b. Desarrolle el método para encontrar la raíz de la ecuación planteada, con al menos tres iteraciones.

c. Estime la cota de error, compare con la tolerancia descrita en el literal a.

Nota: Si el resultado del tema 1 no es satisfactorio, desarrolle el tema con y(x) = 70sin(0.1πx+0.5)

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

2023-02-072023-02-11

3Eva_2022PAOII_T1 Globo meteorológico espía derribado

3ra Evaluación 2022-2023 PAO II. 7/febrero/2023

Tema 1. (25 puntos) En enero del 2023 se detectó un globo aerostático supuestamente espía sobre el territorio soberano de un país, que sobrevoló a 18 Km de altura en la estratosfera y que «no representaba ningún riesgo militar o físico los ciudadanos en la superficie».

Otro país vecino al mismo tiempo hacía seguimiento a otro «posible segundo incidente», se anunció en los medios de comunicación. En el primer caso se decidió no destruir el aparato por el temor de que la caída de sus escombros podría haber sido peligrosa para la superficie y no representaba el globo un peligro inmediato.

Como seguimiento al caso, se requiere describir la trayectoria del globo mediante ecuaciones a partir de las coordenadas de avistamiento reportadas por civiles.

ti    = [11, 12, 14, 15, 17, 19]
x(ti) = [15, 18, 25, 27, 31, 40]
y(ti) = [45, 55, 65, 58, 55, 40]

a. Plantear el ejercicio, describiendo los criterios, método numérico, segmentos a usar en las ecuaciones para realizar la interpolación polinómica de Lagrange.

minimizando oscilaciones del polinomio que puedan resultar en interpretaciones erradas.

b. Realizar el desarrollo analítico de las ecuaciones planteadas y presente el polinomio simplificado.

c. Validar los resultados usando el algoritmo, adjunte los archivos.py, resultados.txt, gráfica.png

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

Referencias: Detectan un globo aerostático espía sobre territorio. Rtve.es/Agencias 03/febrero/2023. https://www.rtve.es/noticias/20230203/eeuu-detecta-globo-aerostatico-espia-china-sobre-su-territorio/2420646.shtml

Derriban globo «espía» sobre la costa del Atlántico. DW 04/febrero/2023. https://www.dw.com/es/eeuu-derriba-globo-esp%C3%ADa-chino-sobre-la-costa-del-atl%C3%A1ntico/a-64613403

EE.UU. derriba el presunto globo espía de China. CNN en Español. 4 feb 2023.

Globos chinos en América desatan preocupación mundial. DW Español
DW Español. 10 feb 2023

2023-01-242023-01-27

2Eva_2022PAOII_T3 EDP Parabólica con coseno 3/4π

2da Evaluación 2022-2023 PAO II. 24/Enero/2023

Tema 3. (35 puntos) Aproxime la solución a la siguiente ecuación diferencial parcial parabólica

\frac{\partial^2 u}{\partial x^2} = b \frac{\partial u}{\partial t}

Con las siguientes condiciones de frontera:
u(0,t)=1
u(1,t)=0

Y las condiciones iniciales
$u(x,0) = \cos \Big( \frac{3π}{2}x\Big)$

Utilice diferencias finitas centradas para x, para t hacia adelante.

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)

Suponga que b = 2, Aproxime la solución con Δx = 0.2, Δt = Δx/100.

d. Realice al menos tres iteraciones en el eje tiempo.
e. Estime el error de u(xi,tj), y presente observaciones sobre la convergencia del método.

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

Referencia: Chapra & R. Canale (2010). Métodos Numéricos para Ingenieros. Ejercicio 30.15 p904,
Solving the heat equation | DE3. 3Blue1Brown 16 Junio 2019.

2023-01-242023-01-27

s2Eva_2022PAOII_T3 EDP Parabólica con coseno 3/4π

Ejercicio: 2Eva_2022PAOII_T3 EDP Parabólica con coseno 3/4π

\frac{\partial^2 u}{\partial x^2} = b \frac{\partial u}{\partial t}

\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} = b\frac{u_{i,j+1}-u_{i,j}}{\Delta t}

agrupando variables

\frac{\Delta t}{b} \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} = \frac{\Delta t}{b}b\frac{u_{i,j+1}-u_{i,j}}{\Delta t}

λ = \frac{\Delta t}{b(\Delta x)^2}

λ = \frac{0.002}{2(0.2)^2} =0.025

como λ<0.5 el método converge.

\lambda \Big[u[i+1,j]-2u[i,j]+u[i-1,j]\Big] = u[i,j+1]-u[i,j]

por el método explícito:

u[i,j+1] =\lambda \Big[u[i+1,j]-2u[i,j]+u[i-1,j]\Big] + u[i,j]

u[i,j+1] =\lambda u[i+1,j]+(1-2\lambda)u[i,j]+\lambda u[i-1,j]

iteración i=1, j=0

u[1,1] =\lambda u[0,0]+(1-2\lambda)u[1,0]+\lambda u[2,0]

u[1,1] =0.025(1) +(1-2(0.025))\cos \Big( \frac{3π}{2}0.2\Big)+0.025 \cos \Big( \frac{3π}{2}0.4\Big)

iteración i=2, j=0

u[2,1] =\lambda u[1,0]+(1-2\lambda)u[2,0]+\lambda u[3,0]

u[2,1] =0.025\cos \Big( \frac{3π}{2}0.2\Big) +(1-2(0.025))\cos \Big( \frac{3π}{2}0.4\Big)

+0.025 \cos \Big( \frac{3π}{2}0.6\Big)

iteración i=3, j=0

u[3,1] =\lambda u[2,0]+(1-2\lambda)u[3,0]+\lambda u[4,0]

u[3,1] =0.025\cos \Big( \frac{3π}{2}0.4\Big) +(1-2(0.025))\cos \Big( \frac{3π}{2}0.6\Big)

+0.025 \cos \Big( \frac{3π}{2}0.8\Big)

iteración i=4, j=0

u[4,1] =\lambda u[3,0]+(1-2\lambda)u[4,0]+\lambda u[5,0]

u[4,1] =0.025\cos \Big( \frac{3π}{2}0.6\Big) +(1-2(0.025))\cos \Big( \frac{3π}{2}0.8\Big)

+0.025 (0)

continuar con las iteraciones en el algoritmo

Resultados con el algoritmo

Tabla de resultados
[[ 1.    1.    1.    1.    1.    1.    1.    1.    1.    1.  ]
 [ 0.59  0.58  0.56  0.55  0.54  0.53  0.53  0.52  0.51  0.5 ]
 [-0.31 -0.3  -0.3  -0.29 -0.28 -0.28 -0.27 -0.27 -0.26 -0.26]
 [-0.95 -0.93 -0.91 -0.89 -0.88 -0.86 -0.84 -0.82 -0.81 -0.79]
 [-0.81 -0.79 -0.78 -0.76 -0.74 -0.73 -0.71 -0.7  -0.68 -0.67]
 [ 0.    0.    0.    0.    0.    0.    0.    0.    0.    0.  ]]

Instrucciones en Python

# EDP parabólicas d2u/dx2  = K du/dt
# método explícito,usando diferencias divididas
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# Valores de frontera
Ta = 1
Tb = 0
#T0 = 25
fx = lambda x: np.cos(3*np.pi/2*x)
# longitud en x
a = 0
b = 1
# Constante K
K = 2
# Tamaño de paso
dx = 0.2
dt = dx/100
# iteraciones en tiempo
n = 10

# PROCEDIMIENTO
# iteraciones en longitud
xi = np.arange(a,b+dx,dx)
fi = fx(xi)
m = len(xi)
ultimox = m-1

# Resultados en tabla u[x,t]
u = np.zeros(shape=(m,n), dtype=float)

# valores iniciales de u[:,j]
j=0
ultimot = n-1
u[0,:]= Ta
u[1:ultimox,j] = fi[1:ultimox]
u[ultimox,:] = Tb

# factores P,Q,R
lamb = dt/(K*dx**2)
P = lamb
Q = 1 - 2*lamb
R = lamb

# Calcula U para cada tiempo + dt
j = 0
while not(j>=ultimot): # igual con lazo for
    for i in range(1,ultimox,1):
        u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]
    j=j+1

# SALIDA
print('Tabla de resultados')
np.set_printoptions(precision=2)
print(u)

# Gráfica
salto = int(n/10)
if (salto == 0):
    salto = 1
for j in range(0,n,salto):
    vector = u[:,j]
    plt.plot(xi,vector)
    plt.plot(xi,vector, '.r')
    
plt.xlabel('x[i]')
plt.ylabel('t[j]')
plt.title('Solución EDP parabólica')
plt.show()

2023-01-242023-01-26

2Eva_2022PAOII_T2 EDO – población de protestantes en una sociedad

2da Evaluación 2022-2023 PAO II. 24/Enero/2023

Tema 1. (35 puntos) En el libro titulado “Looking at History Through Mathematics”, Rashevsky propone un modelo que se puede relacionar con el “protestantismo” en el siglo XVI como una reacción y denuncia de abusos impuestos sobre la sociedad de la época.

En un modelo de Rashevsky modificado con la ecuación logística de Verhulst, la población x(t) de individuos en la sociedad para cada año t, con tasas de natalidad b=0.02 y mortalidad d=0.015, cambia según la ecuación:

\frac{\delta}{\delta t}x(t) = b x(t) - d (x(t))^2

x(0)=1

La cantidad de individuos “protestantes” y(t) en la población se incrementa según la ecuación diferencial compuesta de dos términos.

\frac{\delta}{\delta t}y(t) = b y(t) - d (y(t))^2 +r b (x(t)-y(t))

y(0)=0.01

El primer término supone que todas familias de padre y madre “protestantes” tienen hijos que también se identifican como tales.

El segundo término supone que una porción r = 0.1 de jóvenes descendientes de los “conformistas” al meditar sobre la situación actual, los hechos y los argumentos de protesta se convierten a “protestantes”.

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

b. Desarrolle tres iteraciones para x(t), y(t) con tamaño de paso h=0.5.

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

d. Realice una observación sobre el crecimiento de población y(t) a lo largo del tiempo.

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

Referencia: Burden 5.2 Ejercicio 17 p276, Rashevsky, MIT 1968. pp102-110, Protestantismo https://es.wikipedia.org/wiki/Protestantismo. 3Eva_IIT2014_T2 Crecimiento demográfico. http://blog.espol.edu.ec/analisisnumerico/3eva_iit2014_t2-crecimiento-demografico/

La Reforma protestante y Lutero. Academia Play. 27 agosto 2019

2023-01-242023-01-29

s2Eva_2022PAOII_T2 EDO – población de protestantes en una sociedad

Ejercicio: 2Eva_2022PAOII_T2 EDO – población de protestantes en una sociedad

\frac{\delta}{\delta t}x(t) = b x(t) - d (x(t))^2

\frac{\delta}{\delta t}y(t) = b y(t) - d (y(t))^2 +r b (x(t)-y(t))

literal a

simplificando la nomenclatura

x' = b x - d x^2

y' = b y - d y^2 +r b (x-y)

sustituyendo constantes, y considerando x(0)=1 ; y(0)=0.01 ; h=0.5

x' = 0.02 x - 0.015 x^2

y' = 0.02 y - 0.015 y^2 +0.1(0.02) (x-y)

el planteamiento de Runge Kutta se hace junto a la primera iteración, además de encontrarse en las instrucciones con Python.

literal b

Se describen 3 iteraciones usando los resultados de la tabla con el algoritmo, para mostrar la comprensión del algoritmo.

t = 0
K1x = 0.5 (0.02 (1) – 0.015 (1)² = 0.0025
K1y = 0.5(0.02 (0.01) – 0.015 (0.01)² +0.1(0.02) (1-0.01)= 0.001089

K2x = 0.5 (0.02 (1+0.0025) – 0.015 (1+0.0025)²= 0.00248
K2y = 0.5(0.02 (0.01+0.00108) – 0.015 (0.01+0.00108)²+0.1(0.02) ((1+0.0025)-(0.01+0.00108)) = 0.001101

x₁ = 1 + (1/2)(0.0025+0.00248) = 1.0025
y₁ = 0.01 + (1/2)(0.001089+0.001101) = 0.01109
t₁ = 0 + 0.5 =0.5

t=0.5
K1x = 0.5 (0.02 (1.0025) – 0.015 (1.0025)² = 0.002487
K1y = 0.5(0.02 (0.01109) – 0.015 (0.01109)² +0.1(0.02) (1.0025-0.01109)= 0.001101

K2x = 0.5 (0.02 (1.0025+ 0.002487) – 0.015 (1.0025+ 0.002487)²= 0.002474
K2y = 0.5(0.02 (0.01109+0.001101) – 0.015 (0.01109+0.001101)²+0.1(0.02) ((1.0025+ 0.002487)-(0.01109+0.001101)) = 0.001113

x₂ = 1.0025 + (1/2)(0.002487+0.002474) = 1.0050
y₂ = 0.01109 + (1/2)(0.001101+0.001113) = 0.01220
t₂ = 0.5 + 0.5 = 1

t=1
K1x = 0.5 (0.02 (1.0050) – 0.015 (1.0050)² = 0.002474
K1y = 0.5(0.02 (0.01220) – 0.015 (0.01220)² +0.1(0.02) (1.0050-0.01220)= 0.001113

K2x = 0.5 (0.02 (1.0050+ 0.002474) – 0.015 (1.0050+ 0.002474)²= 0.002462
K2y = 0.5(0.02 (0.01220+0.001113) – 0.015 (0.01220+0.001113)²+0.1(0.02) ((1.0050+ 0.002474)-(0.01220+0.001113)) = 0.001126

x₃ = 1.0050 + (1/2)(0.002474+0.002462) = 1.0074
y₃ = 0.01220 + (1/2)(0.001113+0.001126) = 0.01332
t₃ = 1 + 0.5 = 1.5

Resultado con el algoritmo

Para obtener los datos de las iteraciones, primero se ejecuta el algoritmo para pocas iteraciones.
Para la pregunta sobre 200 años, se incrementa las iteraciones a 2 por año y las condiciones iniciales, es decir 401 muestras.

 [ ti, xi, yi]
 [ ti, xi, yi]
[[0.0000e+00 1.0000e+00 1.0000e-02 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00]
 [5.0000e-01 1.0025e+00 1.1095e-02 2.5000e-03 1.0892e-03 2.4875e-03 1.1014e-03]
 [1.0000e+00 1.0050e+00 1.2203e-02 2.4875e-03 1.1014e-03 2.4749e-03 1.1136e-03]
 [1.5000e+00 1.0074e+00 1.3323e-02 2.4749e-03 1.1137e-03 2.4623e-03 1.1260e-03]
 [2.0000e+00 1.0099e+00 1.4455e-02 2.4624e-03 1.1260e-03 2.4497e-03 1.1384e-03]
 [2.5000e+00 1.0123e+00 1.5600e-02 2.4498e-03 1.1384e-03 2.4371e-03 1.1509e-03]
 [3.0000e+00 1.0148e+00 1.6757e-02 2.4371e-03 1.1509e-03 2.4245e-03 1.1634e-03]
 [3.5000e+00 1.0172e+00 1.7926e-02 2.4245e-03 1.1635e-03 2.4118e-03 1.1761e-03]
 [4.0000e+00 1.0196e+00 1.9109e-02 2.4118e-03 1.1761e-03 2.3991e-03 1.1888e-03]
...
 [1.9950e+02 1.3252e+00 1.1561e+00 ... 1.7202e-03 8.1217e-05 1.7059e-03]
 [2.0000e+02 1.3252e+00 1.1578e+00 ... 1.7060e-03 8.0418e-05 1.6918e-03]
 [2.0050e+02 1.3253e+00 1.1595e+00 ... 1.6919e-03 7.9628e-05 1.6778e-03]]
>>>

Observación: La población identificada como protestante, continua creciendo, mientras que la proporción de «conformistas» se reduce según los parámetros indicados en el ejercicio. Los valores de natalidad y defunción cambian con el tiempo mucho más en años por otras variables, por lo que se deben realizar ajustes si se pretende extender el modelo.

Instrucciones en Python

# Modelo predador-presa de Lotka-Volterra
# Sistemas EDO con Runge Kutta de 2do Orden
import numpy as np

def rungekutta2_fg(f,g,t0,x0,y0,h,muestras):
    tamano = muestras +1
    tabla = np.zeros(shape=(tamano,7),dtype=float)
    tabla[0] = [t0,x0,y0,0,0,0,0]
    ti = t0
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1x = h * f(ti,xi,yi)
        K1y = h * g(ti,xi,yi)
        
        K2x = h * f(ti+h, xi + K1x, yi+K1y)
        K2y = h * g(ti+h, xi + K1x, yi+K1y)

        xi = xi + (1/2)*(K1x+K2x)
        yi = yi + (1/2)*(K1y+K2y)
        ti = ti + h
        
        tabla[i] = [ti,xi,yi,K1x,K1y,K2x,K2y]
    tabla = np.array(tabla)
    return(tabla)

# PROGRAMA ------------------

# INGRESO
# Parámetros de las ecuaciones
b = 0.02
d = 0.015
r = 0.1

# Ecuaciones
f = lambda t,x,y : (b-d*x)*x
g = lambda t,x,y : (b-d*y)*y + r*b*(x-y)

# Condiciones iniciales
t0 = 0
x0 = 1
y0 = 0.01

# parámetros del algoritmo
h = 0.5
muestras = 401

# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,t0,x0,y0,h,muestras)
ti = tabla[:,0]
xi = tabla[:,1]
yi = tabla[:,2]

# SALIDA
np.set_printoptions(precision=6)
print(' [ ti, xi, yi, K1x, K1y, K2x, K2y]')
print(tabla)

# Grafica tiempos vs población
import matplotlib.pyplot as plt

plt.plot(ti,xi, label='xi poblacion')
plt.plot(ti,yi, label='yi protestante')

plt.title('población y protestantes')
plt.xlabel('t años')
plt.ylabel('población')
plt.legend()
plt.grid()
plt.show()

# gráfica xi vs yi
plt.plot(xi,yi)

plt.title('población y protestantes [xi,yi]')
plt.xlabel('x población')
plt.ylabel('y protestantes')
plt.grid()
plt.show()

2023-01-242023-07-04

2Eva_2022PAOII_T1 Altura de cohete en 30 segundos

2da Evaluación 2022-2023 PAO II. 24/Enero/2023

Tema 1. (30 puntos) La velocidad hacia arriba de un cohete se calcula con la fórmula:

v = u \ln\Big(\frac{m_0}{m_0-qt}\Big) - gt

Donde:
v   = velocidad hacia arriba,
u   = 1800 m/s, velocidad a que se expele el combustible en relación con el cohete,
m₀ = 160 000 kg, masa inicial del cohete en el tiempo t = 0,
q    = 2 500 kg/s, tasa de consumo de combustible y
g    = 9.8 m/s2, aceleración de la gravedad

Para determinar la altura alcanzada por el cohete en un vuelo de 30 segundos desarrolle la parte analítica con los siguientes métodos y compare los resultados.

a. Utilice la regla de Simpson, en el planteamiento incluya la cantidad de tramos o segmentos a usar

b. Use el método de cuadratura de Gauss para la misma cantidad de segmentos que el literal anterior

c. Compare y comente los resultados, sobre los errores entre los métodos.

Rúbrica: Planteamiento de tramos (5 puntos), integral con Simpson (10 puntos), cuadratura de Gauss (10 puntos), literal c (5 puntos).

Referencia: Chapra ejercicio 24.46 p701. NASA y SpaceX realizan con éxito el despegue del primer vuelo de EE. UU. hacia la Estación Espacial Internacional en nueve años. EFE 30 mayo 2020 https://youtu.be/npcgpQUKAbg

2023-01-242023-01-27

s2Eva_2022PAOII_T1 Altura de cohete en 30 segundos

Ejercicio: 2Eva_2022PAOII_T1 Altura de cohete en 30 segundos

literal a

v = u \ln\Big(\frac{m_0}{m_0-qt}\Big) - gt

v = 1800 \ln\Big(\frac{160000}{160000-2500t}\Big) - 9.8t

Seleccionando el método de Simpson de 3/8, se requieren al menos 3 tramos o segmentos para usarlo, que generan 4 muestras. El vector de tiempo se obtiene como:

v = lambda t: 1800*np.log(160000/(160000-2500*t))-9.8*t
a = 0
b = 30
tramos = 3
h = (b-a)/tramos
ti = np.linspace(a,b,tramos+1)
vi = v(ti)

siendo los vectores:

ti = [ 0. 10. 20. 30.]
vi = [ 0. 207.81826623 478.44820899 844.54060574]

la aplicación del método de Simpson de 3/8 es:

I = \frac{3}{8}(10) \Bigg(1800 \ln\Big(\frac{160000}{160000-2500(0)}\Big) - 9.8(0)

+3(1800 \ln\Big(\frac{160000}{160000-2500(10)}\Big) - 9.8(10))

+3(1800 \ln\Big(\frac{160000}{160000-2500(20)}\Big) - 9.8(20))

+1800 \ln\Big(\frac{160000}{160000-2500(30)}\Big) - 9.8(30) \Bigg) =

I = \frac{3}{8}(10) \Big(v(0)+3(v(10))+3(v(20))+v(30) \Big)

I = \frac{3}{8}(10) \Big(0+3(207.81)+3(478.44)+844.54 \Big)

I = 10887.52

literal b

para el primer segmento se usa t entre [0,10]

x_a = \frac{0+10}{2} + \frac{1}{\sqrt{3}}\frac{10-0}{2} = 7.88

x_b = \frac{0+10}{2} - \frac{1}{\sqrt{3}}\frac{10-0}{2} = 2.11

I = \frac{10-0}{2}\Big(v(7.88)+v(2.11)\Big)=995.79

para el 2do segmento se usa t entre [10,20]

x_a = \frac{10+20}{2} + \frac{1}{\sqrt{3}}\frac{20-10}{2} = 17.88

x_b = \frac{10+20}{2} - \frac{1}{\sqrt{3}}\frac{20-10}{2} = 12.11

I = \frac{20-10}{2}\Big(v(17.88)+v(12.11)\Big) =3368.42

para el 3er segmento se usa t entre [20,30]

x_a = \frac{20+30}{2} + \frac{1}{\sqrt{3}}\frac{30-20}{2} = 27.88

x_b = \frac{20+30}{2} - \frac{1}{\sqrt{3}}\frac{30-20}{2} = 22.11

I = \frac{30-20}{2}\Big(v(27.88)+v(22.11)\Big) = 6515.23

Altura = 995.79+ 3368.42 + 6515.23 = 10879.44

literal c

el error es la diferencia entre los métodos
error_entre = |10887.52-10879.44| = 8.079

Resultados con algoritmo

Método de Simpon 3/8
ti
[ 0. 10. 20. 30.]
vi
[ 0. 207.81826623 478.44820899 844.54060574]
Altura con Simpson 3/8 : 10887.52511781406
segmento Cuad_Gauss :    [995.792, 3368.421, 6515.231]
Altura Cuadratura Gauss: 10879.445437288954
diferencia s3/8 y Cuad_Gauss: 8.079680525106596
>>>

Instrucciones en Python

# 2Eva_2022PAOII_T1 Altura de cohete en 30 segundos
import numpy as np

# INGRESO
v = lambda t: 1800*np.log(160000/(160000-2500*t))-9.8*t
a = 0
b = 30
tramos = 3

# PROCEDIMIENTO literal a
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)

h = (b-a)/tramos
ti = np.linspace(a,b,tramos+1)
vi = v(ti)
altura = integrasimpson38_fi(ti,vi)

# SALIDA
print('Método de Simpon 3/8')
print('ti')
print(ti)
print('vi')
print(vi)
print('Altura con Simpson 3/8 :',altura)

# PROCEDIMIENTO literal b
# cuadratura de Gauss de dos puntos
def integraCuadGauss2p(funcionx,a,b):
    x0 = -1/np.sqrt(3)
    x1 = -x0
    xa = (b+a)/2 + (b-a)/2*(x0)
    xb = (b+a)/2 + (b-a)/2*(x1)
    area = ((b-a)/2)*(funcionx(xa) + funcionx(xb))
    return(area)

area = 0
area_i =[]
for i in range(0,tramos,1):
    deltaA = integraCuadGauss2p(v,ti[i],ti[i+1])
    area = area + deltaA
    area_i.append(deltaA)
# SALIDA
print('segmento Cuad_Gauss :   ', area_i)
print('Altura Cuadratura Gauss:', area)

print('diferencia s3/8 y Cuad_Gauss:',altura-area)

import matplotlib.pyplot as plt
plt.plot(ti,vi)
plt.plot(ti,vi,'o')
plt.title('v(t)')
plt.xlabel('t (s)')
plt.ylabel('v (m/s)')
plt.grid()
plt.show()