Categoría: Sol_2Eva 2007-2010

Ejercicio: 2Eva_IT2010_T2 Movimiento angular

Para resolver, se usa Runge-Kutta_fg de segundo orden como ejemplo

$y'' + 10 \sin (y) =0$

se hace

$y' = z = f(t,y,z)$

y se estandariza:

$y'' =z'= -10 \sin (y) = g(t,y,z)$

teniendo como punto de partida t₀=0, y₀=0 y z₀=0.1

$y(0)=0, y'(0)=0.1$

Se desarrolla el algotitmo para obtener los valores:

 [ t, 		 y, 	 dyi/dti=z]
[[ 0.          0.          0.1       ]
 [ 0.2         0.02        0.08000133]
 [ 0.4         0.03200053  0.02401018]
 [ 0.6         0.03040355 -0.04477916]
 [ 0.8         0.01536795 -0.09662411]
 [ 1.         -0.00703034 -0.10803459]]

que permiten generar la gráfica de respuesta:

---

Algoritmo en Python

# 2Eva_IT2010_T2 Movimiento angular
import numpy as np
import matplotlib.pyplot as plt

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,z0]
    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)

# INGRESO theta = y
ft = lambda t,y,z: z
gt = lambda t,y,z: -10*np.sin(y)

t0 = 0
y0 = 0
z0 = 0.1
h=0.2
muestras = 5

# PROCEDIMIENTO
tabla = rungekutta2_fg(ft,gt,t0,y0,z0,h,muestras)

# SALIDA
print(' [ t, \t\t y, \t dyi/dti=z]')
print(tabla)

# Grafica
ti = np.copy(tabla[:,0])
yi = np.copy(tabla[:,1])
zi = np.copy(tabla[:,2])
plt.subplot(121)
plt.plot(ti,yi)
plt.xlabel('ti')
plt.title('yi')
plt.subplot(122)
plt.plot(ti,zi, color='green')
plt.xlabel('ti')
plt.title('dyi/dti')
plt.show()

Ejercicio: 2Eva_IT2010_T3 EDP elíptica, Placa no rectangular

la fórmula a resolver, siendo f(x,y)=20

$\frac{\delta^2 u}{\delta x^2} + \frac{\delta^2 u}{\delta y^2} = f$ $\frac{\delta^2 u}{\delta x^2} + \frac{\delta^2 u}{\delta y^2} = 20$

Siendo las condiciones de frontera en los bordes marcados:

Se convierte a forma discreta la ecuación

$\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} = 20$ $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) = 20 \Delta x^2$

siendo λ = 1, al tener la malla en cuadrícula Δx=Δy

$u_{i+1,j}-4u_{i,j}+u_{i-1,j} + u_{i,j+1}+u_{i,j-1} = 20 \Delta x^2$

despejando u[i,j]

$u_{i,j} =\frac{1}{4}\Big(u_{i+1,j}+u_{i-1,j} + u_{i,j+1}+u_{i,j-1} - 20 \Delta x^2 \Big)$

Para el ejercicio, se debe usar las fronteras para los valores de cálculo que se puede hacer de dos formas:
1. Por posición del valor en la malla [i,j]
2. Por valor xi,yj que es la forma usada cuando la función es contínua.

1. Por posición del valor en la malla

Se busca la posición de la esquina derecha ubicada en xi = 0 y yj=1.5 y se la ubica como variable ‘desde’ como referencia para ubicar los valores de las fronteras usando los índices i, j.

# valores en fronteras
# busca poscición de esquina izquierda
desde = -1 
for j in range(0,m,1):
    if ((yj[j]-1.5)<tolera):
        desde = j
# llena los datos de bordes de placa
for j  in range(0,m,1):
    if (yj[j]>=1):
        u[j-desde,j] = Ta
    u[n-1,j] = Tb(xi[n-1],yj[j])
for i in range(0,n,1):
    if (xi[i]>=1):
        u[i,m-1]=Td
    u[i,desde-i] = Tc(xi[i],yj[i])
