s2Eva_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

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

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_IT2012_T2 Modelo de clima

Se deja de tarea realizar las tres primeras iteraciones en papel.

Se presenta el resultado usando el algoritmo de segundo grado como una variante a la respuesta usada como ejemplo.

 [ ti, xi, yi, zi]
[[ 0.0000e+00  1.0000e+01  7.0000e+00  7.0000e+00]
 [ 2.5000e-03  9.9323e+00  7.5033e+00  7.1335e+00]
 [ 5.0000e-03  9.8786e+00  7.9988e+00  7.2774e+00]
 ...
 [ 2.4995e+01 -8.4276e+00 -2.7491e+00  3.3021e+01]
 [ 2.4998e+01 -8.2860e+00 -2.6392e+00  3.2858e+01]
 [ 2.5000e+01 -8.1453e+00 -2.5346e+00  3.2692e+01]]

Algoritmo en Python

 2Eva_IT2012_T2 Modelo de clima
# MATG1013 Análisis Numérico
import numpy as np
import matplotlib.pyplot as plt

def rungekutta2_fg(f,g,j,t0,x0,y0,z0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,4),dtype=float)
    # incluye el punto [t0,x0,y0,z0]
    estimado[0] = [t0,x0,y0,z0]
    ti = t0
    xi = x0
    yi = y0
    zi = z0
    for i in range(1,tamano,1):
        K1x = h * f(ti,xi,yi,zi)
        K1y = h * g(ti,xi,yi,zi)
        K1z = h * j(ti,xi,yi,zi)
        
        K2x = h * f(ti+h,xi + K1x, yi + K1y, zi + K1z)
        K2y = h * g(ti+h,xi + K1x, yi + K1y, zi + K1z)
        K2z = h * j(ti+h,xi + K1x, yi + K1y, zi + K1z)
        
        xi = xi + (K1x+K2x)/2
        yi = yi + (K1y+K2y)/2
        zi = zi + (K1z+K2z)/2
        ti = ti + h
        
        estimado[i] = [ti,xi,yi,zi]
    return(estimado)


#INGRESO
to = 0
xo = 10
yo = 7
zo = 7
f = lambda t,x,y,z: 10*(y-x)
g = lambda t,x,y,z: x*(28-z) - y
j = lambda t,x,y,z: -(8/3)*z + (x*y)
h = 0.0025
muestras = 10000

#PROCEDIMIENTO
#Rugen-Kutta 2_orden
tabla = rungekutta2_fg(f,g,j,to,xo,yo,zo,h,muestras)

#SALIDA
np.set_printoptions(precision=4)
print(' [ ti, xi, yi, zi]')
print(tabla)

# Gráfica
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(tabla[:,1], tabla[:,2], tabla[:,3])
plt.show()

s2Eva_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ónd el 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 xi,yj

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

s2Eva_IT2012_T3_MN EDO Taylor 2 Contaminación de estanque

La ecuación a resolver con Taylor es:

s'- \frac{26s}{200-t} - \frac{5}{2} = 0

Para lo que se plantea usar la primera derivada:

s'= \frac{26s}{200-t}+\frac{5}{2}

con valores iniciales de s(0) = 0, h=0.1

La fórmula de Taylor para trestérminos es:

s_{i+1}= s_{i}+s'_{i}h + \frac{s''_{i}}{2}h^2 + error

Para el desarrollo se compara la solución con dos términos, tres términos y Runge Kutta.

1. Solución con dos términos de Taylor

Iteraciones

i = 0, t0 = 0, s(0)=0

s'_{0}= \frac{26s_{0}}{200-t_{0}}+\frac{5}{2} = \frac{26(0)}{200-0}+\frac{5}{2} = \frac{5}{2} s_{1}= s_{0}+s'_{0}h = 0+ \frac{5}{2}*0.1= 0.25

t1 =  t0+h = 0+0.1 = 0.1

i=1


s'_{1}= \frac{26s_{1}}{200-t_{1}}+\frac{5}{2} = \frac{26(0.25)}{200-0.1}+\frac{5}{2} = 2.5325 s_{2}= s_{1}+s'_{1}h = 0.25 + (2.5325)*0.1 = 0.5032

t2 =  t1+h = 0.1+0.1 = 0.2

i=2,

resolver como tarea


2. Resolviendo con Python

estimado
 [xi,yi Taylor,yi Runge-Kutta, diferencias]
