s2Eva_2022PAOI_T3 EDP parabólica barra enfriada en centro

Ejercicio: 2Eva_2022PAOI_T3 EDP parabólica barra enfriada en centro

Para la ecuación dada con Δx = 1/3, Δt = 0.02, en una revisíón rápida para cumplir la convergencia dt<dx/10, condición que debe verificarse con la expresión obtenida para λ al desarrollar el ejercicio.

\frac{\partial U}{\partial t} - \frac{1}{9} \frac{\partial ^2 U}{\partial x^2} = 0 0 \leq x \leq 2, t>0

literal a. grafica de malla

literal b. Ecuaciones de diferencias divididas a usar

\frac{\partial U}{\partial t} - \frac{1}{9} \frac{\partial ^2 U}{\partial x^2} = 0 \frac{\partial ^2 U}{\partial x^2} = 9 \frac{\partial U}{\partial t} \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{(\Delta x)^2} = 9 \frac{u_{i,j+1}-u_{i,j}}{\Delta t}

se agrupan las constantes,

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

literal d Determine el valor de λ

\lambda = \frac{\Delta t}{9(\Delta x)^2} =\frac{0.02}{9(1/3)^2} = 0.02

valor de λ que es menor que 1/2, por lo que el método converge.

continuando luego con la ecuación general,

\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 u[i,j] + \lambda u[i+1,j] \Big) = u[i,j+1]-u[i,j]

literal c. Encuentre las ecuaciones considerando las condiciones dadas en el problema.

\lambda u[i-1,j]+(1-2 \lambda ) u[i,j] + \lambda u[i+1,j] = u[i,j+1]

el punto que no se conoce su valor es u[i,j+1] que es la ecuación buscada.

u[i,j+1] = \lambda u[i-1,j]+(1-2 \lambda ) u[i,j] + \lambda u[i+1,j]

literal e iteraciones

iteración  i=1, j=0

u[1,1] = \lambda u[0,0]+(1-2 \lambda ) u[1,0] + \lambda u[2,0] u[1,1] =0.02 \cos \Big( \frac{\pi}{2}(0-3)\Big) + (1-2(0.02) ) \cos \Big( \frac{\pi}{2}\big(\frac{1}{3}-3\big)\Big) + 0.02 \cos \Big( \frac{\pi}{2}\big( \frac{2}{3}-3\big) \Big) u[1,1] =0.02(0)+(0.96)(-0.5)+0.02(-0.8660)=-0.4973

iteración  i=2, j=0

u[2,1] = \lambda u[1,0]+(1-2 \lambda ) u[2,0] + \lambda u[3,0] u[2,1] = 0.02 \cos \Big( \frac{\pi}{2}(\frac{1}{3}-3)\Big) + (1-2(0.02) ) \cos \Big( \frac{\pi}{2}(\frac{2}{3}-3)\Big)+ + 0.02 \cos \Big( \frac{\pi}{2}\big(\frac{3}{3}-3\big)\Big) u[2,1] = 0.02 (-0.5) + (0.96 ) (-0.866025) + 0.02 (-1) =-0.8614

iteración  i=3, j=0

u[3,1] = \lambda u[2,0]+(1-2 \lambda ) u[3,0] + \lambda u[4,0] u[3,1] = 0.02 \cos \Big( \frac{\pi}{2}\big( \frac{2}{3}-3\big)\Big)+(1-2 (0.02) ) \cos \Big( \frac{\pi}{2}(1-3)\Big) + + 0.02 \cos \Big( \frac{\pi}{2}\big(\frac{4}{3}-3\big)\Big) u[3,1] = 0.02 (-0.866025)+(0.96 ) (-1) + 0.02 (-0,866025) = -0,9946

literal f

la cotas de errores de truncamiento en la ecuación corresponden a segunda derivada O(hx2) y el de primera derivada O(ht), al reemplazar los valores será la suma}

O(hx2) + O(ht) = (1/3)2 + 0.02 = 0,1311

literal g

Resultados usando el algoritmo en Python

Tabla de resultados
[[ 0.      0.      0.      0.      0.      0.      0.      0.      0.       0.    ]
 [-0.5    -0.4973 -0.4947 -0.492  -0.4894 -0.4867 -0.4841 -0.4815 -0.479   -0.4764]
 [-0.866  -0.8614 -0.8568 -0.8522 -0.8476 -0.8431 -0.8385 -0.8341 -0.8296  -0.8251]
 [-1.     -0.9946 -0.9893 -0.984  -0.9787 -0.9735 -0.9683 -0.9631 -0.9579  -0.9528]
 [-0.866  -0.8614 -0.8568 -0.8522 -0.8476 -0.8431 -0.8385 -0.8341 -0.8296  -0.8251]
 [-0.5    -0.4973 -0.4947 -0.492  -0.4894 -0.4867 -0.4841 -0.4815 -0.479   -0.4764]
 [ 0.      0.      0.      0.      0.      0.      0.      0.      0.       0.    ]]

Instrucciones en Python

# EDP parabólicas d2u/dx2  = K du/dt
# método explícito, usando diferencias finitas
# 2Eva_2022PAOI_T3 EDP parabólica barra enfriada en centro
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
# Valores de frontera
Ta = 0
Tb = 0
T0 = lambda x: np.cos((np.pi/2)*(x-3))
# longitud en x
a = 0.0
b = 2.0
# Constante K
K = 9
# Tamaño de paso
dx = 1/3
dt = 0.02
tramos = int(np.round((b-a)/dx,0))
muestras = tramos + 1
# iteraciones en tiempo
n = 10

# PROCEDIMIENTO
# iteraciones en longitud
xi = np.linspace(a,b,muestras)
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(xi[1:ultimox])
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('t[j]')
plt.title('Solución EDP parabólica')
plt.show()

s2Eva_2022PAOI_T2 EDO de circuito RLC con interruptor intermedio

Ejercicio: 2Eva_2022PAOI_T2 EDO de circuito RLC con interruptor intermedio

La corriente del inductor y(t) para t≥0 se deriva para tener la expresión solo derivadas:

\frac{\delta}{\delta t}y(t) + 2 y(t) + 5 \int_{-\infty}^t y(\tau) \delta \tau = 10 \mu(t)

Para t>0 que es donde transcurre el experimento, el escalón es una constante, se tiene que:

\frac{\delta ^2}{\delta t^2}y(t) + 2 \frac{\delta}{\delta t}y(t) + 5 y(t) = 0

tomando las condiciones iniciales dadas para t=0, y(0)=2, y'(0)=–4

literal a

EL resultadoes perado es el planteamiento del problema. Se reescribe la ecuación con la nomenclatura simplificada y se resordena segun el modelo del método:

