s3Eva_2020PAOI_T2 Modelo epidemiológico no letal

Ejercicio: 3Eva_2020PAOI_T2 Modelo epidemiológico no letal

El ejercicio representa un sistema de ecuaciones diferenciales ordinarias, que serán resueltas usando Runge-Kutta de 2do Orden.

De compararse con la curva de contagios de Covid-19 se tienen diferencias en la población recuperada, pues el modelo se considera no letal por lo que no se contabiliza el número de fallecidos.

El módelo es el más básico y permite cambiar por ejemplo la tasa de infección, y se ve los cambios en la curva de infectados. Se puede observar lo que se indicaba como objetivo de «aplanar la curva» al disminuir la población expuesta mediante cambiar la tasa de infección al exponer más o menos población al contagio por iteacción entre «suceptibles» e «infectados.


Desarrollo analítico

Las fórmulas para el algoritmo se identifican como:
binfecta = 1.4
grecupera = 1/4

f(t,S,I,R) = -1.4(S)(I) g(t,S,I,R) = 1.4(S)(I) - \frac{1}{4}(I) w(t,S,I,R) = \frac{1}{4}(I)

que luego se usan en cada iteración que se registra en la tabla empezando con las condiciones iniciales

itera Si Ii Ri
0 1 0.001 0
1 0.9977 0.002809 0.0003937
2 0.9916 0.007862 0.0014987
3 tarea

itera = 1
K_{1S} = h. f(t_i,S_i,I_i,R_i) =

= 1(-1.4)( 1)(0.001) = -0.0014 K_{1I} = h. g(t_i,S_i,I_i,R_i) = = 1\Big((1.4)(1)(0.001) - \frac{1}{4}0.001 \Big) = 0.00115 K_{1R} = h. w(t_i,S_i,I_i,R_i) = = 1 \Big( \frac{1}{4} 0.001 \Big) = 0.00025 K_{2S} = h. f(t_i+h, S_i + K_{1S}, I_i+K_{1I}, R_i +K_{1R}) = = 1 \Big( (-1.4)(1-0.0014)(0.001+0.00115)\Big) = = -0.003005786 K_{2I} = h. g(t_i+h, S_i + K_{1S}, I_i+K_{1I}, R_i + K_{1R}) = = 1\Big((1.4)(1-0.0014)(0.001+0.00115) -\frac{1}{4} (0.001+0.00115)\Big) = = 0.002468286 K_{2R} = h w(t_i+h, S_i + K_{1S}, I_i+K_{1I}, R_i +K_{1R}) = 1\Big( \frac{1}{4}(0.001+0.00115) \Big) = 0.0005375 S_i = S_i + \frac{K_{1S}+K_{2S}}{2} = 1 + \frac{-0.0014 -0.003005786}{2} = 0.997797107 I_i = I_i + \frac{K_{1I}+K_{2I}}{2} = 0.001 + \frac{0.00115+0.002468286}{2} = 0.002809143 R_i = R_i + \frac{K_{1R}+K_{2R}}{2} = 0 + \frac{0.00025 + 0.0005375}{2} = 0.00039375 t_i = t_i + h = 1 + 1 = 2

 

itera = 2

K_{1S} = 1(-1.4)(0.9977)(0.002809) = -0.003924 K_{1I} = 1 \Big(1.4(0.9977)(0.002809) - \frac{1}{4}0.002809 \Big) = 0.003221 K_{1R} = 1(\frac{1}{4} 0.002809) = 0.0007022 K_{2S} = 1 \Big( -1.4*(0.9977-0.003924)(0.002809+0.003221) \Big) = = -0.008391 K_{2I} = 1 \Big(1.4(0.9977-0.003924)(0.002809+0.003221) - \frac{1}{4}(0.002809+0.003221) \Big) = 0.006883 K_{2R} = 1 \Big( \frac{1}{4}(0.002809+0.0032218) \Big) = 0.001507

 

S_i = 0.9977 + \frac{-0.003924 -0.008391}{2} = 0.9916

 

I_i = 0.002809 + \frac{0.003221+0.006883}{2} = 0.007862

 

R_i = 0.0003937 + \frac{0.0007022 + 0.001507}{2} = 0.001498

 

t_i = t_i + h = 2 + 1 = 3

 

itera 3 – TAREA


Instrucciones en Python

# 3Eva_2020PAOI_T2 Modelo epidemiológico no letal
# Sistemas EDO con Runge-Kutta de 2do Orden
import numpy as np

