Etiqueta: ejercicios resueltos Python

Ejercicio: 2Eva_IT2018_T2 Deducir Simpson 1/3

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)

\int_a^b f(x) dx \approx \int_a^b p(x) dx p(x) = L_a f(a) + L_c f(c) + L_b f(b) \int_a^b p(x) dx = \int_a^b \Big[ L_a f(a) +L_c f(c) + L_b f(b) \Big] dx = = \int_a^b L_a f(a) dx + \int_a^b L_c f(c) dx + \int_a^b L_b f(b) dx \int_a^b p(x) dx = I_1 + I_2 + I_3

Como referencia se usa la gráfica para relacionar a, b, c y h:

---

Primer Integral

Para el primer integral I_1= \int_a^b L_a f(a) dx se tiene que:

L_a = \frac{(x-b)(x-c)}{(a-b)(a-c)} = \frac{(x-b)(x-c)}{(-2h)(-h)} L_a = \frac{(x-b)(x-c)}{2h^2}

se convierte en:

I_1 = \int_a^b \frac{(x-b)(x-c)}{2h^2} f(a) dx = \frac{f(a)}{2h^2} \int_a^b (x-b)(x-c)dx

Para dejar la parte del integral en función de h, a y b, teniendo que c está en la mitad de [a,b], es decir c = (a+b)/2 , se usa:

u = x-a

por lo que \frac{du}{dx}=1 y du = dx

x-c = (u+a) - \frac{a+b}{2} = u+ \frac{a-b}{2} = u - \frac{b-a}{2} x-c = u-h x-b = (u+a)-b = u-2\Big(\frac{b-a}{2}\Big) = u-2h

Se actualiza el integral de x entre [a,b]  usando u = x-a, que se convierte el rango para u en [0, b-a] que es lo mismo que [0,2h]

\int_a^b (x-b)(x-c)dx = \int_0^{2h} (u-2h)(u-h)du = = \int_0^{2h} \Big( u^2 - 2hu - uh + 2h^2 \Big) du = \int_0^{2h} \Big( u^2 - 3hu + 2h^2 \Big) du = \frac{u^3}{3}- 3h\frac{u^2}{2}+ 2h^2u \Big|_0^{2h} = \frac{(2h)^3}{3}- 3h\frac{(2h)^2}{2} + 2h^2(2h) -(0-0+0) = \frac{8h^3}{3}- 6h^3 + 4h^3 =\frac{8h^3}{3}- 2h^3 = \frac{2h^3}{3}

resultado que se usa en I₁

I_1= \frac{f(a)}{2h^2}\frac{2h^3}{3} =\frac{h}{3} f(a)

que es el primer término de la fórmula general de Simpson 1/3

---

Segundo Integral

Para el Segundo integral I_2= \int_a^b L_c f(c) dx se tiene que:

L_c = \frac{(x-a)(x-b)}{(c-a)(c-b)} = \frac{(x-a)(x-b)}{(h)(-h)} L_c = \frac{(x-a)(x-b)}{-h^2}

se convierte en:

I_2 = \frac{f(c)}{-h^2} \int_a^b (x-a)(x-b) dx = \frac{f(c)}{-h^2} \int_0^{2h} (u)(u-2h) du

siendo:

\int_0^{2h}(u^2-2hu)du=\Big(\frac{u^3}{3}-2h\frac{u^2}{2}\Big)\Big|_0^{2h} =\frac{(2h)^3}{3}-h(2h)^2-(0-0) =\frac{8h^3}{3}-4h^3 = -\frac{4h^3}{3}

usando en I₂

I_2 = \frac{f(c)}{-h^2}\Big(-\frac{4h^3}{3}) = \frac{h}{3}4f(c)

---

Tarea: Continuar las operaciones para y tercer integral para llegar a la fórmula general de Simpson 1/3:

I = \frac{h}{3} \Big( f(a)+4f(c) + f(b) \Big)

Ejercicio: 2Eva_IT2018_T4 Dragado acceso marítimo

a) La matriz para remover sedimentos se determina como la diferencia entre la profundidad y la matriz de batimetría.

Considere el signo de la profundidad para obtener el resultado:

