Ejercicio: 2Eva_IIT2019_T4 Integrar con Cuadratura de Gauss

f(x) = x ln(x)

1 ≤x≤4

se requiere:

I = \int_1^4 x ln(x) dx

literal a. Usando el método de Cuadratura de Gauss con 2 términos

x_a = \frac{b+a}{2} + \frac{b-a}{2}x_0 = \frac{4+1}{2} + \frac{4-1}{2}\Big(\frac{-1}{\sqrt{3}} \Big)

x_a =1.6339745962155612

x_b = \frac{b+a}{2} + \frac{b-a}{2}x_1 = \frac{4+1}{2} + \frac{4-1}{2}\Big(\frac{1}{\sqrt{3}} \Big)

x_b =3.366025403784439

I \cong \frac{b-a}{2}(f(x_a) + f(x_b)) I \cong \frac{4-1}{2}(x_a ln(x_a) + x_b ln(x_b))

I = 7.33164251999249

literal b.  De la fórmula , despejar el valor del error<0.0001

\Big|\frac{(b-a)}{180}h^4 f^{(4)} (\xi)\Big| <0.0001; \xi \in[a,b] h^4 <0.0001\frac{180}{(4-1)}\frac{1}{f^{(4)} (\xi)} h^4 < 0.006\frac{1}{f^{(4)} (\xi)} h <\Big(0.006\frac{1}{f^{(4)} (\xi)}\Big)^{1/4}

obteniendo la 4ta derivada de la función:

f(x) = x ln(x) f'(x) = ln(x) + x\Big(\frac{1}{x} \Big) = ln(x) +1 f''(x) = \frac{1}{x} f'''(x) = -\frac{1}{x^2} f^{(4)}(x) = 2\frac{1}{x^3}

se tiene que:

h <\Big(0.006\frac{1}{f^{(4)} (\xi)}\Big)^{1/4} h <\Big(0.006\frac{1}{2\frac{1}{\xi^3}}\Big)^{1/4} h <\Big(0.003\xi^3\Big)^{1/4} h <(0.003)^{1/4}\xi^{3/4}

en el peor de los casos, se toma el valor menor de ξ =1

h <(0.003)^{1/4} h<0.2340347319320716

 

Ejercicio: 2Eva_IIT2019_T3 EDP elíptica, placa en (1,1)

dada la ecuación del problema:

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

1 <  x < 2
1 <  y < 2

Se convierte a la versión discreta usando diferencias divididas centradas:

---

\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} = \frac{x_i}{y_j} + \frac{y_j}{x_i}

---

Se agrupan los términos Δx, Δy semejante a formar un λ al multiplicar todo por Δy²

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

---

los tamaños de paso en ambos ejes son de igual valor, se simplifica la ecuación

\lambda= \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] = =\Delta y^2\Big( \frac{x_i}{y_j} + \frac{y_j}{x_i}\Big)

---

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

---

Iteraciones

que permite plantear las ecuaciones para cada punto en posición [i,j]

---

i=1, j=1

u[0,1]-4u[1,1]+u[2,1] + + u[1,0]+u[1,2] =(0.25)^2\Big( \frac{0.25}{0.25} + \frac{0.25}{0.25}\Big) 0.25ln(0.25)-4u[1,1]+u[2,1] + + 0.25ln(0.25)+u[1,2] = 0.125 -4u[1,1]+u[2,1] +u[1,2] = 0.125 - 2(0.25)ln(0.25) -4u[1,1]+u[2,1] +u[1,2] = 0.8181

---

i=2, j=1

u[1,1]-4u[2,1]+u[3,1]+ +0.5 ln(0.5)+u[2,2]=0.15625 u[1,1]-4u[2,1]+u[3,1]+u[2,2]=0.15625-0.5 ln(0.5) u[1,1]-4u[2,1]+u[3,1]+u[2,2]=0.8493

---

i=3, j=1