y'' = - 2y' - 5 y

luego se sustituye la variable y se convierte a las ecuaciones:

z =y' = f_x(t,y,z) z' = - 2z - 5 y = g_z(t,y,z)

se usa una tabla para llevar el registro de operaciones:

Se plantea las operaciones:

K1y = h * f(ti,yi,zi)
K1z = h * g(ti,yi,zi)

K2y = h * f(ti+h, yi + K1y, zi + K1z)
K2z = h * g(ti+h, yi + K1y, zi + K1z)

yi = yi + (K1y+K2y)/2
zi = zi + (K1z+K2z)/2
ti = ti + h

literal b

El resultado esperado es la aplicación correcta de los valores en las expresiones para al menos tres iteraciones usando h=0.01

itera = 0

K1y = 0.01 y'(0) = 0.01(-4) = -0.04 K1z = 0.01 (- 2z(0) - 5 y(0)) = 0.01(- 2(-4) - 5 (2)) = -0.02 K2y = 0.01 (-4-0.02) = -0.0402 K2z = 0.01 (-2(-4-0.02)-5(2-0.04)) = -0.0176 yi = yi + \frac{K1y+K2y}{2} = 2+\frac{-0.04-0.0402} {2} = 1.9599 zi = zi + \frac{K1z+K2z}{2} = -4 +\frac{-0.02-0.0176}{2} = -4.0188 ti = ti + h = 0+0.01 = 0.01

itera = 1

K1y = 0.01(-4.0188) = -0.040188 K1z = 0.01(- 2(-4.0188) - 5 (1.9599)) = -0.0176 K2y = 0.01 (-4.0188-0.0176) = -0.0403 K2z = 0.01 (-2(-4.0188-0.0176)-5(1.9599-0.040188)) = -0.0152 yi = 1.9599 +\frac{-0.040188-0.0403} {2} = 1.9196 zi = -4.0188 +\frac{-0.0176-0.0152}{2} = -4.0352 ti = ti + h = 0.01+0.01 = 0.02

itera = 2

K1y = 0.01(-4.0352) = -0.040352 K1z = 0.01(- 2(-4.0352) - 5 (1.9196)) = -0.0152 K2y = 0.01 (-4.0352-0.0152) = -0.0405 K2z = 0.01 (-2(-4.0352-0.0152)-5(1.9196-0.040352)) = -0.0129 yi = 1.9196 +\frac{-0.040352-0.0405} {2} =1.8792 zi = -4.0352 +\frac{-0.0152-0.0129}{2} = -4.0494 ti = ti + h = 0.02+0.01 = 0.03

Resultados con el algoritmo en Python

   ti,   yi,    zi,      K1y,    K1z,    K2y,     K2z
[[ 0.00  2.0000 -4.0000  0.0000  0.0000  0.0000   0.0000]
 [ 0.01  1.9599 -4.0188 -0.0400 -0.0200 -0.0402  -0.0176]
 [ 0.02  1.9196 -4.0352 -0.0401 -0.0176 -0.0403  -0.0152]
 [ 0.03  1.8792 -4.0494 -0.0403 -0.0152 -0.0405  -0.0129]
...

Literal c

Runge-Kutta 2do Orden tiene error de truncamiento O(h3)

por lo que el error está en el orden de (0.01)3 = 0.000001


Literal d

Se requiere presentar el resultado para el intervalo t entre [0,5]. Siendo el tamaño de paso h=0.01 que es pequeño, se requieren realizar (5-0)/0.01=500 iteraciones, que es más práctico realizarlas usando el algoritmo.

Instrucciones en Python

# Respuesta a entrada cero
# solucion para (D^2+ D + 1)y = 0
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,7),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,z0,0,0,0,0]
    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,K1y,K1z,K2y,K2z]
    return(estimado)

# PROGRAMA
f = lambda t,y,z: z
g = lambda t,y,z: -2*z -5*y + 0

t0 = 0
y0 = 2
z0 = -4

h = 0.01
tn = 5
muestras = int((tn-t0)/h)

tabla = rungekutta2_fg(f,g,t0,y0,z0,h,muestras)
ti = tabla[:,0]
yi = tabla[:,1]
zi = tabla[:,2]

# SALIDA
np.set_printoptions(precision=4)
print('ti, yi, zi, K1y, K1z, K2y, K2z')
print(tabla)

# GRAFICA
plt.plot(ti,yi, color = 'orange', label='y_RK(t)')
plt.ylabel('y(t)')
plt.xlabel('t')
plt.title('y(t) con Runge-Kutta 2do Orden d2y/dx2 ')
plt.legend()
plt.grid()
plt.show()

Nota: En el curso TELG1001 Señales y Sistemas, la solución se realiza con Transformadas de Laplace

s2Eva_2022PAOI_T1 Comparar integrales numéricos Simpson y Cuadratura de Gauss

Ejercicio: 2Eva_2022PAOI_T1 Comparar integrales numéricos Simpson y Cuadratura de Gauss

Literal a. Integral con Simpson 1/3

Para la ecuación en el intervalo x entre [0,3] aplicando dos veces la fórmula en el intervalo se requieren al menos dos tramos cada Simpson de 1/3. Por lo que la cantidad de tramos es 4 (dos por cada formula, y dos fórmulas) que corresponden a 5 muestras.

El tamaño de paso se calcula como:

h = \frac{b-a}{tramos}=\frac{3-0}{4} = \frac{3}{4} = 0.75

representados en una gráfica como

A = \int_0^3 \frac{e^x \sin(x)}{1+x^2} \delta x

con lo que se define la función del integral f(x)

f(x) = \frac{e^x \sin(x)}{1+x^2}

Con lo que aplicando la fórmula se puede obtener los valores de las muestras:

xi= [0. 0.75   1.5    2.25   3.    ]
fi= [0. 0.9235 1.3755 1.2176 0.2834]

Nota: realizar las expresiones completas para las fórmulas si no adjunta el algoritmo en Python

Aplicando Simpson de 1/3 en cada tramo se tiene:

A_s= \frac{1}{3} \Big( \frac{3}{4} \Big ) \Big( 0 + 4(0.9235) + 1.3755 \Big) + + \frac{1}{3} \Big( \frac{3}{4} \Big ) \Big( 1.3755 + 4(1.2176) +0.2834 \Big) A_s = 2.8998

Literal b. Integral con Cuadratura de Gauss

Para usar dos veces el método de Cuadratura de Gauss se usan dos intervalos, con lo que las muestras en x serán:

xj= [0. 1.5 3. ]

se calculan los valores para el tramo [0, 1.5]:

