s2Eva_IIT2018_T3 EDP

Tema 3.  Se indica en el enunciado que b = 0

\frac{\delta u}{\delta t} = \frac{\delta ^2 u}{\delta x^2} + b\frac{\delta u}{\delta x}

simplificando la ecuación a:

\frac{\delta u}{\delta t} = \frac{\delta ^2 u}{\delta x^2}

Reordenando la ecuación a la forma estandarizada:

\frac{\delta ^2 u}{\delta x^2} = \frac{\delta u}{\delta t}

Seleccione un método: explícito o implícito.
Si el método es explícito, las diferencias finitas a usar son hacia adelante y centrada:

U'(x_i,t_j) = \frac{U(x_i,t_{j+1})-U(x_i,t_j)}{\Delta t} + O(\Delta t) U''(x_i,t_j) = \frac{U(x_{i+1},t_j)-2U(x_{i},t_j)+U(x_{i-1},t_j)}{\Delta x^2} + O(\Delta x^2)

como referencia se usa la gráfica.

Se selecciona la esquina inferior derecha como 0,  por la segunda ecuación de condiciones y facilidad de cálculo. (No hubo indicación durante el examen que muestre lo contrario)

condiciones de frontera U(0,t)=0, U(1,t)=1
condiciones de inicio U(x,0)=0, 0≤x≤1

aunque lo más recomendable sería cambiar la condición de inicio a:

condiciones de inicio U(x,0)=0, 0<x<1

Siguiendo con el tema de la ecuación, al reemplazar las diferencias finitas en la ecuación:


\frac{U(x_{i+1},t_j)-2U(x_{i},t_j)+U(x_{i-1},t_j)}{\Delta x^2} = = \frac{U(x_i,t_{j+1})-U(x_i,t_j)}{\Delta t}

se reagrupan los términos que son constantes y los términos de error se acumulan:

\frac{\Delta t}{\Delta x^2} \Big[U(x_{i+1},t_j)-2U(x_i,t_j)+U(x_{i-1},t_j) \Big] = U(x_i,t_{j+1})-U(x_i,t_j)

siendo,

\lambda= \frac{\Delta t}{\Delta x^2} error \cong O(\Delta t) + O(\Delta x^2)

continuando con la ecuación, se simplifica la escritura usando sólo los índices i,j y se reordena de izquierda a derecha como en la gráfica

\lambda \Big[U[i-1,j]-2U[i,j]+U[i+1,j] \Big] = U[i,j+1]-U]i,j] \lambda U[i-1,j]+(-2\lambda+1)U[i,j]+\lambda U[i+1,j] = U[i,j+1] U[i,j+1] = \lambda U[i-1,j]+(-2\lambda+1)U[i,j]+\lambda U[i+1,j] U[i,j+1] = P U[i-1,j]+QU[i,j]+R U[i+1,j] P=R = \lambda Q = -2\lambda+1

En las iteraciones, el valor de P,Q y R se calculan a partir de λ ≤ 1/2

iteraciones: j=0, i=1

U[1,1] = P*0+Q*0+R*0 = 0

j=0, i=2

U[2,1] = P*0+Q*0+R*0=0

j=0, i=3

U[3,1] = P*0+Q*0+R*1=R=\lambda=\frac{1}{2}

iteraciones: j=1, i=1

U[1,2] = P*0+Q*0+R*0 = 0

j=1, i=2

U[2,2] = P*0+Q*0+R*\lambda = \lambda ^2 = \frac{1}{4}

j=1, i=3

U[3,2] = P*0+Q*\frac{1}{4}+R (\lambda) U[3,2] = (-2\lambda +1) \frac{1}{4}+\lambda^2 = \Big(-2\frac{1}{2}+1\Big) \frac{1}{4}+\Big(\frac{1}{2}\Big)^2 U[3,2] =0\frac{1}{4} + \frac{1}{4} = \frac{1}{4}

Literal b. Para el cálculo del error:

\lambda \leq \frac{1}{2} \frac{\Delta t}{\Delta x^2} \leq \frac{1}{2} \Delta t \leq \frac{\Delta x^2}{2}

en el enunciado se indica h = 0.25 = ¼ = Δ x

