Autor: Edison Del Rosario

Ejercicio: 2Eva_IIT2007_T2_AN Lanzamiento vertical proyectil

la ecuación del problema:

m \frac{\delta v}{\delta t} = -mg - kv|v|

se despeja:

\frac{\delta v}{\delta t} = -g - \frac{k}{m}v|v|

y usando los valores indicados en el enunciado:

\frac{\delta v}{\delta t} = -9,8 - \frac{0.002}{0.11}v|v|

con valores iniciales de:

t₀ = 0 , v₀ = 8 , h=0.2

Como muestra inicial, se usa Runge-Kutta de 2do Orden

iteración 1

K_1 = h\frac{\delta v}{\delta t}(0, 8) = 0.2[-9,8 - \frac{0.002}{0.11}8|8|] = -2.1927 K_2 = h\frac{\delta v}{\delta t}(0+0.2, 8 -2.1927 ) = 0.2[-9,8 - \frac{0.002}{0.11}(8 -2.1927)|8 -2.1927|] =-2.0826 v_1 = -9,8 +\frac{-2.1927-2.0826 }{2} = 5.8623 t_1 = t_0 + h = 0 + 0.2 = 0.2 error = O(h^3) = O(0.2^3) = O(0.008)

iteración 2

K_1 = h\frac{\delta v}{\delta t}(0.2, 5.8623) = 0.2[-9,8 - \frac{0.002}{0.11}(5.8623)|5.8623|] = -2.085 K_2 = h\frac{\delta v}{\delta t}(0+0.2, 5.8623 -2.085) = 0.2[-9,8 - \frac{0.002}{0.11}(5.8623 -2.085)|5.8623 -2.085|] =-2.0119 v_2 = -9,8 +\frac{-2.085-2.0119}{2} = 3.8139 t_2 = t_1 + h = 0.2 + 0.2 = 0.4

iteración 3

K_1 = h\frac{\delta v}{\delta t}(0.4, 3.8139) = 0.2[-9,8 - \frac{0.002}{0.11}( 3.8139)| 3.8139|] = -2.0129 K_2 = h\frac{\delta v}{\delta t}(0+0.2, 3.8139 -2.0129) = 0.2[-9,8 - \frac{0.002}{0.11}(3.8139 -2.0129)|3.8139 -2.0129|] =-1.9718 v_3 = -9,8 +\frac{-2.0129-1.9718}{2} = 1.8215 t_3 = t_2 + h = 0.4 + 0.2 = 0.6

---

Tabla y gráfica del ejercicio para todo el intervalo:

 [xi,      yi,     K1,     K2    ]
[[ 0.      8.      0.      0.    ]
 [ 0.2     5.8623 -2.1927 -2.0826]
 [ 0.4     3.8139 -2.085  -2.0119]
 [ 0.6     1.8215 -2.0129 -1.9718]
 [ 0.8    -0.1444 -1.9721 -1.9599]
 [ 1.     -2.0964 -1.9599 -1.9439]]

Algoritmo en Python

# 3Eva_IT2009_T2 EDO Taylor Seno(x)
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
d1y = lambda t,v: -9.8-(0.002/0.11)*v*np.abs(v)
x0 = 0
y0 = 8
h = 0.2
a = 0
b = 1

# PROCEDIMIENTO
muestras = int((b -a)/h)+1
tabla = np.zeros(shape=(muestras,4),dtype=float)

i = 0
xi = x0
yi = y0
tabla[i,:] = [xi,yi,0,0]

i = i+1
while not(i>=muestras):
    K1 = h*d1y(xi,yi)
    K2 = h*d1y(xi+h,yi+K1)
    yi = yi + (K1+K2)/2
    xi = xi +h
    tabla[i,:] = [xi,yi,K1,K2]
    i = i+1
# vector para gráfica
xg = tabla[:,0]
yg = tabla[:,1]

# SALIDA
# muestra 4 decimales
np.set_printoptions(precision=4)
print(' [xi, yi, K1, K2]')
print(tabla)
# Gráfica
plt.plot(xg,yg)
plt.xlabel('ti')
plt.ylabel('yi')
plt.grid()
plt.show()

Tarea: Realizar iteraciones para Runge-Kutta de 4to Orden

Ejercicio: 3Eva_IIT2007_T3 EDO Taylor orden 2

La ecuación del problema es:

y'= 1 +\frac{y}{t} + \Big(\frac{y}{t}\Big) ^2 1\leq t\leq 2 y(1)=0, h=0.2