x_{a1} = \frac{0+1.5}{2} - \frac{1}{\sqrt{3}}\frac{1.5-0}{2} = 0.3169 x_{b1} = \frac{0+1.5}{2} + \frac{1}{\sqrt{3}}\frac{1.5-0}{2} = 1.1830 A_{g1} =\frac{1.5-0}{2} \Big( \frac{f(0.3169)+f(1.1830)}{2} \Big) = 1.2361

se calculan los valores para el tramo [1.5, 3]:

x_{a2} = \frac{1.5+3}{2} - \frac{1}{\sqrt{3}}\frac{3-1.5}{2} = 1.8169 x_{b2} = \frac{1.5+3}{2} + \frac{1}{\sqrt{3}}\frac{3-1.5}{2} = 2.6830 A_{g2} =\frac{3-1.5}{2} \Big( \frac{f(1.8169)+f(2.6830)}{2} \Big) = 1.6329

El total del integral para el intervalo  [0,3]

A_g = A_{g1} + A_{g2} = 2.8691

Al comparar los resultados entre los métodos del literal a y b

errado = |A_s - A_g| = 2.8998 - 2.8691 = 0.0307

Instrucciones integradas en Python

# 2Eva_2022PAOI_T1
# Comparar integrales numéricos Simpson
#   y Cuadratura de Gauss
import numpy as np
import matplotlib.pyplot as plt

# INGRESO
fx = lambda x: (np.exp(x)*np.sin(x))/(1+x**2)
a = 0
b = 3

# PROCEDIMIENTO
# Aplicando Simpson 1/3
tramos = 4
muestras = tramos+1
xi = np.linspace(a,b,muestras)
fi = fx(xi)
hs = xi[1]-xi[0]
As = (1/3)*(3/4)*(fi[0]+4*fi[1]+fi[2])
As = As + (1/3)*(3/4)*(fi[2]+4*fi[3]+fi[4])
erradoS = 2*(hs**5)

# Aplicando Cuadratura de Gauss
tramosG = 2
muestrasG = tramosG+1
xj = np.linspace(a,b,muestrasG)
hg = xj[1]-xj[0]

xa1 = (xj[0]+xj[1])/2 - (1/np.sqrt(3))*(xj[1]-xj[0])/2
xb1 = (xj[0]+xj[1])/2 + (1/np.sqrt(3))*(xj[1]-xj[0])/2
Ag1 = (hg/2)*(fx(xa1)+fx(xb1))

xa2 = (xj[1]+xj[2])/2 - (1/np.sqrt(3))*(xj[2]-xj[1])/2
xb2 = (xj[1]+xj[2])/2 + (1/np.sqrt(3))*(xj[2]-xj[1])/2
Ag2 = (hg/2)*(fx(xa2)+fx(xb2))

Ag = Ag1 + Ag2

# error entre métodos
errado = np.abs(As-Ag)

# SALIDA
print('xi=',xi)
print('fi=',fi)
print('A Simpson =', As)
print('error Truncamiento Simpson 2*(h^5):',
      erradoS)
print('Cuadratura de Gauss xa,xb por tramos:',
      [xa1,xb1,xa2,xb2])
print('  fx(xa),fx(xb) por tramos:',
      [fx(xa1),fx(xb1),fx(xa2),fx(xb2)])
print('  integral de cada tramo:', [Ag1,Ag2])
print('A Gauss =', Ag)
print('errado entre Simpson y Gauss',errado)

# Grafica con mejor resolucion
xk = np.linspace(a,b,5*tramos+1)
fk = fx(xk)

plt.plot(xk,fk)
plt.plot(xi,fi,'o')
for unx in xi:
    plt.axvline(unx,color='red',
                linestyle='dotted')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('(np.exp(x)*np.sin(x))/(1+x**2)')
plt.grid()
plt.show()

s2Eva_2021PAOII_T2 EDO – Embudos cónicos para llenar botellas

Ejercicio: 2Eva_2021PAOII_T2 EDO – Embudos cónicos para llenar botellas

literal a

La expresión dada en el enunciado para EDO, se reordena para definir la funcion a usar con Runge-Kutta:

\frac{\delta y(t)}{\delta t} + \frac{d^2}{4}\sqrt{2 g \text{ }y(t)}\Bigg[\frac{tan \theta}{y(t)} \Bigg]^2 = 0 \frac{\delta y(t)}{\delta t} = - \frac{d^2}{4}\sqrt{2 g \text{ }y(t)}\Bigg[\frac{tan \theta}{y(t)} \Bigg]^2

siendo h = 0.5,  con y(0) = 0.15 m y d= 0.01 m ajustando las unidades de medida.

\frac{\delta y(t)}{\delta t} = - \frac{0.01^2}{4}\sqrt{2 (9.8) \text{ }y(t)}\Bigg[\frac{tan (\pi/4)}{y(t)} \Bigg]^2 \frac{\delta y(t)}{\delta t} = - (1.1068e-4) \sqrt{ y(t)}\Bigg[\frac{1}{y(t)} \Bigg]^2 \frac{\delta y(t)}{\delta t}= - (1.1068e-4) \frac{\sqrt{ y(t)}}{y(t)^2}

literal b

se inicia el cálculo del siguiente punto de la tabla

i t y
0 0 0.15
1 0.5 0.1490
2 1 0.1480
3 1.5 0.1471

i = 0

K_1 = h\Bigg(- (1.1068e-4) \frac{\sqrt{ y(t)}}{y(t)^2} \Bigg) K_1 = 0.5\Bigg(- (1.1068e-4) \frac{\sqrt{0.15}}{0.15^2}\Bigg) = -9.5258e-04 K_2 = h\Bigg(- (1.1068e-4) \frac{\sqrt{ y(t)+K_1}}{(y(t)+K_1)^2} \Bigg) K_2 = 0.5\Bigg(- (1.1068e-4) \frac{\sqrt{ 0.15+-9.5258e-04}}{(0.15-9.5258e-04)^2} \Bigg) K_2 = -9.6173e-04 y_1 = y_0 + \frac{K_1 + K_2}{2} y_1 = 0.15 + \frac{-9.5258e-04 -9.6173e-04}{2} = 0.149

i = 1

K_1 = 0.5\Bigg(- (1.1068e-4) \frac{\sqrt{0.149}}{0.149^2}\Bigg) =-9.6177e-04 K_2 = 0.5\Bigg(- (1.1068e-4) \frac{\sqrt{ 0.149-9.7120e-04}}{(0.149-9.7120e-04)^2} \Bigg) K_2 = -9.7116e-04 y_2 = y_1 + \frac{-9.6177e-04 + -9.7116e-04}{2} = 0.1480