# valor inicial de iteración

2. Por valor x_i, y_j

Se establecen las ecuaciones que forman la frontera:

ymin = lambda x,y: 1.5-(1.5/1.5)*x
ymax = lambda x,y: 1.5+(1/1)*x

que se usan como condición para hacer el cálculo en cada nodo

if ((yj[j]-yjmin)>tolera and (yj[j]-yjmax)<-tolera):
                u[i,j] = (u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1]-20*(dx**2))/4

Para minimizar errores por redondeo o truncamiento al seleccionar los puntos, la referencia hacia cero se toma como «tolerancia»; en lugar de más que cero o menos que cero.

Con lo que se obtiene el resultado mostrado en la gráfica aumentando la resolución con Δx=Δy=0.05:

converge =  1
xi=
[0.   0.05 0.1  0.15 0.2  0.25 0.3  0.35 0.4  0.45 0.5  0.55 0.6  0.65
 0.7  0.75 0.8  0.85 0.9  0.95 1.   1.05 1.1  1.15 1.2  1.25 1.3  1.35
 1.4  1.45 1.5 ]
yj=
[0.   0.05 0.1  0.15 0.2  0.25 0.3  0.35 0.4  0.45 0.5  0.55 0.6  0.65
 0.7  0.75 0.8  0.85 0.9  0.95 1.   1.05 1.1  1.15 1.2  1.25 1.3  1.35
 1.4  1.45 1.5  1.55 1.6  1.65 1.7  1.75 1.8  1.85 1.9  1.95 2.   2.05
 2.1  2.15 2.2  2.25 2.3  2.35 2.4  2.45 2.5 ]
matriz u
[[  0.     0.     0.   ...   0.     0.     0.  ]
 [  0.     0.     0.   ...   0.     0.     0.  ]
 [  0.     0.     0.   ...   0.     0.     0.  ]
 ...
 [  0.     0.     6.67 ...  96.1   98.03 100.  ]
 [  0.     3.33   5.31 ...  96.03  98.   100.  ]
 [  0.     2.     4.   ...  96.    98.   100.  ]]
>>>

---

Algoritmo en Python

# Ecuaciones Diferenciales Parciales
# Elipticas. Método iterativo para placa NO rectangular
import numpy as np

# INGRESO
# Condiciones iniciales en los bordes
Ta = 100
Tb = lambda x,y:(100/2.5)*y
Tc = lambda x,y: 100-(100/1.5)*x
Td = 100
# dimensiones de la placa no rectangular
x0 = 0
xn = 1.5
y0 = 0
yn = 2.5
ymin = lambda x,y: 1.5-(1.5/1.5)*x
ymax = lambda x,y: 1.5+(1/1)*x
# discretiza, supone dx=dy
dx = 0.05 
dy = 0.05 
maxitera = 100
tolera = 0.0001

# PROCEDIMIENTO
xi = np.arange(x0,xn+dx,dx)
yj = np.arange(y0,yn+dy,dy)
n = len(xi)
m = len(yj)
# Matriz u
u = np.zeros(shape=(n,m),dtype = float)

# valores en fronteras
# busca posición de esquina izquierda
desde = -1 
for j in range(0,m,1):
    if ((yj[j]-1.5)<tolera): desde = j 

# llena los datos de bordes de placa 

for j in range(0,m,1): if (yj[j]>=1):
        u[j-desde,j] = Ta
    u[n-1,j] = Tb(xi[n-1],yj[j])
for i in range(0,n,1):
    if (xi[i]>=1):
        u[i,m-1]=Td
    u[i,desde-i] = Tc(xi[i],yj[i])
# valor inicial de iteración
# se usan los valores cero