u[2,1]-4u[3,1]+u[4,1] + + u[3,0]+u[3,2] =(0.25)^2\Big( \frac{0.75}{0.25} + \frac{0.25}{0.75}\Big) u[2,1]-4u[3,1]+u[4,1]+u[3,2] = 0.20833333-0.75 ln(0.75) u[2,1]-4u[3,1]+u[4,1]+u[3,2] =0.4240

---

Tarea: continuar con el ejercicio hasta plantear todo el sistema de ecuaciones.

A = np.array([
[-4, 1, 0, 1, 0, 0, 0, 0, 0],
[ 1,-4, 1, 0, 1, 0, 0, 0, 0],
[ 0, 1,-4, 0, 0, 1, 0, 0, 0],
[ 1, 0, 0,-4, 1, 0, 1, 0, 0],
[ 0, 1, 0, 1,-4, 1, 0, 1, 0],
[ 0, 0, 1, 0, 1,-4, 0, 0, 1],
[ 0, 0, 0, 1, 0, 0,-4, 1, 0],
[ 0, 0, 0, 0, 1, 0, 1,-4, 1],
[ 0, 0, 0, 0, 0, 1, 0, 1,-4]])
B = np.array(
[0.125 - 2(0.25)ln(0.25),
 0.15625 - 0.5ln(0.5),
 0.20833 - 0.75 ln(0.75),
 ...])

B = [0.8181,
     0.8493,
     0.4240,
     ...]

---

Algoritmo con Python

Con valores para la matriz solución:

iteraciones:  15
error entre iteraciones:  6.772297286980838e-05
solución para u: 
[[0.        0.2789294 0.6081976 0.9793276 1.3862943]
 [0.2789294 0.6978116 1.1792239 1.7127402 2.2907268]
 [0.6081976 1.1792239 1.8252746 2.5338403 3.2958368]
 [0.9793276 1.7127402 2.5338403 3.4280053 4.3846703]
 [1.3862943 2.2907268 3.2958368 4.3846703 5.5451774]]
>>>

y algoritmo detallado:

# 2Eva_IIT2019_T2 EDO, problema de valor inicial
# método iterativo
import numpy as np

# INGRESO
# longitud en x
a = 1
b = 2
# longitud en y
c = 1
d = 2
# tamaño de paso
dx = 0.25
dy = 0.25
# funciones en los bordes de la placa
abajo     = lambda x,y: x*np.log(x)
arriba    = lambda x,y: x*np.log(4*(x**2))
izquierda = lambda x,y: y*np.log(y)
derecha   = lambda x,y: 2*y*np.log(2*y)
# función de la ecuación
fxy = lambda x,y: x/y + y/x

# control de iteraciones
maxitera = 100
tolera = 0.0001

# PROCEDIMIENTO
# tamaño de la matriz
n = int((b-a)/dx)+1
m = int((d-c)/dy)+1
# vectores con valore de ejes
xi = np.linspace(a,b,n)
yj = np.linspace(c,d,m)
# matriz de puntos muestra
u = np.zeros(shape=(n,m),dtype=float)

# valores en los bordes
u[:,0]   = abajo(xi,yj[0])
u[:,m-1] = arriba(xi,yj[m-1])
u[0,:]   = izquierda(xi[0],yj)
u[n-1,:] = derecha(xi[n-1],yj)

# valores interiores la mitad en intervalo x,
# mitad en intervalo y, para menos iteraciones
mx = int(n/2)
my = int(m/2)
promedio = (u[mx,0]+u[mx,m-1]+u[0,my]+u[n-1,my])/4
u[1:n-1,1:m-1] = promedio

# método iterativo
itera = 0
converge = 0
while not(itera>=maxitera or converge==1):
    itera = itera +1
    # copia u para calcular errores entre iteraciones
    nueva = np.copy(u)
    for i in range(1,n-1):
        for j in range(1,m-1):
            # usar fórmula desarrollada para algoritmo
            fij = (dy**2)*fxy(xi[i],yj[j])
            u[i,j]=(u[i-1,j]+u[i+1,j]+u[i,j-1]+u[i,j+1]-fij)/4
    diferencia = nueva-u
    erroru = np.linalg.norm(np.abs(diferencia))
    if (erroru<tolera):
        converge=1

