s2Eva_IIT2019_T1 Canteras y urbanizaciones

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:

Icantera = 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:

Iu = Ius – Iui = 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()

s2Eva_IT2019_T3 EDP Elíptica Placa 6×5

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 hx =2, hy= 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)

s2Eva_IT2019_T2 Péndulo vertical

Ejercicio: 2Eva_IT2019_T2 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(h3)
Runge-Kutta 4do Orden tiene error de truncamiento O(h5)


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 Movimiento angular

s2Eva_IT2019_T1 Esfuerzo en pulso cardiaco

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

s2Eva_IIT2018_T3 EDP

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

s2Eva_IIT2018_T2 Kunge Kutta 2do Orden x»

Ejercicio: 2Eva_IIT2018_T2 Kunge Kutta 2do Orden x»

\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 (\pi t) = 0

Tarea: Desarrollar 3 iteraciones en Papel.

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

   t,              x,              z
[[ 0.          6.          1.5       ]
 [ 0.2         6.3         0.92679462]
 [ 0.4         6.38218195 -0.27187703]
 [ 0.6         6.19792527 -1.17287944]
 [ 0.8         5.88916155 -1.23638799]
 [ 1.          5.6491005  -0.61819399]
 [ 1.2         5.5872811   0.17288691]
 [ 1.4         5.69750883  0.69945284]
 [ 1.6         5.8992535   0.77223688]
 [ 1.8         6.09372469  0.43437943]
 [ 2.          6.20586248 -0.12630953]]

Instrucciones en 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

Ejercicio: 2Eva_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

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 t0=0, y0=0 y z0=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()

s2Eva_IT2018_T3 EDP Eliptica

Ejercicio: 2Eva_IT2018_T3 EDP Eliptica

Generar las ecuaciones a resolver usando diferencias finitas divididas centradas:

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

por facilidad se sustituye también en la forma discreta de la ecuación:

f[i,j] = f(x_i,y_j) = 2(x_{i} ^2+y_{j} ^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] = = \Delta y^2 f[i,j]
u[i-1,j]-4u[i,j]+u[i+1,j] + + u[i,j-1]+u[i,j+1] = \Delta y^2 f[i,j]
4u[i,j] = u[i-1,j]+u[i+1,j] + + u[i,j-1]+u[i,j+1]-\Delta y^2 f[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]-\Delta y^2 f[i,j] \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+hy,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
erroru = 2*tolera # para calcular al menos una matriz
while not(erroru=maxitera):
    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]-(hy**2)*fij[i,j])/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('iteraciones = ', itera)
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)
figura = plt.figure()
ax = Axes3D(figura)
U = np.transpose(u) # ajuste de índices fila es x
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

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