Tomando como referencia Taylor de 3 términos más el término de error O(h³)

y_{i+1} = y_i +\frac{h}{1!}y'_i + \frac{h^2}{2!}y''_i + \frac{h^3}{3!}y'''_i

Se usa hasta el tercer término para el algoritmo.

y_{i+1} = y_i +\frac{h}{1!}y'_i + \frac{h^2}{2!}y''_i

Se determina que se requiere la segunda derivada para completar la aproximación. A partir de la ecuación del problema se aplica en cada término:

\Big(\frac{u}{v}\Big)' = \frac{u'v-uv' }{v^2} y'= 1 +\frac{y}{t} + \frac{y^2}{t^2} y''= 0+\frac{y't-y}{t^2} + \frac{2yy't^2-2y^2t}{t^4} y''= \frac{y't-y}{t^2} + \frac{2y\Big(1 +\frac{y}{t} + \frac{y^2}{t^2}\Big) t^2-2y^2t}{t^4}

Con lo que se realizan las iteraciones para llenar la tabla

iteración 1

t₀ = 1

y₀= 0

y'= 1 +\frac{0}{1} + \Big(\frac{0}{1}\Big)^2 = 1 y''= \frac{(1)(1)-0}{(1)^2} + \frac{2(0)(1)(1)^2-2(0)^2(1)}{(1)^4} = 1 y_{1} = 0 +\frac{0.2}{1!}(1) + \frac{0.2^2}{2!}(1) = 0.22

t₁ = t₀ + h = 1 + 0.2 = 1.2

iteración 2

t₁ = 1.2

y₁= 0.22

y'= 1 +\frac{0.22}{1.2} + \Big(\frac{0.22}{1.2}\Big) ^2 = 1.21694444 y''= \frac{y'(1.2)-0.22}{1.2^2} + \frac{2(0.22)y'(1.2)^2-2(0.22)^2(1.2)}{(1.2)^4} = 1.17716821 y_{2} = 0.22 +\frac{0.2}{1!}(1.21694444) + \frac{0.2^2}{2!}(1.17716821) = 0.48693225

t₂ = t₁ + h = 1.2 + 0.2 = 1.4

iteración 3

t₂ = 1.4

y₂= 0.48693225

y'= 1 +\frac{0.48693225}{1.4} + \Big(\frac{0.48693225}{1.4}\Big) ^2 = 1.46877968 y''= \frac{(1.46877968)(1.4)-0.48693225}{1.4^2} + \frac{2(0.48693225)(1.46877968)(1.4)^2-2(0.48693225)^2(1.4)}{1.4^4} = 1.35766995 y_{3} = 0.48693225 +\frac{0.2}{1!}(1.46877968) + \frac{0.2^2}{2!}(1.35766995) = 0.80784159

t₃ = t₂ + h = 1.4 + 0.2 = 1.6

Con lo que se puede realizar el algoritmo, obteniendo la siguiente tabla:

estimado
 [xi,        yi,        d1yi,      d2yi]
[[1.         0.         1.         1.        ]
 [1.2        0.22       1.21694444 1.17716821]
 [1.4        0.48693225 1.46877968 1.35766995]
 [1.6        0.80784159 1.759826   1.57634424]
 [1.8        1.19133367 2.09990017 1.8564435 ]
 [2.         1.64844258 0.         0.        ]]

Ejercicio: 3Eva_IIT2007_T1 EDP Eliptica, problema de frontera

Con los datos del ejercicio se plantean de la siguiente forma en los bordes:

La ecuación a resolver es:

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

Que en su forma discreta, con diferencias divididas centradas es:

\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} = 4

Al agrupar constantes se convierte en:

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

Siendo Δx= 1/3 y Δy =2/3, se mantendrá la relación λ = (Δx/Δy)²

u[i-1,j]-2u[i,j]+u[i+1,j] + + \lambda u[i,j-1]-2\lambda u[i,j]+\lambda u[i,j+1] = = 4(\Delta x)^2

---

u[i-1,j]-2(1+\lambda)u[i,j]+u[i+1,j] + + \lambda u[i,j-1]+\lambda u[i,j+1] = 4(\Delta x)^2

---

Se pueden realizar las iteraciones para los nodos y reemplaza los valores de frontera:

i=1 y j=1