# SALIDA
print('iteraciones: ',itera)
print('error entre iteraciones: ',erroru)
print('solución para u: ')
print(u)

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

# matrices de ejes para la gráfica 3D
X, Y = np.meshgrid(xi, yj)
U = np.transpose(u) # ajuste de índices fila es x
figura = plt.figure()

grafica = Axes3D(figura)
grafica.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_IIT2019_T2 EDO, problema de valor inicial

la ecuación del problema planteado es:

y'(t) = f(t,y) = \frac{y}{2t^3}

0 ≤ t ≤ 1
y(0.5) = 1.5

---

literal a

La solución empieza usando la Serie de Taylor por ejemplo para tres términos:

y_{i+1} = y_{i} + h y'_i + \frac{h^2}{2!} y''_i x_{i+1} = x_{i} + h E = \frac{h^3}{3!} y'''(z) = O(h^3)

Se observa que se tiene el valor inicial y la primera derivada, si usamos tres términos se puede usar la segunda derivada.

y'(t) = f(t,y) = \frac{y}{2t^3} y''(t) = f'(t,y) = \frac{y'}{2t^3} + y \Big(\frac{-3}{2t^4}\Big) y''(t) = \frac{1}{2t^3}\Big( \frac{y}{2t^3} \Big) - \frac{3y}{2t^4} y''(t) = \frac{y}{4t^6} -\frac{3y}{2t^4}

por lo que la ecuación de Taylor a usar queda de la siguiente forma:

y_{i+1} = y_{i} + h y'_i + \frac{h^2}{2!} y''_i y_{i+1} = y_{i} + h\frac{y}{2t^3}+\frac{h^2}{2} \Big( \frac{y}{4t^6} -\frac{3y}{2t^4}\Big)

que es la ecuacion que se usará con un error de O(h³)

Reemplazando los valores en la fórmula se obtiene la siguiente tabla:

estimado
 [ti,      yi,      d1yi,    d2yi]
[[  0.5    1.5      6.      -12.        ]
 [  0.6    2.04     4.7222  -12.68004115]
 [  0.7    2.4488   3.5697  -10.09510219]
 [  0.8    2.7553   2.6907  -7.46259891]
 [  0.9    2.9870   2.0487  -5.42399041]
 [  1.     3.1648   0.       0.        ]]

cuya gráfica es:

---

literal b

Para desarrollar Runge-Kutta de 2do orden se dispone de los siguientes datos:

y'(t) = f(t,y) = \frac{y}{2t^3}

t₀ = 0.5, y₀ = 1.5, h = 0.1

pasos del algoritmo,

K_1 = h * y'(t_i) K_2 = h * y'(t_i+h, y_i + K_1) y_{i+1} = y_i + \frac{K_1+K_2}{2} t_i = t_i + h</p> <p>

iteración 1: i = 0

K_1 = 0.1 * y'(0.5) = 0.6 K_2 = 0.1 * y'(0.5+0.1, 1.5 + 0.6) = 0.4861 y_{1} = 1.5 + \frac{0.6+0.4861}{2} = 2.0430 t_1 = 0.5 + 0.1 = 0.6</p> <p>

iteración 2: i = 1

K_1 = 0.1 * y'(0.6) = 0.4729 K_2 = 0.1 * y'(0.6+0.1, 2.0430 + 0.4729) = 0.3667 y_{1} = 2.0430 + \frac{0.4729+0.3667}{2} = 2.4629 t_1 = 0.6 + 0.1 = 0.7</p> <p>

iteración 3: i = 2

K_1 = 0.1 * y'(0.7) = 0.3590 K_2 = 0.1 * y'(0.7+0.1, 2.4629 + 0.3590) = 0.3667 y_{1} = 2.4629 + \frac{0.3590+0.3667}{2} = 2.7802 t_1 = 0.7 + 0.1 = 0.8</p> <p>

obteniendo la siguiente tabla:

estimado
 [ti,   yi,     K1,     K2]
[[0.5   1.5     0.      0.    ]
 [0.6   2.0430  0.6     0.4861]
 [0.7   2.4629  0.4729  0.3667]
 [0.8   2.7802  0.3590  0.2755]
 [0.9   3.0206  0.2715  0.2093]
 [1.    3.2048  0.2071  0.1613]]
Diferencias entre Taylor y Runge-Kutta2
[ 0.   -0.0030 -0.0140 -0.0248 -0.0335 -0.0400]

---

Algoritmo en Python

# EDO. Método de Taylor 3 términos
# Runge-Kutta de 2 Orden
# 2Eva_IIT2019_T2 EDO, problema de valor inicial
import numpy as np

# Funciones desarrolladas en clase
def edo_taylor3t(d1y,d2y,x0,y0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,4),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,0,0]
    x = x0
    y = y0
    for i in range(1,tamano,1):
        estimado[i-1,2:]= [d1y(x,y),d2y(x,y)]
        y = y + h*d1y(x,y) + ((h**2)/2)*d2y(x,y)
        x = x+h
        estimado[i,0:2] = [x,y]
    return(estimado)
