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)
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.
total de fuente de generación para la tabla es: [1500,400,600,200]
literal a Planteamiento
Considerando que el factor de la tabla es el factor de consumo por tipo de cliente, para encontrar las cantidades por tipo cliente que se puedan atender con la capacidad de cada «fuente de generación», se plantea:
Al observar la tercera columna j=2, desde la diagonal i=2, el valor mayor ya se encuentra en la posición de la diagonal, no hay cambio de filas. Para la última fila no se tiene otra fila para comparar, con lo que el pivoteo se terminó.
literal c. Método de Jacobi
A partir de la matriz del literal anterior, las expresiones quedan:
literal e. convergencia y revisión de resultados obtenidos
Para el asunto de convergencia, el error disminuye entre iteraciones, por lo que el método es convergente.
Analizando los resultados desde la segunda iteración, se obtienen valores negativos, lo que no es acorde al concepto del problema. Se continuarán con las iteraciones siguientes y se observará el resultado para dar mayor «opinión» sobre lo obtenido.
literal f Número de condición
Número de condición = 5.77 que al ser «bajo y cercano a 1» se considera que el método converge.
El resultado muestra que los clientes residenciales serán -257 según las ecuaciones planteadas. Lo que se interpreta según lo presentado por los estudiantes:
1. energía insuficiente: causa por los apagones diarios en el país en éstos días.
2. Se requiere mejorar el planteamiento, pues la respuesta no debe contener valores negativos según el concepto planteado en el ejercicio.
Se considerarán ambas válidas para la evaluación al observar que el resultado numérico es correcto, pero no es acorde al ejercicio.
Instrucciones en Python
# 1Eva_2024PAOII_T3 matriz de distribución de energía por sectoresimport numpy as np
#INGRESO
A= [[0.7,0.19,0.1,0.01],
[0.12,0.18,0.68,0.02],
[0.2,0.38,0.4,0.02],
[0.11,0.23,0.15,0.51]]
B = [1500,400,600,200]
# PROCEDIMIENTO
A = np.array(A,dtype=float)
B = np.transpose([np.array(B)])
X0 = [100,100,100,100]
tolera = 0.1
iteramax = 100
verdecimal = 5
defjacobi(A,B,X0, tolera, iteramax=100, vertabla=False, precision=4):
''' Método de Jacobi, tolerancia, vector inicial X0
para mostrar iteraciones: vertabla=True
'''
A = np.array(A,dtype=float)
B = np.array(B,dtype=float)
X0 = np.array(X0,dtype=float)
tamano = np.shape(A)
n = tamano[0]
m = tamano[1]
diferencia = np.ones(n, dtype=float)
errado = np.max(diferencia)
X = np.copy(X0)
xnuevo = np.copy(X0)
tabla = [np.copy(X0)]
itera = 0
if vertabla==True:
print('Iteraciones Jacobi')
print('itera,[X],errado')
print(itera, xnuevo, errado)
np.set_printoptions(precision)
whilenot(errado<=tolera or itera>iteramax):
for i inrange(0,n,1):
nuevo = B[i]
for j inrange(0,m,1):
if (i!=j): # excepto diagonal de A
nuevo = nuevo-A[i,j]*X[j]
nuevo = nuevo/A[i,i]
xnuevo[i] = nuevo
diferencia = np.abs(xnuevo-X)
errado = np.max(diferencia)
X = np.copy(xnuevo)
tabla = np.concatenate((tabla,[X]),axis = 0)
itera = itera + 1
if vertabla==True:
print(itera, xnuevo, errado)
# No convergeif (itera>iteramax):
X=itera
print('iteramax superado, No converge')
return(X,tabla)
defpivoteafila(A,B,vertabla=False):
'''
Pivotea parcial por filas, entrega matriz aumentada AB
Si hay ceros en diagonal es matriz singular,
Tarea: Revisar si diagonal tiene ceros
'''
A = np.array(A,dtype=float)
B = np.array(B,dtype=float)
# Matriz aumentada
nB = len(np.shape(B))
if nB == 1:
B = np.transpose([B])
AB = np.concatenate((A,B),axis=1)
if vertabla==True:
print('Matriz aumentada')
print(AB)
print('Pivoteo parcial:')
# Pivoteo por filas AB
tamano = np.shape(AB)
n = tamano[0]
m = tamano[1]
# Para cada fila en AB
pivoteado = 0
for i inrange(0,n-1,1):
# columna desde diagonal i en adelante
columna = np.abs(AB[i:,i])
dondemax = np.argmax(columna)
# dondemax no es en diagonalif (dondemax != 0):
# intercambia filas
temporal = np.copy(AB[i,:])
AB[i,:] = AB[dondemax+i,:]
AB[dondemax+i,:] = temporal
pivoteado = pivoteado + 1
if vertabla==True:
print(' ',pivoteado, 'intercambiar filas: ',i,'y', dondemax+i)
if vertabla==True:
if pivoteado==0:
print(' Pivoteo por filas NO requerido')
else:
print('AB')
return(AB)
# PROGRAMA Búsqueda de solucion --------# PROCEDIMIENTO
AB = pivoteafila(A,B,vertabla=True)
# separa matriz aumentada en A y B
[n,m] = np.shape(AB)
A = AB[:,:m-1]
B = AB[:,m-1]
[X, puntos] = jacobi(A,B,X0,tolera,vertabla=True,precision=verdecimal)
iterado = len(puntos)
# numero de condicion
ncond = np.linalg.cond(A)
# SALIDAprint('numero de condición:', ncond)
print('respuesta X: ')
print(X)
print('iterado, incluyendo X0: ', iterado)
De la tabla del ejercicio, solo se toman las dos primeras filas ti,xi
ti
1
1.2
1.4
1.8
2
xi
-80.0108
-45.9965
3.1946
69.5413
61.1849
yi
-8.3002
-22.6765
-20.9677
15.8771
33.8999
zi
113.8356
112.2475
110.5523
106.7938
104.71
literal a. Planteamiento
En el planteamiento del problema para un polinomio de grado 3 indicado en el enunciado, se requieren al menos 4 puntos (grado=puntos-1). Al requerir cubrir todo el intervalo, los puntos en los extremos son obligatorios, quedando solo dos puntos por seleccionar, que se sugiere usar alrededor de 1.63 que es el valor buscado. quedando de ésta forma la selección de los puntos:
xi = [1. , 1.4, 1.8, 2. ]
fi = [-80.0108, 3.1946, 69.5413, 61.1849]
En el parcial, se revisaron dos métodos: polinomio de interpolación con sistema de ecuaciones en la Unidad 3, y el conceptual con diferencias divididas de Newton. Por la extensión del desarrollo se realiza con diferencias divididas de Newton, que por facilidad de espacios se muestra en dos tablas, la de operaciones y resultados.
Tabla de operaciones:
i
ti
f[ti]
Primero
Segundo
Tercero
0
1
-80.0108
1.4−13.1946−(−80.0108)
1.8−1165.8667−208.0135
2−1−346.0812−52.6834
1
1.4
3.1946
1.8−1.469.5413−3.1946
2−1.4−41.782−165.8667
2
1.8
69.5413
2−1.861.1849−69.5413
3
2
61.1849
Tabla de resultados
i
ti
f[ti]
Primero
Segundo
Tercero
0
1
-80.0108
208.0135
-52.6834
-293.3978
1
1.4
3.1946
165.8667
-346.0812
2
1.8
69.5413
-41.782
3
2
61.1849
Se construye el polinomio con los datos de la primera fila:
Para encontrar el valor t cuando la cofia alcanza 30 mts o 0.030 Km se usa la fórmula del eje z(t)=0.030 dado que se indica que las fórmulas se encuentran en Km.
también puede verificarse usando las gráficas de las trayectorias de cada eje vs el tiempo en el intervalo. Para el eje z, se puede marcar el punto de interés al trazar una línea horizontal en cero o en la altura z(t) = 0.030.
literal c. Iteraciones
Se indica usar el método del punto fijo, donde es necesario definir g(t) como un lado de la balanza, mientra que del otro lado es la identidad con t. de la fórmula se despeja por ejemplo del término t2
el punto inicial t0 se puede escoger por ejemplo la mitad del intervalo. Según se observa la gráfica.
La tolerancia puede ser de 10 cm, o 1 m, depende de la precisión a proponer.
Para la parte escrita se selecciona 1m, como las unidades se encuentran en Km: tolera = 0.001
la gráfica del intervalo muestra que es necesario ampliar(zoom) el eje f(x) para observar mejor.
Instrucciones en Python
# Algoritmo de punto fijo# [a,b] son seleccionados desde la gráfica# error = toleraimport numpy as np
# INGRESO
g = 35.28
alt = 0.030
fx = lambda t: -alt + 120 + 450/60*t*(1-np.exp(-0.35*t))-0.5*g*t**2 + 0.012*(7.5 - g*t)**2
gx = lambda t: np.sqrt((2/g)*(-alt + 120 + 450/60*t*(1-np.exp(-0.35*t)) + 0.012*(7.5 - g*t)**2))
a = 5
b = 10
tolera = 0.001 # tolera 1m, probar con 10 cm
iteramax = 100
# PROCEDIMIENTO
i = 0 # iteración
b = gx(a)
tramo = abs(b-a)
print('i,[a,b,tramo]')
np.set_printoptions(precision=4)
print(0,np.array([a,b,tramo]))
whilenot(tramo<=tolera or i>iteramax):
a = b
b = gx(a)
tramo = abs(b-a)
i = i + 1
print(i,np.array([a,b,tramo]))
respuesta = b
# Valida convergenciaif (i>=iteramax):
respuesta = np.nan
# SALIDAprint('raíz en: ', respuesta)
print('errado : ', tramo)
print('itera : ', i)
# GRAFICAimport matplotlib.pyplot as plt
# calcula los puntos para fx y gx
a = 5
b = 10
muestras = 21
xi = np.linspace(a,b,muestras)
fi = fx(xi)
gi = gx(xi)
yi = xi
plt.plot(xi,fi, label='f(t)',
linestyle='dashed')
plt.plot(xi,gi, label='g(t)')
plt.plot(xi,yi, label='y=t')
plt.axvline(respuesta, color='magenta',
linestyle='dotted')
plt.axhline(0)
plt.xlabel('t')
plt.ylabel('f(t)')
plt.title('Punto Fijo')
plt.grid()
plt.legend()
plt.show()
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.
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.