# iterar puntos interiores
itera = 0
converge = 0
while not(itera>=maxitera and converge==1):
    itera = itera + 10000
    nueva = np.copy(u)
    for i in range(1,n-1):
        for j in range(1,m-1):
            yjmin = ymin(xi[i],yj[j])
            yjmax = ymax(xi[i],yj[j])
            if ((yj[j]-yjmin)>tolera and (yj[j]-yjmax)<-tolera):
                u[i,j] = (u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1]-20*(dx**2))/4
    diferencia = nueva-u
    erroru = np.linalg.norm(np.abs(diferencia))

    if (erroru<tolera):
        converge=1

# SALIDA
np.set_printoptions(precision=2)
print('converge = ', converge)
print('xi=')
print(xi)
print('yj=')
print(yj)
print('matriz u')
print(u)

Para incorporar la gráfica en forma de superficie.

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

Ejercicio: 2Eva_IT2008_T1_MN Producción petroleo

Literal a

Para el cálculo de las derivadas se hace uso de las fórmulas de diferenciación presentadas en la unidad 6, y basadas en el polinomio de Taylor:

$f'(x_i) = \frac{f(x_{i+1})-f(x_i)}{h} + O(h)$ $f''(x_i) = \frac{f(x_{i+2})-2f(x_{i+1})+f(x_i)}{h^2} + O(h)$

[ dia, prod, dprod, d2prod]
[[ 1.000e+00  3.345e+03 -1.000e+02  6.600e+01]
 [ 2.000e+00  3.245e+03 -3.400e+01  1.320e+02]
 [ 3.000e+00  3.211e+03  9.800e+01 -5.600e+01]
 [ 4.000e+00  3.309e+03  4.200e+01  1.900e+01]
 [ 5.000e+00  3.351e+03  6.100e+01 -2.430e+02]
 [ 6.000e+00  3.412e+03 -1.820e+02  8.700e+01]
 [ 7.000e+00  3.230e+03 -9.500e+01  9.200e+01]
 [ 8.000e+00  3.135e+03 -3.000e+00  0.000e+00]
 [ 9.000e+00  3.132e+03 -3.000e+00  0.000e+00]
 [ 1.000e+01  3.129e+03  0.000e+00  0.000e+00]]

representados en las siguientes gráfica:

literal b

Dado que las fórmlas de error usadas tienen error del orden h: O(h), el error de las fórmulas es del orden de:

h= dia[1]-dia[0] = 2-1 = 1

literal c

Para el día dos se observa un decrecimiento en la producción, tal como lo refleja el valor negativo de la primera derivada.
Sin embargo para el día siguiente, la producción no mantiene la tasa de decrecimiento, se observa la segunda derivada positiva, Empieza a «acelerar».

---

Las instrucciones en Python para la tabla presentada son:

# 2Eva_IT2008_T1_MN Producción petroleo
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
dia = np.array( [1., 2, 3, 4, 5, 6, 7, 8, 9, 10])
produccion = [3345., 3245, 3211, 3309, 3351, 3412, 3230, 3135, 3132, 3129]
produccion = np.array(produccion, dtype = float)

# PROCEDIMIENTO
n = len(dia)

# primera derivada
dp = np.zeros(n,dtype=float)
for i in range(0,n-1,1):
    dp[i] = (produccion[i+1]-produccion[i])/(dia[i+1]-dia[i])

# segunda derivada
d2p = np.zeros(n,dtype=float)
h = dia[1]-dia[0]
for i in range(0,n-2,1):
    d2p[i] = (produccion[i]-2*produccion[i+1]+ produccion[i+2])/(h**2)

tabla = np.concatenate(([dia],[produccion],[dp],[d2p]),axis=0)
tabla = np.transpose(tabla)

# SALIDA
print(" [ dia, prod, dprod, d2prod]")
print(tabla)

# gráfica
plt.subplot(121)
plt.plot(dia,produccion)
plt.xlabel('dia')
plt.ylabel('producción')
plt.grid()
plt.subplot(122)
plt.plot(dia,dp,color='green',label='dp')
plt.xlabel('dia')
plt.plot(dia,d2p,color='orange',label='d2p')
plt.axhline(0)
plt.legend()
plt.grid()
plt.show()

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