i = 2

K_1 = 0.5\Bigg(- (1.1068e-4) \frac{\sqrt{0.1480}}{0.1480^2}\Bigg) = -9.7120e-04 K_2 = 0.5\Bigg(- (1.1068e-4) \frac{\sqrt{ 0.1480-9.7120e-04}}{(0.1480-9.7120e-04)^2} \Bigg) K_2= -9.8084e-04 y_3 = y_2 + \frac{-9.7120e-04 + -9.8084e-04}{2} = 0.1471

literal c

Resultados usando Algoritmo, se encuentra que el embudo se vacia entre 3.15 y 3.20 segundos

 [ t , y , K1 , K2 ]
[[ 0.0000e+00  1.5000e-01  0.0000e+00  0.0000e+00]
 [ 5.0000e-01  1.4904e-01 -9.5258e-04 -9.6173e-04]
 [ 1.0000e+00  1.4808e-01 -9.6177e-04 -9.7116e-04]
 [ 1.5000e+00  1.4710e-01 -9.7120e-04 -9.8084e-04]
 [ 2.0000e+00  1.4611e-01 -9.8088e-04 -9.9078e-04]
 [ 2.5000e+00  1.4512e-01 -9.9083e-04 -1.0010e-03]
...
[ 3.1000e+01  2.8617e-02 -7.5631e-03 -1.0583e-02]
 [ 3.1500e+01  1.0620e-02 -1.1431e-02 -2.4563e-02]
 [ 3.2000e+01         nan -5.0566e-02         nan]
 [ 3.2500e+01         nan         nan         nan]

Instrucciones Python

# 2Eva_2021PAOII_T2 EDO – Embudos cónicos para llenar botellas
import numpy as np

def rungekutta2(d1y,x0,y0,h,muestras):
    tamano   = muestras + 1
    estimado = np.zeros(shape=(tamano,4),dtype=float)
    # incluye el punto [x0,y0,K1,K2]
    estimado[0] = [x0,y0,0,0]
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1 = h * d1y(xi,yi)
        K2 = h * d1y(xi+h, yi + K1)

        yi = yi + (K1+K2)/2
        xi = xi + h
        
        estimado[i] = [xi,yi,K1,K2]
    return(estimado)

# INGRESO
d = 0.01
theta = np.pi/4
g = 9.8
d1y = lambda t,y: -(d**2)/4*np.sqrt(2*g*y)*(np.tan(theta)/y)**2

t0 = 0
y0 = 0.15
h  = 0.5
muestras = 70

# PROCEDIMIENTO
tabla = rungekutta2(d1y,t0,y0,h,muestras)

# SALIDA
np.set_printoptions(precision=4)
print('[ t , y , K1 , K2 ]')
print(tabla)

s2Eva_2021PAOII_T1 Promedio de precipitación de lluvia en área

Ejercicio: 2Eva_2021PAOII_T1 Promedio de precipitación de lluvia en área

Los datos de la tabla corresponden a las medidas de precipitación lluviosa en cada cuadrícula. El área observada es de 300×250, los tramos horizontales 6 y los verticales 4.

Los tamaños de paso son:
h = 300/6 = 50
k = 250/4 = 62.5

f(xi,yj)
i \ j 1 2 3 4 5 6
1 0.02 0.36 0.82 0.65 1.7 1.52
2 3.15 3.57 6.25 5 3.88 1.8
3 0.98 0.98 2.4 1.83 0.04 0.01
4 0.4 0.04 0.03 0.03 0.01 0.08

Para obtener el integral doble, usando la forma compuesta de Simpson, primero desarrolla por filas.

i=1

I_1 = \frac{3}{8}(50) ( 0.02+3(0.36)+3(0.82)+0.65) + \frac{1}{3}(50) (0.65+4(1.7)+1.52) = 228.43

i=2

I_2 = \frac{3}{8}(50) ( 3.15+3(3.57)+3(6.25)+5) + \frac{1}{3}(50) (5+4(3.88)+1.8) = 1077.18

i=3

I_3 = \frac{3}{8}(50) ( 0.98+3(0.98)+3(2.4)+1.83) + \frac{1}{3}(50) (1.83+4(0.04)+0.01) = 276.14

i=4

I_4 = \frac{3}{8}(50) ( 0.4+3(0.04)+3(0.03)+0.03) + \frac{1}{3}(50) (0.03+4(0.01)+0.08) = 14.5

con los resultados, luego se desarrolla por columnas:

columnas = [ 228.43,  1077.18,  276.14,  14.5 ]

I_{total} = \frac{3}{8}(62.5) ( 228.43+3(1077.18)+3(276.14)+14.5) I_{total} = 100850.09

para el promedio se divide el integral total para el área de rectángulo.

Precipitacion_{promedio} = \frac{100850.09}{(300)(250)} = 1.34

Solución con algoritmo

Resultados usando algoritmos con Python:

integral por fila: [ 228.4375     1077.1875      276.14583333   14.5       ]
integral total: 100850.09765625
promedio precipitación:  1.34466796875
>>>

Instrucciones con Python

# 2Eva_2021PAOII_T1 Promedio de precipitación de 
#  lluvia en área
import numpy as np

def simpson_compuesto(vector,h):
    m = len(vector)
    suma = 0
    j = 0
    while (j+3)<m:  # Simpson 3/8
        f = vector[j:j+4]
        s = (3*h/8)*(f[0]+3*f[1]+3*f[2]+f[3])
        suma = suma + s
        j = j + 3
    while (j+2)<m: # Simpson 1/3
        f = vector[j:j+3]
        s = (h/3)*(f[0]+4*f[1]+f[2])
        suma = suma + s
        j = j + 2
    while (j+1)<m: # Trapecio
        f = vector[j:j+2]
        s = (h/2)*(f[0]+f[1])
        suma = suma + s
        j = j + 1
    return(suma)
    

# INGRESO
A = [[0.02, 0.36, 0.82, 0.65, 1.7 , 1.52],
     [3.15, 3.57, 6.25, 5.  , 3.88, 1.8 ],
     [0.98, 0.98, 2.4 , 1.83, 0.04, 0.01],
     [0.4 , 0.04, 0.03, 0.03, 0.01, 0.08]]
base = 300
altura = 250

h = base/6
k = altura/4

# PROCEDIMIENTO
A = np.array(A)
[n,m] = np.shape(A)
columna = np.zeros(n,dtype=float)
for i in range(0,n,1):
    vector = A[i]
    integralfila = simpson_compuesto(vector,h)
    columna[i] = integralfila

total = simpson_compuesto(columna,k)
promedio = total/(base*altura)
    
# SALIDA
print('integral por fila:', columna)
print('integral total:', total)
print('promedio precipitación: ', promedio)

s2Eva_2021PAOII_T3 EDP – Línea de transmisión sin pérdidas

Ejercicio: 2Eva_2021PAOII_T3 EDP – Línea de transmisión sin pérdidas

Desarrollo para determinar el Voltaje  v(t),

\frac{\partial ^2 V}{\partial x^2} =LC \frac{\partial ^2 V}{\partial t^2}

0 < x < 200
t>0

Seleccionando diferencias divididas centradas para x y t

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

agrupando las constantes como lambda,

\frac{(\Delta t)^2}{LC}\frac{v_{i,j+1}-2v_{i,j}+v_{i,j-1}}{(\Delta x)^2} =LC \frac{v_{i+1,j}-2v_{i,j}+v_{i-1,j}}{(\Delta t)^2} \frac{(\Delta t)^2}{LC} \frac{(\Delta t)^2}{LC(\Delta x)^2} \Big[ v_{i,j+1} -2v_{i,j}+v_{i,j-1} \Big] =v_{i+1,j}-2v_{i,j}+v_{i-1,j} \lambda = \frac{(\Delta t)^2}{LC(\Delta x)^2} \lambda = \frac{(0.1)^2}{(0.1)(0.3)(10)^2} = 0.003333

con lo que se comprueba que λ ≤ 1, por lo que el algoritmo es convergente.

Para reducir la cantidad de muestras, se iguala λ = 1, manteniendo el valor de Δx y calculando Δt

\Delta t = \Delta x \sqrt{LC} \Delta t = 10 \sqrt{0.1(0.3)} = 1.7320

continuando con el valor de λ en la ecuacion general,

\lambda \Big[ v_{i,j+1} -2v_{i,j}+v_{i,j-1} \Big] =v_{i+1,j} -2v_{i,j}+v_{i-1,j} \lambda v_{i,j+1} +(-2\lambda +2) v_{i,j} + \lambda v_{i,j-1} = v_{i+1,j}+v_{i-1,j}

usando la condición inicial, para j = 0, y diferencias divididas centradas:

\frac{\partial V}{\partial t}(x,0) = 0 \frac{v_{i,1}-v_{i,-1}}{2\Delta t} = 0 v_{i,1}-v_{i,-1} = 0 v_{i,1} = v_{i,-1}

y en la ecuacion del problema, con j = 0

\lambda v_{i,1} +(-2\lambda +2) v_{i,0} + \lambda v_{i,-1} = v_{i+1,0}+v_{i-1,0} 2 \lambda v_{i,1} +(-2\lambda +2) v_{i,0} = v_{i+1,0}+v_{i-1,0} 2 \lambda v_{i,1} = (2\lambda -2) v_{i,0} + v_{i+1,0}+v_{i-1,0} v_{i,1} = \frac{1}{2 \lambda}\Bigg[ (2\lambda -2) v_{i,0} + v_{i+1,0}+v_{i-1,0} \Bigg]

si λ=1,

v_{i,1} = \frac{v_{i+1,0}+v_{i-1,0}}{2 }

con j =0, i =1

v_{1,1} = \frac{v_{2,0}+v_{0,0}}{2 } =\frac{1}{2}\Big( 110 \sin \frac{\pi(0+2\Delta x) }{200}+0 \Big) =17

con j =0, i =2

v_{2,1} = \frac{v_{3,0}+v_{1,0}}{2 } v_{2,1} = \frac{1}{2}\Big( 110 \sin \frac{\pi(30) }{200}+110 \sin \frac{\pi(10) }{200} \Big)=33.57

con j =0, i =3

v_{3,1} = \frac{v_{4,0}+v_{2,0}}{2} v_{2,1} = \frac{1}{2}\Big( 110 \sin \frac{\pi(40) }{200}+110 \sin \frac{\pi(20) }{200} \Big)=49.32

para j>0, se dispone de los valores de los puntos dentro del rombo, excepto v[i,j+1]

\lambda v_{i,j+1} +(-2\lambda +2) v_{i,j} + \lambda v_{i,j-1} = v_{i+1,j}+v_{i-1,j} \lambda v_{i,j+1} = (2\lambda -2) v_{i,j} - \lambda v_{i,j-1} + v_{i+1,j}+v_{i-1,j} v_{i,j+1} = \frac{1}{\lambda}\Big( (2\lambda -2) v_{i,j} - \lambda v_{i,j-1} + v_{i+1,j}+v_{i-1,j} \Big)

teniendo como valores de las filas:

u[:,1] = [  0.  ,  17.  ,  33.57,  49.32,  63.86,  76.82,  87.9 , ...  ])
u[:,0] = [  0.  ,  17.21,  33.99,  49.94,  64.66,  77.78,  88.99, ...  ])