matriz remover sedimentos: 
[[ 4.21  0.    0.96  0.    3.76  3.09]
 [ 2.15  0.11  2.05  3.77  0.    3.07]
 [ 0.    1.14  1.65  0.    1.62  1.35]
 [ 3.7   0.    0.59  2.33  0.    4.23]
 [ 0.    1.38  3.53  4.49  1.98  1.4 ]
 [ 0.    0.77  0.32  1.06  4.24  3.54]]

se obtiene con la instrucciones:

# 2da Evaluación I Término 2018
# Tema 4. canal acceso a Puertos de Guayaquil
import 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 in range(0,n,1):
    for j in range(0,m,1):
        if remover[i,j]<0:
            remover[i,j]=0
# SALIDA
print('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:

# literal b) ---------------------------
def simpson13(fi,h):
    s13 = (h/3)*(fi[0] + 4*fi[1] + fi[2])
    return(s13)
def simpson38(fi,h):
    s38 = (3*h/8)*(fi[0] + 3*fi[1] + 3*fi[2]+ fi[3])
    return(s38)

Integralx = np.zeros(n,dtype = float)

for i in range(0,n,1):
    hx = xi[1]-xi[0]
    fi = remover[i, 0:(0+4)]
    s38 = simpson38(fi,hx)
    fi = remover[i, 3:(3+3)]
    s13 = simpson13(fi,hx)
    Integralx[i] = s38 + s13

hy = yi[1] - yi[0]
fj = Integralx[0:(0+4)]
s38 = simpson38(fj,hy)
fj = Integralx[3:(3+3)]
s13 = simpson13(fj,hy)
volumen = s38 + s13

# Salida
np.set_printoptions(precision=2)
print('Integral en eje x: ')
print(Integralx)
print('Volumen: ', volumen)

Para el examen escrito, se requieren realizar al menos 3 iteraciones/ filas para encontrar el integral.

Ejercicio: 2Eva_IT2018_T1 Paracaidista wingsuit

El problema para un tiempo de observación t>0, se puede dividir en dos partes: velocidad y altura.

1. Determinar velocidad: Se aplica Runge-Kutta con los valores iniciales dados, descartar calcular las alturas.
2. Determinar las altura:  Con los valores de velocidades y la altura inicial de 1 km = 1000 m puede integrar cada tramo.
   Observe las unidades de medida y que la  velocidad es contraria  al eje de altura dy/dt = -v.
3. En el ejercicio se presentan los resultados integrando por trapecios y como otra forma usando la fórmula de Taylor, tomando los puntos encontrados para velocidad.
4. Otra opción para las alturas, en caso de requerir determinar el tiempo exacto en que sucede el resultado es sustituir v = dy/ty y resolver nuevamente la ecuación inicial. (Tarea).

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/s, el objeto se estrella…!

altura1 con trapecio, altura2 con Taylor
 [tiempo, velocidad, altura1, altura2]
[[    0.       0.    1000.    1000.  ]
 [    2.      18.64   981.36   980.4 ]
 [    4.      34.04   928.68   925.26]
 [    6.      45.02   849.62   843.37]
 [    8.      52.13   752.47   743.87]
 [   10.      56.49   643.85   633.59]
 [   12.      59.07   528.3    516.97]
 [   14.      60.58   408.65   396.68]
 [   16.      61.45   286.63   274.28]
 [   18.      61.94   163.24   150.66]
 [   20.      62.23    39.07    26.36]
 [   22.      62.39   -85.55   -98.34]]

Los cálculos se realizaron usando las instrucciones en Python:

# 2da Evaluación I Término 2018
# Tema 1. Paracaidista wingsuit
import numpy as np

def rungekutta2(d1y,x0,y0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,2),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0]
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1 = h * d1y(xi,yi)
        K2 = h * d1y(xi+h, yi + K1)

        yi = yi + (K1+K2)/2
        xi = xi + h
        
        estimado[i] = [xi,yi]
    return(estimado)