def rungekutta2(d1y,x0,y0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,4),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,0,0]
    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,K1,K2]
    return(estimado)

# PROGRAMA PRUEBA
# INGRESO
# d1y = y', d2y = y''
d1y = lambda t,y: y/(2*t**3)
d2y = lambda t,y: y/(4*t**6)-(3/2)*y/(t**4)
t0 = 0.5
y0 = 1.5
h = 0.1
muestras = 5

# PROCEDIMIENTO
# Edo con Taylor
puntos = edo_taylor3t(d1y,d2y,t0,y0,h,muestras)
ti = puntos[:,0]
yi = puntos[:,1]

# Runge-Kutta
puntosRK2 = rungekutta2(d1y,t0,y0,h,muestras)
# ti = puntosRK2[:,0] # lo mismo del anterior
yiRK2 = puntosRK2[:,1]

# diferencias
diferencia = yi-yiRK2

# SALIDA
print('estimado[ti, yi, d1yi, d2yi]')
print(puntos)

print('estimado[ti, yi, K1, K2]')
print(puntosRK2)

print('Diferencias entre Taylor y Runge-Kutta2')
print(diferencia)

# Gráfica
import matplotlib.pyplot as plt
plt.plot(ti[0],yi[0],'o', color='r',
         label ='[t0,y0]')
plt.plot(ti[1:],yi[1:],'o', color='g',
         label ='y con Taylor 3 términos')
plt.plot(ti[1:],yiRK2[1:],'o', color='m',
         label ='y Runge-Kutta 2 Orden')
plt.title('EDO: Solución numérica')
plt.xlabel('t')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show())

Ejercicio: 2Eva_IIT2019_T1 Canteras y urbanizaciones

Literal a. área de cantera

Canteras- frontera superior

xi	55	85	195	305	390	780	1170
f(xi)	752	825	886	1130	1086	1391	1219

Para proceder se calculan los tamaños de paso, h, en cada intervalo:

dx - tamaños de paso

dxi

30

110

110

85

390

390

___

Ics = \frac{30}{2}(752+825) + \frac{110}{2}(825+886) + + \frac{110}{2}(1130+886) + \frac{85}{2}(1130+1086) + \frac{390}{2}(1086+1391) +\frac{390}{2}(1391+1219)

que tiene como resultado: Ics = 1342435.0

Canteras- frontera inferior

xi	55	...	705	705	850	850	1010	1170
f(xi)	260	...	260	550	741	855	855	1055

Para proceder se calculan los tamaños de paso, h, en cada intervalo:

dx - tamaños de paso

dxi

650

...

0.0

145

0.0

160

160

____

Se observa que existen rectángulos en los intervalos, por lo que se simplifica la fórmula.

Ici = (650)(260) + \frac{145}{2}(741+550) + + (160)(855) + \frac{160}{2}(1055+855)

cuyo resultado es: Ici =552197.5