\Delta t \leq \frac{(1/4)^2}{2} = \frac{1}{32} error \cong O(\Delta t) + O(\Delta x^2) error \cong \frac{\Delta x^2}{2}+ \Delta x^2 error \cong \frac{3}{2}\Delta x^2 error \cong \frac{3}{2}( \frac{1}{4})^2 error \cong \frac{3}{32} = 0.09375

s2Eva_IIT2018_T2 Kunge Kutta 2do Orden x»

Tema 2.

\frac{\delta ^2 x}{\delta t^2} + 5t\frac{\delta x}{\delta t} +(t+7)\sin (\pi t) = 0 x'' + 5tx' +(t+7)\sin (\pi t) = 0 x'' = -5tx' +(t+7)\sin (\pi t) = 0

si se usa z=x’

z' = -5tz +(t+7)\sin (\pi t) = 0

se convierte en:
f(t,x,z) = x’ = z
g(t,x,z) = x» = z’ = -5tz +(t+7)sin (π t) = 0

Donde se aplica el algoritmo de Runge Kutta
http://blog.espol.edu.ec/matg1013/8-2-2-runge-kutta-d2y-dx2/

   t,              x,              z
[[ 0.00000000e+00  6.00000000e+00  1.50000000e+00]
 [ 2.00000000e-01  6.30000000e+00  1.77320538e+00]
 [ 4.00000000e-01  6.70381805e+00  2.26987703e+00]
 [ 6.00000000e-01  7.20775473e+00  2.41163944e+00]
 [ 8.00000000e-01  7.68994485e+00  1.90531839e+00]
 [ 1.00000000e+00  8.01027755e+00  9.52659193e-01]
 [ 1.20000000e+00  8.10554347e+00 -5.65431040e-03]
 [ 1.40000000e+00  8.00869435e+00 -6.09147239e-01]
 [ 1.60000000e+00  7.81236802e+00 -7.16247408e-01]
 [ 1.80000000e+00  7.62013640e+00 -3.92947221e-01]
 [ 2.00000000e+00  7.50882725e+00  1.63598524e-01]]

instrucciones Python:

# 2Eva_IIT2018_T2 Kunge Kutta 2do Orden x''
import numpy as np

def rungekutta2_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 in range(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
# INGRESO
f = lambda t,x,z: z
g = lambda t,x,z: -5*t*z+(t+7)*np.sin(np.pi*t)
t0 = 0
x0 = 6
z0 = 1.5
h = 0.2
muestras = 10

# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,t0,x0,z0,h,muestras)

# SALIDA
print(tabla)
# GRAFICA
import matplotlib.pyplot as plt
plt.plot(tabla[:,0],tabla[:,1])
plt.xlabel('t')
plt.ylabel('x(t)')
plt.show()

s2Eva_IIT2018_T1 Masa entra o sale de un reactor

a) Se pueden combinar los métodos para realizar la integral. Se usa el método de Simpson 1/3 para los primeros dos tramos y Simpson 3/8 para los 3 tramos siguientes.  Siendo f(x) equivalente a Q(t)C(t). El tamaño de paso h es constante para todo el ejercicio con valor 5.

a.1 Simpson 1/3, tramos 2, puntos 3:

I_1 \cong \frac{h}{3}[f(x_0)+4f(x_1) + f(x_2)] I_1 \cong \frac{5}{3}[(10)(4)+4(18)(6) + (27)(7)] I_1 \cong 1101,66

a.2 Simpson de 3/8, tramos 3, puntos 4:

I_2 \cong \frac{3h}{8}[f(x_0)+3f(x_1) +3 f(x_2)+f(x_3)] I_2 \cong \frac{3(5)}{8}[(27)(7)+3(35)(6) +3(40)(5)+(30)(5)] I_2 \cong 2941,88 I_1 + I_2 \cong 4043,54

b) El error se calcula por tramo y se acumula.

b.1 se puede estimar como la diferencia entre la parábola del primer tramo y simpson 1/3
b.2 siguiendo el ejemplo anterior, como la diferencia entre la interpolación de los tramos restantes y simpson 3/8.