con j = 1 ; i = 1 ; (λ=1)

v_{1,2} = \frac{1}{\lambda}\Big( (2\lambda -2) v_{1,1} - \lambda v_{1,0} + v_{2,1}+v_{0,1} \Big) v_{1,2} = \frac{1}{\lambda}\Big( (2\lambda-2) 17-\lambda 17.21+ 33.57+0\Big)=16.37

con j = 1 ; i = 2 ; (λ=1)

v_{2,2} = \frac{1}{\lambda}\Big( (2\lambda -2) v_{2,1} - \lambda v_{2,0} + v_{3,1}+v_{1,1} \Big) v_{2,2} = \frac{1}{\lambda}\Big( (2\lambda-2) 33.57-\lambda 33.99+ 49.32+17\Big)=32.33

con j = 1 ; i = 3 ; (λ=1)

v_{2,2} = \frac{1}{\lambda}\Big( (2\lambda -2) v_{3,1} - \lambda v_{3,0} + v_{4,1}+v_{2,1} \Big) v_{2,2} = \frac{1}{\lambda}\Big( (2\lambda-2) 49.32-\lambda 49.94 + 63.86+33.57\Big)=47.49

usando el algoritmo con los valores y ecuaciones encontradas se obtiene:

21 41
xi =
[  0  10  20  30  40  50  60  70  80  90 100 110 120 130 140 150 160 170
 180 190 200]
tj =
[ 0.    1.73  3.46  5.2   6.93  8.66 10.39 12.12 13.86 15.59 17.32 19.05
 20.78 22.52 24.25 25.98 27.71 29.44 31.18 32.91 34.64 36.37 38.11 39.84
 41.57 43.3  45.03 46.77 48.5  50.23 51.96 53.69 55.43 57.16 58.89 60.62
 62.35 64.09 65.82 67.55 69.28]