El área correspondiente a la cantera es:

I_cantera = Ics -Ici =1342435.0 - 552197.5 = 790237.5

---

Literal b. área de urbanización
la frontera inferior está referenciada a la eje x con g(x)=0, por lo que solo es necesario realizar el integral para la frontera superior. El valor de la integral de la frontera inferior de la urbanización es cero.

Urbanización - frontera superior

xi	720	800	890	890	1170	1220
g(xi)	527	630	630	760	760	533

dx - tamaños de paso

dxi

80

90

0.0

280

50

____

Ius = \frac{80}{2}(527+630) + (90)(630) + + (280)(760) + \frac{50}{2}(760+533)

El valor del área de la urbanización es:

I_u = I_us - I_ui = 348105.0 - 0 = 348105.0

---

literal c

Se pude mejora la precisión para los intervalos donde el tamaño de paso es igual, sin necesidad de aumentar o quitar puntos.

Observando los tramaños de paso en cada sección se sugiere usar el método de Simpson de 1/3 donde existen dos tamaños de paso iguales y de forma consecutiva.

Cantera - frontera superior: en el intervalo xi= [85,195,305] donde h es= 110

Cantera - frontera inferior: en el intervalo xi = [850,110,1170] donde h es= 160

---

Algoritmo con Python

Para trapecios en todos los intervalos. Considera que si es un rectángulo, la fórmula del trapecio también funciona.

# 2Eva_IIT2019_T1 Canteras y urbanizaciones
import numpy as np
import matplotlib.pyplot as plt

# Funciones para integrar realizadas en clase
def itrapecio (xi,fi):
    n=len(fi)
    integral=0
    for i in range(0,n-1,1):
        h = xi[i+1]-xi[i]
        darea = (h/2)*(fi[i]+fi[i+1])
        integral = integral + darea 
    return(integral)

# INGRESO
# Canteras - frontera superior
xcs = [  55.,  85, 195,  305,  390,  780, 1170]
ycs = [ 752., 825, 886, 1130, 1086, 1391, 1219]
# Canteras - frontera inferior
xci = [ 55., 705, 705, 850, 850, 1010, 1170]
yci = [260., 260, 550, 741, 855,  855, 1055]

# Urbanización - frontera superior
xus = [720., 800, 890, 890, 1170, 1220]
yus = [527., 630, 630, 760,  760,  533]
# Urbanización - frontera inferior
xui = [720., 1220]
yui = [  0.,    0]

# PROCEDIMIENTO

# Area de cantera
Ics = itrapecio(xcs,ycs)
Ici = itrapecio(xci,yci)
Icantera = Ics-Ici

# Area de urbanización
Iurb = itrapecio(xus,yus)

# SALIDA
print('Area canteras: ',Icantera)
print('Area urbanización: ', Iurb)

# Gráfica canteras
plt.plot(xcs,ycs,color='brown')
plt.plot(xci,yci,color='brown')
plt.plot([xci[0],xcs[0]],[yci[0],ycs[0]],color='brown')
plt.plot([xci[-1],xcs[-1]],[yci[-1],ycs[-1]],color='brown')

# Gráfica urbanizaciones
plt.plot(xus,yus, color='green')
plt.plot(xui,yui, color='green')
plt.plot([xui[0],xus[0]],[yui[0],yus[0]], color='green')
plt.plot([xui[-1],xus[-1]],[yui[-1],yus[-1]], color='green')
plt.show()

Ejercicio: 2Eva_IT2019_T3 EDP Elíptica Placa 6×5

La ecuación se discretiza con diferencias divididas centradas

\frac{\partial^2 u}{\partial x^2} (x,y)+\frac{\partial ^2 u}{\partial y^2 } (x,y) = -\frac{q}{K} \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}= -\frac{q}{K}

Para procesar los datos se toma como referencia la malla

se agrupan los términos conocidos de las diferencias divididas

\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}=-\frac{q}{K}(\Delta y)^2

Considerando que h_x =2, h_y= 2.5 son diferentes, se mantiene el valor de λ para el desarrollo.