s2Eva_IT2010_T2 Movimiento angular

Propuesta de solución

# 2da Evaluación I Término 2010
# Tema 2. Movimiento Angular
import numpy as np
import matplotlib.pyplot as plt

f = lambda x, y: -10*np.sin(y)

x0 = 0
y0 = 0
dy0 = 0.1
h = 0.1
muestras = 11

# RUNGE KUTTA
estimado = np.zeros(shape=(muestras+1,3),dtype=float)
estimado[0] = [x0,dy0,0]
xi = x0
vi = dy0
i = 0
while not(i>=muestras):
    K1 = h*f(xi,vi)

    K2 = h*f(xi+h,vi + K1)

    vi = vi + (K1+K2)/2
    xi = xi + h
    i = i+1

    estimado[i] = [xi,vi,0]
    print(estimado[i])

# INTEGRANDO PARA OBTENER y
i = 0
suma = 0
while not(i>=(muestras)):
    dI = h*(estimado[i,1]+estimado[i+1,1])/2
    suma = suma +dI
    i = i+1
    estimado[i,2] = suma

# Grafica
plt.subplot(211)
plt.plot(estimado[:,0],estimado[:,1], label='velocidad')
plt.ylabel('Velocidad')
plt.subplot(212)
plt.plot(estimado[:,0],estimado[:,2], label='altura')
plt.ylabel('Altura')
plt.legend()
plt.show()

s2Eva_IT2018_T3 EDP Eliptica

Generar las ecuaciones a resolver usando diferencias finitas divididas centradas:

\frac{\delta ^2 u}{\delta x^2} + \frac{\delta ^2 u}{\delta y^2} = 2(x^2+y^2)

por facilidad se sustituye también:

f[i,j] = f(x,y) = 2(x^2+y^2)
\frac{u[i-1,j]-2u[i,j]+u[i+1,j]}{\Delta x^2} + + \frac{u[i,j-1]-2u[i,j]+u[i,j+1]}{\Delta y^2} = f[i,j]
\frac{\Delta y^2}{\Delta x^2}\Big(u[i-1,j]-2u[i,j]+u[i+1,j]\Big) + + u[i,j-1]-2u[i,j]+u[i,j+1] = \Delta y^2 f[i,j]

dado que hx = hy = 1/3

\frac{\Delta y^2}{\Delta x^2} = 1
(u[i-1,j]-2u[i,j]+u[i+1,j]) + + u[i,j-1]-2u[i,j]+u[i,j+1] = = \Big(\frac{1}{3}\Big)^2 f[i,j]
u[i-1,j]-4u[i,j]+u[i+1,j] + + u[i,j-1]+u[i,j+1] = \frac{f[i,j]}{9}
4u[i,j] = u[i-1,j]+u[i+1,j] + + u[i,j-1]+u[i,j+1]-\frac{f[i,j]}{9}
u[i,j] = \frac{1}{4} \Big(u[i-1,j]+u[i+1,j] + + u[i,j-1]+u[i,j+1]-\frac{f[i,j]}{9} \Big)

La última ecuación puede ser usada de forma iterativa, para lo cual hay que definir los valores iniciales de la matriz u.

Al conocer el rango de operación para los ejes x, y, hx, hy se realizan los cálculos para:

1. Evaluar los valores para cada eje x[i], y[j]

x[i]:
[ 0.    0.33  0.67  1.  ]
y[j]:
[ 0.    0.33  0.67  1.  ]

2. Evaluar en cada punto generando una matriz f(i,j):

f[i,j]:
[[ 0.    0.22  0.89  2.  ]
 [ 0.22  0.44  1.11  2.22]
 [ 0.89  1.11  1.78  2.89]
 [ 2.    2.22  2.89  4.  ]]

3. Se evaluan las funciones indicadas para la frontera y se tiene la matriz inicial para u:

matriz inicial u[i,j]:
[[ 1.    1.33  1.67  2.  ]
 [ 1.33  0.    0.    2.44]
 [ 1.67  0.    0.    3.11]
 [ 2.    2.44  3.11  4.  ]]

con lo que se puede trabajar cada punto i,j de forma iterativa, teniendo como resultado para la matriz u:

resultado para u, iterando: 
converge =  1
[[ 1.    1.33  1.67  2.  ]
 [ 1.33  1.68  2.05  2.44]
 [ 1.67  2.05  2.53  3.11]
 [ 2.    2.44  3.11  4.  ]]

La gráfica usando una mayor resolución para tener una idea de la solución:


Los resultados se obtienen usando las siguientes instrucciones:

# 2da Evaluación I Término 2018
# Tema 3. EDP Eliptica
import numpy as np

# INGRESO
# ejes x,y
x0 = 0 ; xn = 1 ; hx = 1/3 # (1/3)/10
y0 = 0 ; yn = 1 ; hy = 1/3 # (1/3)/10
# Fronteras
fux0 = lambda x: x+1
fu0y = lambda y: y+1
fux1 = lambda x: x**2 + x + 2
fu1y = lambda y: y**2 + y + 2

fxy = lambda x,y: 2*(x**2+y**2)

# PROCEDIMIENTO
xi = np.arange(x0,xn+hx,hx)
yj = np.arange(y0,yn+hx,hy)
n = len(xi)
m = len(yj)
# funcion f[xi,yi]
fij = np.zeros(shape=(n,m), dtype = float)
for i in range(0,n,1):
    for j in range(0,m,1):
        fij[i,j]=fxy(xi[i],yj[j])
# matriz inicial u[i,j]
u = np.zeros(shape=(n,m), dtype = float)
u[:,0] = fux0(xi)
u[0,:] = fu0y(yj)
u[:,m-1] = fux1(xi)
u[n-1,:] = fu1y(yj)

uinicial = u.copy()

# Calcular de forma iterativa
maxitera = 100
tolera = 0.0001
# valor inicial de iteración
promedio = (np.max(u)+np.min(u))/2
u[1:n-1,1:m-1] = promedio
# iterar puntos interiores
itera = 0
converge = 0
while not(itera>=maxitera and converge==1):
    itera = itera +1
    nueva = np.copy(u)
    for i in range(1,n-1):
        for j in range(1,m-1):
            u[i,j] = (u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1]-fij[i,j]/9)/4
    diferencia = nueva-u
    erroru = np.linalg.norm(np.abs(diferencia))
    if (erroru<tolera):
        converge=1

# SALIDA
np.set_printoptions(precision=2)
print('x[i]:')
print(xi)
print('y[j]:')
print(yj)
print('f[i,j]:')
print(fij)
print('matriz inicial u[i,j]:')
print(uinicial)
print('resultado para u, iterando: ')
print('converge = ', converge)
print(u)

para obtener la gráfica se debe añadir:

# Gráfica
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

X, Y = np.meshgrid(xi, yj)
U = np.transpose(u) # ajuste de índices fila es x
figura = plt.figure()
ax = Axes3D(figura)
ax.plot_surface(X, Y, U, rstride=1, cstride=1, cmap=cm.Reds)
plt.title('EDP elíptica')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

s2Eva_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 I1

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 I2

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)

s2Eva_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)
    f = remover[i, 3:(3+3)]
    s13 = simpson13(fi,hx)
    Integralx[i] = s38 + s13

