Revista Ingeniería y Región

ISSN 1657 - 6985 | e-ISSN 2216 - 1325



Revista Ingeniería y Región, Volumen 26. 29-44 pp
Julio - Diciembre de 2021/Universidad Surcolombiana



Artículo de revisión



Formación de patrones en un modelo difusivo bidimensional depredador - presa tipo Holling II


Pattern formation in a two-dimensional diffusive predator model - Holling II type prey


Allison María Ramírez Fierro

https://orcid.org/0000-0002-0935-1654

Matemático, Universidad Surcolombiana

E-mail: maria21907@hotmail.com


Ingrid Tatiana Cumbe-Morales

https://orcid.org/0000-0002-4274-3705

Matemático, Universidad Surcolombiana

E-mail: ingri.b.5656@hotmail.com


Christian Camilo Cortes García

https://orcid.org/0000-0002-8955-4530

PhD(c) en Ingeniería Matemática, Universidad Carlos III de Madrid. Docente Universidad Surcolombiana. Investigador predoctoral Centro Nacional de biotecnología, Madrid – España.

E-mail: christian.cortes@usco.edu.co, chcortes@math.uc3m.es, cc.cortes@cnb.csic.es




Fecha de envío: 08/03/2021

Fecha de Revisión: 16/04/2021

Fecha de Aprobación: 29/09/2021


DOI: https://doi.org/10.25054/22161325.2972






Resumen


En este trabajo se presenta un método numérico para observar la distribución y el comportamiento entre la interacción de presas y depredadores bajo un modelo difusivo bidimensional, con crecimiento logístico para las presas y funcional de depredación tipo Holling II. Al realizar algunas perturbaciones en los parámetros del modelo, determinar condiciones de contornos apropiadas y establecer intervalos de tiempo para la convergencia del método, se presentan diversos patrones sobre las soluciones simuladas en el modelo. En vista que el modelo matemático sin difusión presenta ciclo límite, un equilibrio que puede ser localmente un nodo o una espiral estable, las soluciones numéricas del modelo difusivo reflejan dichos comportamientos.


Palabras clave: Método de diferencias finitas; operador laplaciano; bifurcación de Hopf; coeficiente de difusión.


Abstract


This paper presents a numerical method to observe the distribution and behavior between prey and predator interaction under a two-dimensional diffusive model, with logistic growth for prey and Holling II type predation functional. By performing some perturbations in the model parameters, determining appropriate boundary conditions and establishing time intervals for the convergence of the method, several patterns are presented on the simulated solutions in the model. Since the mathematical model without diffusion presents a limit cycle, an equilibrium that can be locally a node or a stable spiral, the numerical solutions of the diffusive model reflect these behaviors.


Keywords: Finite difference method; Laplacian operator; Hopf bifurcation; diffusion coefficient.


1. Introducción


La formulación matemática que describe un fenómeno biológico determinista puede ser descrito por un sistema dinámico, el cual permite prever cualquier estado futuro si el operador de evolución y su estado inicial son conocidos. Si el tiempo para ese operador se define en la recta real, se habla de un sistema dinámico continuo.


En particular, las ecuaciones diferenciales de la forma (Ecuación 1):