matriz u =
....

Instrucciones en Python

# 2Eva_2021PAOII_T2 EDP – Línea de transmisión sin pérdidas
# Ecuaciones Diferenciales Parciales
# Hiperbólica. Método explicito
import numpy as np

# INGRESO
L = 0.1
C = 0.3

x0 = 0
xn = 200 # Longitud de cable

vx0 = lambda x: 110*np.sin(np.pi*x/200)

y0 = 0
yn = 0 # Puntos de amarre

t0 = 0
tn = 68

# discretiza
dx = 10
dt = dx*np.sqrt(L*C)

# PROCEDIMIENTO
xi = np.arange(x0,xn+dx,dx)
tj = np.arange(0,tn+dt,dt)
n = len(xi)
m = len(tj)
print(n,m)

u = np.zeros(shape=(n,m),dtype=float)
u[:,0] = vx0(xi)
u[0,:] = y0
u[n-1,:] = yn

# calcula lamb
lamb = (dt**2)/(L*C*(dx**2))

# Aplicando condición inicial
j = 0
for i in range(1,n-1,1):
    u[i,j+1] = ((2*lamb-2)*u[i,j]+u[i+1,j]+u[i-1,j])/(2*lamb)
# Para los otros puntos
for j in range(1,m-1,1):
    for i in range(1,n-1,1):
        u[i,j+1] = (1/lamb)*((2*lamb-2)*u[i,j]-lamb*u[i,j-1]+u[i+1,j]+u[i-1,j])

# SALIDA
np.set_printoptions(precision=2)
print('xi =')
print(xi)
print('tj =')
print(tj)
print('matriz u =')
print(u)

# GRAFICA
import matplotlib.pyplot as plt

for j in range(0,m,1):
    y = u[:,j]
    plt.plot(xi,y)

plt.title('EDP hiperbólica')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

# **** GRÁFICO CON ANIMACION ***********
import matplotlib.animation as animation

# Inicializa parametros de trama/foto
retardo = 70   # milisegundos entre tramas
tramas  = m
maximoy = np.max(np.abs(u))
figura, ejes = plt.subplots()
plt.xlim([x0,xn])
plt.ylim([-maximoy,maximoy])

# lineas de función y polinomio en gráfico
linea_poli, = ejes.plot(xi,u[:,0], '-')
plt.axhline(0, color='k')  # Eje en x=0

plt.title('EDP hiperbólica')
# plt.legend()
# txt_x = (x0+xn)/2
txt_y = maximoy*(1-0.1)
texto = ejes.text(x0,txt_y,'tiempo:',
                  horizontalalignment='left')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()

# Nueva Trama
def unatrama(i,xi,u):
    # actualiza cada linea
    linea_poli.set_ydata(u[:,i])
    linea_poli.set_xdata(xi)
    linea_poli.set_label('tiempo linea: '+str(i))
    texto.set_text('tiempo['+str(i)+']')
    # color de la línea
    if (i<=9):
        lineacolor = 'C'+str(i)
    else:
        numcolor = i%10
        lineacolor = 'C'+str(numcolor)
    linea_poli.set_color(lineacolor)
    return linea_poli, texto

# Limpia Trama anterior
def limpiatrama():
    linea_poli.set_ydata(np.ma.array(xi, mask=True))
    linea_poli.set_label('')
    texto.set_text('')
    return linea_poli, texto

# Trama contador
i = np.arange(0,tramas,1)
ani = animation.FuncAnimation(figura,
                              unatrama,
                              i ,
                              fargs=(xi,u),
                              init_func=limpiatrama,
                              interval=retardo,
                              blit=True)
# Graba Archivo video y GIFAnimado :
# ani.save('EDP_hiperbólica.mp4')
ani.save('EDP_LineaSinPerdidas.gif', writer='imagemagick')
plt.draw()
plt.show()

s2Eva_2021PAOI_T3 EDP Elíptica con valores en la frontera f(x) g(y)

Ejercicio: 2Eva_2021PAOI_T3 EDP Elíptica con valores en la frontera f(x) g(y)

Dada la EDP elíptica

\frac{\partial ^2 u}{\partial x^2} +\frac{\partial^2 u}{\partial y^2} = 0 0 \lt x \lt \frac{1}{2}, 0 \lt y\lt \frac{1}{2}

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

 


\frac{u[i-1,j]-2u[i,j]+u[i+1,j]}{\Delta x^2} + + \frac{u[i,j-1]-2u[i,j]+u[i,j+1]}{\Delta y^2} =

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

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

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

\lambda= \frac{\Delta y^2}{\Delta x^2} = 1
u[i-1,j]-2u[i,j]+u[i+1,j] + + u[i,j-1]-2u[i,j]+u[i,j+1] = 0
u[i-1,j]-4u[i,j]+u[i+1,j] + + u[i,j-1]+u[i,j+1] = 0

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

En cada iteración se requiere el uso de los valores en la frontera

u(x,0)=0, 0 \leq x \leq \frac{1}{2} u(0,y)=0 , 0\leq y \leq \frac{1}{2} u\Big(x,\frac{1}{2} \Big) = 200 x , 0 \leq x \leq \frac{1}{2} u\Big(\frac{1}{2} ,y \Big) = 200 y , 0 \leq y \leq \frac{1}{2}

Iteraciones

i=1, j=1

u[0,1]-4u[1,1]+u[2,1] + + u[1,0]+u[1,2] = 0

usando los valores en la frontera,

0-4u[1,1]+u[2,1] + 0+u[1,2] = 0 -4u[1,1]+u[2,1] + u[1,2] = 0

i=2, j=1

u[1,1]-4u[2,1]+u[3,1] + + u[2,0]+u[2,2] = 0

usando los valores en la frontera,

u[1,1]-4u[2,1]+200(1/6) + 0+u[2,2] = 0 u[1,1]-4u[2,1] + u[2,2] = -200(1/6)

i=1, j=2

u[0,2]-4u[1,2]+u[2,2] + + u[1,1]+u[1,3] = 0 0 - 4u[1,2]+u[2,2] + u[1,1]+200\frac{1}{6} = 0 - 4u[1,2] + u[2,2]+u[1,1] = -200\frac{1}{6}

i=2, j=2

u[1,2]-4u[2,2]+u[3,2] + + u[2,1]+u[2,3] = 0 u[1,2]-4u[2,2]+200\frac{2}{6} + u[2,1]+200\frac{2}{6} = 0 u[1,2]-4u[2,2]+ u[2,1] = -(2)200\frac{2}{6}

Sistema de ecuaciones a resolver:

\begin{bmatrix} -4 & 1 &1&0\\1 &-4&0&1\\1&0&-4&1 \\0&1&1&-4\end{bmatrix} \begin{bmatrix} u[1,1]\\u[2,1] \\u[1,2]\\u[2,2] \end{bmatrix} = \begin{bmatrix} 0\\-200(1/6)\\-200(1/6)\\-200(4/6) \end{bmatrix}

Resolviendo el sistema se tiene:

[11.11111111 22.22222222 22.22222222 44.44444444]

Instrucciones Python

import numpy as np

A = np.array([[-4, 1, 1, 0],
              [ 1,-4, 0, 1],
              [ 1, 0,-4, 1],
              [ 0, 1, 1,-4]])

B = np.array([0,-200*(1/6),-200*(1/6),-200*(4/6)])

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

print(x)

s2Eva_2021PAOI_T1 Masa transportada por tubo

Ejercicio: 2Eva_2021PAOI_T1 Masa transportada por tubo

Las expresiones siguientes se usan dentro de la expresión del integral

Q(t)=9+4 \cos ^2 (0.4t) c(t)=5e^{-0.5t}+2 e^{-0.15 t} M = \int_{t_1}^{t_2} Q(t)c(t) dt

literal a

Usando los valores dados para el intervalo [2,8] con 6 tramos h = (8-2)/6 =1

Se usa los valores de cada ti en Se puede obtener una tabla de valores muestreados para integrar f(t) = Q(i)c(t)

[ti,	 Qi,	 Ci,	 fi]
[ 2.     10.9416  3.321  36.3374]
[ 3.      9.5252  2.3909 22.7739]
[ 4.      9.0034  1.7743 15.9747]
[ 5.      9.6927  1.3552 13.1352]
[ 6.     11.175   1.0621 11.8687]
[ 7.     12.5511  0.8509 10.6793]
[ 8.     12.9864  0.694   9.0121]

Para el integral se usan los valores por cada dos tramos

I\cong \frac{1}{3}[36.3374+4(22.7739) + 15.9747] + \frac{1}{3}[15.9747+4(13.1352) + 11.8687] + \frac{1}{3}[11.8687+4(10.6793) + 9.0121] I = 95.7965

literal b

L acota de error de truncamiento por cada fórmula usada, se estima como O(h5),

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

para un valor de z entre [a,b]

por lo que al usar 3 veces la formula de Simpson se podría estimar en:

error_{trunca} = 3(1^5/90) = 0.033

literal c

El resultado se puede mejorar de dos formas:

1. Dado que el número de tramos es múltiplo de 3, se puede cambiar la fórmula a Simpon de 3/8, que tendría una cota de error menor

2. Aumentar el número de tramos disminuyendo el valor de h para que el error disminuya. Por ejemplo si se reduce a 0.5, el error disminuye en el orden de 0.55

Podría recomendar la segunda opión, pues a pesar que se aumenta la cota de error por cada vez que se usa la fórmula, el error de cada una disminuye en ordenes de magnitud 0,03125


La gráfica del ejercicio es:

Instrucciones en Python

import numpy as np
import matplotlib.pyplot as plt

# INGRESO
Q = lambda t: 9+4*(np.cos(0.4*t)**2)
C = lambda t: 5*np.exp(-0.5*t)+2*np.exp(-0.15*t)

t1 = 2
t2 = 8
n  = 6

# PROCEDIMIENTO
muestras = n+1 
dt = (t2-t1)/n
ti = np.arange(t1,t2+dt,dt)
Qi = Q(ti)
Ci = C(ti)
fi = Qi*Ci

# integración con Simpson 1/3
h= dt
I13 = 0
for i in range(0,6,2):
    S13 = (h/3)*(fi[i]+4*fi[i+1]+fi[i+2])
    I13 = I13 + S13

# SALIDA
np.set_printoptions(precision=4)
print("[ti,\t Qi,\t Ci,\t fi]")
for i in range(0,muestras,1):
    print(np.array([ti[i],Qi[i],Ci[i],fi[i]]))
print("Integral S13: ",I13)
# grafica
plt.plot(ti,Qi, label = "Q(t)")
plt.plot(ti,Ci, label = "c(t)")
plt.plot(ti,fi, label = "f(t)")
plt.plot(ti,Qi,'.b')
plt.plot(ti,Ci,'.r')
plt.plot(ti,fi,'.g')
plt.xlabel("t")
plt.ylabel("f(t)=Q(t)*c(t)")
plt.show()

s2Eva_2021PAOI_T2 EDO para cultivo de peces

Ejercicio: 2Eva_2021PAOI_T2 EDO para cultivo de peces

Siendo la captura una constante mas una función periódica,

h(t) = a + b \sin (2 \pi t)

La ecuación EDO del ejercicio, junto a las constantes a=0.9 y b=0.75, r=1

\frac{\delta y(t)}{\delta t} = r y(t)-h(t)

se convierte en:

\frac{\delta y(t)}{\delta t} = (1) y(t)- \Big( 0.9 + .75 \sin (2 \pi t)\Big) \frac{\delta y(t)}{\delta t} = y(t)- 0.9 - .75 \sin (2 \pi t)

Considerando que la población inicial de peces es 1 o 100%, y(0)=1

literal a

h=1/12
tamano = muestras + 1
estimado = np.zeros(shape=(tamano,2),dtype=float)
estimado[0] = [0,1]
ti = 0
yi = 1
for i in range(1,tamano,1):
    K1 = 1/12 * d1y(ti,yi)
    K2 = 1/12 * d1y(ti+1/24, yi + K1/2)
    K3 = 1/12 * d1y(ti+1/24, yi + K2/2)
    K4 = 1/12 * d1y(ti+1/12, yi + K3)

    yi = yi + (1/6)*(K1+2*K2+2*K3 +K4)
    ti = ti + 1/12
        
    estimado[i] = [ti,yi]

literal b

iteración i=0

t(0) = 0

y(0) = 1

K1 = \frac{1}{12} \Big(1- 0.9 - .75 \sin (2 \pi 0)\Big) = 0,008333 K2 = \frac{1}{12} \Big(1- 0.9 - .75 \sin \Big(2 \pi (0+\frac{1}{12})\Big)\Big) = -0.02222 y(1) = 0 + \frac{0.008333+(-0.02222)}{2} = 0.9930 t(1) = 0 + \frac{1}{12} = \frac{1}{12}

iteración i=1

t(1) = \frac{1}{12}

y(1) = 0.9930