def rungekutta2_fgw(f,g,w,t0,x0,y0,z0,h,muestras):
    tamano = muestras +1
    tabla = np.zeros(shape=(tamano,4),dtype=float)
    tabla[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 * w(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 * w(ti+h, xi + K1x, yi+K1y, zi +K1z)

        xi = xi + (1/2)*(K1x+K2x)
        yi = yi + (1/2)*(K1y+K2y)
        zi = zi + (1/2)*(K1z+K2z)
        ti = ti + h
        
        tabla[i] = [ti,xi,yi,zi]
    tabla = np.array(tabla)
    return(tabla)

# Programa
# Parámetros de las ecuaciones

binfecta = 1.4
grecupera = 1/4
# Ecuaciones
f = lambda t,S,I,R : -binfecta*S*I
g = lambda t,S,I,R : binfecta*S*I - grecupera*I
w = lambda t,S,I,R : grecupera*I
# Condiciones iniciales
t0 = 0
S0 = 1.0
I0 = 0.001
R0 = 0.0

# parámetros del algoritmo
h = 1.0
muestras = 51

# PROCEDIMIENTO
tabla = rungekutta2_fgw(f,g,w,t0,S0,I0,R0,h,muestras)
ti = tabla[:,0]
Si = tabla[:,1]
Ii = tabla[:,2]
Ri = tabla[:,3]
# SALIDA
np.set_printoptions(precision=6)
print(' [ ti, Si, Ii, Ri]')
print(tabla)

# Grafica tiempos vs población
import matplotlib.pyplot as plt
plt.plot(ti,Si, label='S')
plt.plot(ti,Ii, label='I')
plt.plot(ti,Ri, label='R')
plt.title('Modelo SIR')
plt.xlabel('t tiempo')
plt.ylabel('población')
plt.legend()
plt.grid()
plt.show()

s3Eva_IIT2019_T4 completar polinomio de interpolación

Ejercicio: 3Eva_IIT2019_T4 completar polinomio de interpolación

Para visualizar la solución, se plantea graficar los polinomios que están completos S0(x) y S2(x).

S_0(x) = 1 + 1.1186x + 0.6938 x^3

 0.0 ≤ x ≤ 0.4

S_1(x) = 1.4918 + 1.4516(x-0.4) + c(x-0.4)^2 +d(x-0.4)^3

0.4 ≤ x ≤ 0.6

S_2(x) = 1.8221 + 1.8848(x-0.6) + +1.3336(x-0.6)^2 - 1.1113(x-0.6)^3

0.6 ≤ x ≤ 1.0

Como en trazadores cúbicos, lo que se usa es un polinomio por cada tramo muestreado para una curva contínua, etc. se tiene que los polinomios deben tener valores iguales en los puntos el eje x = 0.4 y 0.6

Por lo que se evalua con  los polinomios completos:

S_0(0.4) = 1 + 1.1186(0.4) + 0.6938 (0.4)^3 = 1.4918432 S_2(0.6) = 1.8221 + 1.8848(0.6-0.6) + +1.3336(0.6-0.6)^2-1.1113(0.6-0.6)^3=1.8221

Opción 1

Valores que se usan en los extremos del polinomio S1(x) para crear un sistema de dos ecuaciones y determinar los valores de c y d, completando el polinomio.

S_1(0.4) = 1.4918 + 1.4516(0.4-0.4) + c(0.4-0.4)^2 +d(0.4-0.4)^3

como los términos de c y d se hacen cero, hace falta una ecuación.

S_1(0.6) = 1.4918 + 1.4516(0.6-0.4) + + c(0.6-0.4)^2 +d(0.6-0.4)^3 = 1.8221 1.8221 = 1.4918 + 1.4516(0.2) + c(0.2)^2 +d(0.2)^3 0,03998 = 0.04c +0,008d

la otra ecuación se podría obtener usando la propiedad que las primeras derivadas de los polinomios deben ser iguales en los puntos x=0,4 y x= 0.6

S’0(0.2) =S’1(0.2)

Tarea: Desarrollar la siguiente ecuación y resolver


Opción 2

Si no recuerda la propiedad anterior, puede optar por usar otros conceptos para aproximar el resultado.

Si para el tramo en que se busca el polinomio se puede retroceder un tamaño de paso x = 0.2 y evualuar usando S0(0.2), se obtiene otropunto de referencia para crear un polinomio que pase por los mismos puntos.

S_0(0.2) = 1 + 1.1186*(0.2) + 0.6938 (0.2)^3

Se aplica lo mismo para un tamaño de paso más adelante de x = 0.6 es x = 0.8m se evalua S2(0.8) y se tienen suficientes puntos para usar cualquier método de interpolación y determinar el polinomio para el tramo faltante.

S_2(0.8) = 1.8221 + 1.8848(0.8-0.6) + +1.3336(0.8-0.6)^2 - 1.1113(0.8-0.6)^3
xi     = [0.        0.2        0.4      ]
S0(xi) = [1.        1.2292704  1.4918432]
xi     = [0.6       0.8        1.       ]
S2(xi) = [1.8221    2.2435136  2.7182728]

que permite hacer una tabla de puntos, y usando por ejemplo el método de interpolación de Lagrange  con x entre [0.2, 0.8] se obtiene otra forma del polinomio buscado:

p(x)=1.2292704 \frac{(x-0.4)(x-0.6)(x-0.8)}{(0.2-0.4)(0.2-0.6)(0.2-0.8)} + +1.4918432\frac{(x-0.2)(x-0.6)(x-0.8)}{0.4-0.2)(0.4-0.6)(0.4-0.8)}+ +1.8221\frac{(x-0.2)(x-0.4)(x-0.8)}{(0.6-0.2)(0.6-0.4)(0.6-0.8)} + + 2.2435136\frac{(x-0.2)(x-0.4)(x-0.6)}{(0.8-0.2)(0.8-0.4)(0.8-0.6)}

Tarea: continuar con el desarrollo


El literal b, requiere usar un metodo de búsqueda de raíces, para el cual se puede usar incluso bisección.

Tarea: continuar con el desarrollo

s3Eva_IIT2019_T3 Preparación de terreno en refineria

Ejercicio: 3Eva_IIT2019_T3 Preparación de terreno en refineria

Se requiere usar el nivel inicial en la matriz, para restar del nivel requerido que es constante 220

Nivel inicio (m) 0 50 100 150 200
0 241 239 238 236 234
25 241 239 237 235 233
50 241 239 236 234 231
75 242 239 236 232 229
100 243 239 235 231 227

lo que genera la matriz de diferencias. El valor es positivo indica remoción, el valor negativo indica por rellenar.

Diferencia (m) 0 50 100 150 200
0 21 19 18 16 14
25 21 19 17 15 13
50 21 19 16 14 11
75 22 19 16 12 9
100 23 19 15 11 7

El volumen se puede calcular por un método en cada fila, y luego los resultados por columnas por otro método o el mismo.
Por ejemplo Simpson de 1/3

I= \frac{hx}{3}(f(x_0) +4f(x_1)+f(x_2))

con lo que se obtiene:

I_{fila}(0) = \frac{50}{3}(21 +4(19)+18) +\frac{50}{3}(18 +4(16)+14) =24750 I_{fila}(25) = = \frac{50}{3}(21 +4(19)+17) + \frac{50}{3}(17 +4(15)+13) = 23383,33 I_{fila}(50) = \frac{50}{3}(21 +4(19)+16) + \frac{50}{3}(16 +4(14)+11) = 22000 I_{fila}(75) = \frac{50}{3}(22 +4(19)+16) + \frac{50}{3}(16 +4(12)+5) =21850 I_{fila}(100) = \frac{50}{3}(23 +4(19)+15) + \frac{50}{3}(15 +4(11)+7) = 20483,33

y usando el otro eje, se completa el volumen usando dos veces simpson:

Volumen = \frac{h_y}{3}(f(x_0) +4f(x_1)+f(x_2)) Remover = \frac{25}{3}(24750 +4(23383,33)+22000) + + \frac{25}{3}(22000 +4(21850)+20483,33)=2251388,89

El signo lo trae desde la diferencia, y muestra el sentido del desnivel.

Se adjunta la gráfica de superficie en azul como referencia del signo,  respecto al nivel requerido en color verde.

Error de truncamiento

la cota del error de truncamiento se estima como O(h5)

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

para un valor de z entre [a,b]

para cuantificar el valor, se puede usar la diferencia finita Δ4f, pues con la derivada sería muy laborioso.

s3Eva_IIT2019_T2 Diferenciación, valor en frontera

Ejercicio: 3Eva_IIT2019_T2 Diferenciación, valor en frontera

La ecuación de problema de valor de frontera, con h = 1/4 = 0.25:

y'' = -(x+1)y' + 2y + (1-x^2) e^{-x}

0 ≤ x ≤ 1
y(0) = -1
y(1) = 0

Se interpreta en la tabla como los valores que faltan por encontrar:

i 0 1 2 3 4
xi 0 1/4 1/2 3/4 1
yi -1 0

Por ejemplo, se usan las derivadas en diferencias finitas divididas centradas para segunda derivada y hacia adelante para primera derivada.

Semejante al procedimiento aplicado para métodos con EDO.

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

Se formula entonces la ecuación en forma discreta, usando solo los índices para los puntos yi:

\frac{y[i+1]-2y[i]+y[i-1]}{h^2} = -(x_i +1) \frac{y[i+1]-y[i]}{h} + 2y[i] + (1-x_i ^2) e^{-x_i}

se multiplica todo por h2

y[i+1]-2y[i]+y[i-1] = -h(x_i +1)(y[i+1]-y[i]) + + 2 h^2 y[i] + h^2 (1-x_i ^2) e^{-x_i}
y[i+1]-2y[i]+y[i-1] = -h(x_i +1)y[i+1] + + h(x_i +1)y[i] + 2 h^2 y[i] + h^2 (1-x_i ^2) e^{-x_i}
y[i+1](1 +h(x_i +1)) + y[i](-2 - h(x_i +1)-2 h^2)+y[i-1] = = h^2 (1-x_i ^2) e^{-x_i}

Ecuación que se aplica en cada uno de los puntos desconocidos con i =1,2,3


i = 1

y[2](1 +h(x_1 +1)) + y[1](-2 - h(x_1 +1)-2 h^2)+y[0] = = h^2 (1-x_1 ^2) e^{-x_1} y[2]\Big(1 +\frac{1}{4}\Big(\frac{1}{4} +1\Big)\Big) + y[1]\Big(-2 - \frac{1}{4}\Big(\frac{1}{4} +1\Big)-2 \Big(\frac{1}{4}\Big)^2\Big) -1 = = \Big(\frac{1}{4}\Big)^2 \Big(1-\Big(\frac{1}{4}\Big) ^2\Big) e^{-1/4} y[2]\Big(1 +\frac{5}{16} \Big) + y[1]\Big(-2 - \frac{5}{16}- \frac{2}{16}\Big) = = 1+ \frac{1}{16} \Big(1-\frac{1}{16}\Big) e^{-1/4} \frac{21}{16} y[2]- \frac{39}{16}y[1] = 1+ \frac{15}{16^2} e^{-1/4} - \frac{39}{16}y[1] + \frac{21}{16} y[2] = 1+ \frac{15}{16^2} e^{-1/4}

multiplicando ambos lados de la ecuacion por 16 y reordenando

- 39 y[1] + 21 y[2] = 16+ \frac{15}{16} e^{-1/4}

i = 2

y[3](1 +h(x_2 +1)) + y[2](-2 - h(x_2 +1)-2 h^2)+y[1] = = h^2 (1-x_2 ^2) e^{-x_2} y[3]\Big(1 +\frac{1}{4}\Big(\frac{1}{2} +1\Big)\Big) + y[2]\Big(-2 - \frac{1}{4}\Big(\frac{1}{2} +1\Big)-2 \Big(\frac{1}{4}\Big)^2\Big)+y[1] = = \Big(\frac{1}{4}\Big)^2 \Big(1-\Big(\frac{1}{2}\Big) ^2\Big) e^{-\frac{1}{2}} y[3]\Big(1 +\frac{3}{8}\Big) + y[2]\Big(-2 - \frac{1}{4}\Big(\frac{3}{2}\Big)- \frac{1}{8}\Big)+y[1] = = \Big(\frac{1}{16}\Big) \Big(1-\frac{1}{4}\Big) e^{-\frac{1}{2}} \frac{11}{8} y[3] - \frac{21}{8} y[2]+y[1] = \frac{1}{16} \Big(\frac{3}{4}\Big) e^{-\frac{1}{2}}

multiplicando ambos lados por 8 y reordenando,

8y[1] - 21 y[2] + 11 y[3] = \frac{3}{8} e^{-\frac{1}{2}}

i = 3

y[4](1 +h(x_3 +1)) + y[3](-2 - h(x_3 +1)-2 h^2)+y[2] = = h^2 (1-x_3 ^2) e^{-x_3} (0) \Big(1 +h\Big(x_3 +1\Big)\Big) + y[3]\Big(-2 - \frac{1}{4}\Big(\frac{3}{4} +1\Big)-2 \frac{1}{16}\Big)+y[2] = = \frac{1}{16}\Big (1-\Big(\frac{3}{4}\Big) ^2 \Big) e^{-3/4} y[3]\Big(-2 - \frac{7}{16} - \frac{2}{16}\Big)+y[2] = \frac{1}{16}\Big (1-\frac{9}{16}\Big) e^{-3/4}

multiplicando todo por 16 y reordenando:

16 y[2] - 41 y[3] = \frac{7}{16} e^{-3/4}

Con lo que se puede crear la forma matricial del sistema de ecuaciones Ax=B

\begin{bmatrix} -39 && 21 && 0\\8 && 21 && 11 \\ 0 && 16 && 41 \end{bmatrix}\begin{bmatrix} y[1]\\ y[2] \\y[3] \end{bmatrix} =\begin{bmatrix} 16+ \frac{15}{16} e^{-1/4} \\ \frac{3}{8} e^{-\frac{1}{2}} \\ \frac{7}{16} e^{-3/4} \end{bmatrix}

con lo que resolviendo la matriz se obtienen los valroes para y[1], y[2], y[3]

solución de matriz: 
[-0.59029143 -0.29958287 -0.12195088]

con lo que se completan los puntos de la tabla,

Solución de ecuación
x[i] =  [ 0.   0.25        0.5         0.75        1.  ]
y[i] =  [-1.  -0.59029143 -0.29958287 -0.12195088  0.  ]

con la siguiente gráfica:


Algoritmo para solución de matriz con Python

tarea: modificar para cambiar el valor del tamaño de paso.

# Problema de Frontera
import numpy as np
import matplotlib.pyplot as plt

# INGRESO

h = 1/4
y0 = -1
y1 = 0

xi = np.arange(0,1+h,h)
n = len(xi)
yi = np.zeros(n,dtype=float)
yi[0] = y0
yi[n-1] = y1
    
A = np.array([[-39.0, 21,  0],
              [  8.0, -21, 11],
              [  0.0, 16, -41]])
B = np.array([16+(15/16)*np.exp(-1/4),
              (3/8)*np.exp(-1/2),
              (7/16)*np.exp(-3/4)])

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

for i in range(1,n-1,1):
    yi[i] = x[i-1]

# SALIDA
print('solución de matriz: ')
print(x)
print('Solución de ecuación')
print('x[i] = ',xi)
print('y[i] = ',yi)

# Grafica
plt.plot(xi,yi,'ro')
plt.plot(xi,yi)
plt.show()

s3Eva_IIT2019_T1 Lanzamiento de Cohete

Ejercicio: 3Eva_IIT2019_T1 Lanzamiento de Cohete

A partir de la tabla del enunciado  se realiza la tabla de diferencias finitas.

i ti fi Δfi Δ2fi Δ3fi Δ4fi Δ5fi
1 0 0 32 -6 0 0 0
2 25 32 26 -6 0 0
3 50 58 20 -6 0
4 75 78 14 -6
5 100 92 8
6 125 100

Observando que a partir de la tercera diferencia finita  los valores son cero, por lo que se usa la fórmula general de diferencias finitas divididas hasta el polinomio de grado 2.

p_2 (x) = f_0 + \frac{\Delta f_0}{h} (x - x_0) + + \frac{\Delta^2 f_0}{2!h^2} (x - x_0)(x - x_1)

al sustituir los valores conocidos, se convierte en,

p_2 (t) =0 + \frac{32}{25} (t -0) + + \frac{-6}{2(25)^2} (t -0)(t - 25) =\frac{32}{25}t + \frac{-3}{(25)^2} (t^2 - 25t) =\frac{32}{25}t + \frac{-3}{(25)^2} t^2 - \frac{-3}{(25)^2}25t =\frac{7}{5}t - \frac{3}{625} t^2 y(t) =p_2 (t) =1.4 t - 0.0048 t^2

Con lo que se puede obtener la velocidad:

y'(t) = 1.4 - 0.0096 t

y luego la aceleración:

y''(t) = - 0.0096

Si el error es el próximo término del polinomio Δ3fi  entonces se estima en cero.

Tarea:  Evaluar la velocidad y aceleración para cada punto de la tabla

La gráfica del polinomio encontrado es:

Algoritmo en Python

El algoritmo realizado en Python entrega los siguientes resultados:

[[  i,  ti,  fi, df1, df2, df3, df4, df5,  df6]]
[[  1.   0.   0.  32.  -6.   0.   0.   0.   0.]
 [  2.  25.  32.  26.  -6.   0.   0.   0.   0.]
 [  3.  50.  58.  20.  -6.   0.   0.   0.   0.]
 [  4.  75.  78.  14.  -6.   0.   0.   0.   0.]
 [  5. 100.  92.   8.   0.   0.   0.   0.   0.]
 [  6. 125. 100.   0.   0.   0.   0.   0.   0.]]
polinomio:
-0.0048*t**2 + 1.4*t

las instrucciones en Python son:

# 3Eva_IIT2019_T1 Lanzamiento de Cohete
# Tarea: Verificar tamaño de vectores
#        considerar puntos no equidistantes en eje t
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

# INGRESO , Datos de prueba
ti = np.array([0.0, 25, 50, 75, 100, 125])
fi = np.array([0.0, 32, 58, 78, 92, 100])

# PROCEDIMIENTO
# Tabla de diferencias finitas
titulo = ['i','ti','fi']
n = len(ti)

# cambia a forma de columnas
i = np.arange(1,n+1,1)
i = np.transpose([i])
ti = np.transpose([ti])
fi = np.transpose([fi])

# Añade matriz de diferencias
dfinita = np.zeros(shape=(n,n),dtype=float)
tabla = np.concatenate((i,ti,fi,dfinita), axis=1)

# Sobre matriz de diferencias, por columnas
[n,m] = np.shape(tabla)
c = 3
diagonal = n-1
while (c<m):
    # Aumenta el título para cada columna
    titulo.append('df'+str(c-2))
    # calcula cada diferencia por fila
    f = 0
    while (f < diagonal):
        tabla[f,c] = tabla[f+1,c-1]-tabla[f,c-1]
        f = f+1
    
    diagonal = diagonal - 1
    c = c+1

# POLINOMIO con diferencias finitas
# caso: puntos en eje t equidistantes
dfinita = tabla[:,3:]
n = len(dfinita)
t = sym.Symbol('t')
h = ti[1,0]-ti[0,0]
polinomio = fi[0,0]
for c in range(1,n,1):
    denominador = np.math.factorial(c)*(h**c)
    factor = dfinita[0,c-1]/denominador
    termino=1
    for f  in range(0,c,1):
        termino = termino*(t-ti[f])
    polinomio = polinomio + termino*factor

# simplifica polinomio, multiplica los (t-ti)
polinomio = polinomio.expand()

# para evaluacion numérica
pt = sym.lambdify(t,polinomio)

# Puntos para la gráfica
a = np.min(ti)
b = np.max(ti)
muestras = 101
ti_p = np.linspace(a,b,muestras)
fi_p = pt(ti_p)

# SALIDA
print([titulo])
print(tabla)
print('polinomio:')
print(polinomio)

# Gráfica
plt.title('Interpolación polinómica')
plt.plot(ti,fi,'o', label = 'Puntos')
plt.plot(ti_p,fi_p, label = 'Polinomio')
plt.legend()
plt.show()

s3Eva_IT2019_T3 Difusión en sólidos

Ejercicio: 3Eva_IT2019_T3 Difusión en sólidos

Siguiendo el procedimiento planteado en la sección EDP parabólicas, se plantea la malla del ejercicio:

Para plantear la ecuación en forma discreta:

\frac{\phi_{i,j+1}-\phi_{i,j}}{\Delta t}=D\frac{\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}}{(\Delta x)^2}

