Para el problema, se usan los puntos: [a,f(a)], [b,f(b)] y [c,f(c)]
por donde pasa la curva f(x) aproximada a un polinomio de grado 2, f(x) \approx p(x)
# 2da Evaluación I Término 2018# Tema 4. canal acceso a Puertos de Guayaquilimport numpy as np
# INGRESO
profcanal = 11
xi = np.array([ 0., 50., 100., 150., 200., 250.])
yi = np.array([ 0., 100., 200., 300., 400., 500.])
batimetria = [[ -6.79,-12.03,-10.04,-11.60, -7.24,-7.91],
[ -8.85,-10.89, -8.95, -7.23,-11.42,-7.93],
[-11.90, -9.86, -9.35,-12.05, -9.38,-9.65],
[ -7.30,-11.55,-10.41, -8.67,-11.84,-6.77],
[-12.17, -9.62, -7.47, -6.51, -9.02,-9.60],
[-11.90,-10.23,-10.68, -9.94, -6.76,-7.46]]
batimetria = np.array(batimetria)
# PROCEDIMIENTO
[n,m] = np.shape(batimetria)
# Matriz remover sedimentos
remover = batimetria + profcanal
for i inrange(0,n,1):
for j inrange(0,m,1):
if remover[i,j]<0:
remover[i,j]=0
# SALIDAprint('matriz remover sedimentos: ')
print(remover)
b) el volumen se calcula usando el algoritmo de Simpson primero por un eje, y luego con el resultado se continúa con el otro eje,
Considere que existen 6 puntos, o 5 tramos integrar en cada eje.
Al usar Simpson de 1/3 que usan tramos pares, faltaría integrar el último tramo.
En el caso de Simpson de 3/8 se requieren tramos multiplos de 3, porl que faltaría un tramo para volver a usar la fórmula.
La solución por filas se implementa usando una combinación de Simpson 3/8 para los puntos entre remover[i, 0:3] y Simpson 1/3 para los puntos entre remover[i, 3:5].
Luego se completa el integral del otro eje con el resultado anterior, aplicando el mismo método.
Se obtuvieron los siguientes resultados:
Integral en eje x:
[ 219.1 309.83 260.44 217.75 511.21 137.85]
Volumen: 160552.083333
que se obtiene usando las instrucciones a continuación de las anteriores:
a. Planteamiento con Runge-Kutta 2do Orden para Segunda derivada
La expresión:
\frac{dv}{dt} = g - \frac{cd}{m} v^2
se puede plantear sustituir la variable con v = -\frac{dy}{dt} al considerar el sentido contrario entre la velocidad al caer y la referencia de altura hacia arriba. Ver figura.
\frac{dy^2}{dt^2} = g - \frac{cd}{m} \Big( \frac{dy}{dt} \Big) ^2
con las condiciones iniciales del ejercicio t0 = 0 , y0 = 1000, v0 = 0
la velocidad se inicia con cero, si el paracaidista se deja caer desde el borde el risco, como en el video adjunto al enunciado.
Para las iteraciones, recuerde que
t se corrige con t+h (en el algoritmo era la posición para x)
y se corrige con y+K1y
v se corrige con v+K1v (en el algoritmo era la posición para z)
itera = 0
K1y = h f(t,y,v) = 2(-(0)) = 0 K1v = h g(t,y,v) = 2(9.8 - \frac{0.225}{90} (0)^2) = 19.6
.. K2y = h f(t+h, y+K1y, v + K1v) = 2(-(0 + 19.6)) = -39.2
K2v = h g(t+h, y+K1y, v + K1v) = 2(9.8 - \frac{0.225}{90} (0+19.6)^2) =17.6792
Resultados
La velocidad máxima, si no hubiese límite en la altura, se encuentra en alrededor de 62.39 m/s.
Sin embargo, luego de 20 segundos se observa que la altura alcanza el valor de cero, es decir se alcanzó tierra con velocidad de 62 m/s, que son algo mas de 223 Km/h, el objeto se estrella…!
# 2Eva_IT2018_T1 Paracaidista wingsuitimport numpy as np
import matplotlib.pyplot as plt
defrungekutta2_fg(f,g,x0,y0,z0,h,muestras,
vertabla=False, precision = 6):
''' solucion a EDO con Runge-Kutta 2do Orden Segunda derivada,
x0,y0 son valores iniciales, h es tamano de paso,
muestras es la cantidad de puntos a calcular.
'''
tamano = muestras + 1
tabla = np.zeros(shape=(tamano,3+4),dtype=float)
# incluye el punto [x0,y0,z0]
tabla[0] = [x0,y0,z0,0,0,0,0]
xi = x0
yi = y0
zi = z0
i=0
if vertabla==True:
print('Runge-Kutta Segunda derivada')
print('i ','[ xi, yi, zi',']')
print(' [ K1y, K1z, K2y, K2z ]')
np.set_printoptions(precision)
print(i,tabla[i,0:3])
print(' ',tabla[i,3:])
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
tabla[i] = [xi,yi,zi,K1y,K1z,K2y,K2z]
if vertabla==True:
txt = ' 'if i>=10:
txt=' 'print(str(i)+'',tabla[i,0:3])
print(txt,tabla[i,3:])
return(tabla)
# PROGRAMA PRUEBA# INGRESO
f = lambda t,y,v: -v # el signo, revisar diagrama cuerpo libre
g = lambda t,y,v: 9.8 - (0.225/90)*(v**2)
t0 = 0
y0 = 1000
v0 = 0
h = 2
muestras = 15+1
# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,t0,y0,v0,h,muestras, vertabla=True)
ti = tabla[:,0]
yi = tabla[:,1]
vi = tabla[:,2]
# SALIDA# print('tabla de resultados')# print(tabla)# GRAFICA
plt.subplot(211)
plt.plot(ti,vi,label='velocidad v(t)', color='green')
plt.plot(ti,vi,'o')
plt.ylabel('velocidad (m/s)')
plt.title('paracaidista Wingsuit con Runge-Kutta')
plt.legend()
plt.grid()
plt.subplot(212)
plt.plot(ti,yi,label='Altura y(t)',)
plt.plot(ti,yi,'o',)
plt.axhline(0, color='red')
plt.xlabel('tiempo (s)')
plt.ylabel('Altura (m)')
plt.legend()
plt.grid()
plt.show()
b. Usando Runge-Kutta 2do Orden para Primera Derivada o velocidad,
El problema para un tiempo de observación t>0, se puede dividir en dos partes: velocidad y altura.
Determinar velocidad: Se aplica Runge-Kutta a la expresión con primera derivada o velocidad. Use los valores iniciales dados, descarte calcular las alturas.
Determinar las altura: Con los valores de velocidades y la altura inicial de 1 km = 1000 m puede integrar con trapecio para obtener la tabla de alturas. Se integra tramo a tramo.
Observe las unidades de medida y que la velocidad es contraria al eje de altura dy/dt = -v.
La expresión:
\frac{dv}{dt} = g - \frac{cd}{m} v^2 f(t,v) = g - \frac{cd}{m} v^2
las siguientes iteraciones se completan con el algoritmo.
Resultados
La velocidad máxima, si no hubiese límite en la altura, se encuentra en alrededor de 62.39 m/s.
Sin embargo, luego de 20 segundos se observa que la altura alcanza el valor de cero, es decir se alcanzó tierra con velocidad de 62 m/s, que son algo mas de 223 Km/h, el objeto se estrella…!
Los cálculos se realizaron usando las instrucciones en Python:
# 2da Evaluación I Término 2018# Tema 1. Paracaidista wingsuitimport numpy as np
defrungekutta2(d1y,x0,y0,h,muestras):
# Runge Kutta de 2do orden
tamano = muestras + 1
tabla = np.zeros(shape=(tamano,2+2),dtype=float)
# incluye el punto [x0,y0]
tabla[0] = [x0,y0,0,0]
xi = x0
yi = y0
for i inrange(1,tamano,1):
K1 = h * d1y(xi,yi)
K2 = h * d1y(xi+h, yi + K1)
yi = yi + (1/2)*(K1+K2)
xi = xi + h
tabla[i] = [xi,yi,K1,K2]
return(tabla)
defintegratrapecio_fi_tabla(xi,fi,y0):
tamano = len(xi)
yi = np.zeros(tamano,dtype = float)
yi[0] = y0
for i inrange(1,tamano,1):
h = xi[i]-xi[i-1]
trapecio = h*(fi[i]+fi[i-1])/2
yi[i]= yi[i-1] + trapecio
return(yi)
# PROGRAMA -------------------------# INGRESO
g = 9.8
cd = 0.225
m = 90
d1v = lambda t,v: g - (cd/m)*(v**2)
t0 = 0
v0 = 0
h = 2
y0 = 1000
muestras = 11
# PROCEDIMIENTO
velocidad = rungekutta2(d1v,t0,v0,h,muestras)
ti = velocidad[:,0]
vi = velocidad[:,1]
# Altura, velocidad es contraria altura,# integrar en cada tramo por trapecios o Cuadratura de Gauss
altura = integratrapecio_fi_tabla(ti,-vi,y0)
# Tabla de resultados de tiempo, velocidad, altura
altura = np.transpose([altura])
tabla = np.concatenate((velocidad[:,0:2],altura), axis = 1)
# SALIDA
np.set_printoptions(precision=4)
print('velocidad con Runge-Kutta primera derivada')
print(' [tiempo, velocidad, K1,K2]')
print(velocidad)
np.set_printoptions(precision=2)
print('velocidad con RK2 y altura con trapecio')
print(' [tiempo, velocidad, altura]')
print(tabla)
En valor práctico 2.028 m usando flexómetro, a menos que use medidor laser con precisión 10-6 usará más dígitos con un tanque de 6 metros de altura ganará una precisión de una gota de agua para usar en duchas o regar el césped .
c) El orden de convergencia del error observando las tres primeras iteraciones es cuadrático
Tarea: Realizar los cálculos con Python, luego aplique otro método. Añada sus observaciones y conclusiones.
Observación: la matriz A ya es diagonal dominante, no requiere pivotear por filas. Se aumentó el punto decimal a los valores de la matriz A y el vector B para que sean considerados como números reales.
El número de condición es: np.linalg.cond(A) = 3.0
que es cercano a 1 en un orden de magnitud, por lo que la solución matricial es «estable» y los cambios en los coeficientes afectan proporcionalmente a los resultados. Se puede aplicar métodos iterativos sin mayores inconvenientes.
b y c) método de Jacobi para sistemas de ecuaciones, con vector inicial
X(0) = [[60.0],
[40],
[70],
[50]]
reemplazando los valores iniciales en cada ecuación sin cambios.
ti = [0.00, 0.15, 0.30, 0.45, 0.60, 0.75, 0.90, 1.05, 1.20]
xi = [0.00, 0.50, 1.00, 1.50, 1.80, 2.00, 1.90, 1.10, 0.30]
yi = [0.00, 4.44, 8.88,13.31,17.75,22.19,26.63,31.06,35.50]
zi = [0.00, 0.81, 1.40, 1.77, 1.91, 1.84, 1.55, 1.03, 0.30]
Observe que, un gol simplificado con física básica es un tiro parabólico cuya trayectoria se compone de movimientos en los ejes, Y y Z. Sin embargo, lo «imposible» del gol mostrado es añadir el movimiento en X. (Referencia de la imagen en el enunciado)
El movimiento de «profundidad» o dirección hacia el arco y(t) es semejante a un polinomio de primer grado, y el movimiento de «altura» z(t) es un polinomio de segundo grado. El movimiento «imposible» en el eje X, podría ser descrito por un polinomio de segundo o mayor grado.
a) Encontrar t para altura máxima, que se encuentra al igualar la derivada dz/dt a cero. Para interpolar el polinomio z(t), de segundo grado, se puede usar tres puntos de los sugeridos: 0, 0.3 y 0.6, que además son equidistantes en t (facilita usar cualquier método de interpolación).
Por ejemplo, con diferencias finitas avanzadas:
t
z(t)
d1
d2
d3
0.00
0.00
1.40
-0.89
0.30
1.40
0.51
0.60
1.91
z(t) = 0 + 1.40\frac{(t-0)}{1!(0.3)} - 0.89 \frac{(t-0)(t-0.3)}{2!(0.3)^2} = 0 + 4.66 t - 4.94(t^2-0.3t) = 4.66 t + 1.48 t - 4.94 t^2 z(t) = 6.42 t - 4.94 t^2
para encontrar el valor máximo se encuentra \frac{dz}{dt} = 0
\frac{dz}{dt} = 6.42 - 2(4.94) t 6.42 - 2(4.94) t = 0 t = \frac{6.42}{2(4.94)} t = 0.649
Observe que el resultado tiene sentido, pues según la tabla, el máximo debería estar entre 0.60 y 0.75
b) El cruce por la «barrera», corresponde al desplazamiento del balón en el eje Y a 9 metros: y(t) = 9.
El polinomio modelo de primer grado usa dos puntos para la interpolación, de los sugeridos se escogen dos, por ejemplo: 0.0 y 0.3.
29.6 t = 9
t = 0.30
z(0.30) = 1.40
(de la tabla de datos)
cuya respuesta con solo dos dígitos decimales es coherente, al estar cercano el valor a una distancia y=8.88 o aproximado a 9.
La respuesta puede ser mejorada usando más digitos significativos en las operaciones.
c) La desviación máxima en eje X.
Determine un polinomio para la trayectoria en el eje X y obtenga el máximo igualando la derivada del polinomio x(t) a cero.
Por simplicidad, se usa un polinomio de segundo grado, alrededor de los valores máximos en el eje X
lo que es coherente con la tabla para el eje x, pues el máximo registrado es 2, y el próximo valor es menor, la curva será decreciente.
Desarrollo extra para observar y verificar resultados:
Usando los puntos y un graficador 3D se puede observar la trayectoria:
Tarea: Realice el ejercicio usando los algoritmos en Python, muestre los polinomios obtenidos y compare.
Nota: La gráfica 3D presentada usa interpolación de Lagrange con todos los puntos. Realice las observaciones y recomendaciones necesarias y presente su propuesta como tarea.
Si el polinomio de Taylor fuera de grado 0, sería una constante, que si se evalúa en x0 = 0 para eliminar los otros términos, se encuentra que sería igual a 1
Como se pide el polinomio de grado 2, se tiene la forma:
p(x) = a + bx + c x ^2 p(x) = 1 + bx + c x^2
Se disponen de dos puntos adicionales que se pueden usar para determinar b y c:
p(0.2) = 1 + 0.2 b + (0.2)^2 c = 1.6 p(0.4) = 1 + 0.4 b + (0.4)^2 c = 2.0
simplificando:
0.2 b + (0.2)^2 c = 1.6 - 1 = 0.60.4 b + (0.4)^2 c = 2.0 - 1 = 1
multiplicando la primera ecuación por 2 y restando la segunda ecuación:
0 - 0.08 c = 1.2-1 = 0.2 c = - 0.2/0.08 = -2.5
sustituyendo el valor de c obtenido en la primera ecuación
0.2 b + 0.04(-2.5) = 0.6 0.2 b = 0.6 - 0.04(-2.5) = 0.6 + 0.1 = 0.7 b = 0.7/0.2 = 3.5
con lo que el polinomio queda: p(x) = 1 + 3.5 x - 2.5 x^2
validando con python:
tomando los puntos de prueba:
xi = [ 0, 0.2, 0.4]
fi = [ 1, 1.6, 2 ]
>>>
se obtiene la gráfica:
se adjunta las instrucciones usadas para validar que el polinomio pase por los puntos requeridos.
# 1Eva_IIT2017_T1 Aproximar a polinomio usando puntosimport numpy as np
import matplotlib.pyplot as plt
# INGRESO
xi = [ 0, 0.2, 0.4]
fi = [ 1, 1.6, 2 ]
px = lambda x: 1 + 3.5*x - 2.5*(x**2)
a = -0.5
b = 1
muestras = 21
# PROCEDIMIENTO
xj = np.linspace(a,b,muestras)
pxj = px(xj)
# SALIDAprint(xj)
print(pxj)
# Gráfica
plt.plot(xj,pxj,label='p(x)')
plt.plot(xi,fi,'o', label='datos')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.grid()
plt.legend()
plt.show()
Nota: Se puede intentar realizar los polinomios aumentando el grado, sin embargo cada término agrega un componente adicional a los términos anteriores por la forma (x – x0)k
literal b
se requiere el integral aproximado de f(x) en el intervalo limitado por los 3 puntos de la tabla:
\int_{0}^{0.4}f(x) dx
Esta aproximación con un polinomio es el concepto de integración numérica con la regla de Simpson de 1/3, tema desarrollado en la unidad 5
Runge kutta de 2do Orden
f: y' = z
g: z' = .....
K1y = h f(xi, yi, zi)
K1z = h g(xi, y1, zi)
K2y = h f(xi+h, yi+K1y, zi+K1z)
K2z = h g(xi+h, yi+K1y, zi+K1z)
y(i+1) = yi + (1/2)(K1y + K2y)
z(i+1) = zi + (1/2)(K1z + K2z)
x(i+1) = xi + h
E = O(h3)
xi ≤ z ≤ x(i+1)
Instruccciones en python, usando el algoritmo desarrollado en clase
# Runge Kutta de 2do# EDO de 2do orden con condiciones de inicioimport numpy as np
import matplotlib.pyplot as plt
defrungekutta2_fg(f,g,v0,h,m):
tabla = [v0]
xi = v0[0]
yi = v0[1]
zi = v0[2]
for i inrange(0,m,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)
yi1 = yi + (1/2)*(K1y+K2y)
zi1 = zi + (1/2)*(K1z+K2z)
xi1 = xi + h
vector = [xi1,yi1,zi1]
tabla.append(vector)
xi = xi1
yi = yi1
zi = zi1
tabla = np.array(tabla)
return(tabla)
# Programa Prueba# Funciones
f = lambda x,y,z : z
g = lambda x,y,z : (-gr/L)*np.sin(y)
gr = 9.8
L = 2
x0 = 0
y0 = np.pi/6
z0 = 0
v0 = [x0,y0,z0]
h = 0.1
xn = 2
m = int((xn-x0)/h)
# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,v0,h,m)
xi = tabla[:,0]
yi = tabla[:,1]
zi = tabla[:,2]
# SALIDA
np.set_printoptions(precision=6)
print('x, y, z')
print(tabla)
plt.plot(xi,yi, label='y')
plt.plot(xi,zi, label='z')
plt.legend()
plt.grid()
plt.show()
