Sistemas EDO. modelo predador-presa

Referencia: Chapra 28.2 p831 pdf855, Rodriguez 9.2.1 p263
https://es.wikipedia.org/wiki/Ecuaciones_Lotka%E2%80%93Volterra
https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations

Modelos depredador- presa y caos. Ecuaciones Lotka-Volterra. En el sistema de ecuaciones:

\frac{dx}{dt} = ax - bxy \frac{dy}{dt} = -cy + dxy

x =  número de presas
y = depredadores
a = razón de crecimiento de la presa, (0.5)
c = razón de muerte del depredador (0.35)
b = efecto de la interacción depredador-presa sobre la muerte de la presa (0.7)
d = efecto de la interacción depredador-presa sobre el crecimiento del depredador, (0.35)

Los términos que multiplican xy hacen que las ecuaciones sean no lineales.


Para resolver el sistema, se plantean las ecuaciones de forma simplificada para el algoritmo:

f =  a x - b x y
g = -c y + d x y
con valores iniciales:
t = 0,x = 2, y = 1

Planteamiento que se ingresan al algoritmo con el algoritmo rungekutta2_fg(fx,gx,x0,y0,z0,h,muestras), propuesto en:

8.2.2 Runge-Kutta d2y/dx2

http://blog.espol.edu.ec/matg1013/8-2-2-runge-kutta-d2y-dx2/

Al ejecutar el algoritmo se obtienen los siguientes resultados:

 [ ti, xi, yi] 
[[  0.         2.         1.      ]
 [  0.5        1.754875   1.16975 ]
 [  1.         1.457533   1.302069]
 [  1.5        1.167405   1.373599]
 [  2.         0.922773   1.381103]
 [  2.5        0.734853   1.33689 ]
 [  3.         0.598406   1.258434]
 ...
 [ 49.         1.11309    1.389894]
 [ 49.5        0.876914   1.38503 ]
 [ 50.         0.698717   1.331107]
 [ 50.5        0.570884   1.246132]]
>>> 

Los resultados se pueden observar de diferentes formas:
a) Cada variable xi, yi versus ti, es decir cantidad de animales de cada especie durante el tiempo de observación

b) Independiente de la unidad de tiempo, xi vs yi, muestra la relación entre la cantidad de presas y predadores. Relación que es cíclica y da la forma a la gráfica.

Las instrucciones del algoritmo en Python usadas en el problema son:

# Modelo predador-presa de Lotka-Volterra
# Sistemas EDO con Runge Kutta de 2do Orden
import numpy as np

def rungekutta2_fg(f,g,t0,x0,y0,h,muestras):
    tamano = muestras +1
    tabla = np.zeros(shape=(tamano,3),dtype=float)
    tabla[0] = [t0,x0,y0]
    ti = t0
    xi = x0
    yi = y0
    for i in range(1,tamano,1):
        K1x = h * f(ti,xi,yi)
        K1y = h * g(ti,xi,yi)
        
        K2x = h * f(ti+h, xi + K1x, yi+K1y)
        K2y = h * g(ti+h, xi + K1x, yi+K1y)

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

# Programa
# Parámetros de las ecuaciones
a = 0.5
b = 0.7
c = 0.35
d = 0.35
# Ecuaciones
f = lambda t,x,y : a*x -b*x*y
g = lambda t,x,y : -c*y + d*x*y
# Condiciones iniciales
t0 = 0
x0 = 2
y0 = 1
# parámetros del algoritmo
h = 0.5
muestras = 101

# PROCEDIMIENTO
tabla = rungekutta2_fg(f,g,t0,x0,y0,h,muestras)
ti = tabla[:,0]
xi = tabla[:,1]
yi = tabla[:,2]

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

Los resultados numéricos se usan para generar las gráficas presentadas, añadiendo las instrucciones:

# Grafica tiempos vs población
import matplotlib.pyplot as plt
plt.plot(ti,xi, label='xi presa')
plt.plot(ti,yi, label='yi predador')
plt.title('Modelo predador-presa')
plt.xlabel('t tiempo')
plt.ylabel('población')
plt.legend()
plt.grid()
plt.show()
# gráfica xi vs yi
plt.plot(xi,yi)
plt.title('Modelo presa-predador [xi,yi]')
plt.xlabel('x presa')
plt.ylabel('y predador')
plt.grid()
plt.show()

Tarea: Añadir la animación de la gráfica usando la variable tiempo para mostrar un punto que marque el par ordenado[x,y] al variar t.

8.2.2 Runge-Kutta d2y/dx2