\lambda = \frac{(\Delta y)^2}{(\Delta x)^2} = \frac{(2.5)^2}{(2)^2} = 1.5625 \lambda u_{i+1,j}-2\lambda u_{i,j} +\lambda u_{i-1,j} + u_{i,j+1}-2u_{i,j}+u_{i,j-1}=-\frac{q}{K}(\Delta y)^2

Se agrupan términos de u que son iguales

\lambda u_{i+1,j}-2(1+\lambda) u_{i,j} +\lambda u_{i-1,j} + u_{i,j+1}+u_{i,j-1}=-\frac{q}{K}(\Delta y)^2

con lo que se puede generar el sistema de ecuaciones a resolver para los nodos de la malla.

i = 1 , j=1

1.5625 u_{2,1}-2(1+ 1.5625) u_{1,1} + 1.5625 u_{0,1} + +u_{1,2}+u_{1,0}=-\frac{1.5}{1.04}(2.5)^2 1.5625 u_{2,1}-5.125u_{1,1} + 1.5625 y_1(5-y_1) + + 0+x_i(6-x_i)=-9.0144 1.5625 u_{2,1}-5.125u_{1,1} + 1.5625[ 1(5-1)]+ + 0+2(6-2)=-9.0144 1.5625 u_{2,1}-5.125 u_{1,1} = =-9.0144 - 1.5625 [1(5-1)] - 0-2(6-2)

 

1.5625 u_{2,1}-5.125 u_{1,1} =-23.2644

i=2, j=1

1.5625 u_{3,1}-2(1+1.5625) u_{2,1} +1.5625 u_{1,1} + +u_{2,2}+u_{2,0}=-9.0144 1.5625 (0)-5.125 u_{2,1} +1.5625 u_{1,1} + 0 +x_2(6-x_2)=-9.0144 -5.125 u_{2,1} +1.5625 u_{1,1} + 0 +4(6-4)=-9.0144

 

-5.125 u_{2,1} +1.5625 u_{1,1} =-17.0144

las ecuaciones a resolver son:

1.5625 u_{2,1}-5.125 u_{1,1} =-23.2644 -5.125 u_{2,1} +1.5625 u_{1,1} =-17.0144

cuya solución es:

u_{2,1} = 5.1858 u_{1,1} =6.1204

Resuelto usando:

import numpy as np
A = np.array([[1.5625,-5.125],
             [-5.125, 1.5625]])
B = np.array([-23.2644,-17.0144])
x = np.linalg.solve(A,B)
print(x)

Ejercicio: 2Eva_IT2019_T2 EDO Péndulo vertical

Para resolver la ecuación usando Runge-Kutta se estandariza la ecuación a la forma:

\frac{d^2\theta }{dt^2}+\frac{g}{L}\sin (\theta)=0 \frac{d^2\theta }{dt^2}= -\frac{g}{L}\sin (\theta)

Se simplifica su forma a:

\frac{d\theta}{dt}=z = f_t(t,\theta,z) \frac{d^2\theta }{dt^2}= z' =-\frac{g}{L}\sin (\theta) = g_t(t,\theta,z)

y se usan los valores iniciales para iniciar la tabla:

\theta(0) = \frac{\pi}{6} \theta '(0) = 0

ti	θ(ti)	θ'(ti)=z
0	π/6 = 0.5235	0
0.2	0.3602	-1.6333
0.4	-0.0815	-2.2639
...

para las iteraciones se usan los valores (-g/L) = (-9.8/0.6)

Iteración 1:  ti = 0 ; yi = π/6 ; zi = 0

K1y = h * ft(ti,yi,zi) 
    = 0.2*(0) = 0
K1z = h * gt(ti,yi,zi) 
    = 0.2*(-9.8/0.6)*sin(π/6) = -1.6333
        
K2y = h * ft(ti+h, yi + K1y, zi + K1z)
    = 0.2*(0-1.6333)= -0.3266
