Los resultados para x(t) de las primeras iteraciones indican que la población total del país disminuye en el intervalo observado. Se puede comprobar al usar el algoritmo para mas iteraciones como se muestra en la gráfica.
literal d
Los resultados para y(t) muestran que la población clasificada como Rojo aumenta en el intervalo observado. Sin embargo la mayoría sigue siendo Azul para las tres iteraciones realizadas.
literal e
La población y(t) alcanza la mitad de la población cuanto t cambia de 30.5 a 31. Tiempo en el que de realizar elecciones, ganarían los Rojos.
Usando el algoritmo, se añade una columna con la condición que MayoriaY = yi>xi/2, que se convierte a 1 o verdadero cuando se cumple la condición. Con lo que se encuentra el cambio de mayoría a Rojo.
Calcular los tamaños de paso dxi en cada frontera y plantear la integración con fórmulas compuestas.
Usando los datos de las coordenadas de obtiene cada dxi = xi[i+1]-xi[i]. De forma semejante se encuentra cada dxj, Seleccionando los métodos según se disponga de tamaño de paso iguales y consecutivos como se muestra en la tabla ampliada.
Frontera superior
Trapecio
Trapecio
Trapecio
Trapecio
Simpson 1/3
Trapecio
Trapecio
Simpson 1/3
Trapecio
Simpson 1/3
dxi
40
100
-30
66
20
20
79
125
54
50
50
20
20
—
xi
410
450
550
520
586
606
626
705
830
884
934
984
1004
1024
yi
131
194
266
337
402
483
531
535
504
466
408
368
324
288
Frontera Inferior
Trapecio
Simpson 3/8
dxj
190
190
190
44
—
xj
410
600
790
980
1024
yj
131
124
143
231
288
Desarrollando con instrucciones sobre el arreglo en Python con la instrucción np.diff(xi).
Desarrollar las expresiones del área para las coordenadas de la frontera superior, según el literal a. Cuando se tienen dos tamaños de paso iguales se usa Simpson de 1/3.
Como referencia par el ejercicio, solo se realizará la primera parte, acorde a:
Encuentre el tiempo tc y la velocidad de la persona cuando se alcanza la longitud de la cuerda extendida y sin estirar (30 m), es decir y<L, aún se entra cayendo signo(v)=1. (solo primera ecuación)
dt2d2y=g−signo(v)mCd(dtdy)2y≤L
Suponga que las condiciones iniciales son:
y(0) = 0
dtdy(0)=0
literal a
Realice el planteamiento del ejercicio usando Runge-Kutta de 2do Orden
dt2d2y=9.8−signo(v)68.10.25(dtdy)2
Usando los valores de las constantes g, Cd, m, haciendo el cambio de variable dy/dt=v. Adicionalmente, en la caída inicial, signo(v)=1 y se mantiene constante hasta el 2do tramo con la cuerda estirada.
v′=9.8−68.10.25v2
Se plantea el ejercicio como Runge-Kutta para 2da derivada
con lo que se observa que para alcanzar los 30m de la cuerda sin estirar se alcanzan entre los [2.5, 3.0] s.
literal d
Indique el valor de tc, muestre cómo mejorar la precisión y realice sus observaciones sobre los resultados.
Para mejorar la precisión del algoritmo se reduce el valor del tamaño de paso h, de esta forma se puede obtener una mejor lectura del tiempo de la primera fase del ejercicio.
Por ejemplo haciendo el tamaño de paso h=0.5/4, el tiempo entre [2.5, 2.625]. que tiene un error segundos menor al valor encontrado inicialmente y se mejora la precisión.
Para el volumen con f(x) con al menos 3 tramos y un método de Simpson, directamente se puede usar 3/8. Por lo que se Se reemplaza en la fórmula de volumen del sólido de revolución f(x) con:
f(x)=sin(x/2)
obteniendo:
Vfx=∫abπ(sin(x/2))2dx=∫abπsin(x/2)dx
La expresión dentro del integral se denomina como fv:
fv(x)=πsin(x/2)
en el intervalo [0.1, 1.8], con al menos 3 tramos, se requieren 4 muestras con tamaño de paso hf: y truncando a 4 decimales los resultados calculados con Python.
hf=tramosb−a=31.8−0.1=0.5666
los puntos de muestra quedan np.linspace(0.1,1.8,3+1):
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 amortiguadoimport numpy as np
import matplotlib.pyplot as plt
defrungekutta2_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 inrange(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()
para la evaluación numérica de 1 se usa un valor muy cercano desplazado con la tolerancia aplicada.
Igx≅83(0.3333)[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.
# 2Eva_2023PAOI_T1 Material para medalla de academiaimport numpy as np
import matplotlib.pyplot as plt
defintegrasimpson13_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
whilenot(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)
defintegrasimpson38_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
whilenot(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
# SALIDAprint('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()
# EDP parabólicas d2u/dx2 = K du/dt# método explícito,usando diferencias divididasimport 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
whilenot(j>=ultimot): # igual con lazo forfor i inrange(1,ultimox,1):
u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]
j=j+1
# SALIDAprint('Tabla de resultados')
np.set_printoptions(precision=2)
print(u)
# Gráfica
salto = int(n/10)
if (salto == 0):
salto = 1
for j inrange(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()
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.
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 Ordenimport numpy as np
defrungekutta2_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 inrange(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ónimport 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()
