Sistemas de 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 preparan las ecuaciones de la forma:

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

resultados que se usan con el algoritmo rungekutta2fg(fx,gx,x0,y0,z0,h,muestras) y 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]]
>>> 

que generan las gráficas para análisis de resultados:

las intrucciones en python que resuelven el problema

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

cuyos resultados numéricos se convierten en gráficas al añadir:

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

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]
    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 + d_k4

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

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' = f, d2y = y'' = f'
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.00017092  0.00037776  0.00062618  0.00092265  0.0012745 ]

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