K2z = h * gt(ti+h, yi + K1y, zi + K1z)
    = 0.2*(-9.8/0.6)*sin(π/6+0) = -1.6333

yi = yi + (K1y+K2y)/2 
   = π/6+ (0+(-0.3266))/2 = 0.3602
zi = zi + (K1z+K2z)/2 
   = 0+(-1.6333-1.6333)/2 = -1.6333
ti = ti + h = 0 + 0.2 = 0.2

estimado[i] = [0.2,0.3602,-1.6333]

Iteración 2: ti = 0.2 ; yi = 0.3602 ; zi = -1.6333

K1y = 0.2*( -1.6333) = -0.3266
K1z = 0.2*(-9.8/0.6)*sin(0.3602) = -1.1515
        
K2y = 0.2*(-1.6333-0.3266)= -0.5569
K2z = 0.2*(-9.8/0.6)*sin(0.3602-0.3266) = -0.1097

yi = 0.3602 + ( -0.3266 + (-0.3919))/2 = -0.0815
zi = -1.6333+(-1.151-0.1097)/2 = -2.2639
ti = ti + h = 0.2 + 0.2 = 0.4

estimado[i] = [0.4,-0.0815,-2.2639]

Se continúan con las iteraciones, hasta completar la tabla.

Tarea: realizar la Iteración 3

Usando el algoritmo RungeKutta_fg se obtienen los valores y luego la gráfica

 [ t, 		 y, 	 dyi/dti=z]
[[ 0.          0.52359878  0.        ]
 [ 0.2         0.36026544 -1.63333333]
 [ 0.4        -0.08155862 -2.263988  ]
 [ 0.6        -0.50774327 -1.2990876 ]
 [ 0.8        -0.60873334  0.62920692]
 [ 1.         -0.29609456  2.32161986]]

si se mejora la resolución disminuyendo h = 0.05 y muestras = 20 para cubrir el dominio [0,1] se obtiene el siguiente resultado:

Tarea: Para el literal b, se debe considerar que los errores se calculan con

Runge-Kutta 2do Orden tiene error de truncamiento O(h³)
Runge-Kutta 4do Orden tiene error de truncamiento O(h⁵)

---

Algoritmo en Python

# 3Eva_IT2019_T2 Péndulo vertical
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
g = 9.8
L  = 0.6
ft = lambda t,y,z: z
gt = lambda t,y,z: (-g/L)*np.sin(y)

t0 = 0
y0 = np.pi/6
z0 = 0
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()

---

Otra solución propuesta puede seguir el problema semejante:

s2Eva_IT2010_T2 EDO Movimiento angular

Ejercicio: 2Eva_IT2019_T1 Esfuerzo en pulso cardiaco

Para resolver el ejercicio, la función a integrar es el cuadrado de los valores. Para minimizar los errores se usarán TODOS los puntos muestreados, aplicando los métodos adecuados.
Con aproximación de Simpson se requiere que los tamaños de paso sean iguales en cada segmento.
Por lo que primero se revisa el tamaño de paso entre lecturas.

tamaño de paso h:
[0.04 0.04 0.02 0.01 0.01 0.01 0.03 0.04 0.03 0.02 0.  ]
tiempos:
[0.   0.04 0.08 0.1  0.11 0.12 0.13 0.16 0.2  0.23 0.25]
ft:
[ 10.  18.   7.  -8. 110. -25.   9.   8.  25.   9.   9.]

Observando los tamaños de paso se tiene que:
- entre dos tamaños de paso iguales se usa Simpson de 1/3
- entre tres tamaños de paso iguales se usa Simpson de 3/8
- para tamaños de paso variables se usa trapecio.

Se procede a obtener el valor del integral,

Intervalo [0,0.8], h = 0.04

I_{S13} = \frac{0.04}{3}(10^2+4(18^2)+7^2)

Intervalo [0.08,0.1], h = 0.02

I_{Tr1} = \frac{0.02}{2}(7^2+(-8)^2)

Intervalo [0.1,0.13], h = 0.01

I_{S38} = \frac{3}{8}(0.01)((-8)^2+3(110)^2+3(-25)^2+9^2)