def integratrapecio(xi,fi,y0):
    tamano = len(xi)
    yi = np.zeros(tamano,dtype = float)
    yi[0] = y0
    for i in range(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
t0 = 0
v0 = 0

dt = 2
y0 = 1000
muestras = 11

d1v = lambda t,v: g - (cd/m)*(v**2)

# PROCEDIMIENTO
velocidad = rungekutta2(d1v,t0,v0,dt,muestras)
ti = velocidad[:,0]
vi = velocidad[:,1]

# Altura, velocidad es contraria altura,
# integrar en cada tramo por trapecios o Cuadratura de Gauss
altura = integratrapecio(ti,-vi,y0)

# Tabla de resultados
altura = np.transpose([altura])
tabla = np.concatenate((velocidad,altura), axis = 1)

# otra forma: usando Taylor
altura2 = np.zeros(muestras+1,dtype = float)
altura2[0] = y0
for i in range(0, muestras,1):
    acelera = d1v(ti[i],vi[i])
    altura2[i+1] = altura2[i] -(vi[i]*dt + acelera*(dt**2)/2)
# Tabla de resultados
altura2 = np.transpose([altura2])
tabla = np.concatenate((tabla,altura2), axis = 1)

# SALIDA
np.set_printoptions(precision=2)
print('altura1 con trapecio, altura2 con Taylor')
print(' [tiempo, velocidad, altura1,altura2]')
print(tabla)

y para la gráfica:

# Gráfica
import matplotlib.pyplot as plt
plt.subplot(211)
plt.plot(ti,vi,'g')
plt.plot(ti,vi,'bo')
plt.title('paracaidista Wingsuit')
plt.ylabel('velocidad (m/s)')
plt.subplot(212)
plt.plot(ti,altura)
plt.plot(ti,altura,'o')
plt.plot(ti,altura2)
plt.axhline(0, color= 'red')
plt.ylabel('altura (m)')
plt.xlabel('tiempo (s)')
plt.show()

Ejercicio: 2Eva_IIT2017_T4 Problema con valor de frontera

\frac{d^2T}{dx^2} + \frac{1}{x}\frac{dT}{dx} +S =0 0 \leq x \leq 1

Las diferencias finitas a usar son:

\frac{dT}{dx} =\frac{T_{i+1} - T_{i-1}}{2h}+O(h^2) \frac{d^2T}{dx^2}=\frac{T_{i+1} - 2T_i + T_{i-1}}{h^2}+O(h^2)

que al reemplazar el la ecuación:

\frac{T_{i+1} - 2T_i + T_{i-1}}{h^2} + \frac{1}{x_i}\frac{T_{i+1} -T_{i-1}}{2h}+S=0 2x_i (T_{i+1} - 4h x_i T_i + T_{i-1}) + h (T_{i+1} - T_{i-1}) = -2h^2 S x_i T_{i+1}(2x_i + h) - 4x_i T_i + T_{i-1}(2x_i - h) = -2h^2 S x_i T_{i-1}(2x_i - h) - 4x_i T_i + T_{i+1}(2x_i + h)= -2h^2 S x_i

con lo que se puede crear un sistema de ecuaciones para cada valor xi con T₀=2 y T₄ =1 que son parte del enunciado del problema.

Siendo h = 0.25 = 1/4,  y se indica al final que S=1, se crea un sistema de ecuaciones a resolver,

x = [0, 1/4, 1/4, 3/4, 1]

T_{i-1}\Big[2x_i - \frac{1}{4} \Big] - 4x_i T_i + T_{i+1}\Big[2x_i + \frac{1}{4} \Big] = -2\Big(\frac{1}{4}\Big)^2 (1) x_i T_{i-1}\Big[2x_i -\frac{1}{4}\Big] - 4x_i T_i + T_{i+1}\Big[2x_i + \frac{1}{4}\Big] = -\frac{1}{8} x_i

se sustituye con los valores conocidos para cada i:

i=1: 
T0[2(1/4) - 1/4] - 4(1/4)T1 + T2[2(1/4) + 1/4] = -(1/8)(1/4)

     - T1 + (3/4)T2 = -1/32 - (1/4)(2)
     - T1 + (3/4)T2 = -17/32

i=2: 
T1[2(1/2) - 1/4] - 4(1/2)T2 + T3[2(1/2) + 1/4] = -(1/8)(1/2)

     (3/4)T1 - 2T2 + (5/4)T3 = -1/16

i=3: 
T2[2(3/4) - 1/4] - 4(3/4)T3 + T4[2(3/4) + 1/4] = -(1/8)(3/4)

     (5/4)T2 - 3T3 = -3/32 - (7/4)(1)
     (5/4)T2 - 3T3 = -59/32

se ponen las ecuaciones en matrices para resolver con algun metodo numérico:

A = [[ -1, 3/4,   0],
     [3/4,  -2, 5/4],
     [  0, 5/4,  -3]]
B = [-17/32, -1/16, -59/32]
np.linalg.solve(A,B)
array([ 1.54,  1.34,  1.17])

Ejercicio: 2Eva_IIT2017_T1 EDO Runge Kutta 2do Orden d2y/dx2

Tema 1

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)