y resolver usando el método explícito para ecuaciones parabólicas, obteniendo el siguiente resultado:

\phi_{i,j+1}-\phi_{i,j}=D\frac{\Delta t }{\Delta x^2}(\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}) \lambda = D\frac{\Delta t }{\Delta x^2} \phi_{i,j+1}-\phi_{i,j}=\lambda (\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}) \phi_{i,j+1} =\lambda \phi_{i+1,j}-2\lambda\phi_{i,j}+\lambda\phi_{i-1,j}+\phi_{i,j} \phi_{i,j+1} =\lambda \phi_{i+1,j}(1-2\lambda)\phi_{i,j}+\lambda\phi_{i-1,j} \phi_{i,j+1} =P \phi_{i+1,j}+Q\phi_{i,j}+R\phi_{i-1,j}

siendo:
P = λ = 0.16 (Δx/100)/Δx2 = 0.0016/Δx = 0.0016/0.02=0.08
Q = 1-2λ = 1-2*(0.08) = 0.84
R = λ =0.08

\phi_{i,j+1} =0.08 \phi_{i+1,j}+ 0.84\phi_{i,j}+0.08\phi_{i-1,j}

Iteración 1 en tiempo:
i=1, j=0

\phi_{1,1} =0.08 \phi_{2,0}+ 0.84\phi_{1,0}+0.08\phi_{0,0} \phi_{1,1} =0.08 (0)+ 0.84(0)+0.08(5)=0.4