u[0,1]-2(1+\lambda)u[1,1]+u[2,1] + + \lambda u[1,0]+\lambda u[1,2] = 4(\Delta x)^2 \Delta y^2-2(1+\lambda)u[1,1]+u[2,1] + + \Delta x^2+\lambda u[1,2] = 4(\Delta x)^2 -2(1+\lambda)u[1,1]+u[2,1] +\lambda u[1,2] = 4(\Delta x)^2 - \Delta y^2 - \Delta x^2 -2(1+0.25)u[1,1]+u[2,1] +0.25 u[1,2] = 4(1/3)^2 - (2/3)^2- (1/3)^2 -2.5u[1,1]+u[2,1] +0.25 u[1,2] = -0.1111

i=2 , j=1

u[1,1]-2(1+\lambda)u[2,1]+u[3,1] + + \lambda u[2,0]+\lambda u[2,2] = 4(\Delta x)^2 u[1,1]-2(1+\lambda)u[2,1]+(\delta y-1)^2 + + \Delta x^2 +\lambda u[2,2] = 4(\Delta x)^2 u[1,1]-2(1+\lambda)u[2,1] + \lambda u[2,2] = 4(\Delta x)^2 - (\Delta y-1)^2 - \Delta x^2 u[1,1]-2(1+0.25)u[2,1] + 0.25 u[2,2] = 4(1/3)^2 - (2/3-1)^2 - (1/3)^2 u[1,1]-2.5u[2,1] + 0.25 u[2,2] = 0.2222

i=1 , j=2

u[0,2]-2(1+\lambda)u[1,2]+u[2,2] + + \lambda u[1,1]+\lambda u[1,3] = 4(\Delta x)^2 (2\Delta y^2)-2(1+\lambda)u[1,2]+u[2,2] + + \lambda u[1,1]+\lambda (\Delta x-1)^2 = 4(\Delta x)^2 -2(1+\lambda)u[1,2]+u[2,2] + \lambda u[1,1] = 4(\Delta x)^2 - (2\Delta y^2) -\lambda (\Delta x-1)^2 -2(1+0.25)u[1,2]+u[2,2] + 0.25 u[1,1] = 4(1/3)^2 - (2(2/3)^2] -0.25 (1/3-1)^2 -2.5u[1,2]+u[2,2] + 0.25 u[1,1] = -0.5555

i=2, j=2

u[1,2]-2(1+\lambda)u[2,2]+u[3,2] + + \lambda u[2,1]+\lambda u[2,3] = 4(\Delta x)^2 u[1,2]-2(1+\lambda)u[2,2]+(\Delta y-1)^2 + + \lambda u[2,1]+\lambda (\Delta x-1)^2 = 4(\Delta x)^2 u[1,2]-2(1+\lambda)u[2,2] + \lambda u[2,1] =4(\Delta x)^2- (\Delta y-1)^2 -\lambda (\Delta x-1)^2 u[1,2]-2(1+0.25)u[2,2] + 0.25 u[2,1] =4(1/3)^2- (2/3-1)^2 -0.25(1/3-1)^2 u[1,2]-2.5u[2,2] + 0.25 u[2,1] =-0.2222

Obtiene el sistema de ecuaciones a resolver:

-2.5u[1,1]+u[2,1] +0.25 u[1,2] = -0.1111 u[1,1]-2.5u[2,1] + 0.25 u[2,2] = 0.2222 -2.5u[1,2]+u[2,2] + 0.25 u[1,1] = -0.5555 u[1,2]-2.5u[2,2] + 0.25 u[2,1] =-0.2222

Se escribe en forma matricial Ax=B

\begin {bmatrix} -2.5 && 1&&0.25&&0 \\1&&-2.5&&0&&0.25\\0.25&&0&&-2.5&&1\\0&&0.25&&1&&-2.5\end{bmatrix} \begin {bmatrix} u[1,1] \\ u[2,1]\\ u[1,2]\\u[2,2] \end{bmatrix} = \begin {bmatrix}-0.1111\\0.2222\\-0.5555\\-0.2222 \end{bmatrix}

que se resuelve con un algoritmo:

import numpy as np

A = np.array([[-2.5,1,0.25,0],
              [1,-2.5,0,0.25],
              [0.25,0,-2.5,1],
              [0,0.25,1,-2.5]])
B = np.array([-0.1111,0.2222,-0.5555,-0.2222])

x = np.linalg.solve(A,B)

print(x)

y los resultados son:

[ 0.05762549 -0.04492835  0.31156835  0.20901451]