[[ 0.0  0.0000e+00  0.0000e+00  0.0000e+00]
 [ 0.1  2.5000e-01  2.5163e-01 -1.6258e-03]
 [ 0.2  5.0325e-01  5.0655e-01 -3.2957e-03]
 [ 0.3  7.5980e-01  7.6481e-01 -5.0106e-03]
 [ 0.4  1.0197e+00  1.0265e+00 -6.7714e-03]
 [ 0.5  1.2830e+00  1.2916e+00 -8.5792e-03]
 [ 0.6  1.5497e+00  1.5601e+00 -1.0435e-02]
 [ 0.7  1.8199e+00  1.8322e+00 -1.2339e-02]
 [ 0.8  2.0936e+00  2.1079e+00 -1.4294e-02]
 [ 0.9  2.3710e+00  2.3873e+00 -1.6299e-02]
 [ 1.0  2.6519e+00  2.6703e+00 -1.8357e-02]
 [ 1.1  2.9366e+00  2.9570e+00 -2.0467e-02]
 [ 1.2  3.2250e+00  3.2476e+00 -2.2632e-02]
 [ 1.3  3.5171e+00  3.5420e+00 -2.4853e-02]
 [ 1.4  3.8132e+00  3.8403e+00 -2.7129e-02]
 [ 1.5  4.1131e+00  4.1426e+00 -2.9464e-02]
 [ 1.6  4.4170e+00  4.4488e+00 -3.1857e-02]
 [ 1.7  4.7248e+00  4.7592e+00 -3.4310e-02]
 [ 1.8  5.0368e+00  5.0736e+00 -3.6825e-02]
 [ 1.9  5.3529e+00  5.3923e+00 -3.9402e-02]
 [ 2.0  5.6731e+00  5.7152e+00 -4.2043e-02]]
error en rango:  0.04204310894163932


2. Algoritmo en Python

# EDO. Método de Taylor 3 términos 
# estima la solucion para muestras espaciadas h en eje x
# valores iniciales x0,y0
# entrega arreglo [[x,y]]
import numpy as np

def edo_taylor2t(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]
    x = x0
    y = y0
    for i in range(1,tamano,1):
        y = y + h*d1y(x,y) # + ((h**2)/2)*d2y(x,y)
        x = x+h
        estimado[i] = [x,y]
    return(estimado)

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)

# PROGRAMA PRUEBA
# 2Eva_IIT2016_T3_MN EDO Taylor 2, Tanque de agua

# INGRESO.
# d1y = y' = f, d2y = y'' = f'
d1y = lambda x,y: 26*y/(200-x)+5/2
x0 = 0
y0 = 0
h = 0.1
muestras = 20

# PROCEDIMIENTO
puntos = edo_taylor2t(d1y,x0,y0,h,muestras)
xi = puntos[:,0]
yi = puntos[:,1]

# Con Runge Kutta
puntosRK2 = rungekutta2(d1y,x0,y0,h,muestras)
xiRK2 = puntosRK2[:,0]
yiRK2 = puntosRK2[:,1]

# diferencias
diferencias = yi-yiRK2
error = np.max(np.abs(diferencias))
tabla = np.copy(puntos)
tabla = np.concatenate((puntos,np.transpose([yiRK2]),
                        np.transpose([diferencias])),
                       axis = 1)

# SALIDA
np.set_printoptions(precision=4)
print('estimado[xi,yi Taylor,yi Runge-Kutta, diferencias]')
print(tabla)
print('error en rango: ', error)

# Gráfica
import matplotlib.pyplot as plt
plt.plot(xi[0],yi[0],'o',
         color='r', label ='[x0,y0]')
plt.plot(xi[0:],yi[0:],'-',
         color='g',
         label ='y Taylor 2 términos')
plt.plot(xiRK2[0:],yiRK2[0:],'-',
         color='blue',
         label ='y Runge-Kutta 2Orden')
plt.axhline(y0/2)
plt.title('EDO: Taylor 2T vs Runge=Kutta 2Orden')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()

s2Eva_IIT2016_T3_MN EDO Taylor 2, Tanque de agua

La solución obtenida se realiza con h=0.5 y usando dos métodos para comparar resultados.

\frac{dy}{dt} = -k \sqrt{y}

1. EDO con Taylor

Usando una aproximación con dos términos de Taylor:

y_{i+1}=y_{i}+ y'_{i} h+\frac{y"_{i}}{2}h^{2}

Por lo que se obtienen las derivadas necesarias:

y'_i= -k (y_i)^{1/2} y"_i= \frac{-k}{2}(y_i)^{-1/2}
1.1 iteraciones

i=0, y0=3, t0=0

y'_0= -k(y_0)^{1/2} =-0.06(3)^{1/2} = -0.1039 y"_0= \frac{-0.06}{2}(3)^{-1/2} = -0.0173 y_{1}=y_{0}+ y'_{0} (1)+\frac{y"_{0}}{2}(1)^{2} y_{1}=3+ (-0.1039) (0.5)+\frac{-0.0173}{2}(0.5)^{2}= 2.9458

t1=t0+h = 0+0.5= 0.5

i=1, y1=2.9458, t1=0.5

y'_1= -k(y_1)^{1/2} =-0.06(2.887)^{1/2} =-0.1029 y"_1= \frac{-0.06}{2}(2.887)^{-1/2} = -0.0174 y_{2}=y_{1}+ y'_{1} (1)+\frac{y"_{1}}{2}(1)^{2} y_{1}=2.9458+ (-0.1029) (1)+\frac{-0.0174}{2}(1)^{2}= 2.8921

t2=t1+h = 0.5+0.5 = 1.0

i=2, y2=2.8921, t2=1.0

Resolver como Tarea

1.2 Resultados con Python

Realizando una tabla de valores usando Python y una gráfica, encuentra que el valor buscado del tanque a la mitad se obtiene en 16 minutos.

estimado[xi,yi]
[[ 0.          3.        ]
 [ 0.5         2.94587341]
 [ 1.          2.89219791]
 [ 1.5         2.83897347]
 [ 2.          2.7862001 ]
 ...
 [14.          1.65488507]
 [14.5         1.61337731]
 [15.          1.57231935]
 [15.5         1.53171109]
 [16.          1.49155239]
 [16.5         1.45184313]
 [17.          1.41258317]
 [17.5         1.37377234]
 [18.          1.33541049]
 [18.5         1.29749744]
 [19.          1.26003297]
 [19.5         1.22301689]
 [20.          1.18644897]]

Algoritmo en Python para Solución EDO con tres términos:

# EDO. Método de Taylor 3 términos 
# estima la solucion para muestras espaciadas h en eje x
# valores iniciales x0,y0
# entrega arreglo [[x,y]]
import numpy as np

def edo_taylor3t(d1y,d2y,x0,y0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,2),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0]
    x = x0
    y = y0
    for i in range(1,tamano,1):
        y = y + h*d1y(x,y) + ((h**2)/2)*d2y(x,y)
        x = x+h
        estimado[i] = [x,y]
    return(estimado)

# PROGRAMA PRUEBA
# 2Eva_IIT2016_T3_MN EDO Taylor 2, Tanque de agua

# INGRESO.
k=0.06
# d1y = y' = f, d2y = y'' = f'
d1y = lambda x,y: -k*(y**0.5)
d2y = lambda x,y: -(k/2)*(y**(-0.5))
x0 = 0
y0 = 3
h = 1/2
muestras = 40

# PROCEDIMIENTO
puntos = edo_taylor3t(d1y,d2y,x0,y0,h,muestras)
xi = puntos[:,0]
yi = puntos[:,1]

# SALIDA
print('estimado[xi,yi]')
print(puntos)
# Gráfica
import matplotlib.pyplot as plt
plt.plot(xi[0],yi[0],'o', color='r', label ='[x0,y0]')
plt.plot(xi[1:],yi[1:],'o', color='g', label ='y estimada')
plt.axhline(y0/2)
plt.title('EDO: Solución con Taylor 3 términos')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()

2. EDO con Runge-Kutta de 2do Orden dy/dx

Para éste método no se requiere desarrollar la segunda derivada, se usa el mismo h =0.5 con fines de comparación de resultados

2.1 ITeraciones

i = 1, y0=3, t0=0

K_1 = h y'(x_0,y_0) = (0.5)*(-0.06)(3)^{1/2} =-0.05196 [\latex]

K_2 = h y'(x_0+h,y_0+K_1) = (0.5)* y'(0.5,3-0.05196) = -0.05150\latex]

y_1 = y_0+\frac{K1+K2}{2} = 3+\frac{-0.05196-0.05150}{2} = 2.9482\latex]

i = 2, y1=2.9482, t1=0.5

K_1 = h y'(x_1,y_1) = (0.5)*(-0.06)(2.9482)^{1/2} =-0.05149 [\latex]