Para sistemas de ecuaciones diferenciales ordinarias con condiciones de inicio en x0, y0, y’0

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

Forma estandarizada de la ecuación:

y'' = y' + etc

se convierte a:

y' = z = f_x(x,y,z) z' = (y')' = z + etc = g_x(x,y,z)

con las condiciones de inicio en x0, y0, z0

y se pueden reutilizar los métodos para primeras derivadas, por ejemplo Runge-Kutta de 2do y 4to orden.

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

Runge-Kutta 2do Orden d2y/dx2 en Python:

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

Runge-Kutta 4do Orden d2y/dx2 en Python:

def rungekutta4_fg(fx,gx,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 * fx(xi,yi,zi)
        K1z = h * gx(xi,yi,zi)
        
        K2y = h * fx(xi+h/2, yi + K1y/2, zi + K1z/2)
        K2z = h * gx(xi+h/2, yi + K1y/2, zi + K1z/2)
        
        K3y = h * fx(xi+h/2, yi + K2y/2, zi + K2z/2)
        K3z = h * gx(xi+h/2, yi + K2y/2, zi + K2z/2)

        K4y = h * fx(xi+h, yi + K3y, zi + K3z)
        K4z = h * gx(xi+h, yi + K3y, zi + K3z)

        yi = yi + (K1y+2*K2y+2*K3y+K4y)/6
        zi = zi + (K1z+2*K2z+2*K3z+K4z)/6
        xi = xi + h
        
        estimado[i] = [xi,yi,zi]
    return(estimado)

Una aplicación del algoritmo en Señales y Sistemas:

http://blog.espol.edu.ec/telg1001/01-respuesta-entrada-cero-ltic/

8.2.1 Runge-Kutta 4to Orden dy/dx

Referencia: Chapra 25.3.3 p746 pdf 770, Rodriguez 9.1.8 p358

Para una ecuación diferencial de primer orden con una condición de inicio, la fórmula de Runge-Kutta de 4to orden se obtiene de la expresión con cinco términos:

y_{i+1} = y_i + aK_1 + bK_2 + cK_3 + dK_4

siendo:

y'(x) = f(x_i,y_i) y(x_0) = y_0

debe ser equivalente a la serie de Taylor de 5 términos:

y_{i+1} = y_i + h f(x_i,y_i) + + \frac{h^2}{2!} f'(x_i,y_i) + \frac{h^3}{3!} f''(x_i,y_i) + +\frac{h^4}{4!} f'''(x_i,y_i) + O(h^5) x_{i+1} = x_i + h

que usando aproximaciones de derivadas, se obtienen:

def rungekutta4(d1y,x0,y0,h,muestras):
    # Runge Kutta de 4do orden
    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/2, yi + K1/2)
        K3 = h * d1y(xi+h/2, yi + K2/2)
        K4 = h * d1y(xi+h, yi + K3)

        yi = yi + (1/6)*(K1+2*K2+2*K3 +K4)
        xi = xi + h
        
        estimado[i] = [xi,yi]
    return(estimado)

Note que el método de Runge-Kutta de 4to orden es similar a la regla de Simpson 1/3. La ecuación representa un promedio ponderado para establecer la mejor pendiente.

8.1 EDO con Taylor

Referencia: ejemplo Rodriguez 9.1.1 p335. Chapra 25.1.3 p731 pdf755

En los métodos con Taylor, se aproxima el resultado a n términos de la serie, para lo cual se ajusta la expresión del problema a la derivada correspondiente.

Ejemplo:

Se pide encontrar puntos de la solución en la ecuación diferencial usando los tres primeros términos de la serie de Taylor con h=0.1 y punto inicial x0=0, y0=1

y'-y -x +x^2 -1 = 0

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

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

Luego, al despejar de la fórmula planteada el valor de  y’, se puede obtener a combinar las ecuaciones.

y' = y -x^2 +x +1 y'' = y' -2x + 1 = (y -x^2 +x +1) -2x + 1 y'' = y -x^2 -x +2

ecuaciones que permiten estimar nuevos valores yi+1 para nuevos puntos  muestra distanciados en i*h desde el punto inicial.

Se empieza evaluando el nuevo punto a una distancia x1= x0+h del punto de origen con lo que se obtiene y1 , repitiendo el proceso para el siguiente punto en forma sucesiva

Algoritmo con Python

Para simplificar los calculos se crea una función edo_taylor3t() para encontrar  los valores para una cantidad de muestras distanciadas entre si h veces del punto inicial [x0,y0]

# EDO. Método de Taylor 3 términos 
# estima solucion para muestras separadas 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)

La función se puede usar en un programa que plantee resolver el problema propuesto:

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

# INGRESO.
# d1y = y', d2y = y''
d1y = lambda x,y: y -x**2 + x + 1
d2y = lambda x,y: y -x**2 - x + 2
x0 = 0
y0 = 1
h = 0.1
muestras = 5

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

# SALIDA
print('estimado[xi,yi]')
print(puntos)

obteniendo como resultados:

estimado[xi,yi]
[[ 0.          1.        ]
 [ 0.1         1.215     ]
 [ 0.2         1.461025  ]
 [ 0.3         1.73923262]
 [ 0.4         2.05090205]
 [ 0.5         2.39744677]]

Cálculo de Error

La ecuación diferencial ordinaria del ejercicio tiene una solución conocida, lo que permite encontrar el error real en cada punto respecto a la aproximación estimada.

y = e^x + x + x^2

Note que el error crece al distanciarse del punto inicial

# ERROR vs solución conocida
y_sol = lambda x: ((np.e)**x) + x + x**2

yi_psol = y_sol(xi)
errores = yi_psol - yi
errormax = np.max(np.abs(errores))

# SALIDA
print('Error máximo estimado: ',errormax)
print('entre puntos: ')
print(errores)

con los siguientes resultados:

Error máximo estimado:  0.0012745047595
entre puntos: 
[ 0.  0.000170  0.000377  0.000626  0.000922  0.00127 ]

Gráfica

Para visualizar los resultados se usan los puntos estimados y muchos más puntos para la solución conocida en el intervalo de estudio. Con la gráfica se podrán observar las diferencias entre las soluciones encontradas.

# GRAFICA [a,b+2*h]
a = x0
b = h*muestras+2*h
muestreo = 10*muestras+2
xis = np.linspace(a,b,muestreo)
yis = y_sol(xis)

# Gráfica
import matplotlib.pyplot as plt
plt.plot(xis,yis, label='y conocida')
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.title('EDO: Solución con Taylor 3 términos')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()

Tarea: Realizar el ejercicio con más puntos muestra, donde se visualice que el error aumenta al aumentar la distancia del punto inicial [x0,y0]

8.2 Runge-Kutta 2do Orden dy/dx

Referencia: Burden 5.4 p272 pdf282, Chapra 25.3 p740 pdf164, Rodriguez 9.1.7 p354

Video Revisar: 4’39» + 2 minutos

¿ Que trayectoria siguen los proyectiles?
¿ Que trayectoria siguen los aviones?
¿ Cuál es la relación entre las trayectorias?


Los métodos de Runge-Kutta  mejoran la aproximación a la respuesta sin requerir determinar las expresiones de las derivadas de orden superior, usan una corrección a la derivada tomando valores de puntos alrededor.

Por ejemplo, Runge-Kutta de 2do Orden usa el promedio entre los incrementos xi y xi+h, calculados como términos K1 y K2.

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

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)