f: z = Θ’
g: z’ = (-gr/L) sin(Θ)

Θ(0) = π/6
z(0) = 0

h=0.1

i=0, t0 = 0, Θ0 = π/6, z0 = 0
    K1y = 0.1(0) = 0
    K1z = 0.1(-9.8/2)sin(π/6) = -0.245

    K2y = 0.1(0+(-0.245)) = -0.0245
    K2z = 0.1(-9.8/2)sin(π/6+0) = -0.245

    Θ1 = π/6 + (1/2)(0+(-0.0245)) = 0.51139
    z1 = 0 + (1/2)(-0.245-0.245) = -0.245
    t1 = 0 + 0.1 = 0.1

i=1, t1 = 0.1, Θ1 = 0.51139, z1 = -0.245
    K1y = 0.1(-0.245) = -0.0245
    K1z = 0.1(-9.8/2)sin(0.51139) = -0.23978

    K2y = 0.1(-0.245+(-0.0245)) = -0.049
    K2z = 0.1(-9.8/2)sin(0.51139+(-0.0245)) = -0.22924

    Θ2 = 0.51139 + (1/2)(-0.0245+(-0.049)) = 0.47509
    z2 = -0.245 + (1/2)(-0.23978+(-0.22924)) = -0.245
    t2 = 0.1 + 0.1 = 0.2

   t         theta     z
[[ 0.        0.523599  0.      ]
 [ 0.1       0.511349 -0.245   ]
 [ 0.2       0.47486  -0.479513]
 [ 0.3       0.415707 -0.692975]
 [ 0.4       0.336515 -0.875098]
 [ 0.5       0.240915 -1.016375]
 [ 0.6       0.133432 -1.108842]
 [ 0.7       0.019289 -1.14696 ]
 [ 0.8      -0.09588  -1.128346]
 [ 0.9      -0.206369 -1.054127]
 [ 1.       -0.306761 -0.92877 ]
 [ 1.1      -0.39224  -0.759461]
 [ 1.2      -0.458821 -0.555246]
 [ 1.3      -0.503495 -0.326207]
 [ 1.4      -0.524294 -0.082851]
 [ 1.5      -0.520315  0.164197]
 [ 1.6      -0.491715  0.404296]
 [ 1.7      -0.439718  0.62682 ]
 [ 1.8      -0.366606  0.821313]
 [ 1.9      -0.275693  0.977893]
 [ 2.       -0.171235  1.087942]]

Literal b), con h= 0.25, con t = 1 ángulo= -0.352484

   t         theta     z
[[ 0.        0.523599  0.      ]
 [ 0.25      0.447036 -0.6125  ]
 [ 0.5       0.227716 -1.054721]
 [ 0.75     -0.070533 -1.170971]
 [ 1.       -0.352484 -0.910162]
 [ 1.25     -0.527161 -0.363031]
 [ 1.5      -0.540884  0.299952]
 [ 1.75     -0.387053  0.890475]
 [ 2.       -0.106636  1.221932]]

El error de del orden h³

---

Instruccciones en python, usando el algoritmo desarrollado en clase

# Runge Kutta de 2do
# EDO de 2do orden con condiciones de inicio
import numpy as np
import matplotlib.pyplot as plt

def rungekutta2_fg(f,g,v0,h,m):
    tabla = [v0]
    xi = v0[0]
    yi = v0[1]
    zi = v0[2]
    for i in range(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()