© Revista de Matemática: Teoría y Aplicaciones 2021 28(2) : 209–236 cimpa – ucr issn: 1409-2433 (Print), 2215-3373 (Online) doi: 10.15517/rmta.v28i2.38645

solución tipo onda viajera en un modelo

difusivo depredador–presa tipo holling ii

traveling wave type solution in a model

diffusive predator–prey type holling ii

Christian Cortés-García^{∗}  Allison Ramírez-Fierro^{†}

Received: 20/Nov/2020; Revised: 9/Mar/2021; Accepted: 12/Apr/2021

∗Universidad Carlos III de Madrid, Departamento de Matemática, Madrid, España. Centro Nacional de Biotecnología, Departamento de Biología de Sistemas, Madrid, España. E-Mail: chcortes@math.uc3m.es , cc.cortes@cnb.csic.es

†Universidad Surcolombiana, Departamento de Matemática y Estadística, Neiva, Colombia. E-Mail: maria21907@hotmail.com

210 c. cortés-garcía — a. ramírez-fierro

Resumen

En este trabajo se demuestra la existencia de ondas viajeras como soluciones para un modelo depredador–presa, con funcional de depredación Holling II y término difusivo unidimensional para los depredadores. Al realizar un análisis cualitativo al modelo sin difusión, se deduce que el modelo con difusión presenta soluciones periódicas. De igual forma, al asumir solución tipo onda viajera al modelo con difusión, se demuestra que posee una órbita heteroclínica que conecta dos puntos de equilibrio, atractora a uno de ellos, y por tanto presenta frentes de ondas.

Palabras clave: modelo de Gause; ciclo límite; teorema de Hartman Grobman; principio de LaSalle; teorema de la bifurcación de Hopf.

Abstract

This paper demonstrates the existence of traveling waves as solutions for a predator - prey model with a Holling II predation function and a onedimensional diffusive term for predators. When performing a qualitative analysis on the model without diffusion, it follows that the model with diffusion presents periodic solutions. Similarly, by assuming a traveling wave-type solution to the diffusion model, it is shown that it has a heteroclinical orbit that connects two equilibrium points, attracted to one of them, and therefore presents wave fronts.

Keywords: Gause model; limit cycle; Hartman Grobman theorem; LaSalle principle; Hopf bifurcation theorem.

Mathematics Subject Classification: 35B10, 35B32, 37C27, 65N06, 92B05.

1 Introducción

Uno de los temas en el área de la biología matemática es determinar y estudiar la dinámica existente entre la interacción de presas y depredadores bajo diversas hipótesis, representados normalmente por sistemas de ecuaciones diferenciales.

Los sistemas tipo depredador–presa exhiben amplios comportamientos dinámicos, los cuales dependen de los parámetros del modelo, que van desde cambios de estabilidad local en puntos de equilibrio hasta la existencia de ciclos límites y órbitas heteroclínicas u homoclínicas.

Al considerar la difusión espacial, en al menos una de las poblaciones, se hace una contribución significativa a la dinámica de las poblaciones, especialmente a la propagación de especies. Según Kuznetsov [6], la interacción de las presas 0 ≤ u ≤ K y depredadores v ≥ 0 puede ser descrita mediante el modelo:

\{ \begin{matrix}u_{t} = D_{1}u_{xx} + ru (1 − \frac{u}{K} ) − \frac{buv}{1 + eu}, \\ v_{t} = D_{2}v_{xx} + \frac{muv}{1 + eu} − dv, \end{matrix}  (1)

con D_{1} > 0 y D_{2} > 0 las constantes de difusión de las especies, respectivamente, r > 0 la tasa de crecimiento de presas en ausencia de depredación, K > 0 la capacidad de carga de las presas, b > 0 y m > 0 las tasas de encuentro entre las dos especies, e > 0 la constante de saturación de las presas y d > 0 la tasa de mortalidad de los depredadores.

Por otro lado, la existencia de soluciones tipo ondas viajeras para modelos tipo depredador–presa con difusión, han sido estudiadas por diversos investigadores. Por ejemplo, Dunbar [1] prueba la existencia de soluciones tipo onda viajera para un modelo equivalente a (1), con re-escalamiento de variables y reducción de sus parámetros. De igual forma, Wan-Tong [8] analiza la existencia de ondas viajeras para un modelo equivalente a (1) con funcional de respuesta tipo Holling III y, Chufen Wu [11] para un funcional tipo Holling de la forma f(u, v) = \frac{u^{p}v}{a+u^{p}}, p > 0, y reducción de sus parámetros.

Al considerar que las presas se distribuyen uniformemente con el fin de simplificar el análisis al modelo (1), esto es D_{1} = 0, y a diferencia de los trabajos presentados por Chufen Wu [11] y Dunbar [1], en este trabajo se analiza la existencia de soluciones tipo ondas viajeras con re-escalamiento de la variable espacial y sin reducción de los parámetros en el modelo sin difusión. Para ello, se desarrollan los siguientes pasos: en la segunda sección se dan algunos previos para probar los resultados de este trabajo; en la tercera sección se hace un análisis cualitativo local y global del modelo (1) sin difusión con el fin de determinar los equilibrios espaciales y una primera inspección para demostrar la existencia de soluciones tipo ondas viajeras, dadas finalmente, en la sección cuatro.

2 Preliminares

2.1 Sistemas no lineales

Considere un sistema autónomo de la forma

ẋ = f(x),  (2)

con f : Ω ⊂ ℝ^{n} ⟶ ℝ^{n} una función no lineal de clase C^{1}, y x_{∗} ∈ ℝ^{n} un punto de equilibrio hiperbólico, esto es, todos los valores propios de la matriz jacobiana A(x_{∗}) del sistema (2) calculada en el punto x_{∗} poseen parte real diferente de cero.

El siguiente teorema determina la estabilidad local del punto de equilibrio hiperbólico x_{∗} ∈ Ω del sistema (2) a partir del análisis de los valores propios de la matriz jacobiana A(x_{∗}), cuya demostración se puede consultar en el trabajo de Sotomayor [10].

Teorema 2.1 (Hartman Grobman) Sea x_{∗} un punto de equilibrio hiperbólico del sistema (2). Entonces existe un homeomorfismo H : U ⟶ V entre las vecindades U de x_{∗} y V de 0 ∈ ℝ^{n} tal que

H ∘ ϕ(x, t) = e^{tA} ∘ H,

donde A = Df(x_{∗}) denota la matriz jacobiana del sistema evaluada en x_{∗} y ϕ(x, t) la trayectoria del sistema (2).

LaSalle [7] muestra condiciones para probar la estabilidad del equilibrio x_{∗} en el sistema (2).

Teorema 2.2 (Principio de invarianza de LaSalle) Considere un punto de equilibrio x_{∗} ∈ Ω y una función de Lyapunov £(t) del sistema (2) en una región Ω alrededor de x_{∗}. Suponga que,

i. \frac{d£}{dt} ≤ 0 para x ∈ Ω.

ii. Una solución para la cual \frac{d£}{dt} = 0 es x(t) = x_{∗}.

Entonces x_{∗} es asintóticamente estable.

Por otro lado, dos sistemas de ecuaciones diferenciales de tipo (2) son topológicamente equivalentes, si existe un homeomorfismo h : ℝ^{n} → ℝ^{n} el cual realice un mapeo de trayectorias del primer sistema en trayectorias del segundo sistema, preservando la dirección del tiempo.

De esta manera, y como lo plantea Kuznetsov [6], el sistema

\{ \begin{matrix}ẋ = f(x, α), \\ α̇ = 0, \end{matrix}  (3)

con x ∈ ℝ^{n} y α ∈ ℝ un parámetro que no depende del tiempo, presenta una bifurcación en α = α_{0} si el sistema ẋ = f(x, α_{0}) no es topológicamente equivalente en virtud a pequeñas variaciones del parámetro α ≠ α_{0} sobre el sistema (3).

En particular, el sistema (3) presenta una bifurcación de Hopf cuando existen un par de valores propios complejos conjugados con parte real nula en la matriz jacobiana de f(x, α), calculada en el punto de equilibrio x_{∗} para algún α_{∗}, esto es, λ_{1,2} = ±iw_{0} con w_{0} > 0, como lo observado en el siguiente Teorema, probado por Kuznetsov [6].

Teorema 2.3 (Bifurcación de Hopf) Considere el sistema ẋ = f(x, α) con x ∈ ℝ^{n} y α ∈ ℝ. Suponga que existe (x_{∗}, α_{∗}) tal que:

i. f(x_{∗}, α_{∗}) = 0 (condición de equilibrio).

ii. Df(x_{∗}, α_{∗}) posee un único par de valores propios en el eje imaginario λ_{1,2}(α_{∗}) = ±iw_{0} y el resto están fuera de él (condición de bifurcación).

iii. \frac{d}{dα}(Re(λ_{1,2}(α)))|_{α=α_{∗}} ≠ 0 (condición de transversalidad).

Entonces el sistema ẋ = f(x, α) posee soluciones periódicas con pequeña amplitud de periodo T_{0} = \frac{2π}{Imλ_{1}(α_{0})}.

Por otro lado, considere la ecuación diferencial

ẋ = f(t, x),  (4)

donde f(t, x) es una función continua definida en un conjunto abierto (t, x) ∈ Ω, Ω_{0} un subconjunto abierto en Ω con ∂Ω_{0} su frontera y \overline{Ω_{0}} su clausura.

Se dice que (t_{0}, x_{0}) es un punto de egreso de Ω_{0} respecto al sistema (4), si para toda solución x = x(t) del sistema (4) del cual x(t_{0}) = x_{0}, existe ϵ > 0 tal que (t, x(t)) ∈ Ω_{0} para t_{0} − ϵ ≤ t < t_{0}. Un punto de egreso (t_{0}, x_{0}) de Ω_{0} es un punto de egreso estricto si (t, x(t)) ∉ \overline{Ω_{0}} para t_{0} < t < t_{0} + ϵ con ϵ > 0 pequeño. Considere Ω_{0}^{e} como el conjunto de puntos de egreso de Ω_{0}, y Ω_{0}^{se} como el conjunto de puntos de egreso estricto de Ω_{0}. Como condición suficiente, Ω_{0}^{e} se puede determinar al verificar ⟨l, n⟩ > 0, donde l es el campo vectorial del sistema (4) y n el vector normal exterior de ∂Ω_{0}.

Si U es un espacio topológico y V un subconjunto de U, la función continua π : U → V es una retracción de U a V, si la restricción π de U a V es la identidad, esto es, π(u) ∈ U y π(v) = v para todo v ∈ V. Si existe una retracción de U a V, V es llamado un retractor de U.

El resultado, probado por Hartman [4], muestra lo siguiente,

Lema 2.1 Sea ϕ(t, x), con condición inicial x(t_{0}) = x_{0}, una solución única para el sistema (4). Consideren Ω_{0} ⊂ Ω como un subconjunto abierto, con Ω_{0}^{e} = Ω_{0}^{se}, y S ⊂ Ω_{0} ∪ Ω_{0}^{e} no vacio tal que S ∩ Ω_{0}^{e} no es un retractor de S pero si de Ω_{0}^{e}. Entonces, existe al menos (t_{0}, x_{0}) ∈ S ∩Ω_{0}^{e} tal que la solución ϕ(t, x), con condición inicial x(t_{0}) = x_{0}, está contenido en Ω_{0} en un intervalo máximo de existencia.

2.2 Solución tipo onda viajera en una ecuación reacción–difusión

Los sistemas de difusión son modelos matemáticos que describen densidades de partículas, sustancias u organismos cuando se dispersan a través del espacio como resultado del movimiento microscópico irregular de cada una de estas. Por otro lado, un proceso de reacción ocurre cuando las partículas presentan reacciones químicas, o procesos biológicos, debido a sus interacciones con el entorno.

Si u(x, t) representa la densidad de una sustancia o de una población en el tiempo t, un modelo de reacción–difusión para u viene dado por:

\frac{∂u}{∂t} = f(u) + δ \frac{∂^{2}u}{∂x^{2}},  (5)

con δ > 0 el coeficiente de difusión y f(u) la función de reacción.

Los modelos de reacción–difusión pueden presentar soluciones tipo ondas viajeras definidas como u(x, t) = U(x − ct) o u(x, t) = U(x + ct), donde la onda se desplaza hacia la derecha o izquierda, respectivamente, para una función U y velocidad de onda c > 0 constante. Si c = 0, la onda resultante se denomina onda estacionaria.

Las ondas viajeras se caracterizan según sus propiedades. Por ejemplo, para una onda viajera que se aproxima a los estados constantes U(−∞) = u_{l} y U(∞) = u_{r}, con u_{l} ≠ u_{r}, la onda se denomina un frente de onda. Si u_{l} = u_{r}, la onda correspondiente se conoce como onda de pulso. Si una onda muestra una periodicidad con U(z + F) = U(z), con F > 0, la onda se denomina onda periódica.

3 Modelo sin difusión

Sean 0 ≤ u ≤ K y v ≥ 0 la cantidad de presas y depredadores, respectivamente, cuya dinámica es representada por:

\{ \begin{matrix}u̇ = ru (1 − \frac{u}{K} ) − \frac{buv}{1 + eu}, \\ v̇ = \frac{muv}{1 + eu} − dv, \end{matrix}  (6)

donde

• f_{1}(u) = ru (1 − \frac{u}{K} ) describe el crecimiento logístico de las presas en ausencia de los depredadores, con r > 0 la tasa de crecimiento intrínseca de las presas y K > 0 la capacidad de carga.

• f_{2}(u, v) = −\frac{buv}{1+eu} representa la disminución en las presas al interactuar con los depredadores, con θ(u) = \frac{bu}{1+eu} un funcional tipo Holling II, donde b > 0 es la tasa de encuentro entre las dos especies y e > 0 la constante de saturación de las presas.

• f_{3}(u, v) = \frac{muv}{1+eu} es 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.

• f_{4}(v) = −dv describe el decrecimiento de los depredadores en ausencia de las presas, con d > 0 la tasa de mortalidad de los depredadores.

3.1 Análisis cualitativo local

Cuando el campo vectorial se iguala a cero, el sistema (6) presenta dos puntos de equilibrios: E_{1} = (0, 0), E_{2} = (K, 0) y, un punto de equilibrio interior

E_{∗} = (u_{∗}, v_{∗}) = (\frac{d}{m − de}, \frac{rm[mK − d(eK + 1)]}{bK(m − de)^{2}}) ,  (7)

si m > de y mK > d(eK + 1).

Para realizar un análisis cualitativo del sistema (6) se debe analizar la estabilidad local para cada punto de equilibrio. Para ello, al calcular la matriz jacobiana A(E_{1}) del sistema (6) evaluada en el equilibrio E_{1}, esto es,

A(E_{1}) = (\begin{matrix}r & 0 \\ 0 & −d \end{matrix}) ,

se tiene que DetA(E_{1}) = −rd < 0 y, por el Teorema 2.1, el equilibrio E_{1} es localmente un punto silla, con variedad estable W^{s}(E_{1}) en la dirección de v, pues para todo punto (0, v), v > 0, se tiene

\{ \begin{matrix}u̇ = 0, \\ v̇ = −dv < 0, \end{matrix}

y variedad inestable W^{u}(E_{1}) en la dirección de u, pues para todo (u, 0), con 0 < u < K, se tiene

\{ \begin{matrix}u̇ = ru (1 − \frac{u}{K} ) > 0, \\ v̇ = 0, \end{matrix}

De igual forma, la matriz jacobiana A(E_{2}) del sistema (6) asociada al punto E_{2} es dada por

A(E_{2}) = (\begin{matrix}−r & \frac{−bK}{1 + eK} \\ 0 & \frac{mK}{1 + eK} − d \end{matrix}) ,

con valores propios λ_{1} = −r < 0 y λ_{2} = \frac{mK−d(1+eK)}{1+eK}. Si K < \frac{d}{m−de}, entonces λ_{2} < 0 y por el Teorema 2.1, el equilibrio E_{2} es localmente un nodo estable. En cambio, si K > \frac{d}{m−de} se tiene λ_{2} > 0 y el punto E_{2} es localmente un punto silla. Del mismo razonamiento para E_{1}, si E_{2} es localmente un punto silla, la variedad estable W^{s}(E_{2}) se encuentra en el eje u.

Por último, la matriz jacobiana A del sistema (6) evaluada en el punto de equilibrio E_{∗} es dada por:

A = A(E_{∗}) = (\begin{matrix}\frac{rd [(eK + 1)(m − de) − 2m]}{mK(m − de)} & \frac{−bd}{m} \\ \frac{r [mK − d(eK + 1)]}{bK} & 0 \end{matrix}) ,

donde

TrA = \frac{rd[(eK + 1)(m − de) − 2m]}{mK(m − de)},
DetA = \frac{rd[mK − d(eK + 1)]}{mK},
Δ = (TrA)^{2} − 4DetA.

(8)

Como DetA > 0, el punto E_{∗} no es localmente un punto silla. Por otro lado, si K > \frac{de+m}{e(m−ed)} entonces TrA > 0 y caso contrario si K < \frac{de+m}{e(m−ed)}. En consecuencia, dichas conclusiones se resumen en el siguiente resultado.

Lema 3.1 El sistema (6) tiene dos puntos de equilibrios E_{1} = (0, 0) y E_{2} = (K, 0). Si m > de y mK > d(eK + 1), el sistema (6) posee además un punto de equilibrio interior E_{∗} = (u_{∗}, v_{∗}) dado en (7). Por otro lado, las condiciones de estabilidad local de estos puntos de equilibrio son dados por:

• El punto (0, 0) es un punto silla.

• El punto (K, 0) es

– silla si K > \frac{d}{m−de},

– nodo estable si K < \frac{d}{m−de}.

• El punto (u_{∗}, v_{∗}) es

– espiral inestable si Δ < 0 y K > \frac{de+m}{e(m−ed)},

– espiral estable si Δ < 0 y K < \frac{de+m}{e(m−ed)},

– nodo estable si Δ > 0 y K < \frac{de+m}{e(m−ed)},

– nodo inestable si Δ > 0 y K > \frac{de+m}{e(m−ed)}.

3.2 Análisis cualitativo global

Como el sistema (6) es equivalente al modelo de Gause

\{ \begin{matrix}u̇ = ug(u) − vp(u), \\ v̇ = v[−γ + q(u)], \end{matrix}  (9)

el sistema (9) presenta un único ciclo límite globalmente estable si satisface la siguiente condición, probado por Yang Kuang [5].

Teorema 3.1 Considere un modelo de Gause (9) para todo u ≥ 0, v ≥ 0,

γ > 0 y g, p, q funciones de clase C^{1}. Si \frac{d}{du} (^{\frac{ug′(u)+g(u)−ug(u) \frac{p′(u)}{p(u)}}{−γ+q(u)}} ) ≤ 0,

en 0 ≤ u < u_{∗} y u_{∗} < u ≤ K, el sistema (9) tiene un único ciclo límite globalmente asintóticamente estable en {(u, v) : u > 0, v > 0} \ {(u_{∗}, v_{∗})}, donde E_{∗} = (u_{∗}, v_{∗}) es un punto interior del sistema (9).

Por consiguiente, el sistema (6) presenta un ciclo límite globalmente estable si E_{∗} es una espiral inestable, como lo enunciado en el siguiente resultado.

Lema 3.2 Si Δ < 0 y K > \frac{de+m}{e(m−ed)}, entonces E_{∗} es una espiral inestable y existe un único ciclo límite globalmente asintóticamente estable para el modelo (6) en el primer cuadrante.

Demostración. Considere

g(u) = r (1 − \frac{u}{K} ) , g′(u) = −\frac{r}{K} , γ = d, q(u) = \frac{mu}{1 + eu},

con

\begin{aligned}\frac{d}{du}(\frac{ug′(u)+g(u)−ug(u)\frac{p′(u)}{p(u)}}{−γ+q(u)}) & = & −\frac{d}{du}(\frac{\frac{2reu^{2}}{K} + \frac{ru}{K} − reu}{u − u_{∗}} ) \\ & = & − (\frac{2reu^{2}}{K} − \frac{4reuu_{∗}}{K} − \frac{ru_{∗}}{K} +reu_{∗}). \end{aligned}

Por el Teorema 3.1 se debe verificar que

\frac{d}{du} (\frac{ug′(u) + g(u) − ug(u)\frac{p′(u)}{p(u)}}{−γ + q(u)} ) ≤ 0.  (10)

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 28(2): 209–236, Jul–Dec 2021

En efecto, (10) se mantiene siempre que

−\frac{2reu^{2}}{K} + \frac{4reuu_{∗}}{K} + \frac{ru_{∗}}{K} − reu_{∗} ≤ 0.

Sea

H(u) = −\frac{2reu^{2}}{K} + \frac{4reuu_{∗}}{K} + \frac{ru_{∗}}{K} − reu_{∗},

una parábola que abre hacia abajo, con discriminante

ΔH(u) = \frac{8r^{2}eu_{∗}}{K^{2}} (2eu_{∗} + 1 − eK).

Como K > \frac{de+m}{e(m−ed)} se deduce u_{∗} ≤ \frac{eK−1}{2e}, entonces 2eu_{∗} + 1 − eK ≤ 0 y ΔH(u) ≤ 0. Por tanto H(u) tiene como máximo una raíz real, lo que indica que H(u) ≤ 0, manteniéndose la desigualdad. En consecuencia, el sistema (6) posee un único ciclo límite globalmente asintóticamente estable.

Por otro lado, si el punto E_{2} = (K, 0) es un nodo estable, esto es mK < d(eK + 1), el punto E_{∗} no pertenece al sistema (6) y, por consiguiente, E_{2} es un punto globalmente asintóticamente estable. De igual forma, si mK > d(eK +1) y \frac{d}{m−de} < K < \frac{de+m}{e(m−ed)}, los puntos E_{2} y E_{∗} coexisten en el sistema (6), con E_{2} localmente un punto silla y E_{∗} un punto globalmente asintóticamente estable. Por consiguiente, se deduce el siguiente resultado.

Lema 3.3 Si K > d(eK + 1), el equilibrio E_{2} del sistema (6) es globalmente asintóticamente estable en el primer cuadrante. Si

\frac{d}{m − de} < K < \frac{de + m}{e(m − ed)},

el equilibrio E_{∗} del sistema (6) es globalmente asintóticamente estable.

3.3 Análisis de bifurcación

Como el sistema (6) podría presentar tres puntos de equilibrio y un ciclo límite, se debe analizar en qué situaciones estos elementos colisionan entre sí para formar bifurcaciones.

En primera instancia, el sistema (6) presenta una bifurcación de Hopf al colisionar el punto E_{∗} con el ciclo límite, esto es, cuando TrA = 0, equivalente a K = \frac{de+m}{e(m−ed)}, como lo observado en el siguiente resultado.

Lema 3.4 Si K = \frac{de+m}{e(m−ed)}, entonces el sistema (6) presenta una bifurcación de Hopf con periodo

T_{0} = \frac{2π}{\sqrt{\frac{rd(mK−d(eK+1))}{mK}}}

.

Demostración. De acuerdo con el Teorema 2.3, se debe garantizar que el sistema (6) cumpla con las siguientes condiciones:

i. Condición de equilibrio.

ii. Condición de bifurcación.

iii. Condición de transversalidad o de control.

En efecto, sea K = K_{0} := \frac{de+m}{e(m−ed)} y (u_{0}, v_{0}) = (\frac{d}{m−de}, \frac{m^{2}r}{b(de+m)(m−de)}), esto es, el punto E_{∗} al reemplazar K = K_{0}. Por consiguiente, el sistema (6) evaluado en (K_{0}, u_{0}, v_{0}) es igual a cero, lo cual satisface la condición de equilibrio.

Se sabe que los valores propios de la matriz jacobiana A evaluada en el punto E_{∗} son dados por λ_{1,2} = \frac{TrA±\sqrt{TrA^{2}−4DetA}}{2}. Si K = K_{0}, entonces TrA = 0, DetA > 0 y de (8) se tiene que,

λ_{1,2}(K_{0}) = ±i\sqrt{\frac{rd(mK − d(eK + 1))}{mK}}.

Por consiguiente, se presentan un par de valores propios imaginarios puros, lo que satisface la condición de bifurcación.

Si TrA^{2} < 4DetA, entonces Re(λ(K)) = \frac{TrA}{2} y de (8) se tiene que:

\frac{d}{dK} (Re(λ(K))) |_{K=K_{0}} = \frac{dre^{2}(m − de)}{2m(de + m)} > 0,

pues m > de. Por tanto se satisface la condición de transversalidad y por ende el sistema (6) presenta una bifurcación de Hopf con periodo T_{0} = _{\frac{2π}{\sqrt{\frac{rd(mK−d(eK+1))}{mK}}}}.

Por otro lado, al colisionar el punto E_{∗} con el punto E_{2}, esto es, cuando K = \frac{d}{m−de}, el sistema (6) presenta una bifurcación transcrítica.

Nótese que el equilibrio E_{1} no colisiona con E_{2} y E_{∗}, debido a que K > 0 y E_{∗} es punto interior. Por lo tanto, el sistema (6) presenta 3 curvas de bifurcación dadas por K = \frac{d}{m−de}, K = \frac{de+m}{e(m−ed)} y Δ = 0.

220 c. cortés-garcía — a. ramírez-fierro

(a) Plano K − d

(b) Región (1)  (c) Región (2)

(d) Región (3)  (e) Región (4)

Figura 1: Curvas de bifurcación y retratos de fase del modelo (6) con parámetros fijos r = 1, b = 1, m = 0.5 y e = 0.125, y parámetros de bifurcación: (b) K = 1.12 y d = 1.2 ; (c) K = 3 y d = 1; (d) K = 0.8 y d = 0.1; (e) K = 12 y d = 0.39.

En la Figura 1 se observan las curvas de bifurcación del modelo (6) con sus retratos de fase. En la Región (1), el equilibrio (K, 0) es localmente un nodo estable, y al cruzar la curva K = \frac{d}{m−de} se observa la formación de un equilibrio E_{∗} localmente un nodo estable dado en la Región (2). De igual forma, en la Región (2) al cruzar la curva Δ = 0, el equilibrio E_{∗} cambia de nodo estable a un foco estable, y a medida que cruza la curva K = \frac{de+m}{e(m−ed)} se presenta un ciclo límite estable como se muestra en la Región (4). Nótese que no es posible la estabilidad local de E_{∗} en un nodo inestable.

4 Modelo con difusión para los depredadores

Cuando los depredadores v(x, t) se distribuyen en un determinado espacio unidimensional, el modelo (6) se modifica como

\{ \begin{matrix}u_{t} = ru (1 − \frac{u}{K} ) − \frac{buv}{1 + eu}, \\ v_{t} = Dv_{xx} + \frac{muv}{1 + eu} − dv, \end{matrix}  (11)

con v_{t} = \frac{∂v}{∂t}, v_{xx} = \frac{∂^{2}v}{∂x^{2}} y D ≥ 0 el coeficiente de difusión.

Al considerar x̃ = \sqrt{\frac{1}{D}} x, en el sistema (11) se tiene

\{ \begin{matrix}u_{t} = ru (1 − \frac{u}{K} ) − \frac{buv}{1 + eu}, \\ v_{t} = v_{xx} + \frac{muv}{1 + eu} − dv. \end{matrix}  (12)

donde se asume que x := x̃.

Una forma de analizar la dinámica de un sistema de ecuaciones diferenciales parciales es convertirlo a un sistema de ecuaciones diferenciales ordinarias por medio de ondas viajeras dado en la sección 2.2. Este método se conoce cuando el sistema es reemplazado por

(u(x, t), v(x, t)) = (U(s), V (s)),  (13)

donde s = x + ct, y c > 0 es la velocidad de onda.

Al usar (13) en el sistema (12) se tiene

\{ \begin{matrix}cU′(s) = rU (1 − \frac{U}{K} ) − \frac{bUV}{1 + eU} , \\ cV′(s) = V′′(s) + \frac{mUV}{1 + eU} − dV, \end{matrix}  (14)

el cual, al realizar un cambio de variables s = cσ y ϵ = \frac{1}{c^{2}} , el sistema (14) se reescribe de la siguiente forma

\{ \begin{matrix}U′(σ) = rU (1 − \frac{U}{K} ) − \frac{bUV}{1 + eU} , \\ V′(σ) = ϵV′′(σ) + \frac{mUV}{(1 + eU)} − dV. \end{matrix}  (15)

Si se toma ϵ como un parámetro cercano a cero, el sistema (15) se reduce a

\{ \begin{matrix}U′(σ) = rU (1 − \frac{U}{K} ) − \frac{bUV}{1 + eU} , \\ V′(σ) = \frac{mUV}{(1 + eU)} − dV, \end{matrix}  (16)

equivalente al sistema (6).

Si K > \frac{de+m}{e(m−ed)}, el sistema (16) tiene un único ciclo de límite globalmente estable con período T_{0}, como lo dado en el Lema 3.4. Por lo tanto, al considerar c^{2} lo suficientemente grande, el sistema (14) debe tener una solución periódica con período cT_{0}. Esto proporciona una motivación para buscar soluciones de onda viajera periódica en el sistema (12).

4.1 Solución tipo onda viajera periódica

Al reescribir el sistema (14) como un sistema de tres ecuaciones, se tiene

\{ \begin{matrix}U′(s) = \frac{r}{c} U (1 − \frac{U}{K} ) − \frac{bUV}{c(1 + eU)}, \\ V′(s) = Z, \\ Z′(s) = cZ + dV − \frac{mUV}{1 + eU} , \end{matrix}  (17)

con dos puntos de equilibrios: \overline{E}_{1} = (0, 0, 0), \overline{E}_{2} = (K, 0, 0). Si m > de y mK > d(eK + 1), el sistema (17) posee un punto de equilibrio interior \overline{E}_{∗} = (u_{∗}, v_{∗}, 0), con u_{∗} y v_{∗} dados en (7).

Obsérvese que \overline{E}_{1} es localmente un punto silla, pues la matriz jacobiana del sistema (17) evaluada en el punto de equilibrio \overline{E}_{1} está dada por

A(\overline{E}_{1}) = (\begin{matrix}\frac{r}{c} & 0 & 0 \\ 0 & 0 & 1 \\ 0 & d & c \end{matrix}) ,

con valores propios: λ_{1} = \frac{r}{c} > 0, λ_{2} = \frac{c+\sqrt{c^{2}+4d}}{2} > 0 y λ_{3} = \frac{c−\sqrt{c^{2}+4d}}{2} < 0, pues c^{2} < c^{2} + 4d con c, d > 0.

De igual forma, la estabilidad local del punto de equilibrio \overline{E}_{2} está sujeta a la naturaleza de los valores propios de la matriz jacobiana del sistema (17) evaluada en el punto \overline{E}_{2}. En efecto, dicha matriz es dada por

A(\overline{E}_{2}) = (\begin{matrix}−\frac{r}{c} & −\frac{bK}{c(eK + 1)} & 0 \\ 0 & 0 & 1 \\ 0 & d − \frac{mK}{eK + 1} & c \end{matrix}) ,

con valores propios:

λ_{1} = −\frac{r}{c} < 0 y,
λ_{2,3} = \frac{c}{2} ∓ \frac{\sqrt{c^{2} − 4ρ}}{2},  (18)

donde ρ = \frac{mK−d(eK+1)}{eK+1} > 0. Si c ≥ 2\sqrt{ρ}, entonces λ_{2,3} ∈ ℝ y λ_{1} < 0 < λ_{2} < λ_{3}. Por el Teorema 2.1, el equilibrio \overline{E}_{2} es localmente un punto silla, con variedad estable en dirección de U.

Para el punto de equilibrio \overline{E}_{∗}, y a diferencia de los resultados por Chufen Wu [11] y Dunbar [1], se tiene el siguiente resultado:

Teorema 4.1 Consideren m > de, mK > d(eK + 1), K > \frac{de+m}{e(m−ed)},

M := \frac{dr[de(eK + 1) − m(eK − 1)]}{mK(m − de)} y N := \frac{dr[d(eK + 1) − mK]}{mK}.

Cuando el parámetro K cruza la curva c^{2} = M + \frac{N}{M} para algún K_{0} en el plano K − c^{2}, el punto de equilibrio \overline{E}_{∗} del sistema (17) presenta una bifurcación de Hopf, con solución periódica de período,

T (K_{0}) = \frac{2π}{\sqrt{−M}} + O (K − K_{0}) .

Esta solución periódica corresponde a una pequeña amplitud de solución de ondas de (12).

Demostración. La matriz jacobiana del sistema (17) en el punto de equilibrio \overline{E}_{∗} = (u_{∗}, v_{∗}, 0) es dada por,

A(\overline{E}_{∗}) = (\begin{matrix}\frac{r(K − 2u_{∗})}{cK} − \frac{bv_{∗}}{c(eu_{∗} + 1)^{2}} & −\frac{bd}{cm} & 0 \\ 0 & 0 & 1 \\ \frac{r[d(eK + 1) − mK]}{bK} & 0 & c \end{matrix}) ,

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 28(2): 209–236, Jul–Dec 2021

con polinomio característico

F (K, λ) := −λ^{3} + (c − \frac{M}{c} ) λ^{2} + Mλ − \frac{N}{c} .  (19)

La bifurcación de Hopf ocurre cuando hay un par de raíces imaginarias en el polinomio (19). Para esto, al sustituir λ_{0} = iφ en (19),

F (K, λ_{0}) = \frac{φ^{2}(M − c^{2})}{c} − \frac{N}{c} + i(φ^{3} + φM),

y al igualar la parte real y la parte imaginaria a cero se tiene:

φ^{2} = −M, φ^{2} = \frac{N}{M − c^{2}} .

Por lo tanto, existen un par de valores propios imaginarios si los parámetros K y c satisfacen la condición:

− M = \frac{N}{M − c^{2}} ⟺ c^{2} = M + \frac{N}{M} .  (20)

Al considerar λ = λ(K) como una función de K y al diferenciar (19) implícitamente con respecto a K, se obtiene:

\frac{dλ}{dK} = − \frac{\frac{∂F (K, λ)}{∂K}}{\frac{∂F (K, λ)}{∂λ}} = \frac{M′λ^{2} − cM′λ + N′}{−3cλ^{2} + 2 (c^{2} − M) λ + cM} ,

donde

M′ = \frac{dM}{dK} y N′ = \frac{dN}{dK} .

Ahora, al sustituir λ = iφ en la anterior ecuación y al usar el hecho que λ^{2} = −φ^{2} = M dado en (20), se tiene

\frac{dλ}{dK} = \frac{MM′ − cM′iφ + N′}{−2cM + 2(c^{2} − M)iφ}
= \frac{−c [M′φ^{2}(c^{2} − M) + M(MM′ + N′)]}{2[φ^{2}(c^{2} − M)^{2} + c^{2}M^{2}]}
+ i \{\frac{φ[M(MM′ + N′) − c^{2}N]}{2 [φ^{2}(c^{2} − M)^{2} + c^{2}M^{2}]}\} ,

Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 28(2): 209–236, Jul–Dec 2021

donde

Re (\frac{dλ}{dK} ) = \frac{−Mc (MM′ − \frac{N}{M} M′ + N′)}{2(Mc)^{2} − 2M(c^{2} − M)^{2}}.

Como K > \frac{de+m}{e(m−de)} entonces de(eK + 1) − m(eK − 1) < 0 y así M < 0.

Por tanto −Mc > 0 y 2(Mc)^{2} − 2M(c^{2} − M)^{2} > 0. Así, el signo de Re (\frac{dλ}{dK} )

es determinado por el signo de MM′ − \frac{N}{M} M′ + N′.

Se sabe de (20) que

\frac{dc^{2}}{dK} = M′ + \frac{N′}{M} − \frac{N}{M^{2}} M′,

entonces

M (\frac{dc^{2}}{dK} ) = MM′ − \frac{N}{M} M′ + N′,

y así, el signo de Re (\frac{dλ}{dK} ) es el opuesto de \frac{dc^{2}}{dK} para M < 0. Así, se debe ver

que \frac{dc^{2}}{dK} ≠ 0 para algún K_{0} para verificar que Re (\frac{dλ}{dK} ) ≠ 0 y concluir que el sistema (17) presenta una bifurcación de Hopf.

Como m > de, entonces

\begin{aligned}M′ & = & \frac{dr(de + m)}{mK^{2}(de − m)} < 0, \\ \frac{N′}{M} − \frac{N}{M^{2}} M′ & = & −\frac{m(de − m)^{2}}{[de(eK + 1) − m(eK − 1)]^{2}} < 0. \end{aligned}

Así, \frac{dc^{2}}{dK} < 0 para K > \frac{de+m}{e(m−de)} y Re (\frac{dλ}{dK} ) > 0. En la Figura 2 se observa el

comportamiento de la curva c^{2} = M + \frac{N}{M}, donde m^{2} − de(er + m) > 0 pues m > de.

Esto muestra que un par de raíces conjugadas pasan transversalmente desde la mitad izquierda del plano al medio plano derecho a medida que el parámetro K pasa la curva de bifurcación c^{2} = M + \frac{N}{M} en algún K_{0} > \frac{de+m}{e(m−de)}. Ahora, por el Teorema 2.3, el sistema (17) presenta una bifurcación de Hopf.

Por otro lado, el sistema (17) presenta una bifurcación transcrítica cuando colisiona el punto \overline{E}_{∗} con el punto \overline{E}_{2}, es decir, cuando K = \frac{d}{m−de}.

Figura 2: Plano K − c^{2}.

4.2 Existencia de un frente de onda

Debido a la naturaleza del sistema (12), se debe garantizar la existencia de soluciones no negativas que satisfacen las condiciones de frontera

\begin{aligned}lim_{s→−∞}(U(s), V (s)) & = & (K, 0), \\ lim_{s→∞}(U(s), V (s)) & = & (u_{∗}, v_{∗}). \end{aligned}  (21)

Una solución de este tipo, es una órbita heteroclínica punto a punto del sistema (14) que corresponde a una solución de frente de onda del sistema (12), que conecta dos equilibrios, \overline{E}_{2} y \overline{E}_{∗}.

El primer resultado garantiza condiciones c > 0 tal que no satisface soluciones de ondas viajeras para el sistema (14) que satisfacen las condiciones de frontera (21).

Teorema 4.2 Si c ∈ (0, 2\sqrt{ρ}), no existen soluciones no-negativas del sistema (17) que satisfacen las condiciones de frontera (21).

Demostración. Si c < 2\sqrt{ρ}, los valores propios de la matriz jacobiana del sistema (20) evaluada en \overline{E}_{2} están dados por:

\begin{aligned}λ_{1} & = & −\frac{r}{c} < 0, \\ λ_{2,3} & = & \frac{1}{2} (c ∓ i\sqrt{4ρ − c^{2}}) , \end{aligned}

lo cual expresa, localmente, un punto silla, con λ_{1} < 0 y λ_{2,3} con parte real positiva. Por tanto, existe una variedad estable a lo largo del eje u que tiende al punto de equilibrio \overline{E}_{2} y una variedad inestable de dos dimensiones que se alejan de él y todas las otras soluciones son una combinación de estos comportamientos, por lo que se mueven en espiral alrededor de la línea recta de soluciones como se muestra en la Figura 3. Esto muestra que, para condiciones iniciales cerca del punto \overline{E}_{2},

lim_{s→−∞}(U(s), V (s), Z(s)) = (K, 0, 0),  (22)

se debería tener V (s) < 0 para algún s, lo cual contradice (21) para V (s) ≥ 0.

Figura 3: Estabilidad local de \overline{E}_{2} cuando c < 2ρ.

Considere c ≥ 2\sqrt{ρ}. Sea W = cV − V′, entonces por (14), se tiene un sistema equivalente a (14),

\{ \begin{matrix}U′(s) & = \frac{r}{c} U (1 − \frac{U}{K} ) − \frac{bUV}{c(1 + eU)}, \\ V′(s) & = cV − W, \\ W′(s) & = \frac{mUV}{1 + eU} − dV, \end{matrix}  (23)

con puntos de equilibrio dados por: \overline{E}_{1w} = \overline{E}_{1}, \overline{E}_{2w} = \overline{E}_{2} y \overline{E}_{∗w} = (u_{∗}, v_{∗}, cv_{∗}).

Para c ≥ 2\sqrt{ρ}, la matriz jacobiana del sistema (23) calculada en \overline{E}_{2} es dada por:

A(\overline{E}_{2}) = (\begin{matrix}−\frac{r}{c} & −\frac{bK}{c(eK + 1)} & 0 \\ 0 & c & −1 \\ 0 & \frac{mK}{eK + 1} − d & 0 \end{matrix}) ,

con valores propios:

\begin{aligned}λ_{1} & = & −\frac{r}{c} < 0, \\ λ_{2} & = & \frac{1}{2} (c − \sqrt{c^{2} − 4ρ}) , \\ λ_{3} & = & \frac{1}{2} (c + \sqrt{c^{2} − 4ρ}) , \end{aligned}

correspondientes a los valores propios de la linealización del sistema (17) en \overline{E}_{2}. Si c > 2\sqrt{ρ}, entonces λ_{1} < 0 < λ_{2} < λ_{3} y los vectores propios asociados a los valores propios λ_{2,3} son dados por:

X_{2} = X_{2}(λ_{2}) = (−\frac{bK}{(cλ_{2} + r)(eK + 1)}, 1, c − λ_{2})^{T} ,

X_{3} = X_{3}(λ_{3}) = (−\frac{bK}{(cλ_{3} + r)(eK + 1)}, 1, c − λ_{3})^{T} ,

con un variedad inestable local en dos dimensiones W^{u}(\overline{E}_{2}) tangente a Gen\{X_{2}, X_{3}\} en el equilibrio \overline{E}_{2}.

Como 0 < U < K y V > 0, se observa del sistema (23) que

− dV < W_{s} = V (\frac{mU}{1 + eU} − d) < V (\frac{mK}{1 + eK} − d) = ρV,  (24)

esto es, −V < W_{s} < ρV. Esto proporciona una motivación para intercalar el subsistema dado de (23) entre dos sistemas lineales dados por:

\{ \begin{matrix}V′(s) = cV − W, \\ W′(s) = ρV, \end{matrix}  (25)

y

\{ \begin{matrix}V′(s) = cV − W, \\ W′(s) = −dV, \end{matrix}  (26)

con punto de equilibrio (V, W) = (0, 0) y valores propios μ_{2,3} = \frac{1}{2} (c ∓ \sqrt{c^{2} + 4d}) para el subsistema (25) y λ_{2,3} para el subsistema (26).

Como

0 < c − λ_{3} = \frac{c}{2} − \frac{1}{2}\sqrt{c^{2} − 4ρ} < \frac{c}{2} < c < c − μ_{2},

entonces, el plano W = (c − μ_{2})V está por encima del plano W = (c − λ_{3})V. Por tanto, considere el siguiente conjunto ℑ limitado por los siguientes planos, como lo observado en la Figura 4,

\begin{aligned}S_{1} & = & \{(U, V, W) : 0 < U < K, V > 0, W = (c − λ_{3})V \} , \\ S_{2} & = & \{(U, V, W) : 0 < U < K, V > 0, W = (c − μ_{2})V \} , \\ S_{3} & = & \{(U, V, W) : U = 0, V > 0, W > 0\} , \\ S_{4} & = & \{(U, V, W) : U = K, V > 0, W > 0\} . \end{aligned}  (27)

Figura 4: Conjunto ℑ.

De (24) se observa que el equilibrio \overline{E}_{∗w} se encuentra en ℑ y los equilibrios \overline{E}_{1w} y \overline{E}_{2w} en la frontera de ℑ. Por otro lado, el siguiente resultado muestra que toda condición inicial que se encuentra en ℑ tiene trayectorias generadas por el sistema (23) que salen de ℑ en S_{1} o S_{2}.

Lema 4.1 Considere c ≥ 2\sqrt{ρ} y P_{0} = (U(s_{0}), V (s_{0}), W (s_{0})) un punto inicial en ℑ. Entonces, la trayectoria ϕ(s, P_{0}), con s ≥ s_{0}, generado por el sistema (23) solo puede salir de ℑ desde S_{1} o S_{2}.

Demostración. En primera instancia, las trayectorias con punto inicial P_{0} salen de ℑ si ⟨n_{i}, h⟩ > 0, donde n_{i} es el vector normal exterior de S_{i}, i = 1, 2, 3, 4, y h_{i} es el campo vectorial del sistema (23) evaluado en S_{i}. Como n_{1} = (0, c − λ_{3}, −1) y n_{2} = (0, μ_{2} − c, 1), entonces ⟨n_{i}, h(S_{i})⟩ > 0, con i = 1, 2, esto es, las trayectorias salen de S_{1} y S_{2}. Además, como h(S_{3}) = (0, cV − W, −dV) y n_{3} = (−1, 0, 0), entonces ⟨n_{3}, h(S_{3})⟩ = 0. Así, el conjunto S_{3} es invariante.

Por otro lado, como −\frac{bKV}{c(1+eK)} < 0 para todo V > 0 y n_{4} = (1, 0, 0), entonces ⟨n_{4}, h(S_{4})⟩ = −\frac{bKV}{c(1+eK)} < 0. Luego, las trayectorias entran a ℑ.

Para finalizar, se debe ver que las trayectorias no salen en el conjunto

H = S_{1} ∩ S_{2} = {(U, V, W) : 0 < U < K, V = W = 0} .

Como n = (0, 0, −1) es el vector normal de H, entonces ⟨n, h(H)⟩ = 0, por lo cual, el conjunto H es invariante. Por tanto, S_{1} y S_{2} son los únicos conjuntos de los cuales sale la trayectoria cuyo punto inicial pertenece a ℑ.

Del Teorema 4.2 y del Lema 4.1, se evidencia que las soluciones del sistema (23) son positivas bajo la condición inicial dado en (21). Por ende, las soluciones de tipo onda viajera del sistema (11) son positivas.

Por otro lado, como los vectores propios X_{2} y X_{3} están entre S_{1} y S_{2} y apuntan en ℑ, el conjunto inestable W^{u}(\overline{E}_{∗2}) intercepta con ℑ y forma una curva Γ en ℑ que se encuentra en W^{u}(\overline{E}_{∗2}) con dos puntos finales en S_{i}, i = 1, 2. Dichos puntos son llamados C_{1} ∈ S_{1} y C_{2} ∈ S_{2}. Por tanto, se tiene el siguiente resultado.

Lema 4.2 Existe un punto P_{∗} ∈ Γ tal que la trayectoria ϕ(s, P_{∗}) del sistema (23) permanecerá en ℑ para todo s ≥ s_{∗}.

Demostración. Como el sistema (11) es autónomo, los conjuntos en el Lema 2.1 no dependen de t. Sea Ω = ℝ y Ω_{0} = ℑ. Por el Lema 4.3, Ω_{0}^{e} = Ω_{0}^{se} = S_{1} ∪S_{2}. Dado que C_{1} y C_{2} pertenecen al conjunto de egreso S_{1} ∪S_{2}, para algún P_{∗} ∈ Γ, la trayectoria ϕ(s, P_{0}), con s ≥ s_{0}, presenta uno de los tres casos:

1. Toca S_{1} y sale de ℑ.

2. Toca S_{2} y sale fuera ℑ.

3. Permanece en ℑ.

Si P_{0} está cerca de C_{1}, el primer caso es válido. Si P_{0} está cerca de C_{2}, el segundo caso es verdadero. Como C_{1} ∪ C_{2} no es una retracción de Γ sino una retracción de Ω_{0}, por el Lema 2.1, existe un punto P_{∗} ∈ Γ tal que ϕ(s, P_{∗}) permanece en ℑ para todo s ≥ s_{∗}. Esta situación corresponde al caso (3).

Al despejar V en U′(s) = 0, se considera la isóclina para las presas

g(u) = \frac{ru (1 − \frac{u}{K} ) (1 + eu)}{bu} = \frac{r(K − u)(1 + eu)}{bK},

con

g′(u) = \frac{r(e(K − 2u) − 1)}{bK}.

Para \frac{K}{2} ≤ u ≤ K, se tiene que e(K − 2u) − 1 ≤ −1 < 0 y g′(u) < 0, esto es, g(u) es decreciente para u ∈ [\frac{K}{2} , K]. Por otro lado, si eK < 1 entonces

\frac{K− \frac{1}{e}}{2} < 0 y así, g(u) es decreciente para todo u ∈ (0, K].

Al considerar eK < 1, para u ≤ u_{∗} se tiene g(u_{∗}) ≤ g(u) y (u − u_{∗})[g(u) − g(u_{∗})] ≤ 0. Para u ≥ u_{∗} se tiene g(u) ≤ g(u_{∗}) y (u − u_{∗})[g(u) − g(u_{∗})] ≤ 0. Por tanto, se ha verificado el siguiente resultado.

Lema 4.3 Si m > de, mK > d(eK + 1) y eK < 1, entonces (u − u_{∗})[g(u) − g(u_{∗})] ≤ 0 para u ∈ (0, K].

El siguiente resultado muestra la convergencia de trayectorias con punto inicial P_{∗} ∈ Γ.

Lema 4.4 Sea P_{∗} ∈ Γ como lo definido en el Lema 4.2. Entonces ϕ(s, P_{∗}) → \overline{E}_{∗w} cuando s → ∞, con ϕ(s, P_{∗}) la trayectoria generada por el sistema (23).

Demostración. Considere la función de Lyapunov

£(s) = \frac{1}{b(m − de)} ∫_{u_{∗}}^{U} \frac{(m − de)ξ − d}{ξ}dξ + \frac{1}{m − de} [\frac{W}{c} − g(u_{∗})\frac{W}{cV} − g(u_{∗})lnV ] ,

\begin{aligned}\frac{d£}{ds} & = & \frac{∂£}{∂U} U′(s) + \frac{∂£}{∂V} V′(s) + \frac{∂£}{∂W} W′(s) \\ & = & \frac{U − u_{∗}}{c(1 + eU)}[g(U) − g(u_{∗})] − \frac{g(u_{∗})(W − cV)^{2}}{cV^{2}(m − de)} ≤ 0, \end{aligned}  (28)

\begin{aligned}£(0) ≥ £(s) & ≥ & \frac{1}{m − de} [\frac{(c − λ_{3})V}{c} − g(u_{∗})\frac{c − μ_{2}}{c} − g(u_{∗})lnV ] \\ & = & \frac{1}{m − de} [(1 − \frac{λ_{3}}{c} ) V − g(u_{∗}) (1 − \frac{μ_{2}}{c} ) − g(u_{∗}) ln V ] . \end{aligned}

donde U(s) = U(s, P_{∗}), V (s) = V (s, P_{∗}), W (s) = W (s, P_{∗}). Del Lema 4.3, esto es, £(s) es no creciente. Como (U(s), V (s), W (s)) ∈ ℑ, se tiene que (c−λ_{3})V ≤ W ≤ (c − μ_{2})V, esto es, \frac{W}{V} ≤ c − μ_{2} y (c − λ_{3})V ≤ W; de (28) se tiene Rev.Mate.Teor.Aplic. (ISSN print: 1409-2433; online: 2215-3373) Vol. 28(2): 209–236, Jul–Dec 2021

Así (m − de)£(0) + g(u_{∗}) (1 − \frac{μ_{2}}{c} ) ≥ (1 − \frac{λ_{3}}{c} ) V − g(u_{∗}) ln V. Esto implica que existe un K_{1} > 0 lo suficientemente grande tal que, 0 < V (s) ≤ K_{1} y 0 < W (s) ≤ (c − μ_{2})K_{1}. Como 0 < U < K, la solución delimitada (U(s), V (s), W (s)) muestra que \frac{d£}{ds} → 0 cuando s → ∞. Además, \frac{d£}{ds} = 0 si y solo si (U, V, W) ∈ Λ = {U = u_{∗}, W = cV, 0 < V < K_{1}}.

Como el mayor conjunto invariante Λ en ℑ es el punto de equilibrio \overline{E}_{∗w} = (u_{∗}, v_{∗}, cv_{∗}), por el Teorema 2.2 se concluye que, ϕ(s, P_{∗}) → \overline{E}_{∗w} cuando s → ∞.

Para finalizar, al elegir un punto P_{∗} ∈ Γ, por los Lemas 4.2, 4.3 y 4.4, se concluye que la trayectoria del sistema (23) permanecerá en ℑ y se aproxima al equilibrio interior \overline{E}_{∗w} = (u_{∗}, v_{∗}, cv_{∗}) cuando s → ∞. De igual forma, la trayectoria tiende a \overline{E}_{2∗} = (K, 0, 0) cuando s → −∞ cuando P_{∗} ∈ W^{u}(\overline{E}_{∗2}). Por lo tanto, y a diferencia de los resultados planteados por Chufen Wu [11] y Dunbar [1] para un modelo equivalente a (11) con reducción de sus parámetros y re-escalamiento de la variable espacial y temporal, se ha construido una solución de onda viajera para c ≥ 2\sqrt{ρ}, como lo evidenciado en el siguiente Teorema.

Teorema 4.3 Si eK < 1 y c ≥ 2\sqrt{ρ}, existe una solución no-negativa del sistema (23) que satisface las condiciones de contorno (21), que corresponden a soluciones de frente de ondas del sistema (11).

(a) Depredador  (b) Presa

Figura 5: Solución aproximada tipo onda viajera del modelo (11) en diferentes momentos T, Δt = 1/24, h = 0.1, con coeficientes D = 1, r = 1, K = 3, b = 1, e = 0.125, m = 0.5, d = 1 y condiciones iniciales u(x, 0) = 3 + 1/(1 + exp(x − 150)) y v(x, 0) = 1/(1 + exp(x − 150)^{2}).

onda viajera en un modelo difusivo depredador-presa… 233

(a) T = 150

(b) T = 700

Figura 6: Solución aproximada tipo onda viajera del modelo (11) para diferentes momentos T, Δt = 1/24, h = 0.1, con coeficientes D = 1, r = 1, K = 12, b = 1, e = 0.125, m = 0.5, d = 0.39 y condición inicial u(x, 0) = 0.02 + exp(−(x − 100)^{2})/5 y v(x, 0) = 0.02 + exp(−(x − 100)^{2})/360.

En la Figura 5 muestra un frente de onda, que desde E_{2} converge al equilibrio espacialmente homogéneo E_{∗}. Por otro lado, en la Figura 6 se observa la presencia de órbitas periódicas cada vez más notorias cuando el momento T se hace más grande.

(a) T = 150

(b) T = 700

Figura 7: Solución aproximada tipo onda viajera del modelo (11) para diferentes momentos T, Δt = 1/24, h = 0.1, con coeficientes D = 1, r = 1, K = 12, b = 1, e = 0.125, m = 0.5, d = 0.39 y condición inicial u(x, 0) = 0.02 + exp(−(x − 100)^{2})/5 y v(x, 0) = 2/5.

Conclusiones

En este trabajo se demostró que el modelo con difusión (11) posee soluciones tipo onda viajera con oscilaciones periódicas y que conectan dos puntos de equilibrio a partir del análisis cualitativo del sistema sin difusión (6). A diferencia del trabajo de Chufen Wu [11], el análisis para la existencia de solución tipo onda viajera en el modelo (11) se planteó sin reducción de parámetros del sistema sin difusión (6) y se contrastaron los resultados obtenidos de manera numérica.

Por otro lado, se realizó un análisis cualitativo al sistema sin difusión (6) y se demostró la existencia de un único ciclo límite globalmente asintóticamente estable. De igual forma se construyó un diagrama de bifurcación para analizar la existencia y condiciones de estabilidad en sus puntos de equilibrio.

La mayoría de los resultados expuestos en este trabajo fueron inspirados en los trabajos de Chufen Wu [11] y Dunbar [1] pero con algunos resultados y demostraciones más simples, los cuales se basan básicamente en usar el Teorema de la bifurcación de Hopf, para evidenciar soluciones oscilatorias, y el uso de la función de Lyapunov y el principio de LaSalle, para demostrar que dichas soluciones conectan dos puntos de equilibrios dados por el sistema sin difusión.

Agradecimientos

Los autores agradecen tanto a los revisores como a los editores por su esfuerzos para la publicación del presente artículo.

Financiación

Los autores agradecen a la CAPES, a la Universidad Carlos III de Madrid y al CSIC por las becas de posgrado asignadas.

Referencias

[1] S. Dunbar, Traveling waves in diffusive predator-prey equations: Periodic orbits and point-to-periodic heteroclinic orbits, SIAM Journal on Applied Mathematics 46(1986), no. 6, 1057–1078. Doi: 10.1137/0146063

[2] M.R. Garvie, Finite-difference schemes for reaction-diffusion equations modeling predator-prey interactions in Matlab, Bulletin of Mathematical Biology 69(2007), no. 3, 931–956. Doi: 10.1007/s11538-006-9062-3 [3] G. Griffiths, W. Schiesser, Traveling Wave Analysis of Partial Differential

Equations: Numerical and Analytical Methods with Matlab and Maple, Academic Press, New York, 2010. Doi: 10.1016/C2009-0-64536-0

[4] P. Hartman, Ordinary Differential Equations, 2nd ed., Birkhäuser, Boston, 1982. Doi: 10.1137/1.9780898719222

[5] Y. Kuang, H.I. Freedman, Uniqueness of limit cycles in Gause-type models of predator-prey systems, Mathematical Biosciences, 88(1988), no. 1, 67–84. Doi: 10.1016/0025-5564(88)90049-1

[6] Y. Kuznetsov, Elements of Applied Bifurcation Theory, Applied Mathematical Sciences 112, Springer, New York, 1995. Doi: 10.1007/978-1-4757-3978-7

[7] J.P. LaSalle, Stability theory for ordinary differential equations, Journal of Differential Equations 4(1968), no. 1, 57–65. Doi: 10.1016/0022-0396(68)90048-X

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

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

[10] J. Sotomayor, Lições de equações diferenciais ordinárias, Projeto Euclides, Instituto de Matemática Pura e Aplicada, Rio de Janeiro, 11(1979).

[11] C. Wu, Y. Yang, P. Weng, Traveling waves in a diffusive predatorprey system of Holling type: Point-to-point and point-to-periodic heteroclinic orbits, Chaos, Solitons & Fractals 48(2013), 43–53. Doi: 10.1016/j.chaos.2013.01.003

[12] X. Wu, Y. Luo, Y. Hu, Traveling waves in a diffusive predator-prey model incorporating a prey refuge, Abstract and Applied Analysis, 2014, Article ID 679131, Doi: 10.1155/2014/679131