i=2,j=0

\phi_{2,1} =0.08 \phi_{3,0}+ 0.84\phi_{2,0}+0.08\phi_{1,0} = 0

Para los proximos valores i>2, todos los resultados son 0

Iteración 2 en tiempo
i=1, j=1

\phi_{1,2} =0.08 \phi_{2,0}+ 0.84\phi_{1,0}+0.08\phi_{0,0}

\phi_{1,2} =0.08 (0)+ 0.84(0.4)+0.08(5)=0.736
i=2, j=1

\phi_{2,2} =0.08 \phi_{3,1}+ 0.84\phi_{2,1}+0.08\phi_{1,1} \phi_{2,2} =0.08(0)+ 0.84(0)+0.08(0.4) = 0.032

i=3, j=1

\phi_{3,2} =0.08\phi_{4,1}+ 0.84\phi_{3,1}+0.08\phi_{2,1}=0

Para los proximos valores i>3, todos los resultados son 0

Tarea: Desarrollar la iteración 3 en el tiempo.

siguiendo las iteraciones se tiene la siguiente tabla:

[[5.0, 0.000, 0.000, 0.00000, 0.00000 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
 [5.0, 0.400, 0.000, 0.00000, 0.00000 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
 [5.0, 0.736, 0.032, 0.00000, 0.00000 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
 [5.0, 1.021, 0.085, 0.00256, 0.00000 , 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
 [5.0, 1.264, 0.153, 0.00901, 0.00020, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
...
]

Con lo que se obtiene la siguiente gráfica.

El resultado se interpreta mejor con una animación: (tarea)

Tarea: Presentar el orden de error de la ecuación basado en las fórmulas de diferenciación


Algorirmo en Python

# 3Eva_IT2019_T3 Difusión en sólidos
# EDP parabólicas. método explícito,usando diferencias finitas
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# Valores de frontera
Ta = 5
Tb = 0
T0 = 0
# longitud en x
a = 0
b = 0.1
# Constante K
K = 1/(1.6e-1)
# Tamaño de paso
dx = 0.02
dt = dx/100
# iteraciones en tiempo
n = 50

# PROCEDIMIENTO
# iteraciones en longitud
xi = np.arange(a,b+dx,dx)
m = len(xi)
ultimox = m-1

# Resultados en tabla u[x,t]
u = np.zeros(shape=(m,n), dtype=float)

# valores iniciales de u[:,j]
j=0
ultimot = n-1
u[0,j]= Ta
u[1:ultimox,j] = T0
u[ultimox,j] = Tb

# factores P,Q,R
lamb = dt/(K*dx**2)
P = lamb
Q = 1 - 2*lamb
R = lamb

# Calcula U para cada tiempo + dt
j = 0
while not(j>=ultimot):
    u[0,j+1] = Ta
    for i in range(1,ultimox,1):
        u[i,j+1] = P*u[i-1,j] + Q*u[i,j] + R*u[i+1,j]
    u[m-1,j+1] = Tb
    j=j+1

# SALIDA
print('Tabla de resultados')
np.set_printoptions(precision=2)
print(u)

# Gráfica
salto = int(n/10)
if (salto == 0):
    salto = 1
for j in range(0,n,salto):
    vector = u[:,j]
    plt.plot(xi,vector)
    plt.plot(xi,vector, '.r')
plt.xlabel('x[i]')
plt.ylabel('phi[i,j]')
plt.title('Solución EDP parabólica')
plt.show()

La animación se complementa con lo mostrado en la sección de Unidades.

s3Eva_IT2019_T2 Integral con interpolación

Ejercicio: 3Eva_IT2019_T2 Integral con interpolación

El ejercicio considera dos partes: interpolación e integración

a. Interpolación

Se requiere aproximar la función usando tres puntos. Para comprender la razón del método solicitado, se compara la función con dos interpolaciones:

a.1 Lagrange
a.2 Trazador cúbico sujeto

Observando la gráfica se aclara que en éste caso, una mejor aproximación se obtiene con el método  de trazador cúbico sujeto. Motivo por lo que el tema tiene un peso de 40/100 puntos

Los valores a considerar para la evaluación son:

puntos referencia xi,yi: 
[0.         0.78539816 1.57079633]
[ 0.          0.62426595 -0.97536797]
derivadas en los extremos:  
    3.141592653589793 
    0.6929852019184021
Polinomio de Lagrange
-1.80262534301178*x**2 + 2.21061873102778*x
Trazadores cúbicos sujetos
[0.         0.78539816]
-0.548171611756137*x**3 - 2.55744517923506*x**2 + 3.14159265358979*x

[0.78539816 1.57079633]
4.66299098804068*x**3 - 14.8359577843727*x**2 + 12.7851139029174*x - 2.52466795930204

------------------
Valores calculados para Trazadores cúbicos sujetos:
Matriz A: 
[[-0.26179939 -0.13089969  0.        ]
 [ 0.78539816  3.14159265  0.78539816]
 [ 0.          0.13089969  0.26179939]]
Vector B: 
[  2.34675256 -16.9893436    2.72970237]
coeficientes S: 
[-5.11489036 -7.69808822 14.27573913]
coeficientes a,b,c,d
[-0.54817161  4.66299099]
[-2.55744518 -3.84904411]
[ 3.14159265 -1.89005227]
[0.         0.62426595]

b. Integración

Como forma de comparacíon de resultados, se requiere integrar con varios métodos para comparar resultados y errores.

b.1 Integración con Cuadratura de Gauss, usando el resultado de trazador cúbico.

Se integra en cada tramo de cada polinomio:

Trazadores cúbicos sujetos
[0.         0.78539816]
-0.548171611756137*x**3 - 2.55744517923506*x**2 + 3.14159265358979*x

Se obtienen los puntos del método de cuadratura desplazados en el rango:

xa:  0.16597416116944688
xb:  0.6194240022280014
area:  0.5037962958529855

Para el segundo tramo:

[0.78539816 1.57079633]
4.66299098804068*x**3 - 14.8359577843727*x**2 + 12.7851139029174*x - 2.52466795930204
xa:  0.9513723245668951
xb:  1.4048221656254496
area:  -0.2706563884589365

Con lo que el integral total es:

Integral total:  0.23313990739404894

b.2 Integración analítica

\int_0^{\pi /2}sin(\pi x) dx

u = πx
du/dx = π
dx = du/π

se convierte en:

\frac{1}{\pi}\int sin(u) du \frac{1}{\pi}(-cos(u))

volviendo a la variable x:

\frac{1}{\pi}(-cos(\pi x)) \Big\rvert_{0}^{\frac{\pi}{2}} -\frac{1}{\pi}(cos(\pi \frac{\pi}{2})-cos(\pi(0))) = 0.24809580527879377

c. Estimación del error

Se restan los resultados de las secciones b.1 y b.2

error = |0.24809580527879377 – 0.23313990739404894 |

error = 0.014955897884744829


Algoritmo en Python

separado por literales

# 3Eva I T 2019 Interpola e Integra
import numpy as np
import sympy as sym
import matplotlib.pyplot as plt

def interpola_lagrange(xi,yi):
    '''
    Interpolación con método de Lagrange
    resultado: polinomio en forma simbólica
    '''
    # PROCEDIMIENTO
    n = len(xi)
    x = sym.Symbol('x')
    # Polinomio
    polinomio = 0
    for i in range(0,n,1):
        # Termino de Lagrange
        termino = 1
        for j  in range(0,n,1):
            if (j!=i):
                termino = termino*(x-xi[j])/(xi[i]-xi[j])
        polinomio = polinomio + termino*yi[i]
    # Expande el polinomio
    polinomio = polinomio.expand()
    return(polinomio)

def traza3sujeto(xi,yi,u,v):
    '''
    Trazador cúbico sujeto, splines
    resultado: polinomio en forma simbólica
    '''
    n = len(xi)
    # Valores h
    h = np.zeros(n-1, dtype=float)
    # Sistema de ecuaciones
    A = np.zeros(shape=(n,n), dtype=float)
    B = np.zeros(n, dtype=float)
    S = np.zeros(n-1, dtype=float)
    # coeficientes
    a = np.zeros(n-1, dtype=float)
    b = np.zeros(n-1, dtype=float)
    c = np.zeros(n-1, dtype=float)
    d = np.zeros(n-1, dtype=float)
    
    polinomios=[]
    
    if (n>=3):
        for i in range(0,n-1,1):
            h[i]=xi[i+1]-xi[i]
        A[0,0] = -h[0]/3
        A[0,1] = -h[0]/6
        B[0] = u-(yi[1]-yi[0])/h[0]
        for i in range(1,n-1,1):
            A[i,i-1] = h[i-1]
            A[i,i] = 2*(h[i-1]+h[i])
            A[i,i+1] = h[i]
            B[i] = 6*((yi[i+1]-yi[i])/h[i] - (yi[i]-yi[i-1])/h[i-1])
        A[n-1,n-2] = h[n-2]/6
        A[n-1,n-1] = h[n-2]/3
        B[n-1] = v-(yi[n-1]-yi[n-2])/h[n-2]

        # Resolver sistema de ecuaciones
        S = np.linalg.solve(A,B)

        # Coeficientes
        for i in range(0,n-1,1):
            a[i]=(S[i+1]-S[i])/(6*h[i])
            b[i]=S[i]/2
            c[i]=(yi[i+1]-yi[i])/h[i]-(2*h[i]*S[i]+h[i]*S[i+1])/6
            d[i]=yi[i]
      
        # polinomio en forma simbólica
        x=sym.Symbol('x')
        polinomios=[]
        for i in range(0,n-1,1):
            ptramo = a[i]*(x-xi[i])**3 + b[i]*(x-xi[i])**2 + c[i]*(x-xi[i])+ d[i]
            ptramo = ptramo.expand()
            polinomios.append(ptramo)
        parametros = [A,B,S,a,b,c,d]                                                           
    return(polinomios, parametros)

# INGRESO
f = lambda x: np.sin(np.pi*x)
muestrasf = 20
a = 0
b = np.pi/2
# Derivadas en los extremos
u = np.pi*np.cos(np.pi*a)
v = np.pi*np.cos(np.pi*b)
muestras = 3

# literal a
# PROCEDIMIENTO
xif = np.linspace(a,b,muestrasf)
yif = f(xif)

xi = np.linspace(a,b,muestras)
yi = f(xi)

# Usando Lagrange
x = sym.Symbol('x')
pL = interpola_lagrange(xi,yi)
pxL = sym.lambdify(x,pL)
pxiL =  pxL(xif)

# Trazador cúbico sujeto
pS, parametros = traza3sujeto(xi,yi,u,v)
pxiS = np.zeros(muestrasf,dtype=float)

# Evalua trazadores cúbicos sujetos
i=0
ap = xi[i]
bp = xi[i+1]
poli = sym.lambdify(x, pS[i])
for j in range(0,muestrasf,1):
    punto = xif[j]
    if (punto>bp):
        i = i+1
        ap = xi[i]
        bp = xi[i+1]
        poli = sym.lambdify(x,pS[i])
    pxiS[j] = poli(punto)

# SALIDA
print('puntos referencia xi,yi: ')
print(xi)
print(yi)
print('derivadas en los extremos: ',u,v)
print('Polinomio de Lagrange')
print(pL)
print('Trazadores cúbicos sujetos')
n = len(xi)
for i in range(0,n-1,1):
    print(xi[i:i+2])
    print(pS[i])
# Parametros de Trazadores cúbicos sujetos
print('Matriz A: ')
print(parametros[0])
print('Vector B: ')
print(parametros[1])
print('coeficientes S: ')
print(parametros[2])
print('coeficienetes a,b,c,d')
print(parametros[3])
print(parametros[4])
print(parametros[5])
print(parametros[6])

# Gráficas
plt.plot(xif,yif, label='funcion')
plt.plot(xi,yi,'o', label='muestras')
plt.plot(xif,pxiL, label='p(x)_Lagrange')
plt.plot(xif,pxiS, label='p(x)_Traza3Sujeto')
plt.legend()
plt.xlabel('x')
plt.show()

# literal b
# cuadratura de Gauss de dos puntos
def integraCuadGauss2p(funcionx,a,b):
    x0 = -1/np.sqrt(3)
    x1 = -x0
    xa = (b+a)/2 + (b-a)/2*(x0)
    xb = (b+a)/2 + (b-a)/2*(x1)
    area = ((b-a)/2)*(funcionx(xa) + funcionx(xb))
    print('xa: ',xa)
    print('xb: ',xb)
    print('area: ', area)
    return(area)

# INGRESO
f0 = sym.lambdify(x,pS[0])
f1 = sym.lambdify(x,pS[1])
# Procedimiento
I0 = integraCuadGauss2p(f0,xi[0],xi[1])
I1 = integraCuadGauss2p(f1,xi[1],xi[2])
It = I0+I1

# SALIDA
print('Integral 1: ', I0)
print('Integral 2: ', I1)
print('Integral total: ',It)

s3Eva_IT2019_T1 Ecuaciones simultáneas

Ejercicio: 3Eva_IT2019_T1 Ecuaciones simultáneas

Para plantear la intersección de las ecuaciones se pueden simplificar como:

y_1 = -x^2 +x + 0.75 y+5xy=x^3 y(1+5x)=x^3 y_2=\frac{x^3}{1+5x}

Quedando dos ecuaciones simplificadas:

y_1 = -x^2 +x + 0.75 y_2 = \frac{x^3}{1+5x}

cuyas gráficas son:

dónde se puede observar la intersección alrededor de 1.3

Restando ambas ecuaciones, se tiene que encontrar el valor de x para que el resultado sea cero.

y_1(x)-y_2(x)= -x^2 +x + 0.75 -\frac{x^3}{1+5x} f(x) = -x^2 +x + 0.75 -\frac{x^3}{1+5x} = 0

Para encontrarla derivada se procesa la expresión:

(1+5x)(-x^2 +x + 0.75) -x^3 = 0(1+5x) -6x^3+4x^2+4.75x+0.75 = 0 f'(x)= -18x^2 +8x + 4.75

Se usa el punto inicial x0=1 definido en el enunciado y se realizan las iteraciones siguiendo el algoritmo.

Se tiene que la raiz es:

raiz en:  1.3310736382369661
 [  xi, 	 xnuevo,	 fxi,	 dfxi, 	 tramo]
[[ 1.000e+00  1.111e+00  5.833e-01 -5.250e+00  1.111e-01]
 [ 1.111e+00  1.160e+00  4.173e-01 -8.583e+00  4.862e-02]
 [ 1.160e+00  1.193e+00  3.353e-01 -1.018e+01  3.293e-02]
 [ 1.193e+00  1.217e+00  2.766e-01 -1.131e+01  2.445e-02]
 [ 1.217e+00  1.236e+00  2.313e-01 -1.218e+01  1.899e-02]
 [ 1.236e+00  1.251e+00  1.951e-01 -1.286e+01  1.517e-02]
....
]

Algoritmo en Python

# 3Eva I T 2019 ecuaciones simultaneas
import numpy as np
import matplotlib.pyplot as plt

def newton_raphson(funcionx, fxderiva, xi, tolera):
    '''
    funciónx y fxderiva son de forma numérica
    xi es el punto inicial de búsqueda
    '''
    tabla = []
    tramo = abs(2*tolera)
    while (tramo>=tolera):
        fxi = funcionx(xi)
        dfxi = fxderiva(xi)
        xnuevo = xi - fxi/dfxi
        tramo = abs(xnuevo-xi)
        tabla.append([xi,xnuevo,fxi,dfxi,tramo])
        xi = xnuevo
    return(xi,tabla)

# INGRESO
y1 = lambda x: -x**2 +x +0.75
y2 = lambda x: (x**3)/(1+5*x)
a = 0.5
b = 1.5
muestras = 20

f = lambda x: -x**2+x+0.75-x**3/(1+5*x)
df = lambda x: -18*(x**2)+8*x +4.75
tolera = 1e-4
x0 = 1

# PROCEDIMIENTO
# datos para la gráfica
xi = np.linspace(a,b,muestras)
yi1 = y1(xi)
yi2 = y2(xi)
fi = f(xi)
# determina raiz
raiz, tabla = newton_raphson(f, df, x0, tolera)
tabla = np.array(tabla)

# SALIDA
np.set_printoptions(precision=3)
print('raiz en: ',raiz)
print(' [  xi, \t xnuevo,\t fxi,\t dfxi, \t tramo]')
print(tabla)

# Gráfica
plt.plot(xi,yi1, label ='yi1')
plt.plot(xi,yi2, label ='yi2')
plt.plot(xi,fi, label ='fi=yi1-yi2')
plt.axvline(raiz,linestyle='dashed')
plt.axhline(0)
plt.xlabel('x')
plt.legend()
plt.title('ecuaciones simultáneas')
plt.show()

s3Eva_IIT2018_T2 Drenar tanque cilíndrico

Ejercicio: 3Eva_IIT2018_T2 Drenar tanque cilíndrico

La ecuación a desarrollar es:

\frac{\delta y}{\delta t} = -k\sqrt{y}

con valores de k =0.5, y(0)=9


Formula de Taylor con término de error:

P_{n}(x) = \sum_{k=0}^{n} \frac{f^{(k)}(x_0)}{k!} (x-x_0)^k P_{n}(x) = f(x_0)+\frac{f'(x_0)}{1!} (x-x_0) + + \frac{f''(x_0)}{2!}(x-x_0)^2 + + \frac{f'''(x_0)}{3!}(x-x_0)^3 + \text{...}

Se requiere la 2da y 3ra derivadas:

\frac{\delta^2 y}{\delta t^2} = -k\frac{1}{2} y^{(\frac{1}{2}-1)} = -\frac{k}{2} y^{-\frac{1}{2}} \frac{\delta^3 y}{\delta t^3} = -\frac{k}{2}\Big(-\frac{1}{2}\Big) y^{(-\frac{1}{2}-1)} = \frac{k}{4} y^{-\frac{3}{2}}

con lo que inicia las iteraciones y cálculo del error, con avance de 0.5 para t.


t=0 , y(0) = 9


t = 0.5

\frac{\delta y(0)}{\delta t} = -(0.5)\sqrt{9} = -1.5 \frac{\delta^2 y(0)}{\delta t^2} = -\frac{0.5}{2} 9^{-\frac{1}{2}} = - 0.08333 \frac{\delta^3 y(0)}{\delta t^3} = \frac{0.5}{4} 9^{-\frac{3}{2}} = 0.004628 P_{2}(0.5) = 9 - 1.5 (0.5-0) + \frac{-0.08333}{2}(0.5-0)^2 P_{2}(0.5) = 8.2395

Error orden de:

Error = \frac{0.004628}{3!}(0.5-0)^3 = 9.641 . 10^{-5}

t = 1

\frac{\delta y(0.5)}{\delta t} = -(0.5)\sqrt{8.2395} = -1.4352 \frac{\delta^2 y(0.5)}{\delta t^2} = -\frac{0.5}{2} (8.2395)^{-\frac{1}{2}} = - 0.08709 \frac{\delta^3 y(0.5)}{\delta t^3} = \frac{0.5}{4} (8.2395)^{-\frac{3}{2}} = 0.005285 P_{2}(1) = 8.2395 - 1.4352(1-0.5) + \frac{-0.08709}{2}(1-0.5)^2 P_{2}(1) = 7.5110

Error orden de:

Error = \frac{0.005285}{3!}(1-0.5)^3 = 4.404 . 10^{-4}

t = 1.5

\frac{\delta y(1)}{\delta t} = -(0.5)\sqrt{7.5110} = -1.3703 \frac{\delta^2 y(1)}{\delta t^2} = -\frac{0.5}{2} (7.5110)^{-\frac{1}{2}} = - 0.09122 \frac{\delta^3 y(1)}{\delta t^3} = \frac{0.5}{4} (7.5110)^{-\frac{3}{2}} = 0.006072 P_{2}(1.5) = 7.5110 - 1.3703(1.5-1) + \frac{-0.09122}{2}(1.5-1)^2 P_{2}(1.5) = 6.8144

Error orden de:

Error = \frac{0.006072}{3!}(1.5-1)^3 = 1.4 . 10^{-4}

t = 2

\frac{\delta y(1.5)}{\delta t} = -(0.5)\sqrt{6.8144} = -1.3052 \frac{\delta^2 y(1.5)}{\delta t^2} = -\frac{0.5}{2} (6.8144)^{-\frac{1}{2}} = - 0.09576 \frac{\delta^3 y(1.5)}{\delta t^3} = \frac{0.5}{4} (6.8144)^{-\frac{3}{2}} = 0.007026 P_{2}(2) = 6.8144 - 1.3052 (2-1.5) - \frac{0.09576}{2}(2-1.5)^2 P_{2}(2) = 6.1498

Error orden de:

Error = \frac{0.007026}{3!}(2-1.5)^3 = 1.4637 . 10^{-4}

Se estima que el próximo término pasa debajo de 6 pies.
Por lo que estima esperar entre 2 y 2.5 minutos.

resultados usando el algoritmo:

ti, p_i,  error
[[0.00000000e+00 9.00000000e+00 0.00000000e+00]
 [5.00000000e-01 8.23958333e+00 9.64506173e-05]
 [1.00000000e+00 7.51107974e+00 1.10105978e-04]
 [1.50000000e+00 6.81451855e+00 1.26507192e-04]
 [2.00000000e+00 6.14993167e+00 1.46391550e-04]
 [2.50000000e+00 5.51735399e+00 1.70751033e-04]]

Algoritmo en Python

# 3Eva_IIT2018_T2 Drenar tanque cilíndrico
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
y0 = 9
t0 = 0
buscar = 6
k = 0.5
h = 0.5

dy  = lambda t,y: -k*np.sqrt(y)
d2y = lambda t,y: -(k/2)*(y**(-1/2))
d3y = lambda t,y: (k/4)*(y**(-3/2))

# PROCEDIMIENTO
resultado = [[t0,y0,0]]
yi = y0
ti = t0
while not(yi<buscar):
    ti = ti+h
    dyi = dy(ti,yi)
    d2yi = d2y(ti,yi)
    d3yi = d3y(ti,yi)
    p_i = yi +dyi*(h) + (d2yi/2)*(h**2)
    errado = (d3yi/6)*(h**3)
    yi = p_i
    resultado.append([ti,p_i,errado])
resultado = np.array(resultado)

# SALIDA
print('ti, p_i,  error')
print(resultado)

# Grafica
plt.plot(resultado[:,0],resultado[:,1])
plt.ylabel('nivel de agua')
plt.xlabel('tiempo')
plt.grid()
plt.show()