K1 = \frac{1}{12} \Big(0.9930 - 0.9 - .75 \sin \Big( 2 \pi\frac{1}{12}\Big)\Big) = -0.02349 K2 = \frac{1}{12} \Big(0.9930 - 0.9 - .75 \sin \Big(2 \pi (\frac{1}{12}+\frac{1}{12})\Big)\Big) = -0.04832 y(1) = 0.9930 + \frac{-0.02349+(-0.04832)}{2} = 0.9571 t(1) = \frac{1}{12} + \frac{1}{12} = \frac{2}{12}

iteración i=2

t(2) = \frac{2}{12}

y(1) = 0.9571

K1 = \frac{1}{12} \Big(0.9571 - 0.9 - .75 \sin \Big( 2 \pi\frac{2}{12}\Big)\Big) = -0.04936 K2 = \frac{1}{12} \Big(0.9571 - 0.9 - .75 \sin \Big(2 \pi (\frac{2}{12}+\frac{1}{12})\Big)\Big) = -0.06185 y(1) = 0.9571 + \frac{-0.04936+(-0.06185)}{2} = 0.9015 t(3) = \frac{2}{12} + \frac{1}{12} = \frac{3}{12}

literal c

Resultado del algoritmo, muestra que la estragegia de cosecha, en el tiempo no es sostenible, dado que la población de peces en el tiempo decrece.

estimado[xi,yi,K1,K2]
[[ 0.0000e+00  1.0000e+00  8.3333e-03 -2.2222e-02]
 [ 8.3333e-02  9.9306e-01 -2.3495e-02 -4.8330e-02]
 [ 1.6667e-01  9.5714e-01 -4.9365e-02 -6.1852e-02]
 [ 2.5000e-01  9.0153e-01 -6.2372e-02 -5.9196e-02]
 [ 3.3333e-01  8.4075e-01 -5.9064e-02 -4.1109e-02]
 [ 4.1667e-01  7.9066e-01 -4.0361e-02 -1.2475e-02]
 [ 5.0000e-01  7.6425e-01 -1.1313e-02  1.8994e-02]
 [ 5.8333e-01  7.6809e-01  2.0257e-02  4.4822e-02]
 [ 6.6667e-01  8.0063e-01  4.5845e-02  5.8039e-02]
 [ 7.5000e-01  8.5257e-01  5.8547e-02  5.5053e-02]
 [ 8.3333e-01  9.0937e-01  5.4907e-02  3.6606e-02]
 [ 9.1667e-01  9.5513e-01  3.5844e-02  7.5807e-03]
 [ 1.0000e+00  9.7684e-01  6.4031e-03 -2.4313e-02]
 [ 1.0833e+00  9.6788e-01 -2.5593e-02 -5.0602e-02]
 [ 1.1667e+00  9.2978e-01 -5.1645e-02 -6.4322e-02]
 [ 1.2500e+00  8.7180e-01 -6.4850e-02 -6.1881e-02]
 [ 1.3333e+00  8.0844e-01 -6.1757e-02 -4.4027e-02]
 [ 1.4167e+00  7.5554e-01 -4.3288e-02 -1.5645e-02]
 [ 1.5000e+00  7.2608e-01 -1.4494e-02  1.5549e-02]
 [ 1.5833e+00  7.2661e-01  1.6800e-02  4.1077e-02]
 [ 1.6667e+00  7.5554e-01  4.2089e-02  5.3969e-02]
 [ 1.7500e+00  8.0357e-01  5.4464e-02  5.0630e-02]
 [ 1.8333e+00  8.5612e-01  5.0470e-02  3.1799e-02]
 [ 1.9167e+00  8.9725e-01  3.1021e-02  2.3563e-03]
 [ 2.0000e+00  9.1394e-01  1.1619e-03 -2.9991e-02]
 [ 2.0833e+00  8.9953e-01 -3.1289e-02 -5.6773e-02]
 [ 2.1667e+00  8.5550e-01 -5.7835e-02 -7.1028e-02]
 [ 2.2500e+00  7.9107e-01 -7.1578e-02 -6.9169e-02]
 [ 2.3333e+00  7.2069e-01 -6.9069e-02 -5.1948e-02]
 [ 2.4167e+00  6.6018e-01 -5.1235e-02 -2.4254e-02]
 [ 2.5000e+00  6.2244e-01 -2.3130e-02  6.1924e-03]
 [ 2.5833e+00  6.1397e-01  7.4142e-03  3.0909e-02]
 [ 2.6667e+00  6.3313e-01  3.1888e-02  4.2918e-02]
 [ 2.7500e+00  6.7053e-01  4.3378e-02  3.8619e-02]
 [ 2.8333e+00  7.1153e-01  3.8421e-02  1.8746e-02]
 [ 2.9167e+00  7.4012e-01  1.7926e-02 -1.1830e-02]
 [ 3.0000e+00  7.4317e-01  0.0000e+00  0.0000e+00]]

Instrucciones en Python

# EDO. Método de RungeKutta 2do Orden 
# estima la solucion para muestras espaciadas h en eje x
# valores iniciales x0,y0
# entrega arreglo [[x,y]]
import numpy as np

def rungekutta2(d1y,x0,y0,h,muestras):
    tamano   = muestras + 1
    estimado = np.zeros(shape=(tamano,4),dtype=float)
    # incluye el punto [x0,y0]
    estimado[0] = [x0,y0,0,0]
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1 = h * d1y(xi,yi)
        K2 = h * d1y(xi+h, yi + K1)

        yi = yi + (K1+K2)/2
        xi = xi + h
        estimado[i-1,2:]=[K1,K2]
        estimado[i] = [xi,yi,0,0]
    return(estimado)

# PROGRAMA PRUEBA
# Ref Rodriguez 9.1.1 p335 ejemplo.
# prueba y'-y-x+(x**2)-1 =0, y(0)=1

# INGRESO
# d1y = y' = f, d2y = y'' = f'
a =0.9; b=0.75; r=1
d1y = lambda t,y: r*y-(a+b*np.sin(2*np.pi*t))
x0 = 0
y0 = 1
h  = 1/12
muestras = 12*3

# PROCEDIMIENTO
puntosRK2 = rungekutta2(d1y,x0,y0,h,muestras)
xi = puntosRK2[:,0]
yiRK2 = puntosRK2[:,1]

# SALIDA
np.set_printoptions(precision=4)
print('estimado[xi,yi,K1,K2]')
print(puntosRK2)


# Gráfica
import matplotlib.pyplot as plt


plt.plot(xi[0],yiRK2[0],
         'o',color='r', label ='[x0,y0]')
plt.plot(xi[1:],yiRK2[1:],
         'o',color='m',
         label ='y Runge-Kutta 2 Orden')

plt.title('EDO: Solución con Runge-Kutta 2do Orden')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()