{ w ˙ = f ( 𝑤 , α ) α ˙ = 0 (1)

con w ∈ ℝ𝑛 la variable de estados y 𝛼 ∈ ℝ𝑚 un vector de parámetros es un sistema dinámico y del cual pueden ser usados para describir y analizar todas las posibles dinámica de un fenómeno.


Los retratos de fases para el modelo (1) del cual cambia drásticamente cuando se altera el parámetro 𝛼, se denomina bifurcación. Los trabajos publicados por Hirsch y Smale (1974), Sotomayor (1979), Arrowsmith (1992), Berkooz, et al., (1996), Kuznetsov (1998), Blanchard y Devaney (1999), Chicone (2006) y, Renardy y Rogers (2006) y Cortés (2017) describen la forma matemática para analizar cualitativamente el modelo (1) respecto a sus equilibrios, además de todos los posibles casos de bifurcaciones, tanto local como global.


Los modelos biológicos son planteados para analizar la dinámica en la interacción de dos especies bajo un sistema de ecuaciones diferenciales ordinarias, conocidos generalmente como modelo depredador – presa. Un ejemplo para este tipo de modelos es expresado por (Ecuación 2):


{ u ˙ = r u ( 1 u K ) b u v 1 + e u, v ˙ = m u v 1 + e u d v (2)

con 0 ≤ 𝑢 ≤ 𝐾 la población de presas, la cual tiene un crecimiento logístico, y 0 ≤ v la población de depredadores con funcional de depredación Holling II. En este caso, 𝑟 > 0 representa la tasa de nacimientos de las presas, 𝑑 > 0 la tasa de muertes de los depredadores, 𝑏, 𝑚 > 0 las tasas de encuentro entre las dos especies, respectivamente y, 𝑒 > 0 la constante de saturación de las presas al ser consumidas por los depredadores.


En los trabajos publicados por Freedman (1980), Edelstein (1988), Kuang y Freedman (1988), Wolkowicz (1988), Murray (2002), Zhou (2019), Zhu, et al., (2003), Du y Shi (2006), Liu (2010), Hasik (2010), Zu y Mimura (2010), Tian y Xu (2011), Tang, et al., (2014), Lin, et al., (2018), Mukherjee y Maji (2020) y, Cortés (2021a) se plantean y se analizan cualitativamente algunos modelos biológicos de la forma (Ecuación 1) que describen el comportamiento de las presas y depredadores bajo ciertas suposiciones.


Sin embargo, al incluir algunos elementos matemáticos, como la difusión unidimensional de la especie, el modelo (1) se replantea como un sistema de ecuaciones diferenciales parciales de la forma (Ecuación 3):


{ w t = D w x x + f ( w , α ) D ˙ = α ˙ = 0 (3)

con 𝐷 el vector de difusión, con entradas no negativas, del cual permite analizar de diversas maneras el comportamiento y la dinámica de dicha población de una manera más acorde con el ambiente. En particular, al agregar un término difusivo unidimensional a las presas y/o depredadores al modelo (2) se tiene que (Ecuación 4):


{ u t = D 1 u x x + r u ( 1 u K ) b u v 1 + e u, v t = D 2 u x x + m u v 1 + e u d v (4)

con 𝐷1> 0 y 𝐷2> 0 las constantes de difusión para las presas y depredadores, respectivamente.


El estudio analítico para modelos biológicos de la forma (Ecuación 3) están centrados en mostrar soluciones tipo ondas viajeras de la forma u ( x , t ) = U ( x ± c t )   y   v ( x , t ) = V ( x ± c t ) ,   con   c > 0 la velocidad de onda. En los trabajos publicados por Grindrod (1991), Huang (2003), Li y Wu (2008), Hong y Weng (2013), Wu, et al., (2014), Ghazaryan, et al., (2015), Ai, et al., (2017) y, Cortés y Ramírez (2021b), y muestran mecanismos matemáticos para demostrar la existencia de ondas viajeras como soluciones al modelo de la forma (3).


Sin embargo, el análisis matemático hace cada vez más complejo si se al incluir un término difusivo bidimensional al modelo (1), esto es, (Ecuación 5):


{ w t = D Δ w + f ( w , α ), D ˙ = α ˙ = 0 (5)

con Δ w = w x x + w y y un operador Laplaciano de w . En particular, si el modelo (2) es modificado por (Ecuación 6):


{ u t = D 1 Δ u + r u ( 1 u K ) b u v 1 + e u, v t = D 2 Δ v + m u v 1 + e u d v (6)

su estudio se centra en determinar todos los posibles patrones que presentan dichas soluciones de u y 𝑣 en el plano x y para un tiempo definido. Los trabajos publicados por Cohen y Grossberg (1983), Murray y Oster (1984), Doelman (1995), Golubitsky, et al., (1999), Cruywagen, et al., (2000), Hoyle y Hoyle (2006), Mancini (2007), Cross y Greenside (2009), Cárdenas, et al., (2013), Rao y Kang (2016) y, Wu, et al., (2018) y, muestran diversas maneras para encontrar patrones a modelos biológicos de la forma (Ecuación 5).

Para encontrar una solución explicita de un modelo biológico, y en general para cualquier modelo matemático compuesto por ecuaciones diferenciales parciales, no siempre es fácil. De esa forma, se buscan soluciones aproximadas mediante la aplicación de métodos computacionales apropiados para resolver problemas matemáticos difíciles de resolver o que no se pueden resolver analíticamente. En este caso, existen diversos métodos numéricos para aproximar soluciones en modelos de ecuaciones diferenciales y que podrían ser consultados en los trabajos de Mathews, et al., (2000), Iserles (2009) y Sauser (2013).


Por tanto, en este trabajo se determinan los posibles patrones, al ser uso de un método numérico, para un modelo depredador - presa tipo Holling II (6) y término difusivo bidimensional. Para ello, al tomar los trabajos publicados por Garvie (2005ab, 2007), en la segunda sesión se presenta un resumen del método numérico por diferencias finitas, ejecutado por Garvie, para resolver numéricamente un modelo de ecuaciones diferenciales parciales con término difusivo bidimensional.


Por otro lado, y a partir del trabajo publicado por Cortés (2021b), en la tercera sesión se describe un modelo depredador-presa tipo Holling II con termino difusivo bidimensional y, en la cuarta sesión, se presenta un análisis cualitativo al modelo cuando los parámetros de difusión son nulos.


Por último, al usar los lineamientos que se presentaron en las secciones anteriores y al variar algunos parámetros y condiciones iniciales, en la quinta sesión se determinan una serie de patrones que muestran el comportamiento del modelo difusivo bidimensional depredador-presa con respuesta funcional tipo Holling II y se realiza una comparación con la dinámica del modelo sin difusión.


2. Materiales y métodos


Para generar soluciones numéricas al sistema (Ecuación 7):


{ u t = D 1 Δ u + f ( u , v ) v t = D 2 Δ v + g ( u , v ) (7)

con D 1 , D 2 0   y   Δ z = z x x + z y y conocido como un operador Laplaciano de z, se emplea un método numérico de diferencia finita, el cual sustituye las derivadas parciales de una ecuación diferencial por aproximaciones discretas y graficar las soluciones numéricas del modelo (7) sobre una malla con una programación adecuada de sus nodos


Según Garvie (2007), para discretizar el modelo (2) en el intervalo de tiempo [ 0 , T ] y niveles de tiempo t n = n Δ t, para todo n = 1 , , N , con Δ t = T N , considere una malla cuadrática de puntos dada por Ω = [ A , B ] × [ A , B ], con 0 A < B, y una subdivisión uniforme de puntos en la cuadrícula ( x i , y i ) = ( i h + A , j h + A ) , i , j = 0 , , J, con incremento de espacio h = ( B A ) J .


Como lo observado en la Figura 1 y al considerar J ~ = ( J + 1 ) 2 1, los nodos rellenos representan los valores conocidos de la solución aproximada, a partir de las condiciones iniciales y de fronteras, y los nodos ficticios son los puntos de la malla que se deben programar mediante el método, esto es, 2 ( J ~ + 1 ) incógnitas para resolver, donde U i , j n = ( U i , j n , V i , j n ) T representa la solución aproximada de u = ( u , v ) T del modelo (7) en el punto ( x i , y i , t n ).




Los nodos en la malla están numerados naturalmente, de izquierda a derecha, y comienza con la fila inferior, desde k = 0 , 1 , , J. ~ La relación entre la numeración de los nodos y la indexación ( i , j ) está dada por U k n = U i , j, n donde k = i + j ( J + 1 ) para i , j = 0 , , J.


Una forma discreta del modelo (7), con U i , j n = ( U i , j n , V i , j n ) para n = 1 , , N e i , j = 0 , , J , es dado por (Ecuación 8):


{ n U i , j n = Δ h U i , j n + f ^ ( U i , j n , U i , j n 1 ) n V i , j n = Δ h V i , j n + g ^ ( U i , j n , U i , j n 1 ) (8)

donde f ^ y g ^ son funciones discretas, en términos del tiempo anterior t n 1, modificadas de f y g , esto es (Ecuación 9),


f ^ ( U i , j n , U i , j n 1 ) = f ( U i , j n 1 ) g ^ ( U i , j n , U i , j n 1 ) = g ( U i , j n 1 ), (9)

y las aproximaciones de n Z n y del operador Laplaciano discreto en dos dimensiones Δ h z i , j , son representadas por (Ecuación 10):


n Z n = Z n Z n 1 Δ t Δ h Z i , j = Z i , j + 1 + Z i , j 1 + Z i + 1 , j + Z i 1 , j 4 Z i , j h 2 (10)

A partir de las condiciones iniciales del modelo (Ecuación 7), dadas por U i , j 0 := u ( x i , y j , 0 ) , V i , j 0 := v ( x i , y j , 0 ), y las condiciones de contorno a lo largo de los cuatro bordes del cuadrado, es decir, para 1 i , j J 1 , se tiene que (Ecuación 11),


U i , 1 n := U i , 1 n , U i , j + 1 n : = U i , j 1, n U J + 1 , j n := U J 1 , j n , U 1 , j n := U 1 , j, n (11)

y de forma análogamente para las aproximaciones de v .


Como lo observado en la Figura 1, las condiciones de frontera de las esquinas para u son dadas por (Ecuación 12):


U 0 , 1 n = ( U 0 , 0 n + U 0 , 1 n ) 2 , U 1 , 0 n = ( U 0 , 0 n + U 1 , 0 n ) 2 U J , J + 1 n = ( U J , J n + U J , J 1 n ) 2 , U J + 1 , J n = ( U J , J n + U J 1 , J n ) 2 U J , 1 n = 2 U J , 1 n U j , 0 n , U J + 1 , 0 n = 2 U J 1 , 0 n U J , 0 n U 0 , J + 1 n = 2 U 0 , J 1 n U J , 0 n , U 1 , J n = 2 U 1 , J n U 0 , J n (12)

donde se obtiene de igual forma para las aproximaciones de v


Para calcular las 2 ( J ~ + 1 ) ecuaciones lineales que representan los nodos interiores, se usa la siguiente matriz de bloque, con la numeración natural de los nodos (Ecuación 13),


( B 1 0 0 B 2 ) ( U n V n ) = ( U n 1 + Δ t F V n 1 + Δ t G ) , 1 n N, (13)

con U n = ( U 0 n , , U J ~ n ) T , V n = ( V 0 n , , V j n ) T , { F } k = f ( U k n 1 ) y { G } k = g ( U k n 1 ) para k = 0 , , J ~ .


La matriz de coeficientes constantes 𝐵1 y 𝐵2 son calculadas por (Ecuación 14):


B 1 := I + D 1 Δ t L, B 2 := I + D 2 Δ t L (14)

donde la matriz 𝐿, de dimensión ( J + 1 ) 2 × ( J + 1 ) 2 y bloques ( J + 1 ) × ( J + 1 ) , es dada por (Ecuación 15):


L = 1 h 2 [ S T W X W W X W W X W W X W Y Z ] (15)

con T = diag { 3 2 , 2 , 2 , , 2 , 2 , 3 } , W = I , Y = diag { 3 , 2 , 2 , , 2 , 2 , 3 2 },


X = [ 4 2 1 4 1 1 4 1 1 4 1 1 4 1 2 4 ] Z = [ 6 3 1 4 1 1 4 1 1 4 1 1 4 1 3 2 3 ]   y,   S
= [ 3 3 2 1 4 1 1 4 1 1 4 1 1 4 1 3 6 ]

3. Resultados y discusiones


3.1 Presentación del modelo


Sean u ( x , y , t ) 0   y   v ( x , y , t ) 0 la población de presas y depredadores, respectivamente, cuya distribución bidimensional x , y siguen un modelo de la siguiente forma (Ecuación 16):


{ u t = D 1 Δ u + r u ( 1 u K ) b u v 1 + e u, v t = D 2 Δ v + m u v 1 + e u d v (16)

donde


𝐷1, 𝐷2 ≥ 0 son los coeficientes de difusión de las presas y depredadores, respectivamente,

p ( u ) := r u ( 1 u K ) describe el crecimiento logístico de las presas en ausencia de los depredadores, con r > 0 la tasa de crecimiento y K > 0 la capacidad de carga,

q ( u , v ) := b u v 1 + e u es la disminución en las presas al interactuar con los depredadores, donde b > 0 es la tasa de encuentro entre las dos especies y e > 0 la constante de saturación de las presas,

𝑧 ( u , v ) := m u v 1 + e u describe el aumento en la población de los depredadores al interactuar con las presas, donde m > 0 es la tasa de encuentro entre las dos especies,

s ( v ) := d v es el decrecimiento de los depredadores en ausencia de las presas con d > 0 la tasa de muerte de los depredadores.


Nótese que el modelo (16) es equivalente al modelo (7) con f ( u , v ) := r u ( 1 u K ) b u v 1 + e u y g ( u , v ) := m u v 1 + e u d v


3.2 Análisis cualitativo del modelo sin difusión para ambas especies


En el caso que los coeficientes de difusión en el modelo (16) sean nulos, es decir 𝐷1 = 𝐷2 = 0, el modelo (16) toma la forma de un sistema de ecuaciones diferenciales ordinarias, esto es (Ecuación 17),

{ u ˙ = r u ( 1 u K ) b u v 1 + e u, v ˙ = m u v 1 + e u d v, (17)

el cual describe la cantidad de presas y depredadores a lo largo del tiempo.


El modelo (17) tiene dos puntos de equilibrio sobre los ejes x y dados por: E 1 = ( 0 , 0 )   y   E 2 = ( K , 0 ) . Además, si m > d e   y   m K > d ( e K + 1 ), el modelo (11) presenta un punto de equilibrio interior dado por


E = ( d m d e , r m [ m K d ( e K + 1 ) ] b K ( m d e ) 2 ).

Para determinar la estabilidad local de cada equilibrio sobre el modelo (17) se debe analizar el signo de la traza y el determinante de la matriz jacobiana 𝐴 evaluada en cada equilibrio. En efecto, la matriz Jacobiana del modelo (17) evaluada en el punto de equilibrio E 1 es dada por (Ecuación 18).


A ( E 1 ) = ( r 0 0 d ) (18)

con Det A ( E 1 ) = r d < 0. Por consiguiente, el punto E 1 es localmente un punto silla.


De igual forma, la matriz Jacobiana del modelo (17) asociada al punto E 2 esta dada por (Ecuación 19):


A ( E 1 ) = ( r b K 1 + e K 0 m K 1 + e K d ) (19)

con valores propios: λ 1 = r < 0   у   λ 2 = m K d ( 1 + e K ) 1 + e K . Si K < d m d e entonces λ 2 < 0 y así el punto E 2 es localmente un nodo estable, y para K > d m d e se tiene λ 2 > 0, esto es, el punto E 2 es localmente un punto silla.


Por último, de la matriz Jacobiana A = A ( E ) del modelo (17) evaluada en E se tiene que (Ecuación 20):


Tr A = r d [ ( e K + 1 ) ( m d e ) 2 m ] m K ( m d e ) , Det A = r d [ m K d ( e K + 1 ) ] m K , (20)

Δ = ( Tr A ) 2 4 Det A = { r d [ ( e K + 1 ) ( m d e ) 2 m ] m K ( m d e ) } 2 4 r d [ m K d ( e K + 1 ) ] m K


De (15), el Det A > 0 y por lo tanto E no es localmente un punto silla. Asimismo, de (20) se tiene que Tr A > 0 si K > d e + m e ( m e d )   y   Tr A < 0   si   K < d e + m e ( m e d ) . Por consiguiente, la dinámica del modelo (17) puede ser resumida en el siguiente resultado, cuya prueba puede ser vista en Cortés, et al., (2021).


Proposición 1. El modelo (12) tiene dos puntos de equilibrio E 1 = ( 0 , 0 )   y   E 2 = ( K , 0 ) y un punto de equilibrio interior, si m > de   y   m K > d ( e K + 1 ) , dado por

E = ( d m d e , r m [ m K d ( e K + 1 ) ] b K ( m d e ) 2 )

Las condiciones de estabilidad local para los puntos de equilibrios están dadas de la siguiente manera:


El punto E 1 = ( 0 , 0 ) es locamente un punto silla.


El punto E 2 = ( K , 0 )


   una silla si K > d m d e ,


   un nodo estable si K < d m d e .


El punto E es localmente


   una espiral inestable si Δ < 0   y   K > d e + m e ( m e d ) ,


   una espiral estable si Δ < 0   y   K < d e + m e ( m e d ) ,


   un nodo estable si Δ > 0   y   K < d e + m e ( m e d ) ,


   un nodo inestable si Δ > 0   y   K > d e + m e ( m e d ) .


Además, si el punto de equilibrio interior E es localmente inestable, entonces el modelo (6) presenta un único ciclo límite globalmente estable en el primer cuadrante.


Por otro lado, el modelo (17) presenta una bifurcación de Hopf al colisionar el equilibrio E con el ciclo limite estable, esto es, si Tr A = 0, y cuya prueba puede ser vista en Cortés, et al., (2021).


Lema 1. Si K = d e + m e ( m e d ) entonces el modelo (17) presenta una bifurcación de Hopf.


De igual forma, se deduce que la colisión del equilibrio E con el punto de silla E 2 genera un cambio en la estabilidad de E 2. En este caso, se forma una bifurcación transcrítica si la componente sobre el eje 𝑥 del equlibrio E concuerda con la componente sobre el eje 𝑥 de E 2 cómo se observa en el siguiente resultado.


Lema 2. Si K = d m e d entonces el modelo (17) presenta una bifurcación transcrítica.


En la Figura 2 se muestran las curvas de bifurcación del modelo (17) con sus retratos de fase. En la región 1 se observa que el punto E 2 es un nodo estable que se transforma en un punto silla a medida que K cruza la curva K = d m d e y posteriormente se forma un punto de equilibrio interior estable dado en la región 2.


En la región 3, al traspasar los parámetros la curva Δ = 0, el nodo estable se transforma en un foco estable como lo observado en la región 2, y a medida que traspasa la curva K = d e + m e ( m d e ) , el modelo (17) presenta un ciclo límite estable como se muestra en la región 4.



×

(a) Curvas de bifurcación (𝐾, 𝑑)



×

Figura 2. Curvas de bifurcación y retratos de fase del modelo (17) con parámetros fijos 𝑟 = 1, 𝑏 = 1 , 𝑚 = 0.5 y 𝑒 = 0.125, y parámetros de bifurcación: (b) 𝐾 = 1.12 y 𝑑 = 1.2; (c) 𝐾 = 3 y 𝑑 = 1; (d) 𝐾 = 0.8 y 𝑑 = 0.1; (e) 𝐾 = 12 y 𝑑 = 0.39. BT=Bifurcación Transcrítica y BH=Bifurcación de Hopf. Fuente: Cortés (2021b).


3.3 Simulación de patrones


Al aplicar el método numérico, dado en la sección uno, al modelo (16), se procede a realizar simulaciones para detectar posibles patrones en la forma como se distribuyen las especies sobre una malla. El código de programación fue ejecutado en Matlab 2016.


Si 𝐷1 = 1, 𝐷2 = 10, 𝑟 = 1, 𝐾 = 1, 𝑏 = 2.5, 𝑒 = 2.5, 𝑚 = 5 y 𝑑 = 0.6, se observa en la Figura 3 que tanto los depredadores como las presas se distribuyen en forma de espiral a lo largo de la malla [0,500] × [0,500], y se van desvaneciendo a medida que aumenta el valor de 𝑇, formando diversas manchas. Este fenómeno concuerda con la dinámica del modelo (17) presentado en la Figura 2.



×

×

×

Figura 3. Distribución aproximada de presas y depredadores para el modelo (16) con incremento Δ t = 0.1, condición inicial u ( x , y , 0 ) = 6 35 2 × 10 7 ( x 0.1 y 225 ) ( x 0.1 y 675 ) , v ( x , y , 0 ) = 116 245 3 × 10 5 ( x 450 ) 1.2 × 10 4 ( y 150 ) y un tiempo T : ( a , b ) T = 150 ; ( c , d ) T = 300 ; ( e , f) T = 600 . Distribución de las presas u a la izquierda. Distribución de los depredadores v a la derecha. Fuente: Elaboración propia.


Al usar los mismos parámetros que generan la Figura 3 con otras condiciones iniciales u ( x , y , 0 )   у   v ( x , y , 0 ), se observa en la Figura 4 unos cambios drásticos en la distribución de las presas y depredadores para 𝑇 pequeño, el cual, a diferencia de la Figura 3, los depredadores como las presas se distribuyen en forma de nodos a lo largo de la malla [ 0 , 500 ] × [ 0 , 500 ].



×

×

×

Figura 4. Distribución aproximada de presas y depredadores para el modelo (16) con incremento Δ t = 1 24 , condición inicial u ( x , y , 0 ) = 6 35 2 × 10 7 ( x 180 ) ( x 720 ) 6 × 10 7 ( y 90 ) ( y 210 ) , v ( x , y , 0 ) = 116 245 3 × 10 5 ( x 450 ) 6 × 10 5 ( y 135 ), y tiempo T : ( a , b ) T = 150 ; ( c , d ) T = 200 ; ( e , f ) T = 500 . Distribución de las presas u a la izquierda. Distribución de los depredadores v a la derecha. Fuente. Elaboración propia.


Para los parámetros y condiciones dadas en la Figura 5 solo presenta un único patrón. En dicha Figura se puede observar que la distribución de ambas especies presenta ciclos, lo cual concuerda con la dinámica del modelo (11) mostrado en la Figura 2.



×

Figura 5. Densidades de presas y depredadores para el modelo (16) con parámetros 𝐷1 = 1, 𝐷2 = 10, 𝑟 = 1, 𝐾 = 1, 𝑏 = 1, 𝑒 = 1.5, 𝑚 = 1 y 𝑑 = 5 con tiempo 𝑇 = 110, incremento de espacio ℎ = 0.5 , y condición inicial u ( x , y , 0 ) = 1 , v ( x , y , 0 ) = { 0.2 si ( ( x 200 ) 2 + ( y 200 ) 2 ) < 400 0 si ( x 200 ) 2 + ( y 200 ) 2 >= 400 Distribución de presas a la izquierda. Distribución de depredadores a la derecha. Fuente: Elaboración propia.


4. Conclusiones


El método numérico de diferencias finitas resultó ser efectivo para el desarrollo de la simulación, de tal forma que se logró simular las distintas formas de distribución y patrones entre presas y depredadores, cuyo comportamiento es descrito por el modelo difusivo (16). Sin embargo, el método no permitió detectar los posibles comportamientos posteriores al ciclo formado. Para analizar dichos comportamientos cuando se aumenta el valor de 𝑇, se podría aplicar otros métodos numéricos, como el método de elementos finitos.


La formación de patrones dados en el modelo (16) concuerda con las presentadas en la dinámica del modelo sin difusión (17), por lo que se podría pensar que dicho modelo difusivo presenta una bifurcación de Hopf.


5. Referencias bibliográficas


Arrowsmith, D., & Place, C. M. (1992). Dynamical systems: differential equations, maps, and chaotic behaviour (Vol. 5). CRC Press. https://doi.org/10.1017/cbo9780511622700


Berkooz, G., Holmes, P., & Lumley, J. L. (1996). Coherent structures, dynamical systems and symmetry. Cambridge Monographs on Mechanics, Cambridge University Press.


Blanchard, P., & Devaney, R. L. (1999). Ecuaciones diferenciales (No. 517.38 B5).


Cárdenas, M. M., Cuenca, J. V., & Betancourth, G. L. (2013). Bifurcaciones en modelos sobre identificación de patrones. Revista Entornos, 26 (2), 139-159. https://doi.org/10.25054/01247905.480


Chicone, C. (2006). Ordinary differential equations with applications (Vol. 34). Springer Science & Business Media. Cohen, M. A., & Grossberg, S. (1983). Absolute stability of global pattern formation and parallel memory storage by competitive neural networks. IEEE transactions on systems, man, and cybernetics, (5), 815-826. https://doi.org/10.1109/tsmc.1983.6313075


Cortés, C. (2017). Identificación de una Bifurcación de Hopf con o sin Parámetros. Revista de Ciencias, 21 (2), 59-82. https://doi.org/10.25100/rc.v21i2.6699


Cortés, C. (2021a). Bifurcaciones en Modelo Gause Depredador-Presa con discontinuidad. Revista de Matemática: Teoría y Aplicaciones, 28(2), 183-208.


Cortés, C., Ramírez-Fierro, A. (2021b). Solución tipo onda viajera en un Modelo Difusivo Depredador-Presa tipo Holling II. Revista de Matemática: Teoría y Aplicaciones, 28(2), 209-236.


Cross, M., & Greenside, H. (2009). Pattern formation and dynamics in nonequilibrium systems. Cambridge University Press.


Cruywagen, G. C., Murray, J. D., & Maini, P. K. (2000). An envelope method for analyzing sequential pattern formation. SIAM Journal on Applied Mathematics, 61(1), 213-231. https://doi.org/10.1137/s003613999732857x


Doelman, A., & Van Harten, A. (1995). Nonlinear dynamics and pattern formation in the natural environment (Vol. 335). CRC Press.


Du, Y., & Shi, J. (2006). A diffusive predator–prey model with a protection zone. Journal of Differential Equations, 229(1), 63-91. https://doi.org/10.1016/j.jde.2006.01.013


Edelstein, L., (1988). Mathematical Models in Biology. New York. SIAM Edition.


Freedman, H. I. (1980). Deterministic mathematical models in population ecology (Vol. 57). Marcel Dekker Incorporated.


Garvie,M., Trenchea, C., (2005a). Analysis of two generic spatially extended predator–prey models. Nonlinear Anal.


Garvie, M., Trenchea, C.,(2005b). Finite element approximation of spatially extended predator–prey interactions with the Holling type II functional response. https://doi.org/10.1007/s00211-007-0106-x


Garvie, M.,(2007). Finite-Difference Schemes for Reaction–Diffusion Equations Modeling Predator–Prey Interactions in MATLAB. Bulletin of Mathematical Biology (69): 931–956. https://doi.org/10.1007/s11538006-9062-3


Ghazaryan, A., Manukian, V., & Schecter, S. (2015). Travelling waves in the Holling–Tanner model with weak diffusion. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471(2177), 20150045. https://doi.org/10.1098/rspa.2015.0045


Golubitsky, M., Luss, D., & Strogatz, S. H. (1999). Pattern formation in continuous and coupled systems (pp. 65-82). Springer.


Grindrod, P. (1991). Patterns and waves: The theory and applications of reaction-diffusion equations. Oxford University Press.


Hasik, K. (2010). On a predator–prey system of Gause type. Journal of Mathematical Biology, 60(1), 59-74. https://doi.org/10.1007/s00285-009-0257-8


Hirsch, M. W., & Smale, S. (1974). Difierential Equations. Dynamical Systems, and Linear Algebra. San Diego: Academic Press.


Hong, K., & Weng, P. (2013). Stability and traveling waves of a stage-structured predator–prey model with Holling type-II functional response and harvesting. Nonlinear Analysis: Real World Applications, 14(1), 83-103. https://doi.org/10.1016/j.nonrwa.2012.05.004


Huang, J., Lu, G., & Ruan, S. (2003). Existence of traveling wave solutions in a diffusive predator-prey model. Journal of Mathematical Biology, 46(2), 132-152. https://doi.org/10.1007/s00285-002-0171-9


Ai, S., Du, Y., & Peng, R. (2017). Traveling waves for a generalized Holling–Tanner predator–prey model. Journal of Differential Equations, 263(11), 7782-7814. https://doi.org/10.1016/j.jde.2017.08.021


Hoyle, R., & Hoyle, R. B. (2006). Pattern formation: an introduction to methods. Cambridge University Press.


Iserles, A. (2009). A first course in the numerical analysis of differential equations (No. 44). Cambridge university press.


Kuang, Y., & Freedman, H. I. (1988). Uniqueness of limit cycles in Gause-type models of predator-prey systems. Mathematical Biosciences, 88(1), 67-84. https://doi.org/10.1016/0025-5564(88)90049-1


Kuznetsov, Y. A., (1998). Elements of Applied Bifurcation Theory. Applied Mathematical Sciences, (112) New York. Second Edition.


Li, W. T., & Wu, S. L. (2008). Traveling waves in a diffusive predator–prey model with Holling type-III functional response. Chaos, Solitons & Fractals, 37(2), 476-486. https://doi.org/10.1016/j.chaos.2006.09.039


Lin, Q., Xie, X., Chen, F., & Lin, Q. (2018). Dynamical analysis of a logistic model with impulsive Holling type-II harvesting. Advances in Difference Equations, 2018(1), 1-22. https://doi.org/10.1186/s13662-018-1563-5


Liu, P. P. (2010). An analysis of a predator–prey model with both diffusion and migration. Mathematical and Computer Modelling, 51(9-10), 1064-1070. https://doi.org/10.1016/j.mcm.2009.12.010


Mancini, H. L. (Ed.). (2007). New trends, dynamics and scales in pattern formation. EDP Sciences.


Mathews, J. H., Kurtis, D. F., (2000). Métodos Numéricos con MATLAB. Editorial Pearson Prentice Hall. Tercera edición.


Mukherjee, D., & Maji, C. (2020). Bifurcation analysis of a Holling type II predator-prey model with refuge. Chinese Journal of Physics, 65, 153-162. https://doi.org/10.1016/j.cjph.2020.02.012


Murray, J. D. (2002). Mathematical Biology I. An Introduction. JD Murray.


Murray, J. D., & Oster, G. F. (1984). Generation of biological pattern and form. Mathematical Medicine and Biology: A Journal of the IMA, 1(1), 51-75.


Rao, F., & Kang, Y. (2016). The complex dynamics of a diffusive prey–predator model with an Allee effect in prey. Ecological complexity, 28, 123-144. https://doi.org/10.1016/j.ecocom.2016.07.001


Renardy, M., & Rogers, R. C. (2006). An introduction to partial differential equations (Vol. 13). Springer Science & Business Media.


Sauser, T., (2013). Análisis Numérico. Segunda edición. México. p. 664. Editorial Pearson Educación.


Sotomayor, J. (1979). Liçoes de equaçoes diferenciais ordinárias (Vol. 11). Instituto de Matemática Pura e Aplicada, CNPq.


Tang, G., Tang, S., & Cheke, R. A. (2014). Global analysis of a Holling type II predator–prey model with a constant prey refuge. Nonlinear Dynamics, 76(1), 635-647. https://doi.org/10.1007/s11071-013-1157-4


Tian, X., & Xu, R. (2011). Global dynamics of a predator-prey system with Holling type II functional response. Nonlinear Analysis: Modelling and Control, 16(2), 242-253. https://doi.org/10.15388/na.16.2.14109


Wolkowicz, G. S. K. (1988). Bifurcation analysis of a predator-prey system involving group defence. SIAM Journal on Applied Mathematics, 48(3), 592-606. https://doi.org/10.1137/0148033


Wu, X., Luo, Y., & Hu, Y. (2014). Traveling waves in a diffusive predator-prey model incorporating a prey refuge. In Abstract and Applied Analysis (Vol. 2014). Hindawi. https://doi.org/10.1155/2014/679131


Wu, S., Wang, J., & Shi, J. (2018). Dynamics and pattern formation of a diffusive predator–prey model with predator-taxis. Mathematical Models and Methods in Applied Sciences, 28(11), 2275-2312. https://doi.org/10.1142/s0218202518400158


Zhou, Y., Sun, W., Song, Y., Zheng, Z., Lu, J., & Chen, S. (2019). Hopf bifurcation analysis of a predator–prey model with Holling-II type functional response and a prey refuge. Nonlinear Dynamics, 97(2), 1439-1450. https://doi.org/10.1007/s11071-019-05063-w


Zhu, H., Campbell, S. A., & Wolkowicz, G. S. (2003). Bifurcation analysis of a predator-prey system with nonmonotonic functional response. SIAM Journal on Applied Mathematics, 63(2), 636-682. https://doi.org/10.1137/s0036139901397285


Zu, J., & Mimura, M. (2010). The impact of Allee effect on a predator–prey system with Holling type II functional response. Applied Mathematics and Computation, 217(7), 3542-3556. https://doi.org/10.1016/j.amc.2010.09.029