Ejercicio: Para probar el algoritmo se usa la misma ecuación del problema presentado en  ‘EDO con Taylor‘ :

y'-y -x +x^2 -1 = 0

que aplicado con Runge Kutta, se obtiene:

estimado[xi,yi]
[[ 0.          1.        ]
 [ 0.1         1.2145    ]
 [ 0.2         1.4599725 ]
 [ 0.3         1.73756961]
 [ 0.4         2.04856442]
 [ 0.5         2.39436369]]
Error máximo estimado:  0.00435758459732
entre puntos: 
[ 0.          0.00067092  0.00143026  0.0022892   0.00326028  0.00435758]
>>> 

Compare los resultados con Taylor de 2 y 3 términos.

Los resultados se muestran también en la gráfica:

Se adjunta el programa de prueba que usa la función rungekutta2(d1y,x0,y0,h,muestras)  :

# 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'
d1y = lambda x,y: y -x**2 + x + 1
x0 = 0
y0 = 1
h = 0.1
muestras = 5

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

# SALIDA
print('estimado[xi,yi]')
print(puntosRK2)

# ERROR vs solución conocida
y_sol = lambda x: ((np.e)**x) + x + x**2

yi_psol = y_sol(xi)
errores = yi_psol - yiRK2
errormax = np.max(np.abs(errores))

# SALIDA
print('Error máximo estimado: ',errormax)
print('entre puntos: ')
print(errores)

# GRAFICA [a,b+2*h]
a = x0
b = h*muestras+2*h
muestreo = 10*muestras+2
xis = np.linspace(a,b,muestreo)
yis = y_sol(xis)

# Gráfica
import matplotlib.pyplot as plt
plt.plot(xis,yis, label='y conocida')
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()