hy = yi[1] - yi[0]
f = Integralx[0:(0+4)]
s38 = simpson38(f,hy)
f = Integralx[3:(3+3)]
s13 = simpson13(f,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.

s2Eva_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()

s2Eva_IT2015_T2 Deflexión de mástil

\frac{\delta ^2 y}{\delta x^2} = \frac{f}{2EI} (L-x)^2 x= 0, y = 0, \frac{\delta y}{\delta x} = 0

Ecuación en forma estandarizada:

y'' = \frac{f}{2EI} (L-x)^2

Sistema de ecuaciones

y' = z = f(x,y,z) z' = \frac{f}{2EI} (L-x)^2 = g(x,y,z)

Observe la diferencia de escalas entre los ejes.

los datos obtenidos usando el algoritmo son:

Runge Kutta 2do Orden
 [x, y, z]
[[  0.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  1.00000000e+00   2.16000000e-03   4.17840000e-03]
 [  2.00000000e+00   8.35680000e-03   8.07840000e-03]
 [  3.00000000e+00   1.83168000e-02   1.17096000e-02]
 [  4.00000000e+00   3.17760000e-02   1.50816000e-02]
 [  5.00000000e+00   4.84800000e-02   1.82040000e-02]
 ...
 [  2.90000000e+01   9.29856000e-01   4.32216000e-02]
 [  3.00000000e+01   9.73080000e-01   4.32240000e-02]
 [  3.10000000e+01   1.01630400e+00   4.32264000e-02]]
Runge Kutta 4do Orden
 [x, y, z]
[[  0.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  1.00000000e+00   2.11240000e-03   4.17760000e-03]
 [  2.00000000e+00   8.26240000e-03   8.07680000e-03]
 [  3.00000000e+00   1.81764000e-02   1.17072000e-02]
 [  4.00000000e+00   3.15904000e-02   1.50784000e-02]
 [  5.00000000e+00   4.82500000e-02   1.82000000e-02]
 ...
 [  2.90000000e+01   9.28800400e-01   4.31984000e-02]
 [  3.00000000e+01   9.72000000e-01   4.32000000e-02]
 [  3.10000000e+01   1.01520040e+00   4.32016000e-02]]
>>> 

las instrucciones en python obtener los resultados:

# 2da Evaluación I Término 2015
# Tema 2. Deflexión de mástil
import numpy as np

def rungekutta2fg(fx,gx,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 in range(1,tamano,1):
        K1y = h * fx(xi,yi,zi)
        K1z = h * gx(xi,yi,zi)
        
        K2y = h * fx(xi+h, yi + K1y, zi + K1z)
        K2z = h * gx(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)

def rungekutta4fg(fx,gx,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 in range(1,tamano,1):
        K1y = h * fx(xi,yi,zi)
        K1z = h * gx(xi,yi,zi)
        
        K2y = h * fx(xi+h/2, yi + K1y/2, zi + K1z/2)
        K2z = h * gx(xi+h/2, yi + K1y/2, zi + K1z/2)
        
        K3y = h * fx(xi+h/2, yi + K2y/2, zi + K2z/2)
        K3z = h * gx(xi+h/2, yi + K2y/2, zi + K2z/2)

        K4y = h * fx(xi+h, yi + K3y, zi + K3z)
        K4z = h * gx(xi+h, yi + K3y, zi + K3z)

        yi = yi + (K1y+2*K2y+2*K3y+K4y)/6
        zi = zi + (K1z+2*K2z+2*K3z+K4z)/6
        xi = xi + h
        
        estimado[i] = [xi,yi,zi]
    return(estimado)

# INGRESO
f = 60
L = 30
E = 1.25e8
I = 0.05
x0 = 0
y0 = 0
z0 = 0
tramos = 30

fx = lambda x,y,z: z
gx = lambda x,y,z: (f/(2*E*I))*(L-x)**2

# PROCEDIMIENTO
muestras = tramos + 1
h = L/tramos
tabla2 = rungekutta2fg(fx,gx,x0,y0,z0,h,muestras)
xi2 = tabla2[:,0]
yi2 = tabla2[:,1]
zi2 = tabla2[:,2]

tabla4 = rungekutta4fg(fx,gx,x0,y0,z0,h,muestras)
xi4 = tabla4[:,0]
yi4 = tabla4[:,1]
zi4 = tabla4[:,2]

# SALIDA
print('Runge Kutta 2do Orden')
print(' [x, y, z]')
print(tabla2)
print('Runge Kutta 4do Orden')
print(' [x, y, z]')
print(tabla4)

# GRAFICA
import matplotlib.pyplot as plt
plt.plot(xi2,yi2, label='Runge-Kutta 2do Orden')
plt.plot(xi4,yi4, label='Runge-Kutta 4do Orden')
plt.title('Deflexión de mástil')
plt.xlabel('x')
plt.ylabel('y: deflexión')
plt.legend()
plt.grid()
plt.show()

s2Eva_IIT2017_T4 Problema con valor de frontera

2da Evaluación II Término 2017-2018. Febrero 7, 2018

Tema 4.

\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 T0=2 y T4 =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])