Intervalo [0.13,0.25], h = variable

I_{Tr2} = \frac{0.03}{2}(9^2+8^2) I_{Tr3} = \frac{0.04}{2}(8^2+25^2) I_{Tr4} = \frac{0.03}{2}(25^2+9^2) I_{Tr5} = \frac{0.02}{2}(9^2+9^2)

El integral es la suma de los valores parciales, y con el resultado se obtiene el valor X_rms requerido.

I_{total} = \frac{1}{0.08-0}I_{S13}+\frac{1}{0.1-0.08}I_{Tr1} +\frac{1}{0.13-0.1}I_{S38} + \frac{1}{0.16-0.13}I_{Tr2} + \frac{1}{0.2-0.16}I_{Tr3} +\frac{1}{0.23-0.2}I_{Tr4} + \frac{1}{0.25-0.23}I_{Tr5} X_{rms} = \sqrt{I_{total}}

Los valores resultantes son:

Is13:  19.26666666666667
ITr1:  1.1300000000000001
Is38:  143.7
ITr2:  2.175
ITr3:  13.780000000000001
ITr4:  10.59
ITr5:  1.62
Itotal:  5938.333333333333
Xrms:  77.06058222809722

Tarea: literal b

Para Simpson 1/3

error_{trunca} = -\frac{h^5}{90} f^{(4)}(z)

Para Simpson 3/8

error_{truncamiento} = -\frac{3}{80} h^5 f^{(4)} (z)

Para trapecios

error_{truncar} = -\frac{h^3}{12}f''(z)

---

---

Algoritmo en Python

# 3Eva_IT2019_T1 Esfuerzo Cardiaco
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
t  = np.array([0.0,0.04,0.08,0.1,0.11,0.12,0.13,0.16,0.20,0.23,0.25])
ft = np.array([10., 18, 7, -8, 110, -25, 9, 8, 25, 9, 9])

# PROCEDIMIENTO
# Revisar tamaño de paso h
n = len(t)
dt = np.zeros(n, dtype=float)
for i in range(0,n-1,1):
    dt[i]=t[i+1]-t[i]

# Integrales
Is13 = (0.04/3)*((10)**2 + 4*((18)**2) + (7)**2)
ITr1 = (0.02/2)*((7)**2 + (-8)**2)
Is38 = (3/8)*0.01*((-8)**2 + 3*((110)**2) + 3*((-25)**2) + (9)**2)

ITr2 = (0.03/2)*((9)**2 + (8)**2)
ITr3 = (0.04/2)*((8)**2 + (25)**2)
ITr4 = (0.03/2)*((25)**2 + (9)**2)
ITr5 = (0.02/2)*((9)**2 + (9)**2)

Itotal = (1/(0.08-0.0))*Is13 + (1/(0.1-0.08))*ITr1
Itotal = Itotal + (1/(0.13-0.1))*Is38 + (1/(0.16-0.13))*ITr2
Itotal = Itotal + (1/(0.20-0.16))*ITr3 + (1/(0.23-0.20))*ITr4
Itotal = Itotal + (1/(0.25-0.23))*ITr5
Xrms = np.sqrt(Itotal)

# SALIDA
print('tamaño de paso h:')
print(dt)
print('tiempos:')
print(t)
print('ft: ')
print(ft)
print('Is13: ', Is13)
print('ITr1: ', ITr1)
print('Is38: ', Is38)
print('ITr2: ', ITr2)
print('ITr3: ', ITr3)
print('ITr4: ', ITr4)
print('ITr5: ', ITr5)
print('Itotal: ', Itotal)
print('Xrms: ', Xrms)

# Grafica
plt.plot(t,ft)
for i in range(1,n,1):
    plt.axvline(t[i], color='green', linestyle='dashed')
plt.xlabel('tiempo s')
plt.ylabel('valor sensor')
plt.title('Un pulso cardiaco con sensor')
plt.show()

Ejercicio: 2Eva_IIT2018_T3 EDP

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