Para el planteamiento del integral es necesario observar la gráfica de la función a integrar dentro del intervalo.
La función tiene al menos dos «picos» y dos valles en el intervalo. Por lo que un solo tramo del integral podría aumentar el error de integración numérica con una figura trapezoidal equivalente como propone la cuadratura de Gauss.
Se plantea usar al menos dos tramos, y comparar el resultado con tres tramos para observar el error.
Para dos tramos se dispone de los segmentos entre los puntos
[0, 3.5, 7]
Para tres tramos se tiene los segmentos entre los puntos
[ 0, 7/3, 2(7/3), 7]
literal b
Si se usan dos tramos se tienen los segmentos entre los puntos [0,3.5,7]
Para el desarrollo de las iteraciones se requieren valores iniciales. Para la variable independiente tiempo podría usar t = 0.
Para la variable dependiente población, según la descripción se encontraría entre el intervalo [10,50]. Si x0<10 la población se extingue. Si la variable x0>50 se encuentra saturada la capacidad del medio.
Por lo que se propone usar un valor mayor que el mínimo, por ejemplo x0=11 y otro valor que se encuentre en el intervalo.
Según los resultados de las tres iteraciones anteriores, la población crece. El resultado es acorde al concepto descrito en el enunciado para poblaciones mayores al valor mínimo de extinción. En el ejercicio se usó 11 como valor inicial.
Esto se puede comprobar usando el algoritmo y teniendo los resultados presentados en el literal d
xi = [0. , 2.50, 3.75, 5.00, 6.25, 7.5 ]
yi = [3.7 , 3.25, 4.05, 4.33, 2.95, 3.22]
literal a
Según la gráfica presentada, los puntos a considerar para todo el intervalo, deberían ser al menos el primero y el último. Para un polinomio de grado 3 se requieren usar 4 muestras, por lo que faltarían dos muestras dentro del intervalo. Los puntos adicionales estarían entre los puntos intermedios, tratando de mantener una distancia equidistante entre ellos y tratar de mantener los cambios de dirección.
Los puntos seleccionados para el ejercicio serán
xi = [0. , 2.50, 5.00, 7.5 ]
fi = [3.7 , 3.25, 4.33, 3.22]
Se puede construir el polinomio con los cualquiera de los métodos para interpolación dado que tienen tamaños de paso iguales entre tramos.
Desarrollando con el método de Lagrange, el polinomio se construye como:
Se comprueba que los valores obtenidos corresponden a las muestras, por lo que el polinomio cumple con los criterios básicos de interpolación. La gráfica permite verificar también el resultado.
Observando las gráficas de la trayectoria construida junto a las funciones descritas en el tema 1 se tiene que:
Un polinomio de grado 3 es insuficiente para describir la trayectoria, se debe aumentar el grado del polinomio para ajustar mejor la curva.
Por ejemplo, usando todos los puntos, la trayectoria y el polinomio son mas cercanas aunque no iguales.
literal e
Encontrar el error en P(5), como x=5 y es parte de los puntos de muestra, el error debería ser cero. Siempre y cuando x=5 sea parte de los puntos seleccionados.
p3(5)=3.7−0.982(5)+0.42(5)2−0.03968(5)3=4.33
error = | 4.33-4.33| = 0
literal f
Los resultados con el algoritmo de Lagrange se muestran como:
La distancia entre las aeronaves se determina a partir de las ecuaciones proporcionadas en el enunciado.
d=(x2−x1)2+(y2−y1)2+(z2−z1)2
La ecuación también mantiene la forma y el mínimo si se utiliza el cuadrado de la distancia:
D=d2=(x2−x1)2+(y2−y1)2+(z2−z1)2
Las diferencias de distancia por eje pueden ser convergentes hacia el punto de impacto, dado que se conoce que el accidente ya ocurrió. Por lo que también se pueden revisar las distancias entre ejes en lugar de la distancia total en 3D, simplificando un poco el ejercicio. Observe la gráfica proporcionada:
Se realiza la gráfica para la distancia al cuadrado Di y las distancias por cada eje dxi, dyi, dzi :
Por lo que el resultado también se podría determinar usando por ejemplo el eje y. La ecuación a buscar la distancia mínima, punto de choque o cruce por cero podría ser:
f(t)=dy(t)=sin(0.1t)cos(0.7t)+3.7−0.4t=0
literal b
Por la gráfica se puede obtener un intervalo de búsqueda entre [8, 12], que usando solo las distancias en el eje y, se tiene:
dy(8)=1.0563dy(12)=−1.5839
literal c
Se pide usar un método de búsqueda de raíces, pudiendo seleccionar Bisección en el intervalo encontrado en el literal b.
La variable independiente en el ejercicio es tiempo ‘t’, que podría ser segundos. Si la tolerancia se estima en milisegundos, tolera = 10-3 .
El error se ha calculado en cada iteración
literal e
El error disminuye en cada iteración, por lo que el método converge.
La distancia mínima entre las aeronaves no tiene que llegar a cero para que se produzca un accidente. Las dimensiones de las aeronaves muestran que se encuentran entre 70m y 4 m, por lo que se estima que debe existir una separación mayor a 70 metros en cualquiera de los ejes para evitar un accidente. Si las coordenadas se estiman en Km, la tolerancia sería de 0.070 Km al considerar más de 70 metros como la distancia segura.
Instrucciones en Python para Gráfica de distancias
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
# INGRESO# Avión Aterriza
Ax = lambda t: 5.1 + 0*t
Ay = lambda t: 0.4*t
Az = lambda t: 0.5*np.exp(-0.2*t) + 0.3
# helicóptero
Hx = lambda t: 0.5*t + 0*t
Hy = lambda t: np.sin(0.10*t)*np.cos(0.7*t)+ 3.7
Hz = lambda t: 0.36 + 0*t
a = 0
b = 6/0.5 # x entre[0,6] usando gx(t)= 6
muestras = 41
# PROCEDIMIENTO# Distancia por ejes
dx = lambda t: Hx(t) - Ax(t)
dy = lambda t: Hy(t) - Ay(t)
dz = lambda t: Hz(t) - Az(t)
# Distancia 3D
distancia2 = lambda t: dx(t)**2 + dy(t)**2 + dz(t)**2
# Simulacion en segmento t=[a,b]
ti = np.linspace(a,b,muestras)
dxi = dx(ti)
dyi = dy(ti)
dzi = dz(ti)
Di = distancia2(ti)
# SALIDAprint(dy(8))
print(dy(12))
plt.plot(ti,Di, label='Di')
plt.plot(ti,dxi,label='dxi')
plt.plot(ti,dyi,label='dyi')
plt.plot(ti,dzi,label='dzi')
plt.xlabel('ti')
plt.ylabel('distancia2')
plt.grid()
plt.legend()
plt.show()
Instrucciones en Python – Algoritmo Bisección
# 3Eva_2024PAOII_T1 Accidente entre aeronaves# Algoritmo de Bisección# [a,b] se escogen de la gráfica de la función# error = toleraimport numpy as np
import matplotlib.pyplot as plt
# Algoritmo de Bisección# [a,b] se escogen de la gráfica de la función# error = toleraimport numpy as np
defbiseccion(fx,a,b,tolera,iteramax = 20, vertabla=False, precision=4):
'''
Algoritmo de Bisección
Los valores de [a,b] son seleccionados
desde la gráfica de la función
error = tolera
'''
fa = fx(a)
fb = fx(b)
tramo = np.abs(b-a)
itera = 0
cambia = np.sign(fa)*np.sign(fb)
if cambia<0: # existe cambio de signo f(a) vs f(b)if vertabla==True:
print('método de Bisección')
print('i', ['a','c','b'],[ 'f(a)', 'f(c)','f(b)'])
print(' ','tramo')
np.set_printoptions(precision)
while (tramo>=tolera and itera<=iteramax):
c = (a+b)/2
fc = fx(c)
cambia = np.sign(fa)*np.sign(fc)
if vertabla==True:
print(itera,[a,c,b],np.array([fa,fc,fb]))
if (cambia<0):
b = c
fb = fc
else:
a = c
fa = fc
tramo = np.abs(b-a)
if vertabla==True:
print(' ',tramo)
itera = itera + 1
respuesta = c
# Valida respuestaif (itera>=iteramax):
respuesta = np.nan
else:
print(' No existe cambio de signo entre f(a) y f(b)')
print(' f(a) =',fa,', f(b) =',fb)
respuesta=np.nan
return(respuesta)
# INGRESO
fx = lambda t: np.sin(0.10*t)*np.cos(0.7*t)+ 3.7 - 0.4*t
a = 8
b = 12
tolera = 0.001
# PROCEDIMIENTO
respuesta = biseccion(fx,a,b,tolera,vertabla=True)
# SALIDAprint('raíz en: ', respuesta)
El volumen total es la suma de los resultados de cada una de las secciones:
Vtotal=2120.5750+1255.1971+5539.9805+...
Tarea: continuar con los otros intervalos para obtener el volumen total.
literal e
Resultados con el algoritmo
para f(x):
[xa,xb]= [0.6339745962155612, 2.366025403784439]
[f(xa),f(xb)]= [915.4418590078262, 760.106947830205]
Volumenfx: 2513.323210257047
para f(x):
[xa,xb]= [3.4226497308103743, 4.577350269189626]
[f(xa),f(xb)]= [672.4334354361756, 582.7637293668464]
Volumenfx: 1255.1971648030221
para f(x):
[xa,xb]= [7.106908908089714, 12.863091091910285]
[f(xa),f(xb)]= [641.5176069162055, 469.8124936184865]
Volumenfx: 5539.98055116544
Instrucciones en Python
# 3Eva_2023PAOII_T3 Volumen por solido de revolución de un peónimport numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
# INGRESO# intervalos eje x, funciones#a=0 ; b=3 ; f = lambda x: np.pi*(15)**2 + 0*x#a=3 ; b=5 ; f = lambda x: np.pi*(-0.875*x+17.625)**2# con el polinomio del tema 2
a=5 ; b=14.97 ; f = lambda x: np.pi*(-0.1083*x**2+1.8047*x + 6.9341)**2
# con el circulo del tema 1#a=5 ; b=14.97 ; f = lambda x: np.pi*(np.sqrt(6.7**2-(x-8.6)**2) + 7.6)**2
xmuestras = 31
# PROCEDIMIENTO# Volumen para f(x) con Cuadratura de Gauss
xa = (b+a)/2 + (b-a)/2*(-1/np.sqrt(3))
xb = (b+a)/2 + (b-a)/2*(1/np.sqrt(3))
Volumenfx = (b-a)/2*(f(xa)+f(xb))
xab = [xa,xb]
fab = [f(xa),f(xb)]
# SALIDAprint("para f(x):")
print("[xa,xb]=", xab)
print("[f(xa),f(xb)]=", fab)
print("Volumenfx: ",Volumenfx)
print()
Añadir para graficar en 2D
# grafica 2D para el area de corte en eje y,x --------# muestras en x
xi = np.linspace(a, b, xmuestras)
fi = f(xi)
f0 = np.zeros(xmuestras,dtype=float)
# grafica 2D
figura2D = plt.figure()
grafica2D = figura2D.add_subplot(111)
grafica2D.plot(xi,fi,color='blue',label='f(x)')
grafica2D.fill_between(xi,fi,f0,color='lightblue')
grafica2D.grid()
grafica2D.set_title('Area para sólido de revolución')
grafica2D.set_xlabel('x')
grafica2D.set_ylabel('f(x)')
Añadir para graficar en 3D
# Para grafica 3D -----------------------# angulo w de rotación
w_a = 0
w_b = 2*np.pi
w_muestras = 31
# grafica 3D muestras en x y angulo w
wi = np.linspace(w_a, w_b, w_muestras)
X, W = np.meshgrid(xi, wi)
# proyeccion en cada eje
Yf = f(xi)*np.cos(W)
Zf = f(xi)*np.sin(W)
# grafica 3D
figura3D = plt.figure()
grafica = figura3D.add_subplot(111, projection='3d')
grafica.plot_surface(X, Yf, Zf,
color='blue', label='f(x)',
alpha=0.6, rstride=6, cstride=12)
grafica.set_title('Sólido de revolución')
grafica.set_xlabel('x')
grafica.set_ylabel('y')
grafica.set_zlabel('z')
# grafica.legend()
eleva = 30
rota = -45
deltaw = 5
grafica.view_init(eleva, rota)
# rotacion de ejesfor angulo inrange(rota, 360+rota, deltaw ):
grafica.view_init(eleva, angulo)
plt.draw()
plt.pause(.001)
plt.show()
Para el con los puntos xi = [5, 9.985, 14.97] el perfil se describe con tres puntos, por lo que el polinomio de interpolación puede ser de grado 2. Usando el método de Lagrange:
El error para los polinomios acorde a la gráfica proporcionada es cero para los tramos [0,3] y [3,5], dado que el primero es una constante y el segundo es una recta con pendiente.
El error de mayor magnitud se presenta con la interpolación entre el círculo de la ecuación del tema 1 y el polinomio de grado 2 en el intervalo [5, 14.97]. Tomando como ejemplo el punto intermedio del tramo derecho en:
Se puede mejorar la aproximación del polinomio aumentando el grado y número de puntos a usar dentro del intervalo. El resultado se debería acercar mas al perfil superior del círculo de referencia del tema 1.
# 3Eva_2023PAOII_T1 Intersección con círculoimport sympy as sym
import numpy as np
import matplotlib.pyplot as plt
# INGRESO# forma algebraica con sympy
x = sym.Symbol('x')
f = -(1.75/2)*x + (35.25/2)
g = sym.sqrt(6.7**2-(x-8.6)**2) + 7.6
distancia = f - g
x0 = 3
tolera = 0.0001 # 1e-4# PROCEDIMIENTO# literal c
dfg = sym.diff(distancia,x)
# convierte a forma numerica con numpy# Newton-Raphson
fx = sym.lambdify(x,distancia)
dfx = sym.lambdify(x,dfg)
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)
# SALIDAprint('dfg(x) = ')
sym.pprint(dfg)
print()
print(['xi', 'xnuevo', 'tramo'])
np.set_printoptions(precision = 4)
print(tabla)
print('raiz en: ', xi)
print('con error de: ',tramo)
# Respuesta a entrada cero# solucion para (D^2+ D + 1)y = 0import 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,3),dtype=float)
# incluye el punto [x0,y0]
estimado[0] = [x0,y0,z0]
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]
return(estimado)
# PROGRAMA
f = lambda t,y,z: z
g = lambda t,y,z: -3*z -2*y
t0 = 0
y0 = 0
z0 = -5
h = 0.1
tn = 6
muestras = int((tn-t0)/h)
tabla = rungekutta2_fg(f,g,t0,y0,z0,h,muestras)
ti = tabla[:,0]
yi = tabla[:,1]
zi = tabla[:,2]
# SALIDA
np.set_printoptions(precision=6)
print('t, y, z')
print(tabla)
# GRAFICA
plt.plot(ti,yi, label='y(t)')
plt.ylabel('y(t)')
plt.xlabel('t')
plt.title('Runge-Kutta 2do Orden d2y/dx2 ')
plt.legend()
plt.grid()
plt.show()
El ejercicio usa el resultado del tema anterior, planteando una función de Python como la solución para valores dados. Se requiere una función, para disponer de los valores solución en cada llamada para el intervalo de análisis.
Por lo que básicamente lo que se pide es usar algún algoritmo de búsqueda de raíces. Para simplicidad en la explicación se toma el algoritmo de la bisección.
Los resultados se grafican como theta vs Tensiones, y el punto a buscar es cuando las tensiones en los cables son de igual magnitud, es decir:
La función tiene al menos dos «picos» y dos valles en el intervalo. Por lo que un solo tramo del integral podría aumentar el error de integración numérica con una figura trapezoidal equivalente como propone la cuadratura de Gauss.
