Estudio y comparación de diversos modelos de depredador-presa

Facultad de Matem´aticas Departamento de Ecuaciones Diferenciales y An´alisis Num´erico TRABAJO FIN DE GRADO Estudio y comparaci´on de diversos modelos de depredador-presa Autor: Miriam G´alvez Garc´ıa Tutor: Antonio Su´arez Fern´andez 20 de junio de 2018 Summary In this work we study different population models of ordinary differential equations of one and two species, respectively. In particular, we study in detail the prey-predator Lotka-Volterra model. We compare the behavior of the species in the Lotka-Volterra model and in different variants of that system. An outline of this work is as follows: in Chapter 2 we describe the different population models and show its corresponding behavior. In Chapter 3, we study theoretically the original model of Lotka-Volterra, where its stationary points, trajectories and the local stability of such model are analyzed. In Chapter 4, we study the Lotka-Volterra logistic model. We analyze the existence and uniqueness of solutions for system, stationary points and their local and global stability. Finally, in Chapter 5 we make numerical simulations related to the previous chapters and show the difference between the classic Lotka-Volterra model and some of its variants. 3 4 ´ Indice general 1. Introducci´on 7 2. Descripci´on de distintos modelos de depredador-presa 11 2.1. Modelos poblacionales para una sola especie . . . . . . . . . . . . . . . . . . . . 11 2.1.1. El modelo de crecimiento de poblaci´on de Malthus . . . . . . . . . . . . . 11 2.1.2. El modelo de Verhulst o ecuaci´on log´ıstica . . . . . . . . . . . . . . . . . 12 2.1.3. El efecto Allee ................................. 13 2.2. Modelos poblacionales para dos especies . . . . . . . . . . . . . . . . . . . . . . 14 2.2.1. Modelo original de Lotka-Volterra . . . . . . . . . . . . . . . . . . . . . . 14 2.2.2. Modelo log´ıstico de Lotka-Volterra . . . . . . . . . . . . . . . . . . . . . 15 2.2.3. Respuesta funcional de los depredadores . . . . . . . . . . . . . . . . . . 16 3. Estudio te´orico del modelo original de Lotka-Volterra 21 3.1. Ecuaci´on caracter´ıstica y estabilidad local . . . . . . . . . . . . . . . . . . . . . 21 3.2. Trayectorias...................................... 27 3.3. Estabilidad local del modelo original de Lotka-Volterra . . . . . . . . . . . . . . 29 4. Estudio te´orico del modelo log´ıstico de Lotka-Volterra 31 4.1. Ecuaci´onlog´ıstica................................... 31 4.2. Existencia y unicidad de (2.5) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.3. Puntos estacionarios y estabilidad local . . . . . . . . . . . . . . . . . . . . . . . 35 4.4. Estabilidadglobal................................... 41 5. Simulaciones num´ericas y comparaci´on del modelo de presa-depredador 51 5 ´ Indice general 6 Cap´ıtulo 1 Introducci´on El principal objetivo de este trabajo es estudiar en detalle el modelo basado en ecuaciones diferenciales de presa-depredador de Lotka-Volterra. A partir de ´el, se introducen distintas variantes de dicho modelo. El Cap´ıtulo 2 est´a dedicado a describir los distintos modelos poblacionales, tanto para una sola especie como para dos especies. A continuaci´on se muestran las ecuaciones de los modelos descritos en este cap´ıtulo. Empezamos con el modelo de crecimiento de la poblaci´on de Malthus expresado mediante la ecuaci´on: dx dt =rx. (1.1) M´as tarde, P.Verhulst modifica el modelo de Malthus donde la din´amica poblacional responde a la siguiente ecuaci´on diferencial: dx dt =rx(1 −x K) (1.2) Dicha ecuaci´on es denominada como modelo de Verhulst o ecuaci´on log´ıstica. A continuaci´on introducimos el modelo matem´atico que representa el efecto Allee dado por la ecuaci´on diferencial no lineal: dx dt =rx1−x Kx−m(1.3) Por ´ultimo introducimos el modelo original de Lotka-Volterra para sistemas de poblaciones depredador-presa, que es el que cobra mayor importancia en este trabajo: 7 Introducci´on dx dt =αx −axy =x(α−ay) dy dt =−βy +bxy =y(−β+bx) (1.4) Introduciendo un t´ermino log´ıstico en el sistema anterior obtenemos el siguiente sistema denominado modelo log´ıstico de Lotka-Volterra dx dt =x(α−By −Ax) dy dt =y(β+Cx −Dy) (1.5) Al final de este cap´ıtulo mostramos la respuesta funcional de los depredadores que representa la cantidad de presas que puede consumir un depredador en una unidad de tiempo. Esta respuesta funcional se puede escribir de forma general como: dx dt =x(α−Ax)−Bh(x)y dy dt =y(β−Dy) + Ch(x)y, (1.6) y consideramos 4 tipos. La funci´on utilizada en Holling tipo I es la siguiente: h(x) =    δx si 0 ≤x < c δc si c≥x La funci´on utilizada en la respuesta funcional en Holling tipo II es: h(x) = qx x+a La funci´on utilizada en Holling tipo III, es: h(x) = qx2 x2+a2 8 Introducci´on Las funciones utilizadas en Holling tipo IV, vienen dadas por: a) h(x) = qx x2+a2 b) h(x) = qx2 x2−bx +a c) h(x) = qxm xn+an, con 1 ≤m < n Las gr´aficas que representan dichas funciones est´an introducidas en el cap´ıtulo. En el Cap´ıtulo 3, se hace un estudio te´orico del modelo original de Lotka-Volterra, donde se analiza sus puntos estacionarios, trayectorias y la estabilidad local de dicho modelo. En el Cap´ıtulo 4, se estudia de manera te´orica el modelo log´ıstico de Lotka-Volterra. Se resuelve en detalle la ecuaci´on log´ıstica introducida en el primer cap´ıtulo y posteriormente vemos la existencia y unicidad del sistema, los puntos estacionarios y la estabilidad local y global. Finalmente, en el Cap´ıtulo 5 hacemos simulaciones num´ericas con Matlab relacionadas con el cap´ıtulo anterior y mostramos la diferencia que hay entre el cl´asico modelo de Lotka-Volterra y alguna de sus variantes. 9 Cap´ıtulo 2. Descripci´on de distintos modelos de depredador-presa 2.2.3. Respuesta funcional de los depredadores La respuesta funcional, o tasa de consumo, expresa la acci´on de los depredadores en la tasa de crecimiento de la poblaci´on de presas, y representa la cantidad de presas que puede consumir un depredador en una unidad de tiempo. Esta respuesta funcional fue clasificada en tres tipos por C.S Holling en 1959. En el a˜no 1984, Taylor describe un cuarto tipo llamado Holling tipo IV o respuesta funcional no mon´otona. A continuaci´on se presentan estos cuatro tipos de respuesta funcional de forma que el modelo de Lotka-Volterra log´ıstico queda de la siguiente forma: dx dt =x(α−Ax)−Bh(x)y dy dt =y(β−Dy) + Ch(x)y (2.6) siendo h(x) la funci´on definida en cada uno de los siguientes apartados. Respuesta funcional Holling tipo I La respuesta funcional se obtiene suponiendo que el cambio en la densidad de poblaci´on de depredadores es proporcional a la densidad de la poblaci´on de presas disponibles, es decir, existe un aumento lineal de la tasa de ataque de la poblaci´on de depredadores respecto a la densidad de poblaci´on de presas, hasta llegar a un punto a partir del cual la tasa m´axima de ataque se hace constante. La funci´on utilizada en la respuesta funcional Holling tipo I, viene dada por: h(x) =    δx si 0 ≤x < c δc si c≥x donde δes una constante, crepresenta la saturaci´on y x(t) es la densidad de la poblaci´on de presas en el tiempo t≥0. 16 Cap´ıtulo 2. Descripci´on de distintos modelos de depredador-presa Figura 2.3: La respuesta funcional Holling tipo I Respuesta funcional Holling tipo II En esta respuesta funcional aparece un segundo par´ametro denominado tiempo de manipulaci´on. El tiempo de manipulaci´on incluye el tiempo de perseguir, dominar y comer a cada presa, como tambi´en el tiempo que el depredador toma para prepararse e ir a buscar m´as presas. En este caso, el tiempo de b´usqueda disminuye y el tipo de respuesta resulta en un incremento desacelerado del consumo a medida que aumenta la poblaci´on de presas consumidas hasta que la funci´on se estabiliza puesto que el depredador se satura. La funci´on utilizada en esta respuesta funcional, est´a dada por: h(x) = qx x+a Figura 2.4: La respuesta funcional Holling tipo II donde x(t) representa la densidad de la poblaci´on de presas en el tiempo t≥0, qes la tasa 17 Cap´ıtulo 2. Descripci´on de distintos modelos de depredador-presa m´axima de consumo per c´apita y aes la tasa de saturaci´on media, es decir, la cantidad de presas en el que la tasa de depredaci´on alcanza la mitad de su valor m´aximo. Se puede observar que si x→ ∞ entonces h(x)→q, es decir, qes el valor m´aximo de la tasa de depredaci´on de h(x). Respuesta funcional Holling tipo III Esta respuesta expresa una curva de crecimiento sigmoidea, donde la cantidad de presas atacadas por depredador aumenta lentamente a bajos niveles de poblaci´on de presas, luego disminuye a altas densidades hasta que se estabiliza. Aqu´ı, la cantidad de presas consumidas se acelera con el aumento de la densidad de presas hasta que el tiempo de manipulaci´on comienza a limitar su consumo. La funci´on utilizada en Holling tipo III, es h(x) = qx2 x2+a2 Figura 2.5: La respuesta funcional Holling tipo III donde x(t) representa la densidad de la poblaci´on de presas en el tiempo t≥0, qes la tasa m´axima de consumo per c´apita y arepresenta la tasa de saturaci´on media, es decir, la cantidad de presas en el que la tasa de depredaci´on alcanza la mitad de su valor m´aximo. 18 Cap´ıtulo 2. Descripci´on de distintos modelos de depredador-presa Respuesta funcional Holling tipo IV Esta respuesta es utilizada para modelar comportamientos antidepredatorios (APB), llamado formaci´on de grupos de defensa y el fen´omeno de agregaci´on, el cual es empleado por las presas para evitar la depredaci´on. La primera parte de la curva (ver Figura 2.6), puede corresponder a cualquiera de los tipos previos de respuestas funcionales, por lo tanto, el foco de atenci´on surge por la disminuci´on de la cantidad de presas consumidas que algunas veces se observa en altas densidades de presas. Existen tres posibles causas de esta reducci´on. Primero, los depredadores se concentran menos en una presa individual cuando existen muchas presas para escoger. Segundo, algunas presas pueden cooperar, compartiendo la carga de vigilar la presencia de depredadores e informar a otros de la presencia. Tercero, algunas presas podr´ıan intimidar e incluso defenderse contra el ataque de los depredadores. Las funciones utilizadas en Holling tipo IV, vienen dadas por: a) h(x) = qx x2+a2 b) h(x) = qx2 x2−bx +a c) h(x) = qxm xn+an, con 1 ≤m < n Figura 2.6: La respuesta funcional Holling tipo IV en el caso a) 19 Cap´ıtulo 2. Descripci´on de distintos modelos de depredador-presa 20 Cap´ıtulo 3 Estudio te´orico del modelo original de Lotka-Volterra 3.1. Ecuaci´on caracter´ıstica y estabilidad local En esta secci´on introduciremos algunas definiciones de estabilidad, y adem´as, abordaremos el problema de determinar por medio de la ecuaci´on caracter´ıstica la estabilidad local de un sistema. Definici´on 3.1. Consideremos un sistema aut´onomo de ecuaciones diferenciales de primer orden dx dt =f1(x, y) dy dt =f2(x, y) (3.1) 1. Un punto de equilibrio ocurre cuando x(t) y y(t) son constantes, de forma que dx dt = 0 y dy dt = 0. En un punto de equilibrio decimos que el sistema est´a en reposo. Los puntos de equilibrio se conocen tambi´en como puntos estacionarios, o puntos singulares, o puntos cr´ıticos opuntos de reposo. 2. Un punto de equilibro (x∗, y∗) es localmente asint´oticamente estable si la respuesta a una peque˜na perturbaci´on tiende a cero cuando el tiempo tiende a infinito, o m´as precisamente, si l´ım t→∞ x(t) = x∗y l´ım y→∞ y(t) = y∗ 21 Cap´ıtulo 3. Estudio te´orico del modelo original de Lotka-Volterra siempre y cuando (x(t), y(t)) comience en valores iniciales P0≡(x(0), y(0)) ≡(x0, y0) suficientemente cercanos al punto de equilibrio (x∗, y∗). Los puntos de equilibrio localmente asint´oticamente estables se conocen tambi´en como pozos oatractores. Al conjunto de todos los valores iniciales (x0, y0) cuyas trayectorias convergen a un atractor se le llama el dominio de atracci´on. 3. Un punto de equilibrio (x∗, y∗) es estable, pero no localmente asint´oticamente estable, si la respuesta a una peque˜na perturbaci´on es peque˜na pero no tiende a cero cuando el tiempo tiende a infinito. En este caso, el punto de equilibrio se llama centro. En este caso la soluci´on regresar´a a la condici´on inicial en forma peri´odica, de tal manera que su trayectoria es una curva cerrada alrededor del punto de equilibrio, por lo que es llamada ciclo. 4. Un punto de equilibrio es inestable si condiciones iniciales (cerca del punto de equilibrio) producen soluciones que se alejan del punto de equilibrio al transcurrir el tiempo. Cerca de un punto de equilibrio inestable, peque˜nas perturbaciones a condiciones iniciales tienden a afectar dr´asticamente las soluciones. Los puntos de equilibrio inestables se conocen tambi´en como fuentes orepulsores. Observaci´on 3.1. El sistema (2.4) tiene dos puntos de equilibrio, a saber, 0≡(0,0) y P∗≡β b,α a. M´as adelante, en la Secci´on 3.3, deducimos que el punto de equilibrio 0del sistema (2.4) es inestable y que el punto de equilibro P∗es estable (pero no localmente asint´oticamente) y es un centro. Ahora, consideremos nuevamente el modelo aut´onomo (3.1) y supongamos que tiene un punto de equilibrio en (¯x, ¯y), de forma que f1(¯x, ¯y) = 0 y f2(¯x, ¯y) = 0 (3.2) Para analizar la estabilidad lineal del modelo (3.1) en (¯x, ¯y) sustituimos x(t) = ¯x+u1(t) e y(t) = ¯y+u2(t) en (3.1), donde 0 < ui(t)<< 1 (con i= 1,2), considerando as´ı una peque˜na perturbaci´on. Usando la expansi´on en serie de Taylor obtenemos 22 Cap´ıtulo 3. Estudio te´orico del modelo original de Lotka-Volterra du1 dt =f1(¯x+u1,¯y+u2) =f1(¯x, ¯y) + u1 ∂ ∂x +u2 ∂ ∂yf1(x, y)|(¯x,¯y) +1 2!u1 ∂ ∂x +u2 ∂ ∂y2 f1(x, y)|(¯x,¯y)+··· y du2 dt =f2(¯x+u1,¯y+u2) =f2(¯x, ¯y) + u1 ∂ ∂x +u2 ∂ ∂yf2(x, y)|(¯x,¯y) +1 2!u1 ∂ ∂x +u2 ∂ ∂y2 f2(x, y)|(¯x,¯y)+··· Usando (3.2) y linealizando (ya que u1(t) y u2(t) son muy peque˜nas), obtenemos du1 dt =u1 ∂f1 ∂¯x+u2 ∂f1 ∂¯yydu2 dt =u1 ∂f2 ∂¯x+u2 ∂f2 ∂¯y(3.3) donde ∂f1 ∂¯x,∂f1 ∂¯y,∂f2 ∂¯x,∂f2 ∂¯ydenota los valores de ∂f1 ∂x ,∂f1 ∂y ,∂f2 ∂x ,∂f2 ∂y en el punto de equilibro (¯x, ¯y). Consideraremos la soluci´on en la forma u1=A1eλt yu2=A2eλt.(3.4) Observamos que u1yu2no se anulan simult´aneamente ya que hemos realmente generado una perturbaci´on. Sustituyendo u1yu2de (3.4) en (3.3) y simplificando llegamos a λ−∂f1 ∂¯xA1−∂f1 ∂¯yA2= 0 −∂f2 ∂¯xA1+λ−∂f2 ∂¯yA2= 0, 23 Cap´ıtulo 3. Estudio te´orico del modelo original de Lotka-Volterra (a) Nodo estable (b) Nodo inestable Figura 3.1: Cuando las ra´ıces de la ecuaci´on caracter´ıstica son reales y tienen el mismo signo (caso 1), entonces el punto de equilibrio es estable o inestable, ya sea que el signo de las ra´ıces sea negativo o positivo respectivamente. o bien, en forma matricial,      λ−∂f1 ∂¯x−∂f1 ∂¯y −∂f2 ∂¯xλ−∂f2 ∂¯y          A1 A2     =    0 0     . Ya que u1yu2no se anulan simult´aneamente, lo mismo ocurre con A1yA2, de forma que  λ−∂f1 ∂¯x−∂f1 ∂¯y −∂f2 ∂¯xλ−∂f2 ∂¯y  = 0 y por lo tanto λ2−∂f1 ∂¯x+∂f2 ∂¯yλ+∂f1 ∂¯x ∂f2 ∂¯y−∂f1 ∂¯y ∂f2 ∂¯x= 0.(3.5) Esta ecuaci´on se conoce como la ecuaci´on caracter´ıstica. Es una ecuaci´on cuadr´atica en λy por lo tanto tiene dos soluciones λ1yλ2. La soluci´on general la podemos entonces escribir como u1(t) = a1eλ1t+a2eλ2t u2(t) = a3eλ1t+a4eλ2t (3.6) 24 Cap´ıtulo 3. Estudio te´orico del modelo original de Lotka-Volterra (a) Punto silla (b) Centro Figura 3.2: (a) Si las ra´ıces de la ecuaci´on caracter´ıstica son reales y tienen signos contrarios (caso 2), entonces el punto de equilibrio es un punto silla. (b) Si las ra´ıces de la ecuaci´on caracter´ıstica son imaginarias (caso 4), entonces el punto de equilibrio es un centro. donde ai∈R(con i= 1,...,4). Usando x(t) = ¯x+u1(t) y y(t) = ¯y+u2(t) junto con (3.6), obtenemos los siguientes casos: Caso 1: NODO. Las ra´ıces λ1yλ2de la ecuaci´on caracter´ıstica (3.5) son reales y tienen el mismo signo. En este caso el punto de equilibrio (¯x, ¯y) se clasifica como nodo (o punto nodal). Hay dos subcasos: Estable. Si λ1, λ2<0, entonces el punto de equilibrio (¯x, ¯y) es un nodo estable. En este caso u1→0 y u2→0 cuando t→ ∞ y por lo tanto x(t)→¯xy y(t)→¯ycuando t→ ∞ (ver Figura 3.1.a). Inestable. Si λ1, λ2>0, entonces el punto de equilibrio (¯x, ¯y) es un nodo inestable. En este caso u1(t)→ ∞ yu2(t)→ ∞ cuando t→ ∞ (ver Figura 3.1.b). Caso 2: PUNTO SILLA. Las ra´ıces λ1yλ2de la ecuaci´on caracter´ıstica (3.5) son reales y de signos diferentes. En este caso el punto de equilibrio (¯x, ¯y) se clasifica como punto silla. Ahora u1(t)→ ∞ yu2(t)→ ∞ cuanto t→ ∞ y por lo tanto x(t)→ ∞ yy(t)→ ∞ cuando t→ ∞. De forma que (¯x, ¯y) es inestable (ver Figura 3.2.a). Caso 3: FOCO. Las ra´ıces λ1yλ2de la ecuaci´on caracter´ıstica (3.5) son n´umeros complejos conjugados con parte real no nula. En este caso el punto de equilibrio (¯x, ¯y) se clasifica como foco (o punto focal). Cerca del punto focal a las trayectorias que corresponden a las soluciones se comportan como espirales. 25 Cap´ıtulo 4. Estudio te´orico del modelo log´ıstico de Lotka-Volterra cambio de variable, resolvamos la ecuaci´on (4.2). Supongamos que a6= 0, eat(z0+az) = eatb⇐⇒ (eatz)0=beat ⇐⇒ eatz−z0=bZt 0 easds =b a(eat −1) ⇐⇒ eatz(t) = b a(eat −1) + z0⇐⇒ z(t)=(b a(eat −1) + z0)e−at. Deshacemos el cambio de variables y concluimos que la soluci´on de la ecuaci´on (4.1) es: x(t) = eat b a(eat −1) + 1 x0 .(4.3) Como a∈R, puede ocurrir que a= 0, por tanto la ecuaci´on (4.2) ser´ıa de la forma siguiente:      z0=−x(a−bx) x2=b z(0) = z0=1 x0 (4.4) De nuevo resolvemos la ecuaci´on (4.4) y obtenemos: a0=b⇒z(t) = bt +z0 Volvemos a deshacer el cambio de variables y la soluci´on de la ecuaci´on (4.1) cuando a= 0 es: x(t) = 1 bt +1 x0 (4.5) Ya conocidas las soluciones de la ecuaci´on log´ıstica (4.1), vamos a ver el comportamiento de ambas dependiendo de los valores de a, b. Proposici´on 4.1. Sea x(t)la ´unica soluci´on de (4.1) con b > 0. 1. Si a > 0entonces l´ım t→+∞ x(t) = a b>0. 2. Si a≤0entonces l´ım t→+∞ x(t)=0. 32 Cap´ıtulo 4. Estudio te´orico del modelo log´ıstico de Lotka-Volterra Demostraci´on: 1. Vamos a estudiar en primer lugar el caso en el que a, b > 0 y con dato inicial x0>0. Veamos que la soluci´on de la ecuaci´on log´ıstica (4.1) existe y es ´unica. Como a6= 0, estudiaremos la soluci´on (4.3). Para demostrar que la soluci´on existe s´olo tenemos que ver que el denominador no se anula. Sea el denominador: b a(eat −1) + 1 x0 .(4.6) Como t≥0 entonces eat ≥1⇒eat −1≥0. Por tanto, el denominador nunca se anula. Con esto hemos probado que existe una ´unica soluci´on (4.3) para todo t≥0 cuando a, b > 0. Ya visto que existe una soluci´on, vamos a calcular su l´ımite. l´ım t→+∞ x(t) = l´ım t→+∞ 1 b a1−1 eat +1 x0eat =a b>0. Observamos que cuando t→+∞la soluci´on tiende a un n´umero finito positivo. 2. Sea ahora el caso en el que a < 0, b > 0 y condici´on inicial positiva. En primer lugar probaremos que existe una ´unica soluci´on para todo t≤0. Como hemos hecho anteriormente, queremos ver que el denominador de (4.3) no se anula cuando t≤0, siendo dicho denominador (4.6). Como a < 0, obtenemos que eat ≤1 para todo t≥0, por lo que eat −1≤0. Pero a b<0, por tanto a b(eat −1) >0. De nuevo el denominador de la soluci´on (4.3) no se anula. Ya demostrado que la soluci´on (4.3) existe, vamos a calcular su l´ımite. l´ım t→+∞ x(t) = l´ım t→+∞ eat b a(eat −1) + 1 x0 =0 −b a+1 x0 = 0. 33 Cap´ıtulo 4. Estudio te´orico del modelo log´ıstico de Lotka-Volterra 4.2. Existencia y unicidad de (2.5) En esta secci´on veremos la existencia y unicidad del sistema:          x0=x(α−By −Ax)≡f1(x, y) y0=y(β+Cx −Dy)≡f2(x, y) (x(0), y(0)) = (x0, y0)≥(0,0) (4.7) con α,β∈RyA, B, C, D > 0. Nota 4.2.1. La condici´on inicial se puede comportar de la siguiente forma: Si x0≡y0≡0, por unicidad de la soluci´on, entonces la ´unica soluci´on de (4.7) es (x, y)≡(0,0). Si x0≡0, y0>0 entonces la ´unica soluci´on de (4.7) es x=0eysoluci´on de una ecuaci´on log´ıstica, estudiada en la Secci´on 4.1. Si y0≡0, x0>0 ocurre lo mismo que en el apartado anterior. Si x0>0, y0>0 entonces las soluciones del sistema (4.7) son positivas. En efecto, las ecuaciones del sistema (4.7) se pueden escribir de la forma x0=xa(t)x0>0, donde a(t) es una funci´on que viene dada por (4.7). Luego la soluci´on es de la forma x=eA(t)x0>0, donde A(t) es la primitiva de a(t). An´alogo para y. Por tanto observamos que si (x0, y0)>(0,0), no existe t0>0 tal que x(t0) = 0 e y(t0) = 0. Por tanto las soluciones (x, y) correspondientes a (4.7) son positivas. Proposici´on 4.2. Sea el sistema de Lotka-Volterra (4.7). Entonces existe una ´unica soluci´on positiva global en tiempo. Demostraci´on: Vamos a probar que la soluci´on existe y es ´unica. En primer lugar, sabemos por el Teorema de Picard, que la soluci´on existe localmente en tiempo, ya que f1yf2son regulares. 34 Cap´ıtulo 4. Estudio te´orico del modelo log´ıstico de Lotka-Volterra El objetivo de la prueba es acotar las ecuaciones del sistema (4.7), por separadas, ya que conjuntamente, como un sistema, no cumple las condiciones del Teorema de Comparaci´on. Las ecuaciones que usaremos para acotar las de nuestro sistema son las estudiadas en la Secci´on 4.1, es decir, las ecuaciones log´ısticas. Empecemos con la primera ecuaci´on: x0=x(α−Ax −By)≤x(α−Ax)⇒x(t)≤L(t) donde L(t) es la soluci´on de la ecuaci´on log´ıstica:    x0=x(α−Ax) x(0) = x0 ⇒x(t) existe ∀t > 0. Para la segunda ecuaci´on tenemos lo siguiente: y0=y(β−Dy +Cx)≤y(β+CL(t)−Dy)≤y(β+CL0−Dy) ⇒y(t)≤M(t) donde Mes soluci´on de M0=M(β+CL0−DM), donde L0es la cota de L(t). 4.3. Puntos estacionarios y estabilidad local A continuaci´on, estudiaremos los puntos estacionarios y la estabilidad local de cada uno de ellos. Lema 4.1. Sea el sistema de Lotka-Volterra (2.5). Los puntos de equilibrio son los siguientes: P1= (0,0), P2=α A,0, P3=0,β D P4=Dα −Bβ AD +BC ,Aβ +Cα AD +BC  35 Cap´ıtulo 4. Estudio te´orico del modelo log´ıstico de Lotka-Volterra Se tiene: a) Si α, β < 0entonces P1es localmente estable. b) Si α > 0yβ < −αC Aentonces P2es localmente estable. c) Si β > 0yβ > αD Bentonces P3es localmente estable. d) Si α > βB Dyβ > −αC Aentonces P4es localmente estable. Demostraci´on: Vamos a calcular los puntos de equilibrio del sistema (2.5). Para dicho c´alculo igualaremos las ecuaciones del sistema no lineal a cero.    0 = f1(x, y) 0 = f2(x, y) Para obtener los puntos de equilibrio, vamos a diferenciar varios casos. Si x= 0. Sustituimos en el sistema. y(β−Dy) = 0 =⇒y= 0 ´o y=β D Si x= 0 e y= 0 entonces obtenemos uno de los puntos de equilibrio P1= (0,0). Por otro lado si x= 0 e y=β Del punto de equilibrio es P3=0,β D Si y= 0. La ecuaci´on que conseguimos es: x(α−Ax) = 0 =⇒x= 0 ´o x=α A Con estos valores conseguimos de nuevo el punto P1y volviendo a sustituir dicho valores en el sistema obtenemos el punto P2=α A,0. Por ´ultimo supongamos que x, y 6= 0.En este caso, el sistema ser´ıa:    α−Ax −By = 0 β−Dy +Cx = 0 De la primera ecuaci´on obtenemos: x=α−By A 36 Cap´ıtulo 4. Estudio te´orico del modelo log´ıstico de Lotka-Volterra Sustituyendo el resultado en la segunda ecuaci´on: β D+C D(α−By A) = y β D+Cα DA −CBy DA =y β D+Cα DA −CBy DA =y β D+Cα DA =y1+CB DA y=βA +Cα AD +CB Volvemos con x=α−By Ay sustituimos el valor de yobtenido: α A−β AβA +Cα AD +CB =x α A−βBA A(AD +CB)−BCα A(AD +CB)=x α(AD +CB)−βBA −BCα A(AD +CB)=x αAD +αCB −βBA −BCα A(AD +CB)=x x=αD −βB AD +CB Concluimos que el ´ultimo punto de equilibrio del sistema es: P4=αD −βB AD +CB ,βA +Cα AD +CB  Ya conocidos los puntos de equilibrio, vamos a demostrar que estos cuatro puntos son localmente estables bajo ciertas condiciones en los par´ametros del problema. Primeramente, vamos a calcular la matriz jacobiana correspondiente al sistema (2.5): J=   α−By −2Ax −Bx Cy β +Cx −2Dy    37 Cap´ıtulo 4. Estudio te´orico del modelo log´ıstico de Lotka-Volterra A partir de esta matriz, tenemos que calcular los autovalores de cada punto, y si la parte real de cada autovalor es negativa, entonces el punto de equilibrio ser´a localmente estable. Para el punto P1= (0,0), sustituyendo dicho valor en la matriz jacobiana, obtenemos: a) A1=  α0 0β  Como la matriz resultante s´olo tiene diagonal principal, los autovalores correspondientes al punto P1son los elementos de dicha diagonal, es decir, αyβ. Por tanto, si α,β < 0 el punto P1es localmente estable. b) P2= (α A,0). Al igual que con el punto P1, sustituimos en la matriz jacobiana el valor de dicho punto. El resultado es: A2=    −α−Bα A 0β+Cα A     Es un caso similar al del punto P1. Podemos observar que la matriz tiene una diagonal principal y la diagonal inferior est´a formada por cero, luego los autovalores asociados al punto P2, son de nuevo los elementos de la diagonal. γ1=−α γ2=β+Cα A Para que el punto P2sea localmente estable, tiene que ocurrir que los autovalores sean negativos. Entonces: −α < 0 =⇒α > 0 β+Cα A<0 =⇒β < −Cα A Con ´esto hemos demostrado lo que quer´ıamos. 38 Cap´ıtulo 4. Estudio te´orico del modelo log´ıstico de Lotka-Volterra c) Este caso es similar al anterior. De nuevo sustituimos en la matriz jacobiana el valor del punto P3= (0,β D) y la matriz resultante es: A3=    α−Bβ D0 Cβ D−β     Luego, vemos que la matriz es triangular, por ello los autovalores son los elementos de la diagonal principal. γ1=−β γ2=α−Bβ D Para que sea localmente estable tiene que ocurrir, que sean negativos los autovalores, luego si −β < 0⇒β > 0. Y con el otro autovalor obtenemos que α−Bβ D<0 =⇒α < Bβ D Con esto llegamos a lo que quer´ıamos demostrar. d) Por ´ultimo veamos el comportamiento de los autovalores del punto P4.Observemos que como α > βB Dyβ > −αC Alas componentes de P4son positivas. Sea la matriz A4que consiste en sustituir el punto P4en la matriz jacobiana y denotamos el punto P4= (x∗, y∗), entonces: A4=     −AαD −βB AD +CB −BαD −βB AD +CB  CβA +αC AD +CB −DβA +αC AD +CB          −Ax∗−Bx∗ Cy∗−Dy∗   = diag(x∗, y∗)   −A−B C−D    Para calcular los autovalores, no podemos hacer como en los apartados anteriores, ya que la matriz obtenida por el jacobiano no es triangular. Por lo tanto para obtenerlos tenemos 39 Cap´ıtulo 4. Estudio te´orico del modelo log´ıstico de Lotka-Volterra que calcular el determinante siguiente: 0 = |γId −A4|=   γ+Ax∗Bx∗ −Cy∗γ+Dy∗    Por tanto la ecuaci´on resultante de calcular dicho determinante es: (γ+Ax∗)(γ+Dy∗) + BCx∗y∗=γ2+γ(Ax∗+Dy∗)+(AD +BC)x∗y∗= 0 ⇒γ=−(Ax∗+Dy∗)±p(Ax∗+Dy∗)2−4(AD +BC)x∗y∗ 2 (4.8) Podemos observar los distintos casos de la ecuaci´on (4.8). Si Disc ≥0⇒γ < 0 Si Disc < 0⇒Re(γ)<0 denotando Disc como el discriminante de la ecuaci´on (4.8) y Re como la parte real de γ. Por tanto los autovalores resultantes de la ecuaci´on (4.8) tienen la parte real negativa, que es lo que quer´ıamos probar. Nota 4.3.1. Si las condiciones anteriores no se dan, tenemos que los puntos son inestables. Es decir, a) Si α > 0 o β > 0⇒P1es inestable. b) Si α > 0 y β > −αC A⇒P2es inestable. c) Si β > 0 y β < α D B⇒P3es inestable. 40 Cap´ıtulo 4. Estudio te´orico del modelo log´ıstico de Lotka-Volterra Figura 4.1 A continuaci´on hemos dividido el plano α−βen 4 regiones: a) En A1≡ {(α, β) : α < 0 y β < 0}el punto P1es localmente estable. b) En A2≡ {(α, β) : α > 0 y β < −αC A}el punto P2es localmente estable. c) En A3≡ {(α, β) : β > 0 y β > αD B}el punto P3es localmente estable. d) En A4≡ {(α, β) : α > β B Dyβ > −αC A}el punto P4es localmente inestable. 4.4. Estabilidad global En la secci´on anterior hemos probado como los puntos estacionarios del sistema Lotka-Volterra (2.5) son localmente estables bajo ciertas condiciones en los par´ametros. En este apartado queremos demostrar que tambi´en son globalmente estable. Para que un punto de equilibrio de un sistema de ecuaciones diferenciales ordinarias sea globalmente estable, una de las cosas que tiene que ocurrir es que dicho sistema no posea ´orbitas peri´odicas. Para ello, vamos a exponer, en primer lugar, el teorema de Bendixson-Dulac. Teorema 4.1. (Bendixson - Dulac) Sea U⊆R2un conjunto abierto, simplemente conexo y f:U−→ R2una aplicaci´on continua y diferenciable. Sea w:U−→ Rdiferenciable tal que ∇·[w(x)f(x)] 6= 0. 41 Cap´ıtulo 4. Estudio te´orico del modelo log´ıstico de Lotka-Volterra Para demostrar este apartado, utilizaremos la Proposici´on 4.2. En dicha proposici´on llegamos, por un lado a que x0≤x(α−Ax)⇒x(t)≤L(t) como α < 0⇒L(t)→0 y por tanto x(t)→0. Por otro lado, como x(t)→0, dado β < 0 existe ε > 0 tal que β+Cε < 0.Dado que ε > 0, tenemos que x(t)< ε ∀t≥t0. Entonces y0=y(β+Cx −Dy)≤y(β+Cε −Dy) como β+Cε < 0, por la Secci´on 4.1, obtenemos que las soluciones de las ecuaciones log´ısticas tienden a cero, cuando t→+∞. Por tanto, (x, y)→(0,0), cuando t→+∞. 2. Veamos ahora el caso en el que α > 0 y β < −αC A. Para ello igualaremos los segundos miembros del sistema (4.9) a cero,    h1(x, y) = α−By −Ax = 0 h2(x, y) = β+Cx −Dy = 0 Despejamos en ambas ecuaciones la inc´ognita yy obtenemos:      y=α−Ax B y=β+Cx D A continuaci´on, dibujaremos ambas rectas y veremos el comportamiento que tienen. 48 Cap´ıtulo 4. Estudio te´orico del modelo log´ıstico de Lotka-Volterra Figura 4.3 Observamos que la Figura 4.3, se divide en cuatro ´areas (relacionado con la Secci´on 3.2). Vamos a ver detalladamente cada una de ellas. Comencemos con el ´area A1. En ella h1(x, y)>0, h2(x, y)<0. Por tanto si existe una soluci´on (x∗, y∗) en el recinto A1, dicha soluci´on tiende a ir al punto α A,0. Sin embargo en el recinto A2, las dos rectas son positivas, h1(x, y), h2(x, y)>0, por ello, dado una soluci´on (x∗, y∗) en ese recinto, ´esta tiende a introducirse en A1. En el recinto A3, las dos rectas son negativas, h1(x, y), h2(x, y)<0, por ello, dado una soluci´on (x∗, y∗) en ese recinto, ´esta tiende a ir al punto α A,0. 49 Cap´ıtulo 4. Estudio te´orico del modelo log´ıstico de Lotka-Volterra Por ´ultimo, el ´area A4, las rectas tienen signos distintos, es decir, h1(x, y)<0, h2(x, y)>0, por ello, dado una soluci´on (x∗, y∗) en ese recinto, ´esta tiende a introducirse en A3. Por tanto, vemos en la Figura 4.3 que la soluci´on del sistema (4.9) tiende al punto P2, es decir, cuando t→+∞ (x, y)→α A,0 3. La estabilidad global del punto P3, se hace de forma similar. 50 Cap´ıtulo 5 Simulaciones num´ericas y comparaci´on del modelo de presa-depredador En este cap´ıtulo presentamos simulaciones num´ericas relacionadas con el Cap´ıtulo 4 y mostramos la diferencia que hay entre el cl´asico modelo de Lotka-Volterra y alguna de sus variantes. Para simular el sistema de Lotka-Volterra log´ıstico asignamos los valores {A=C=D= 1, B = 3}. Con estos valores obtenemos 4 zonas en el plano α−β(ver Figura 4.1). Figura 5.1 51 Cap´ıtulo 5. Simulaciones num´ericas y comparaci´on del modelo de presa-depredador Si (α, β)∈A1⇒(x, y)→(0,0) cuando t→+∞. Si (α, β)∈A2⇒(x, y)→(α, 0) cuando t→+∞. Si (α, β)∈A3⇒(x, y)→(0, β) cuando t→+∞. Si (α, β)∈A4⇒(x, y)→α−3β 4,β+α 4cuando t→+∞. Veamos esto con Matlab con, por ejemplo, (x0, y0) = (1,1) Tomando (α, β) = (−1,−1) vemos que efectivamente (α, β)∈A1. La siguiente simulaci´on muestra el comportamiento en este caso de la presa y el depredador. Figura 5.2 Tanto las presas como los depredadores tienden a cero, por lo que en un futuro ambas especies se extinguir´an. Tomando (α, β) = (2,−3) tenemos que (α, β)∈A2. La siguiente gr´afica muestra el comportamiento de la presa y el depredador en este caso. 52 Cap´ıtulo 5. Simulaciones num´ericas y comparaci´on del modelo de presa-depredador Figura 5.3 En este caso, las presas sobrevivir´an aunque haya presencia de depredadores. Sin embargo, los depredadores aunque hay presas, en un futuro terminar´an por extinguirse, lo que es debido a que su tasa de crecimiento es peque˜na. Si tomamos (α, β) = (2,2) efectivamente se tiene (α, β)∈A3. Se muestra el comportamiento de ambas especies para este caso en la siguiente simulaci´on. Figura 5.4 Cuando se da esta situaci´on, podemos ver en la Figura 5.4, que las presas tienden a extinguirse puesto que son deboradas por los depredadores, y aunque en alg´un momento 53 Cap´ıtulo 5. Simulaciones num´ericas y comparaci´on del modelo de presa-depredador los depredadores acaban con todas las presas, ´estos sobreviven. Tomando (α, β) = (2,1 2) tenemos que (α, β)∈A4. La siguiente simulaci´on muestra el comportamiento en este caso de la presa y el depredador. Figura 5.5 Por ´ultimo, ambas especies sobrevivir´an con el paso del tiempo, como se puede apreciar en la Figura 5.5. Veamos c´omo para los mismos valores de los par´ametros, al introducir una modificaci´on en el modelo de presa-depredador, la din´amica cambia. Por ejemplo, lo hacemos con (α, β)∈A3, luego tomamos (α, β) = (2,2) y cambiando el sistema de Lotka-Volterra log´ıstico dx dt =x(α−3y−x) dy dt =y(β+x−y) por el Holling tipo II dx dt =x(α−Ax)−Bqx x+ay dy dt =y(β−Dy) + Cqx x+ay 54 Cap´ıtulo 5. Simulaciones num´ericas y comparaci´on del modelo de presa-depredador tomando q=a= 1, obtenemos la siguiente gr´afica: Figura 5.6 donde se ve que a diferencia de la Figura 5.4, en esta, ambas especies sobreviven en el tiempo. Sin embargo en el sistema de Lotka-Volterra con los mismos valores hemos visto que las presas se extingu´ıan. Esto es debido a que en el sistema de Lotka-Volterra, las presas eran deboradas por los depredadores sin control alguno. En cambio, al comparar dicho sistema con Holling tipo II, existe un control en la depredaci´on, lo que trae consigo que las presas no sean deboradas completamente por los depredadores. Esto mismo podr´ıamos haber hecho con las otras variantes del sistema Lotka-Volterra. 55 Cap´ıtulo 5. Simulaciones num´ericas y comparaci´on del modelo de presa-depredador 56 Bibliograf´ıa [1] G. Samanta, R. G´ omez A´ ıza,Modelos din´amicos de poblaciones simples y de sistemas depredador-presa. India, M´exico, (2014). [2] F.J. Reyes Baham´ on,Sobre la din´amica de algunos modelos depredador-presa tipo Leslie con respuesta funcional no mon´otona y efecto Allee en las presas. Universidad Nacional de Colombia Sede Manizales, (2017). [3] E.M. Zapata Pozo,Estudios de redes complejas cooperativas en Ecolog´ıa con sistemas de ecuaciones diferenciales. Universidad de Sevilla, (2014). [4] J.D. Murray,Mathematical Biology, Springer, New York, (1993). 57