K_2 = h y'(x_1+h,y_1+K_1) = (0.5)* y'(0.5,2.9482-0.05149) = -0.05103\latex]

y_1 = y_0+\frac{K1+K2}{2} = 3+\frac{-0.05149-0.05103}{2} = -2.8946\latex]

i = 3,  y1=2.8946, t1=1.0

Resolver como Tarea

2.2 Resultados con Python

Si comparamos con los resultados anteriores en una tabla, y obteniendo las diferencias entre cada iteración se tiene que:

estimado[xi,yi Taylor, yi Runge-Kutta, diferencias]
[[ 0.0  3.00000000  3.00000000  0.00000000e+00]
 [ 0.5  2.94587341  2.94826446 -2.39104632e-03]
 [ 1.0  2.89219791  2.89697892 -4.78100868e-03]
 [ 1.5  2.83897347  2.84614338 -7.16990106e-03]
 [ 2.0  2.78620010  2.79575783 -9.55773860e-03]
...
 [ 14.0  1.65488507  1.72150488 -6.66198112e-02]
 [ 14.5  1.61337731  1.68236934 -6.89920328e-02]
 [ 15.0  1.57231935  1.64368380 -7.13644510e-02]
 [ 15.5  1.53171109  1.60544826 -7.37371784e-02]
 [ 16.0  1.49155239  1.56766273 -7.61103370e-02]
 [ 16.5  1.45184313  1.53032719 -7.84840585e-02]
 [ 17.0  1.41258317  1.49344165 -8.08584854e-02]
 [ 17.5  1.37377234  1.45700611 -8.32337718e-02]
 [ 18.0  1.33541049  1.42102058 -8.56100848e-02]
 [ 18.5  1.29749744  1.38548504 -8.79876055e-02]
 [ 19.0  1.26003297  1.35039950 -9.03665304e-02]
 [ 19.5  1.22301689  1.31576397 -9.27470733e-02]
 [ 20.0  1.18644897  1.28157843 -9.51294661e-02]]
error en rango:  0.09512946613018003

 EDO. Método de Taylor 3 términos 
# estima la solucion para muestras espaciadas h en eje x
# valores iniciales x0,y0
# entrega arreglo [[x,y]]
import numpy as np

def edo_taylor3t(d1y,d2y,x0,y0,h,muestras):
    tamano = muestras + 1
    estimado = np.zeros(shape=(tamano,2),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0]
    x = x0
    y = y0
    for i in range(1,tamano,1):
        y = y + h*d1y(x,y) + ((h**2)/2)*d2y(x,y)
        x = x+h
        estimado[i] = [x,y]
    return(estimado)

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)

# PROGRAMA PRUEBA
# 2Eva_IIT2016_T3_MN EDO Taylor 2, Tanque de agua

# INGRESO.
k=0.06
# d1y = y' = f, d2y = y'' = f'
d1y = lambda x,y: -k*(y**0.5)
d2y = lambda x,y: -(k/2)*(y**(-0.5))
x0 = 0
y0 = 3
h = 1/2
muestras = 40

# PROCEDIMIENTO
puntos = edo_taylor3t(d1y,d2y,x0,y0,h,muestras)
xi = puntos[:,0]
yi = puntos[:,1]

# Con Runge Kutta
puntosRK2 = rungekutta2(d1y,x0,y0,h,muestras)
xiRK2 = puntosRK2[:,0]
yiRK2 = puntosRK2[:,1]

# diferencias
diferencias = yi-yiRK2
error = np.max(np.abs(diferencias))
tabla = np.copy(puntos)
tabla = np.concatenate((puntos,np.transpose([yiRK2]),
                        np.transpose([diferencias])),
                       axis = 1)

# SALIDA
print('estimado[xi,yi Taylor,yi Runge-Kutta,diferencias]')
print(tabla)
print('error en rango: ', error)

# Gráfica
import matplotlib.pyplot as plt
plt.plot(xi[0],yi[0],'o',
         color='r', label ='[x0,y0]')
plt.plot(xi[1:],yi[1:],'o',
         color='g',
         label ='y Taylor 3 terminos')
plt.plot(xiRK2[1:],yiRK2[1:],'o',
         color='blue',
         label ='y Runge-Kutta 2Orden')
plt.axhline(y0/2)
plt.title('EDO: Taylor 3T vs Runge=Kutta 2Orden')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()

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.