Ecología matemática: Dinámica estructural en el modelo Lotka-Volterra

Trabajo Fin de Grado Ecolog´ıa matem´atica: Din´amica Estructural en el Modelo de Lotka-Volterra Jos´e David Guti´errez de Alba Ecolog´ıa matem´atica: Din´amica Estructural en el Modelo de Lotka-Volterra Jos´e David Guti´errez de Alba Memoria presentada como parte de los requisitos para la obtenci´on del t´ıtulo de Grado en Matem´aticas por la Universidad de Sevilla. Tutorizada por Prof. Jos´e Antonio Langa Rosado Prof. Antonio Su´arez Fern´andez Agradecimientos A mis tutores D. Jos´e Antonio Langa Rosado y D. Antonio Su´arez Fenr´andez por su tiempo y dedicaci´on, no solo durante el Trabajo de Fin de Grado, sino tambi´en por todos los a˜nos que me han guiado como alumno interno del Departamento de Ecuaciones Diferenciales y An´alisis Num´erico de la Universidad de Sevilla. Sin ellos no hubiera sido posible la realizaci´on de este trabajo. A mis familiares y amigos, por su apoyo incondicional y por haber compartido conmigo todos estos a˜nos llenos de valiosos momentos. Gracias por todo. Abstract This work studies the dynamic of mutualistic ecosystems via Lotka-Volterra model. After exploring the existence and uniqueness of positive solutions, several properties of globally stable equilibrium are shown using the Linear Complementarity Theory. In addition, the structural stability of real ecological networks is analyzed in a set of simulations. The evidence arisen from these experiments suggests that cooperation plays a fundamental role in the coexistence of species in large ecosystems. Finally, a new Lotka-Volterra model that includes the indirect interactions between species is proposed, and it is compared with the traditional Lotka-Volterra model. ´ Indice general Introducci´on 7 1. Modelo de Lotka-Volterra 19 1.1. Definiciones y resultados previos . . . . . . . . . . . . . . . . . . . . 19 1.2. Descripci´on del modelo L-V . . . . . . . . . . . . . . . . . . . . . . 22 1.3. Red de un ecosistema . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2. Estabilidad global 27 2.1. Introducci´on al Problema de Complementariedad Lineal . . . . . . . 27 2.2. Estabilidad global en el modelo L-V . . . . . . . . . . . . . . . . . . 31 2.3. Estabilidad global en el modelo mutualista . . . . . . . . . . . . . . 33 3. Estabilidad estructural 41 3.1. Definiciones previas . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.2. F´ormula para el volumen de un cono esf´erico . . . . . . . . . . . . . 43 3.3. An´alisis estructural de redes reales . . . . . . . . . . . . . . . . . . 51 3.3.1. Resultados............................ 53 4. Efectos indirectos 59 4.1. Introducci´on............................... 59 4.2. Modelo de Lotka-Volterra con efectos indirectos . . . . . . . . . . . 61 4.3. Comparaci´on de la estabilidad estructural . . . . . . . . . . . . . . 65 4.3.1. Resultados............................ 66 5. Conclusiones y problemas abiertos 71 Programas de MatLab 72 Bibliograf´ıa 89 5 Introducci´on La ecolog´ıa es la ciencia que estudia las relaciones entre los seres vivos en su h´abitat natural. Es una disciplina cient´ıfica muy diversa que abarca distintas ramas como la ecolog´ıa de poblaciones, la ecolog´ıa evolutiva o la etoecolog´ıa entre otras. El enfoque de este trabajo ser´a m´as limitado: nos centraremos en la ecolog´ıa de comunidades, donde estudiaremos la din´amica de poblaciones, equilibrios, extinci´on y supervivencia de las distintas especies de un ecosistema. Para estudiar la din´amica de poblaciones usaremos el conocido modelo de Lotka-Volterra. Propuesto inicialmente por Alfred J. Lotka a principios del siglo XX, el modelo de Lotka-Volterra es considerado una base para estudiar los procesos ecol´ogicos y ha sido ampliamente estudiado en la literatura. Un sistema Lotka- Volterra determinista con Nespecies viene dado por las ecuaciones dui dt =ui bi+ N X j=1 aijuj!i= 1 . . . N (0.1) donde Nes el n´umero de especies en el ecosistema, ui(t) es el n´umero de individuos de la especie ien el tiempo t,bies la tasa de crecimiento intr´ınseco yaij son par´ametros que describen las interacciones de la especie icon el resto. El vector formado por los bicon i= 1 . . . N se denomina vector de crecimiento intr´ınseco. Seg´un el efecto neto sobre el individuo, las relaciones se clasifican en ocho tipos [20]: Competici´on: La competici´on ocurre cuando dos organismos se ven afectados negativamente al competir por los mismos recursos. Depredaci´on: un organismo (el depredador) consume a otro (la presa). La depredaci´on tambi´en incluye todo consumo de nutrientes de otro organismo, incluidos la herbivor´ıa y el parasitismo. 7 Introducci´on 14 bola unidad que ocupa CAse denomina volumen normalizado del cono esf´erico de A, denotado por V∗(CA), y viene dado por V∗(CA) = µ(pos(−A)∩B(0,1)) µ(B(0,1)) (0.9) donde µes la medida de Lebesgue N-dimensional en RNyB(0,1) es la bola centrada en el origen con radio 1. Uno de los resultados te´oricos descubiertos en este trabajo es la siguiente proposici´on, y nos indica que si las interacciones entre las especies son peque˜nas, el volumen del cono tiende a 0 cuando el n´umero de especies se hace grande: Proposici´on 0.0.11. Sea A∈RN×Ncon A=P−I, siendo Pde la forma P=     0p1,2· · · p1,N p2,10. . . . . ....pN−1,N pN,1· · · pN,N−10     . Denotemos por k·k2a la norma espectral, es decir, kAk2=pλmax(ATA), y supongamos que kPk2≤1 N. Entonces V∗(CA)≤N+ 1 1−1 NN 2N .(0.10) Por tanto, es necesario que los coeficientes de relaci´on interespec´ıfica sean suficientemente grandes para garantizar que el sistema sea estructuralmente estable. Finalizamos el cap´ıtulo analizando la estabilidad estructural de ecosistema reales. Para ello vamos a implementar el siguiente modelo [9]: γpij =γp yij kδγ pi , γaij =γa yij kδγ ai , βpij =βp sign(lij) 1 + lδβ ij , βaij =βa sign(lij) 1 + lδβ ij . Donde γpyγarepresentan el nivel basal de la intensidad de mutualismo para las plantas y los animales. Introducci´on 15 kpiykaies el n´umero de especies con las que coopera la planta o el animal i. yij = 1 si las especies iyjinteract´uan y 0 si no lo hacen. βpyβason el nivel basal de la intensidad de competici´on para las plantas y los animales. lij es el n´umero de cooperadores comunes entre la especie iy la j.sign(lij) es el signo de lij, es decir, 1 si lij >0 y 0 si lij = 0 Los par´ametros δγyδβmodulan el intercambio mutualista y la afectaci´on a la competencia, respectivamente. Los resultados obtenidos muestran que el volumen del cono de A es creciente respecto a la cooperaci´on y decreciente respecto a la competici´on. Adem´as, un aumento en el intercambio mutualista y la afectaci´on a la competencia provoca una disminuci´on del volumen del cono, aunque manteniendo las relaciones de crecimiento y decrecimiento descritas anteriormente. Tambi´en puede observarse que la estabilidad estructural disminuye al aumentar el n´umero de especies. (a) Volumen del cono de Aen funci´on de la cooperaci´on para distintos valores de competici´on, tomando δβ=δγ= 0. (b) ´ Idem que la gr´afica anterior, pero tomando δβ= δγ= 1,5. Figura 3: Estabilidad estructural de un ecosistema mutualista con 4plantas y 5animales (Figura 1). Podemos ver c´omo la cooperaci´on aumenta el volumen del cono, mientras que la competici´on lo disminuye. Cap´ıtulo 4 Cuando no se observa interacci´on entre dos especies de un ecosistema, se asume que entre ´estas existe una relaci´on de neutralismo. Sin embargo, como dijimos Introducci´on 16 anteriormente, el neutralismo rara vez ocurre en la naturaleza. Esto es porque las especies pueden afectarse indirectamente unas a otras. Por ejemplo, entre un p´ajaro insect´ıvoro y una planta, puede que no se observe ninguna interacci´on, y por tanto pensar que existe neutralismo entre ellos. Pero si el p´ajaro se alimenta de una especie de abeja que poliniza esa planta, entonces el aumento o disminuci´on de la poblaci´on de p´ajaros va a afectar a la poblaci´on de plantas. Otro objetivo de este trabajo es implementar un modelo que incluya estas relaciones indirectas en el ecosistema. Aunque existen muchas formas de estudiar estos efectos indirectos, nosotros propondremos un modelo nuevo: Definici´on 0.0.12. Sea τ > 0y consideramos el sistema de Lotka-Volterra u0= u(b+Au)siendo Auna matriz de la forma A=   −1aij ... aij −1   . Definimos el modelo de Lotka-Volterra con efectos indirectos, o simplemente modelo de efectos indirectos, como el s.d.o: (LV I)τ     u0=u(b+e Au), e A=A+τ ∞ X n=2 (I+A)n,(0.11) siempre que exista la serie, y donde τ > 0es un par´ametro que modula la intensidad de los efectos indirectos. Por analog´ıa, llamaremos modelo de efectos directos al modelo de Lotka-Volterra usual. Esta definici´on incluye en la matriz de coeficientes todas las interacciones indirectas que pueden darse en el ecosistema. Como es de esperar, el desarrollo te´orico de este cap´ıtulo se centra en proporcionar condiciones suficientes para la existencia de e A, as´ı como garantizar que el sistema (LV I)τtenga un punto de equilibro globalmente estable en RN ≥0. Proposici´on 0.0.13. Si Aes una matriz como en (0.0.12) y adem´as es diagonal dominante negativa, entonces e A=A+τ ∞ X n=2 (I+A)n est´a bien definida para todo τ > 0. Introducci´on 17 Proposici´on 0.0.14. Sea Auna matriz como en la Definici´on 0.0.12 con A∈ NDD. Entonces existe  > 0tal que si τ <  entonces e A∈NDD. Por ´ultimo, comparamos la estabilidad estructural del modelo de efectos indirectos con la del modelo de Lotka-Volterra usual. Uno de los descubrimientos m´as sorprendentes que hemos obtenido en las simulaciones es que el papel que juegan la cooperaci´on y la competici´on puede llegar a revertirse al considerar los efectos indirectos, es decir, el volumen del cono de e Ase vuelve decreciente respecto a la cooperaci´on y creciente respecto a la competici´on. Sin embargo, los ecosistemas con efectos indirectos son estructuralmente m´as estables que los de efectos directos. (a) Volumen del cono de e Aen funci´on de la cooperaci´on para distintos valores de τ, tomando βp= βa= 0,05 yδβ=δγ= 0. (b) ´ Idem que la gr´afica anterior, pero tomando βp= βa= 0,2yδβ=δγ= 0. Figura 4: Estabilidad estructural del modelo de efectos indirectos para el mismo ecosistema que en la figura 1. Cuando τ= 0 la matriz se corresponde con la del modelo de Lotka-Volterra usual. Cuanto mayor sea el valor de τ, mayor ser´a la intensidad de los efectos indirectos. Cap´ıtulo 5 En este ´ultimo cap´ıtulo hacemos un recorrido por toda la memoria, resumiendo los resultados y conclusiones m´as importantes de ´esta. Adem´as, exponemos los principales problemas abiertos que deja este trabajo Introducci´on 18 Cap´ıtulo 1 Modelo de Lotka-Volterra En este primer cap´ıtulo presentamos algunos resultados b´asicos de ecuaciones diferenciales ordinarias, que posteriormente aplicaremos al modelo de Lotka- Volterra para deducir propiedades de existencia de soluci´on y estabilidad de los equilibrios. Adem´as, definimos las redes de ecosistemas, las cuales aportan informaci´on sobre las interacciones entre especies. 1.1. Definiciones y resultados previos Comenzaremos el cap´ıtulo dando algunas definiciones y teoremas que usaremos en el transcurso de la memoria. Todos estos resultados, as´ı como sus demostraciones, se encuentran en los apuntes de las asignaturas de Ecuaciones Diferenciales Ordinarias y de Amplicaci´on de Ecuaciones Diferenciales del Grado en Matem´aticas de la Universidad de Sevilla [16]. Definici´on 1.1.1. Sea Ω⊂RN+1 abierto, f: Ω →RNy(t0, u0)∈Ω. Se llama problema de Cauchy para un sistema diferencial ordinario (s.d.o.) de primer orden y dimensi´on Nal problema de encontrar soluci´on a (P C)u0=f(t, u), u(t0) = u0.(1.1) Definici´on 1.1.2. Dado un intervalo no degenerado Iy una funci´on ϕ:I⊆R→RN, diremos que (I, ϕ)es soluci´on local del (PC) (1.1) si se cumple: 1. ϕ∈C1(I, RN), 2. (t, ϕ(t)) ∈Ω∀t∈I, 3. ϕ0(t) = f(t, ϕ(t)) ∀t∈I, 19 Cap´ıtulo 1 — Modelo de Lotka-Volterra 20 4. t0∈Iyϕ(t0) = u0. Definici´on 1.1.3. Supongamos que Ω⊂RN+1 es un abierto no vac´ıo, que f∈ CΩ; RN∩Liploc (u, Ω), y que (t0, u0)∈Ω. Definimos el conjunto S(t0, u0) = (I, ϕ):(I, ϕ)es soluci´on local del (PC) (1.1). Definici´on 1.1.4. Dada (I, ϕ)∈S(t0, u0),se dice que es una Soluci´on prolongable por la derecha (resp. izquierda) si existe (J, ψ)∈S(t0, u0) tal que: I⊆J, I 6=J, sup I∈◦ J(resp. ´ınf I∈◦ J). Soluci´on prolongable si bien lo es por la derecha o lo es por la izquierda o por ambos extremos. Soluci´on maximal o global si no es prolongable. Teorema 1.1.5. Si Ω⊆RN+1 es un abierto no vac´ıo y f∈CΩ,RN∩Liploc (u, Ω), entonces para cada (t0, u0)∈Ωexiste una ´unica soluci´on maximal de (1.1). Adem´as, el intervalo Idonde est´a definida la soluci´on maximal es abierto. Nota 1.1.6. En adelante denotaremos por ϕ(·, t0, u0)a la soluci´on maximal de (1.1) e I(t0, u0)a su intervalo de definici´on maximal. Definici´on 1.1.7. Dada una soluci´on local del (PC), (I, ϕ)∈S(t0, u0), con I= (α, β), se denomina: 1. Trayectoria de (I, ϕ)a τϕ={(t, ϕ (t)) : t∈I}. 2. Semitrayectoria derecha de (I, ϕ) τ+ ϕ={(t, ϕ (t)) : t∈I, t ≥t0}. 3. Semitrayectoria izquierda de (I, ϕ) τ− ϕ={(t, ϕ (t)) : t∈I, t ≤t0}. 4. Terminal derecho de (I, ϕ) T+ ϕ=(∅si β= +∞, τ+ ϕ∩(β, u) : u∈RNsi β < +∞.. Cap´ıtulo 1 — Modelo de Lotka-Volterra 21 5. Terminal izquierdo de (I, ϕ) T− ϕ=(∅si α=−∞, τ− ϕ∩(α, u) : u∈RNsi α > −∞.. Teorema 1.1.8. Sea (I, ϕ)∈S(t0, u0)con I= (α, β). Son equivalentes: 1. ϕes prolongable por la derecha. 2. τ+ ϕes acotada y d(τ+ ϕ, ∂Ω) >0. 3. T+ ϕ∩Ω6=∅. 4. T+ ϕ∩Ω≡ {un punto de Ω}. Teorema 1.1.9. Sea (I, ϕ)∈S(t0, u0). Entonces se tienen las siguientes equivalencias: 1. (I, ϕ)es soluci´on maximal. 2. O bien τ+ ϕno est´a acotada o bien dτ+ ϕ, ∂Ω= 0, y, adem´as: o bien τ− ϕno est´a acotada o bien dτ− ϕ, ∂Ω= 0. 3. T+ ϕ∪T− ϕ∈∂Ω. Definici´on 1.1.10. Sean 0< ρ ≤+∞,u∗∈RN. Denotemos Bρ=B(y∗, ρ)⊆RN, I⊆Run intervalo y Ω = I×Bρ. Supongamos dado un s.d.o. u0=f(t, u)con f∈C0Ω; RN∩Liploc (u, Ω). Diremos que u∗es un punto de equilibrio del sistema si se verifica que: f(t, u∗) = 0 ∀t∈I. Observaci´on 1.1.11. Supongamos que I= (0,+∞). Si u∗es un punto de equilibrio de un sistema entonces la funci´on eϕ(t)≡u∗es soluci´on de (PC) u0=f(t, u), u(0) = u∗. M´as a´un: ϕ(t, t0, u∗) = eϕ(t)∀t0>0,∀t > 0, y adem´as I(t0, u∗) = (0,+∞) Definici´on 1.1.12. Se dice que u∗es un equilibrio estable en el sentido de Lyapunov si para cada t0>0yε > 0dados, existe δ > 0tal que si |u0| ≤ δentonces: 1. I(t0, u0)⊃[t0,+∞), 2. |ϕ(t, t0, u0)−u∗|< ε ∀t∈[t0,+∞). Cap´ıtulo 1 — Modelo de Lotka-Volterra 22 Definici´on 1.1.13. Se dice que u∗es un punto de equilibrio inestable si no es estable. Definici´on 1.1.14. Se dice que u∗es un equilibrio atractivo para u0=f(t, u)si para cada t0>0existe γ > 0tal que si |u0| ≤ γentonces I(t0, u0)⊇[t0,+∞)y l´ımt→+∞ϕ(t, t0, u0) = u∗. Definici´on 1.1.15. Se dice que u∗es asint´oticamente estable si es estable y atractivo. Teorema 1.1.16. Consideramos el sistema de ecuaciones diferenciales como en (1.1) y sea u∗un punto de equilibrio. Denotemos por λi, i = 1, . . . , N los autovalores del Jacobiano de fen u∗. Se tiene: 1. Si todos los λi, i = 1, . . . , N tienen parte real negativa, entonces u∗es asint´oticamente estable. 2. Si para alg´un i= 1, . . . , N, la parte real de λies positiva, entonces u∗es inestable. 1.2. Descripci´on del modelo L-V El modelo de Lotka-Volterra es un sistema de ecuaciones diferenciales no lineales que se utiliza para describir la din´amica de ecosistemas [9]. La abundancia de cada especie cambia a lo largo del tiempo conforme al conjunto de ecuaciones:          dui dt =ui bi+ N X j=1 aijuj!i= 1 . . . N, ui(0) = u0i (1.2) donde ui(t) es el n´umero de individuos de la especie ien el tiempo t,Nes el n´umero de especies en el ecosistema y bi, aij son par´ametros que describen las interacciones entre las especies. Tambi´en se puede expresar matricialmente este sistema como: (u0=u(b+Au), u(0) = u0.(1.3) Al vector b∈RNse le denomina vector de crecimiento intr´ınseco, y a la matriz A∈RN×N, matriz de interacciones. Cap´ıtulo 1 — Modelo de Lotka-Volterra 23 Un tipo de modelo de Lotka-Volterra es el modelo mutualista. En ´este, las especies est´an divididas en dos grupos (generalmente plantas y animales) existitendo una relaci´on simbi´otica entre ellos, mientras que las especies de cada grupo est´an en competici´on entre s´ı.                    dui dt =ui bpi− P X j=1 βpij uj+ A X k=1 γpik vk!i= 1, . . . , P, dvi dt =vi bai− A X j=1 βaij vj+ P X k=1 γaik uk!i= 1, . . . , A, ui(0) = ui0i= 1, . . . , P, vi(0) = vi0i= 1, . . . , A (1.4) donde PyAson el n´umero de especies plantas y animales respectivamente; uiyvirepresentan las biomasas de las plantas y los animales; bpi, bai∈Rson las tasas de natalidad; βpij , βaij ≥0 las tasas de competencia de cada subsistema; y γpik , γaik ≥0 la intensidad de la interacci´on mutualista entre las especies de los dos subsistemas. Haciendo un cambio de variable, podemos tomar βpii =βajj = 1, i = 1, . . . , P;j= 1, . . . , A . El sistema de Lotka-Volterra mutualista puede escribirse matricialmente como: u0=u(b+Mu),M=B1Γ1 Γ2B2N×N(1.5) y cada bloque se escribe B1=     −1−βp12 · · · −βp1P −βp21 −1· · · −βp2P . . .. . ..... . . −βpP1−βpP2· · · −1     , B2=     −1−βa12 · · · −βa1A −βa21 −1· · · −βa2A . . .. . ..... . . −βaA1−βaA2· · · −1      Γ1=   γp11 · · · γp1A . . ..... . . γpP1· · · γpP A   ,Γ2=   γa11 · · · γa1P . . ..... . . γaA1· · · γaAP   . A continuaci´on vamos a aplicar al modelo de Lotka-Volterra los resultados te´oricos de existencia de soluci´on y estabilidad vistos en la secci´on anterior [9]. Lema 1.2.1. Sean u1, . . . , uN∈Ryβ < 1. Entonces: N X i=1 u2 i+ 2β N X i<j uiuj≥1 + β(N−1) N N X i=1 ui!2 .(1.6) Cap´ıtulo 2 — Estabilidad global 30 Antes de probar este teorema vamos a enunciar dos lemas que usaremos en la demostraci´on. Lema 2.1.11. Sean A∈RN×N,x= (xi)un vector columna e y= (yi) = Ax. Entonces Aes una P-matriz si y solo si para todo xverificando xiyi≤0∀ise tiene que x= 0 Lema 2.1.12. Sea Muna matriz con todos los todos los menores principales positivos. Entonces existen (w, z)soluci´on de LCP(q, M)para todo q∈RN Las demostraciones de estos lemas pueden encontrarse en [7] y [5]. Demostraci´on. (Teorema 2.1.10).Supongamos que Mno es una P-matriz. Entonces por el Lema 2.1.11 existe x= (xi)∈RN,x6= 0, tal que y= (yi) = Mx y, para cada i= 1,· · · , N,xieyitienen signos contrarios (i.e. xiyi≤0). Sean y+ i= m´ax {yi,0}, y− i= m´ın {yi,0},x+ i= m´ax {xi,0}yx− i= m´ın {xi,0}. Entonces yi=y+ i−y− i, y+ i, y− i≥0, y+ iy− i= 0 ∀i= 1, . . . , N, xi=x+ i−x− i, x+ i, x− i≥0, x+ ix− i= 0 ∀i= 1, . . . , N. Como xiyi≤0 para cada i= 1, . . . , N, se verifica que x+ iy+ i=x− iy− i= 0 para cada i= 1, . . . , N. As´ı que (y+)Tx+= (y−)Tx−= 0. Como y=Mx, tenemos y+−Mx+=y−−Mx−=: bq. Y adem´as x6= 0, x+6=x−. Concluimos que, cuando q=bq, el problema LCP(q, M) tiene dos soluciones distintas: (w, z) = (y+, x+) y (w, z) = (y−, x−). As´ı, si Mno es una P-matriz, existe q∈RNpara el cual el problema LCP(q, M) tiene dos soluciones distintas. Cap´ıtulo 2 — Estabilidad global 31 Ahora supongamos que Mes una P-matriz. Entonces por el Lema 2.1.12, el problema LCP(q, M) tiene al menos una soluci´on para cada q∈RN. Supongamos que existe un q∈RNpara el cual LCP(q, M) tiene dos soluciones distintas, llamadas (w, z) y (bx, bz). Entonces ( ¯w−ˆw) = M(¯z−ˆz) y, como estas dos soluciones son distintas, ¯z−ˆz6= 0. Por la condici´on de complementariedad ( ¯w)T¯z= ( ˆw)Tˆz= 0 se verifica que ( ¯wi−ˆwi) (¯zi−ˆzi)⩽0 para todo i= 1, . . . , N y, como ¯z−ˆz6= 0, esto implica por el Lema 2.1.11 que Mno es una P-matriz, lo cual es una contradicci´on. As´ı que la soluci´on de LCP(q, M) debe ser ´unica para cada q∈RN Definici´on 2.1.13. Una clase de conos convexos en RNse dice que es una partici´on si: 1. Cada cono tiene interior no vac´ıo. 2. La uni´on de los conos es RN. 3. Los interiores de cada par de conos son disjuntos. Teorema 2.1.14. Sea M∈RN×N. La clase de conos complementarios C(M)es una partici´on de RNsi y solo si Mes una P-matriz. 2.2. Estabilidad global en el modelo L-V El siguiente lema es clave para probar la equivalencia entre resolver el LCP y hallar un punto estacionario no negativo de (L−V). Lema 2.2.1. Consideremos el sistema (1.3). El Problema de Complementariedad Lineal (−b, −A)es equivalente al problema de encontrar un punto de equilibrio no negativo u∗del sistema (L−V)que satisfaga bi+ N X j=1 aiju∗ j≤0para i= 1, . . . , N. (2.2) Demostraci´on. Basta tomar z=u∗yw=−b−Au∗ por lo que tomando M=−Ayq=−bse tiene que w=Mz +q, w ≥0, z ≥0 y wizi= 0 Cap´ıtulo 2 — Estabilidad global 32 Definici´on 2.2.2. Sean u∗= (u∗ 1, . . . , u∗ N)un punto estacionario no negativo de (L−V),I⊆ {1, . . . , N}tal que u∗ i= 0 ∀i∈IyJ={1, . . . , N} \ I. Definimos el conjunto: RN I=u= (u1, u2, . . . , uN)∈RN +|ui≥0para i∈Iyuj>0para j∈J. En adelante consideraremos que un sistema es globalmente estable si lo es en RN I, es decir, si existe un punto estacionario asint´oticamente estable cuya cuenca de atracci´on sea RN I. Teorema 2.2.3. Supongamos que A∈Sw.Entonces el sistema (L-V) (1.3) tiene un punto de equilibrio u∗satisfaciendo (2.2) para cada b∈RNque, adem´as, es globalmente estable en RN I. Demostraci´on. Por el Lema 2.1.7 como A∈Sw,se sigue que −Aes una Pmatriz. Ahora, por el Teorema 2.1.10 y el Lema 2.2.1 se tiene la existencia y unicidad de un punto saturado u∗, esto es, existe Iun subconjunto de {1, . . . , N}satisfaciendo que u∗ i= 0 para cualquier i∈Iy bi+ N X j=1 aiju∗ j≤0 para i= 1, . . . , N. Denotemos por J={1, . . . , N}\I. Se considera la funci´on V(u) = X j∈J wjuj−u∗ j−u∗ jlog uj u∗ j+X i∈I wiui donde W= diag (wi), wi>0,es la matriz diagonal tal que WA+ATWes definida negativa, y se define el conjunto Ω = Ω(L) := u∈RN I:V(u)≤Lu(0) donde Lu(0) es una constante positiva que depende del dato inicial u(0),tal que Lu(0) ≥V(u(0)). Se observa en primer lugar que: 1. V(u)≥0 en Ω. 2. V(u) = 0 s´olo en u=u∗. Cap´ıtulo 2 — Estabilidad global 33 Calculamos su derivada a lo largo de las soluciones de (1.2) ˙ V(u(t)) = X j∈J wj1−u∗ j/uju0 j+X i∈I wiu0 i =X j∈J wjuj−u∗ jN X k=1 ajk (uk−u∗ k) + X i∈I wiui N X k=1 aik (uk−u∗ k) +X i∈I wiui bi+ N X k=1 aiku∗ k! =1 2(u−u∗)TWA +ATW(u−u∗) + X i∈I wiui bi+ N X k=1 aiku∗ k!. Ya que A∈Sw,el primer t´ermino es negativo, y puesto que u∗es un punto de equilibrio saturado, el segundo t´ermino tambi´en es negativo. As´ı, cada soluci´on permanece en Ω para todo t≥0 si u(0) ∈Ω, y todas las soluciones que comienzan en Ω cumplen que u→u∗cuando t→ ∞ por el Teorema de Estabilidad de Lyapunov. Por otro lado, es trivial que u∗es estable en Ω. De hecho, Ω es positivamente invariante. La uni´on de los conjuntos Ω(L) cuando L→ ∞ converge a RN I. Por lo tanto, u∗es estable en RN Iy cada soluci´on converge a u∗cuando t→ ∞,si u(0) ∈RN I. Esto completa la prueba. Corolario 2.2.4. Si A∈Sw, entonces el sistema (L-V) y cada uno de sus subsistemas tienen un punto de equilibrio no negativo y globalmente estable para cada b∈RN Definici´on 2.2.5. Sea u∗= (u∗ 1, . . . , u∗ N)un punto estacionario de (1.3). Diremos que u∗es factible si u∗es globalmente estable en RN +y adem´as u∗ i≥0∀i= 1, . . . , N. Llamaremos biodiversidad del ecosistema al n´umero de componentes positivas de u∗. Aqu´ı, la interpretaci´on biol´ogica de la biodiversidad es el n´umero de especies que sobreviven en el ecosistema cuando el tiempo tiende a infinito. En el Cap´ıtulo 3 veremos c´omo podemos utilizar esta definici´on para comparar ecosistemas. 2.3. Estabilidad global en el modelo mutualista Para poder aplicar los resultados de estabilidad global vamos a definir tres tipos de matrices. Estas matrices van a ser como en el modelo mutualista (1.5), pero cumpliendo adem´as unas restricciones que garantizan la existencia de un punto estacionario factible. Cap´ıtulo 2 — Estabilidad global 34 Definici´on 2.3.1. Diremos que una matriz Mdefinida como en (1.5) es de: Tipo I, si es diagonal dominante negativa, es decir A X k=1 γpik + P X j6=i βpij <1,∀i= 1, . . . , P, P X k=1 γaik + A X j6=i βaij <1,∀i= 1, . . . , A. Tipo II, si satisface las condiciones: 2− P X j6=iβpij +βpji − A X k=1 γpik >0,∀i= 1, . . . , P, (2.3) 2− A X j6=iβaij +βaji − P X k=1 γaik >0,∀i= 1, . . . , A, (2.4) sup i=1,...,A P X k=1 γpki 2− A X j6=iβaij +βaji − P X k=1 γaik <´ınf i=1,...,P 2− P X j6=iβpij +βpji − A X k=1 γpik A X k=1 γaki . (2.5) Tipo III, si β1=βpij <1, β2=βaij <1, γ1=γpij , γ2=γaij ∀i= 1, . . . , P, ∀j= 1, . . . , A, (2.6) γ1γ2<(1 + β1(P−1)) (1 + β2(A−1)) PA .(2.7) Nota 2.3.2. A las matrices que satisfacen (2.6) se les denomina matriz de campo medio. Concluiremos esta secci´on con un teorema sobre la estabilidad de estos tipos de matrices. Pero antes introduciremos algunos resultados ´utiles para su demostraci´on. Todos estos resultados pueden encontrarse en [9]. Cap´ıtulo 2 — Estabilidad global 35 Teorema 2.3.3. (Teorema de Gershgorin). Sea A= (aij)∈RN×Nyλun autovalor cualquiera de A. Entonces existe i∈ {1, . . . , N}tal que |λ−aii| ≤ n X j6=i |aij|. Definici´on 2.3.4. Dados N, M ∈N, y a, b ∈R, denotamos por D(a, b;M) =        a b b · · · b bab· · · b b b a · · · b . . .. . .. . ..... . . bbb· · · a       M×M e IM×N= (aij), aij = 1, i = 1, . . . , M;j= 1, . . . , N. Proposici´on 2.3.5. Consideramos la matriz C=DPC1 C2DA(P+A)×(P+A) donde DP=D(a1, d1;P), DA=D(a2, d2;A), C1=c1IP×A, C2=c2IA×P con a1, d1, a2, d2, c1, c2∈R. Se cumple que |C|= (a1−d1)P−1(a2−d2)A−1((a1+ (P−1)d1) (a2+ (A−1)d2)−c1c2AP). Teorema 2.3.6. Sea Muna matriz definida como en (1.5). Si Mes de Tipo I, Tipo II o Tipo III, entonces M∈Sw Demostraci´on. Si Mes de Tipo I, por el Teorema 2.1.8 se tiene que M∈Sw. Sea Muna matriz de Tipo II de la forma MT=BT 1ΓT 2 ΓT 1BT 2(P+A)×(P+A) . Cap´ıtulo 2 — Estabilidad global 36 Escribimos la matriz Wde la siguiente manera W=W1∅ ∅W2(P+A)×(P+A) (2.8) donde las matrices W1yW2se consideran como W1=     w10· · · 0 0w1· · · 0 . . .. . ..... . . 0 0 · · · w1     P×P , W2=     w20· · · 0 0w2· · · 0 . . .. . ..... . . 0 0 · · · w2     A×A con w1>0, w2>0. Se obtiene C=WM+MTW=D1C1 C2D2(P+A)×(P+A) con D1=w1     −2−(βp12 +βp21 )· · · − (βp1P+βpP1) −(βp12 +βp21 )−2· · · − (βp2P+βpP2) . . .. . .. . .. . . −(βp1P+βpP1)−(βp2P+βpP2)· · · −2      D2=w2     −2−(βa12 +βa21 )· · · − (βa1A+βaA1) −(βa12 +βa21 )−2· · · − (βa2A+βaA2) . . .. . .. . .. . . −(βa1A+βaA1)−(βa2A+βaA2)· · · −2      C1=     w1γp11 +w2γa11 w1γp12 +w2γa21 · · · w1γp1A+w2γaA1 w1γp21 +w2γa12 w1γp22 +w2γa22 · · · w1γp2A+w2γaA2 . . .. . .. . .. . . w1γpP1+w2γa1Pw1γpP2+w2γa2P· · · w1γpP A +w2γaAP      C2=     w1γp11 +w2γa11 w1γp21 +w2γa12 · · · w1γpP1+w2γa1P w1γp12 +w2γa21 w1γp22 +w2γa22 · · · w1γp22 +w2γa2P . . .. . .. . .. . . w1γp1A+w2γaA1w1γp2A+w2γaA2· · · w1γpP A +w2γaAP     . Luego, por el Teorema 2.3.3 (Teorema de Gershgorin), de la matriz Cse generan Pdesigualdades de la forma |λ+ 2w1| ≤ w1 P X j6=iβpij +βpji +w1 A X k=1 γpik +w2 A X k=1 γaki ,con i= 1, . . . , P Cap´ıtulo 2 — Estabilidad global 37 yAdesigualdades |λ+ 2w2| ≤ w2 A X j6=iβaij +βaji +w2 P X k=1 γaik +w1 P X k=1 γpki ,con i= 1, . . . , A. Como se desea que los autovalores λsean negativos, reescribimos las PyA desigualdades anteriores como λ≤w1PP j6=iβpij +βpji +PA k=1 γpik −2+w2PA k=1 γaki <0,∀i= 1, . . . , P, λ≤w1PP k=1 γpki +w2PA j6=iβaij +βaji +PP k=1 γaik −2<0,∀i= 1, . . . , A. Obtenemos el siguiente sistema de P+Adesigualdades con variables w1y w2. PP j6=iβpij +βpji +PA k=1 γpik −2w1+PA k=1 γaki w2<0,∀i= 1, . . . , P, PP k=1 γpki w1+PA j6=iβaij +βaji +PP k=1 γaik −2w2<0,∀i= 1, . . . , A. Resolviendo las desigualdades para w2y considerando que w2>0,se tiene w2 w1 <−PP j6=iβpij +βpji +PA k=1 γpik −2 PA k=1 γaki ,para toda i= 1, . . . , P. Ahora resolvemos la otra desigualdad tambi´en para w2>0, w2 w1 >−PP k=1 γpki PA j6=iβaij +βaji +PP k=1 γaik −2,para toda i= 1, . . . , A. Por tanto, podemos encontrar la matriz Wsi fijamos un valor para w1>0 yw2que cumplan estas desigualdades, si suponemos (2.5). Sea Muna matriz de Tipo III de orden (P+A), escrita por bloques como en (1.5), con 0 < β1, β2<1, γ1, γ2>0. La matriz W W=W1∅ ∅W2(P+A)×(P+A) donde las matrices W1yW2son como en (2.8). Como antes C=WM+MTW=D1C1 C2D2(P+A)×(P+A) Cap´ıtulo 2 — Estabilidad global 38 donde D1=D(−2w1,−2w1β1;P), D2=D(−2w2,−2w2β2;A), C1= (w1γ1+w2γ2)IP×A, C2= (w1γ1+w2γ2)IA×P. Nos planteamos encontrar la condici´on para que los autovalores de WM+ MTWsean negativos. Para eso, calculamos las soluciones del polinomio caracter´ıstico de C p(λ) = |C(λ)|=|C−λIA×P|=WM+MTW−λIA×P= 0 donde IA×Pdenota la matriz identidad de orden A×P. La matriz C(λ) se puede escribir de la siguiente forma C(λ) = F1C1 C2F2(P+A)×(P+A) donde F1=D(−2w1−λ, −2w1β1;P), F2=D(−2w2−λ, −2w2β2;A) Por tanto, aplicando la Proposici´on 2.3.5, se tiene p(λ) = (−2w1−λ+ 2w1β1)P−1(−2w2−λ+ 2w2β2)A−1·q(λ) donde q(λ) = (−2w1−λ−2w1β1(P−1)) (−2w2−λ−2w2β2(A−1)) −PA (w1γ1+w2γ2)2 = 4"w1+λ 2+w1β1(P−1)w2+λ 2+w2β2(A−1) −PA 4(w1γ1+w2γ2)2#. Los autovalores de la matriz Cson:    λ= 2w1(β1−1) con multiplicidad P−1, λ= 2w2(β2−1) con multiplicidad A−1, las ra´ıces de q(λ). Cap´ıtulo 2 — Estabilidad global 39 Las ra´ıces de q(λ) son las ra´ıces del siguiente polinomio de grado 2 , m(λ) = λ2 4+R1(w1, w2)λ 2+R2(w1, w2) donde R1(w1, w2) = w1(1 + β1(P−1)) + w2(1 + β2(A−1)) , R2(w1, w2) = w1w2(1 + β1(P−1)) (1 + β2(A−1)) −PA 4(w1γ1+w2γ2)2. Teniendo en cuenta que R1>0,para que las ra´ıces de m(λ) sean negativas, se necesita que R2(w1, w2)>0. Esto es equivalente a w1w2(1 + β1(P−1)) (1 + β2(A−1)) >PA 4w2 1γ2 1+w2 2γ2 2+ 2w1w2γ1γ2 de donde, 0> w2 1γ2 1+w2 2γ2 2+ 2w1w2γ1γ2−2 PA (1 + β1(P−1)) (1 + β2(A−1)) Dividiendo por w2 1se tiene que 0>w2 w12 γ2 2+2 w2 w1γ1γ2−2 PA (1 + β1(P−1)) (1 + β2(A−1))+γ2 1. Teniendo en cuenta (1.7), para que existan w1, w2>0 que verifiquen la desigualdad anterior, es suficiente que γ1γ2−2 PA (1 + β1(P−1)) (1 + β2(A−1))2 > γ2 1γ2 2 y, por tanto, se debe cumplir que γ1γ2<2 PA (1 + β1(P−1)) (1 + β2(A−1)) −γ1γ2 o, equivalentemente, γ1γ2<(1 + β1(P−1)) (1 + β2(A−1)) PA . Cap´ıtulo 3 — Estabilidad estructural 46 Sustituyendo (3.3) en (3.2) obtenemos: V(CA) = |det(A)|ZSA rN−1 N−2 Y i=1 sinN−i−1φidr dφ1· · · dφN−1. Integrando en ry aplicando el l´ımite de integraci´on sobre r, 0 < r < 1 kAΦk, llegamos a la f´ormula: V(CA) = |det(A)| NZπ/2 0 · · · Zπ/2 0QN−2 i=1 sinN−i−1(φi) kAΦ (φ1, . . . , φN−1)kNdφ1· · · dφN−1 con Φ(φ1, . . . , φN−1) =      cos φ1 sin φ1cos φ2 . . . sin φ1. . . sin φN−1     . Recordemos que la estabilidad estructural se define trav´es del volumen normalizado del cono. Usando que el volumen de la bola unidad en RNviene dado por [19] V(BN) = 2bN+1 2cπbN 2c N!! (3.4) donde b·c es la parte entera y N!! = N(N−2)(N−4) . . . 1 es el doble factorial, podemos dar el siguiente corolario para el volumen normalizado del cono de A: Corolario 3.2.3. Sea A∈RN×Nuna matriz invertible. Entonces V∗(CA) = |det(A)|(N−2)!! 2bN+1 2cπbN 2cZπ/2 0 · · · Zπ/2 0QN−2 i=1 sinN−i−1(φi) kAΦkNdφ1· · · dφN−1. (3.5) En general, calcular de manera exacta el volumen del cono de una matriz cualquiera es imposible, salvo en casos muy particulares. Uno de esos casos es el del cono de la matriz identidad. Como en RNhay 2Nortantes, y el cono de la identidad es precisamente el primer ortante, la porci´on que ocupa el cono en la bola unidad deber´ıa ser 1 2N. Vamos a probar rigurosamente este resultado, pero antes enunciaremos un lema que usaremos en la demostraci´on. Cap´ıtulo 3 — Estabilidad estructural 47 Lema 3.2.4. Sea n∈N. Se tienen las siguientes identidades: 1. Rπ/2 0sinn(x)dx =(n−1)!! n!! ×(1si nes impar π 2si nes par 2. n=n+1 2+n 2. Omitiremos la demostraci´on de este lema, ya que la primera identidad puede probarse f´acilmente usando integraci´on por partes, mientras que para la segunda solo hay que usar las propiedades de la funci´on parte entera. Proposici´on 3.2.5. Sea Ila matriz identidad N×N. Entonces V∗(CI) = 1 2N.(3.6) Demostraci´on. Primero, calcularemos la integral Zπ/2 0 · · · Zπ/2 0QN−2 i=1 sinN−i−1(φi) kIΦkNdφ1· · · dφN−1. Observemos que, como kΦ(φ1, . . . , φN−1)k= 1 para todo φ1, . . . , φN−1en 0,π 2, la integral nos queda Zπ/2 0 · · · Zπ/2 0 N−2 Y i=1 sinN−i−1(φi)dφ1· · · dφN−1. Esta integral la podemos expresar como el producto de integrales: N−2 Y i=1 Zπ/2 0 sinN−i−1(φi)dφi!·Zπ/2 0 dφN−1. Usando el lema anterior nos queda que Zπ/2 0 · · · Zπ/2 0 N−2 Y i=1 sinN−i−1(φi)dφ1· · · dφN−1= N−2 Y i=1 (N−i−2)!! (N−i−1)!! π 2απ 2 siendo αla cantidad de n´umeros pares menores o iguales que N−2. Puede comprobarse que α=N−2 2y que α+1 = N 2. Adem´as, el producto del segundo miembro de la igualdad anterior es telesc´opico: N−2 Y i=1 (N−i−2)!! (N−i−1)!! =(N−3)!! (N−2)!! (N−4)!! (N−3)!! · · · 1!! 2!! =1 (N−2)!!. Cap´ıtulo 3 — Estabilidad estructural 48 Por tanto Zπ/2 0 · · · Zπ/2 0 N−2 Y i=1 sinN−i−1(φi)dφ1· · · dφN−1=1 (N−2)!! πbN/2c 2bN/2c. Sustituyendo en el volumen de la identidad nos queda: V∗(CI) = (N−2)!! 2bN+1 2cπbN 2cZπ/2 0 · · · Zπ/2 0 N−2 Y i=1 sinN−i−1(φi)dφ1· · · dφN−1 =(N−2)!! 2bN+1 2cπbN 2c 1 (N−2)!! πbN 2c 2bN 2c =1 2bN+1 2c2bN 2c=1 2N. Es decir, el volumen del cono de la identidad tiende a 0 cuando la dimensi´on tiende a infinito. Cabe esperar que si una matriz Ase parece a la identidad, entonces V∗(CA) tambi´en tienda a 0: Proposici´on 3.2.6. Sea A∈RN×Ncon A=P−I, siendo Pde la forma P=     0p1,2· · · p1,N p2,10. . . . . ....pN−1,N pN,1· · · pN,N−10     . Denotemos por k·k2a la norma espectral, es decir, kAk2=pλmax(ATA), y supongamos que kPk2≤1 N. Entonces V∗(CA)≤N+ 1 1−1 NN 2N .(3.7) Demostraci´on. Tenemos que el volumen normalizado de Aviene dado por V∗(CA) = |det(A)|(N−2)!! 2bN+1 2cπbN 2cZπ/2 0 · · · Zπ/2 0QN−2 i=1 sinN−i−1(φi) kAΦkNdφ1· · · dφN−1. Para acotar superiormente esta expresi´on, vamos a: 1. Acotar superiormente |det(A)|por una constante. Cap´ıtulo 3 — Estabilidad estructural 49 2. Acotar inferiormente kAΦkNpor otra constante Ambas constantes dependientes de N. Para acotar kAΦkNusaremos que kx−yk ≥ kxk − kykpara cualquier x, y en RN: kAΦk=kPΦ−Φk ≥ kPΦk−kΦk. Recordemos que kΦk= 1. Como la norma espectral es la norma subordinada de la norma eucl´ıdea, se tiene que kPΦk≤kPk2kΦk=kPk2≤1 N. De esta forma kAΦk ≥ kPΦk − 1= 1 − kPΦk ≥ 1−1 N y, consecuentemente, kAΦkN≥1−1 NN . Ahora probaremos que |det(A)| ≤ N+ 1. Denotemos por m= m´ax 1≤i,j≤N|pi,j|. Puede verse en [12] (p´ag. 365) que m≤ kPk2. Por otro lado, se tiene una f´ormula cerrada para el determinante de una matriz, que viene dada por: det(A) = X σ∈SN sign(σ) N Y i=1 ai,σ(i)! donde SNes el conjunto de todas las permutaciones de Nelementos, y se tiene que |SN|=N! Aplicando la desigualdad triangular, obtenemos que: |det(A)| ≤ X σ∈SN N Y i=1 ai,σ(i). Pero por la forma en la que hemos definido A, se tiene que para cualquier σ∈SN y para todo i= 1, . . . , N ai,σ(i)= 1, si ies un punto fijo de σ(es decir, si σ(i) = i), Cap´ıtulo 3 — Estabilidad estructural 50 ai,σ(i)≤m, si σ(i)6=i. As´ı que podemos acotar superiormente Pσ∈SNQN i=1 ai,σ(i)por la siguiente cantidad: 1·n º de permutaciones con Npuntos fijos + m·n º de permutaciones con N−1 puntos fijos + m2·n º de permutaciones con N−2 puntos fijos + . . . + mN−1·n º de permutaciones con 1 punto fijos + mN·n º de permutaciones con 0 puntos fijos. Las permutaciones que no contienen puntos fijos se denominan desarreglos. Un resultado conocido de la combinatoria es que el n´umero de desarreglos de kelementos es k! k X r=0 (−1)r r!. De esta forma, el n´umero de permutaciones de SNcon N−kpuntos fijos viene dado por N N−kk! k X r=0 (−1)r r!. Como Pk r=0 (−1)r r!≤1 para todo k, y m≤1 N, tenemos que |det(A)| ≤ N X k=0 mkN N−kk! k X r=0 (−1)r r! ≤ N X k=0 1 Nk N! k!(N−k)!k!·1 = N X k=0 1 Nk N! (N−k)! =1+1+N(N−1) N2+N(N−1)(N−2) N3+· · · +N(N−1) · · · 3·2 NN−1+N! NN ≤1 + 1 + N·N N2+N·N·N N3+· · · +NN−1 NN−1+NN NN =N+ 1. Cap´ıtulo 3 — Estabilidad estructural 51 Uniendo las acotaciones de |det(A)|ykAΦkNobtenemos finalmente que V∗(CA) = |det(A)|(N−2)!! 2bN+1 2cπbN 2cZπ/2 0 · · · Zπ/2 0QN−2 i=1 sinN−i−1(φi) kAΦkNdφ1· · · dφN−1 ≤N+ 1 1−1 NN·(N−2)!! 2bN+1 2cπbN 2cZπ/2 0 · · · Zπ/2 0 N−2 Y i=1 sinN−i−1(φi)dφ1· · · dφN−1 =N+ 1 1−1 NNV∗(CI) = N+ 1 1−1 NN 2N . Nota 3.2.7. Una interpretaci´on de este resultado es que los ecosistemas con un gran n´umero de especies no son estructuralmente estables cuando las interacciones entre estas especies son d´ebiles. Esto es porque el volumen del cono de la matriz es tan peque˜no que cualquier perturbaci´on del vector de crecimiento intr´ınseco hace que al menos una especie se extinga. 3.3. An´alisis estructural de redes reales Para finalizar el cap´ıtulo, vamos a calcular el volumen de los conos de matrices de ecosistemas reales. Los m´etodos num´ericos deterministas resultan ser demasiado lentos para este tipo de c´alculos. Es por ello que optaremos por desrrollar un m´etodo de Montecarlo para calcular la porci´on de la bola unidad que ocupa el cono. Adem´as del cono de una matriz, existen otros conos complementarios a los que puede pertenecer el vector de crecimiento intr´ınseco. En funci´on de a qu´e cono pertenezca este vector, la biodiversidad ser´a mayor o menor. Dada una matriz A de un sistema Lotka-Volterra, llamaremos cono de biodiversidad ka un cono complementario de A, tal que si bpertenece a este cono entonces la biodiversidad del sistema es k. De esta forma, el conjunto de conos complementarios de una matriz divide RNen 2Nconos, existiendo N kconos de biodiversidad k. Llamaremos regi´on de biodiversidad ka la uni´on de todos los conos de biodiversidad k. Cap´ıtulo 3 — Estabilidad estructural 52 Para profundizar en el estudio de los ecosistemas reales vamos a medir tanto el volumen del cono como el del resto de conos complementarios. La idea detr´as del algoritmo que utilizaremos es la siguiente: 1. Tomamos A∈RN×Ncomo dato de entrada. 2. Generamos Mpuntos aleatorios uniformemente distribuidos en la bola unidad de RN. Llamaremos Ma esta nube de puntos. 3. Definimos el vector X= 0 ∈RN+1. 4. Tomamos un punto b∈My calculamos la biodiversidad del sistema u0= u(b+Au), es decir, el n´umero de componentes positivas del punto estacionario globalmente estable. 5. Si esta biodiversidad es k, entonces aumentamos en 1 la componente k+1- ´esima de X:Xk+1 ←Xk+1 + 1. 6. Haciendo los pasos 4) y 5) para cada b∈M, nos queda que el vector X almacena en la componente k+ 1 el n´umero de puntos de Mque generan biodiversidad k. 7. Dividiendo Xpor Mobtenemos en Xk+1 una aproximaci´on de la porci´on de la bola unidad que ocupan los conos de biodiversidad k. 8. Por ´ultimo, dividimos cada Xk+1 por N k, de forma que cada Xk+1 representa el volumen promedio del cono de biodiversidad k. En resumen: X1yXNdan una aproximaci´on del volumen de los conos de biodiversidad 0 y Nrespectivamente. X2,...XN−1dan una aproximaci´on de la media del volumen de los conos de biodiversidad 1, . . . , N −1, respectivamente. A la hora de implementar el modelo, vamos a utilizar los siguientes coeficientes para las matrices mutualistas [9] γpij =γp yij kδγ pi , γaij =γa yij kδγ ai , βpij =βp sign(lij) 1 + lδβ ij , βaij =βa sign(lij) 1 + lδβ ij . Donde Cap´ıtulo 3 — Estabilidad estructural 53 γpyγarepresentan el nivel basal de la intensidad de mutualismo para las plantas y los animales. kpiykaies el n´umero de especies con las que coopera la planta o el animal i. yij = 1 si las especies iyjinteract´uan y 0 si no lo hacen. βpyβason el nivel basal de la intensidad de competici´on para las plantas y los animales. lij es el n´umero de cooperadores comunes entre la especie iy la j.sign(lij) es el signo de lij, es decir, 1 si lij >0 y 0 si lij = 0 Los par´ametros δγyδβmodulan el intercambio mutualista y la afectaci´on a la competencia, respectivamente. 3.3.1. Resultados A continuaci´on vamos a estudiar la estabilidad estructural de tres ecosistemas: M SD 029: Ecosistema del Parque Natural de Do˜nana de tipo plantas - dispersores de semilla con 4 plantas y 5 animales. M PL 042: Ecosistema de las islas Gal´apagos de tipo plantas - polinizadores a con 12 plantas y 6 animales. M AF 001: Ecosistema del mar Ar´abigo de tipo an´emonas - peces con 10 an´emonas y 26 peces. Aunque las an´emonas sean animales, juegan el papel de plantas, luego nos referiremos a ellas como plantas. Las redes de ´estos ecosistemas pueden encontrarse en [38]. Lo primero que vamos a estudiar es el efecto de la cooperaci´on en el volumen del cono. En la figura 3.2 podemos observar c´omo aumenta el volumen del cono de cada ecosistema conforme aumenta la cooperaci´on (figuras 3.2 b, d y f). Cuanto mayor es la competici´on, menor es el volumen del cono. Sin embargo, en todas las gr´aficas el volumen es creciente con respecto a la cooperaci´on. Cuando aumentamos βpyβa, las gr´aficas est´an definidas en un intervalo m´as peque˜no. Esto es debido a que disminuye el rango de valores de cooperaci´on para los cuales la matriz es Tipo I. Es posible que la matriz siga siendo de clase Swaunque ya no sea Tipo I, sin embargo, por simplificar el estudio, solo vamos a considerar las matrices NDD. Al aumentar el intercambio mutualista o la afectaci´on a la Cap´ıtulo 3 — Estabilidad estructural 54 competencia (δγyδβ), disminuimos los coeficientes de cooperaci´on y competici´on, lo que aumenta el rango de valores de γp,γa,βpyβppara los cuales la matriz es digonal dominante. En la figura 3.3 podemos ver c´omo cambia el volumen del cono para distintos valores de δγyδβ. Si nos fijamos en las gr´aficas de 3.2 y 3.3 podemos observar que, fijando los los valores de γp,γa,βpyβp, el volumen del cono disminuye conforme aumentan δγyδβ. Sin embargo, se mantiene la relaci´on de crecimiento - decrecimiento respecto a la cooperaci´on y la competici´on. Es posible conseguir una mejor descripci´on de c´omo afectan la cooperaci´on y la competici´on a la estabilidad del ecosistema midiendo el volumen de los conos complementarios. En la figura 3.4 podemos ver algunos ejemplos del volumen de los conos complementarios para distintos valores de cooperaci´on y competici´on. Los conos complementarios que tienen menor biodiversidad ocupan, de media, un volumen proporcionalmente mayor que los de biodiversidad alta cuando la cooperaci´on es peque˜na respecto a la competici´on. Sin embargo, este efecto se revierte por completo cuando aumentamos la cooperaci´on. Si repetimos el experimento aumentando los valores de δγyδβ, vemos que la cooperaci´on pierde esta capacidad de recuperar biodiversidad, es decir, de aumentar el volumen de los conos de biodiversidad alta (figura 3.5). Cap´ıtulo 3 — Estabilidad estructural 55 (a) Red del ecosistema M SD 029. (b) Volumen del cono de M SD 029 (c) Red del ecosistema M PL 042. (d) Volumen del cono de M PL 042. (e) Red del ecosistema M AF 001. (f) Volumen del cono de M AF 001. Figura 3.2: Para cada ecosistema representamos la red mutualista (izquierda) y el volumen del cono respecto a la cooperaci´on para distintos valores de competici´on (derecha). Podemos ver c´omo el volumen del cono aumenta respecto a la cooperaci´on. Adem´as, cuanto mayor sea la competici´on, menor es este volumen. En todos los c´alculos hemos tomado δγ=δβ= 0. Cap´ıtulo 4 — Efectos indirectos 62 Caso n →n+1. Supongamos que Bn ij =PN k1,...,kn−1=1 bik1· · · bkn−1j Bn+1 ij =Bi·Bn ·j= N X r=1 bir N X k1,...,kn−1=1 brk1bk1k2· · · bkn−1j = N X r,k1,...,kn−1=1 birbrk1bk1k2· · · bkn−1j. Como los ´ındices en el sumatorio son mudos y todos var´ıan entre 1 y N, podemos renombrarlos de forma que nos queda BN+1 ij = N X k1,...,kn=1 bik1· · · bknj. Definici´on 4.2.2. Sea τ > 0y consideramos el sistema de Lotka-Volterra u0= u(b+Au)siendo Auna matriz de la forma A=   −1aij ... aij −1   . Definimos el modelo de Lotka-Volterra con efectos indirectos, o simplemente modelo de efectos indirectos, como el s.d.o: (LV I)τ     u0=u(b+e Au), e A=A+τ ∞ X n=2 (I+A)n,(4.3) siempre que exista la serie. Por analog´ıa, llamaremos modelo de efectos directos al modelo de Lotka-Volterra usual. El par´ametro τnos permite controlar la intensidad de los efectos indirectos. Cuanto mayor sea τ, mayor ser´a la diferencia entre Aye A. Por la Proposici´on 4.2.1 tenemos que la entrada (i, j) de la matriz P∞ n=2(I+A)ncontiene la suma de los productos de todos los coeficientes de los caminos que conectan a la especie icon la j(sin incluir el camino de una especie a s´ı misma), luego e Aes una matriz cuyas entradas son como en (4.1). Ahora probaremos la convergencia de P∞ n=2(I+A)n, pero antes enunciaremos un teorema que nos servir´a en la demostraci´on [6]. Cap´ıtulo 4 — Efectos indirectos 63 Teorema 4.2.3. (Lema de Neumann) Sea P∈CN×Ny supongamos que l´ımk→∞ Pk= 0. Entonces I−Pes regular y (I−P)−1= ∞ X k=0 Pk.(4.4) Una condici´on suficiente para que Pk→0es que kPk<1para alguna norma consistente. En tal caso se tiene I+P+P2+· · · +Pk−(I−P)−1≤kPkk+1 1− kPk.(4.5) Proposici´on 4.2.4. Si Aes una matriz como en (4.2.2) y adem´as es diagonal dominante negativa, entonces e A=A+τ ∞ X n=2 (I+A)n est´a bien definida para todo τ > 0. Demostraci´on. Solo hay que probar que P∞ n=0(I+A)nes convergente. Applicando el Teorema 4.2.3 con P=I+Ay considerando la norma matricial k·kFtenemos que kI+AkF=   0aij ... aij 0   F = m´ax 1≤i≤N N X j=1 |aij|<1 por ser Adiagonal dominante. Luego l´ımn→∞(I+A)n= 0 y en consecuencia P∞ n=0(I+A)nconverge. Adem´as se tiene ∞ X n=0 (I+A)n=−A−1.(4.6) Esta proposici´on no solo prueba que el modelo de efectos indirectos est´a bien definido, tambi´en nos dice el valor de P∞ n=0(I+A)n. Corolario 4.2.5. Si Aes diagonal dominante negativa con diagonal −1entonces e A= (1 −τ)A−2τ I −τA−1.(4.7) Cap´ıtulo 4 — Efectos indirectos 64 Solo nos queda ver que la matriz e Asigue siendo Sw. A diferencia de la convergencia, la estabilidad s´ı depende de τ. Sin embargo, para valores de τsuficientemente peque˜nos se cumple que e Aes una matriz Tipo I: Proposici´on 4.2.6. Sea Auna matriz como en la Definici´on 4.2.2 con A∈NDD. Entonces existe  > 0tal que si τ <  entonces e A∈NDD. Demostraci´on. Denotemos A= (aij), e A= (eaij) y A−1= (a0 ij). Veamos las condiciones que debe cumplir τpara que e Asea diagonal dominante negativa. Diagonal negativa. Sea i∈ {1, . . . , N}. Tenemos que eaii =−(1−τ)−2τ− τa0 ii =−1−(1 + a0 ii)τ. Distinguimos tres casos: 1. Si a0 ii =−1, entonces eaii <0 para cualquier τ > 0. 2. Si a0 ii >−1, entonces 1 + a0 ii >0 y por tanto eaii <0 para todo τ > 0. 3. Si a0 ii <−1, entoces −(1 + a0 ii)τ > 0 y se tiene eaii <0⇔τ < −1 1 + a0 ii . Tomando η > 0 tal que η < m´ın n−1 1+a0 ii :a0 ii <−1otenemos que si τ∈(0, η) entonces la diagonal de e Aes negativa. Diagonal dominante. Definimos la funci´on F: [0,+∞)→Rdada por F(τ) = m´ın 1≤i≤N(|eaii| − N X j=1,j6=i |eaij|)(4.8) = m´ın 1≤i≤N(|−1−(1 + a0 ii)τ| − "(1 −τ) N X j6=i |aij|+τ N X j6=ia0 ij#). (4.9) Fes continua por ser composici´on de funciones continuas. Adem´as, e A∈ NDD si y solo si F(τ)>0. Como A∈NDD,F(0) >0, y por la continuidad de Fsabemos que existe δ > 0 tal que si τ∈(0, δ) entonces e Aes diagonal dominante. Por ´ultimo, tomamos = m´ın{η, δ}y obtenemos el resultado. Cap´ıtulo 4 — Efectos indirectos 65 4.3. Comparaci´on de la estabilidad estructural Para estudiar la influencia de los efectos indirectos en el modelo de Lotka- Volterra, vamos a comparar los ecosistemas analizados en el Cap´ıtulo 3 con sus respectivos modelos de efectos indirectos. A continuaci´on presentamos una nueva forma de comparar la estabilidad estructural de los ecosistemas. Adem´as de medir los conos complementarios de Ay e A, vamos a estudiar tambi´en las intersecciones de dichos conos. Las intersecciones entre dos regiones de biodiversidad de Aye Arepresentan el conjunto de vectores de crecimiento intr´ınseco para los cuales se produce un cambio en la biodiversidad al considerar los efectos indirectos. Por ejemplo, si la regi´on de biodiversidad 3 de Ainterseca la regi´on de biodiversidad 4 de e A, significa que existen ciertos vectores de crecimiento intr´ınseco para los cuales la biodiversidad del sistema aumenta en 1 al considerar los efectos indirectos. En general, cuando regiones de biodiversidad superior de e Aintersecan regiones de biodiversidad inferior de A, se produce una mejora de la biodiversidad, mientras que, por el contrario, cuando las regiones de biodiversidad superior de Aintersecan regiones de biodiversidad inferior de e A, se produce una p´erdida de biodiversidad. Cuanto mayor sea el tama˜no de estas intersecciones, mayor ser´a la p´erdida o recuperaci´on de especies que se produce al considerar el modelo de efectos indirectos. El m´etodo que usaremos para medir las intersecciones es una variante del m´etodo de Montecarlo presentado en el cap´ıtulo anterior. 1. Tomamos Aye Acomo datos de entrada. 2. Generamos Mpuntos aleatorios uniformemente distribuidos en la bola unidad de RN. Llamaremos Ma esta nube de puntos. 3. Definimos la matriz Y= 0 ∈R(N+1)×(N+1). 4. Para cada punto ben Mcalculamos la biodiversidad del sistema de e AyA. 5. Denotando por i, j a las biodiversidades de e AyArespectivamente, sumamos 1 a la componente (i+ 1, j + 1) de Y. 6. Al dividir por el n´umero de puntos generados en la nube, obtenemos en la entrada (i+ 1, j + 1) de Yla proporci´on de la bola unidad que ocupa la intersecci´on de los conos de biodiversidad ide e A, y biodiversidad jde A. Cap´ıtulo 4 — Efectos indirectos 66 Figura 4.3: Ejemplo de una matriz Ypara 5especies y utilizando 1000 puntos. En cada celda se almacena el n´umero de puntos que cae en la intersecci´on de los conos de distintas biodiversidades para Aye A 4.3.1. Resultados En el cap´ıtulo anterior vimos que en el modelo de efectos directos, en general, la cooperaci´on aumenta la estabilidad estructural, mientras que la competici´on la disminuye. El papel que desempe˜na cada uno de estos par´ametros puede cambiar cuando consideramos los efectos indirectos. En la figura 4.4 podemos ver que cuando la competici´on es baja, el volumen del cono sigue siendo creciente respecto a γpyγa. Sin embargo, cuando la competici´on es mayor, el volumen del cono se vuelve decreciente respecto a la cooperaci´on. Adem´as, si fijamos valores de cooperaci´on peque˜nos, la competici´on afecta positivamente a la estabilidad estructural. Estas nuevas relaciones de crecimiento - decrecimiento se mantienen cuando aumentamos δγyδβ(figura 4.5). Respecto al volumen de los conos complementarios, podemos ver en la figura 4.6 que se produce un aumento en el volumen de los conos de biodiversidad mayor, y una disminuci´on en los conos de biodiversidad inferior. Este cambio se hace m´as notable cuando la cooperaci´on es alta. Las intersecciones de los conos complementarios de Aye Arepresentan el conjunto de vectores bpara los cuales se produce un aumento (o p´erdida) de la biodiversidad. En la figura 4.7 observamos que la medida de estas intersecciones es mayor cuando los conos de e Aintersecan conos de biodiversidad inferior de A, indicando que se produce una recuperaci´on de especies. Los conos de Atambi´en intersecan conos de biodiversidad inferior de e A, aunque el volumen es menor. Pese a todos estos cambios en el comportamiento del volumen del cono, podemos afirmar que en general los efectos indirectos mejoran la estabilidad estructural, salvo en algunos casos concretos. Estas excepciones ocurren cuando tanto la competici´on como la cooperaci´on son altas. Cap´ıtulo 4 — Efectos indirectos 67 (a) Ecosistema M SD 029 con βp=βa= 0,05.(b) Ecosistema M SD 029 con βp=βa= 0,2. (c) Ecosistema M PL 042 con βp=βa= 0,02.(d) Ecosistema M PL 042 con βp=βa= 0,1. (e) Ecosistema M AF 001 con βp=βa= 0,005.(f) Ecosistema M AF 001 con βp=βa= 0,05. Figura 4.4: Volumen del cono de la matriz de efectos indirectos para distintos valores de τ. Podemos ver que cuando la competici´on es peque˜na, el volumen del cono sigue siendo creciente respecto a la cooperaci´on. Sin embargo, al aumentar la competici´on, el volumen se vuelve decreciente respecto a γpyγapara valores de τaltos. En todos los c´alculos hemos tomado δγ=δβ= 0. Cap´ıtulo 4 — Efectos indirectos 68 (a) Ecosistema M SD 029 con βp=βa= 0,05.(b) Ecosistema M SD 029 con βp=βa= 0,2. (c) Ecosistema M PL 042 con βp=βa= 0,02.(d) Ecosistema M PL 042 con βp=βa= 0,1. (e) Ecosistema M AF 001 con βp=βa= 0,005.(f) Ecosistema M AF 001 con βp=βa= 0,05. Figura 4.5: Volumen del cono de la matriz de efectos indirectos para distintos valores de τ, tomando δγ=δβ= 1,5. Observamos que el crecimiento y decrecimiento del volumen es similar a la gr´afica anterior. Adem´as, para valores fijados de cooperaci´on y competici´on, el volumen del cono de e Aes mayor cuando aumentamos δγyδβ. Cap´ıtulo 4 — Efectos indirectos 69 (a) Ecosistema M SD 029 con γp=γa= 0,07, βp= βa= 0,3yτ= 5. (b) Ecosistema M SD 029 con γp=γa= 0,1, βp= βa= 0,05 yτ= 5. (c) Ecosistema M PL 042 con γp=γa= 0,05, βp= βa= 0,15 yτ= 4. (d) Ecosistema M PL 042 con γp=γa= 0,06, βp= βa= 0,01 yτ= 4. (e) Ecosistema M AF 001 con γp=γa= 0,01, βp= βa= 0,07 yτ= 2. (f) Ecosistema M AF 001 con γp=γa= 0,04, βp= βa= 0,005 yτ= 2. Figura 4.6: Volumenes de los conos complementarios de las matrices de efectos directos (azul) e indirectos (naranja). Vemos que los conos complementarios de biodiversidad alta de e Atienen mayor volumen que los de A. Este efecto se hace m´as notable cuando la cooperaci´on es mayor que la competici´on. En todos los c´alculos hemos tomado δγ=δβ= 0. Cap´ıtulo 4 — Efectos indirectos 70 (a) Ecosistema M SD 029 con γp=γa= 0,07, βp= βa= 0,3yτ= 5. (b) Ecosistema M SD 029 con γp=γa= 0,1, βp= βa= 0,05 yτ= 5. (c) Ecosistema M PL 042 con γp=γa= 0,05, βp= βa= 0,15 yτ= 4. (d) Ecosistema M PL 042 con γp=γa= 0,06, βp= βa= 0,01 yτ= 4. (e) Ecosistema M AF 001 con γp=γa= 0,01, βp= βa= 0,07 yτ= 2. (f) Ecosistema M AF 001 con γp=γa= 0,04, βp= βa= 0,005 yτ= 2. Figura 4.7: Volumenes de las intersecciones de Aye A. Podemos ver que las regiones de biodiversidad superior de e Aintersecan las regiones de biodiversidad inferior de A, indicando que se produce un aumento en la biodiversidad para los vectores de crecimiento intr´ınseco que est´en en dichas intersecciones. En todos los c´alculos hemos tomado δγ=δβ= 0. Cap´ıtulo 5 Conclusiones y problemas abiertos En esta memoria hemos estudiado las propiedades del modelo de Lotka-Volterra mutualista. A trav´es de la Teor´ıa LCP hemos dado condiciones suficientes para la existencia de un punto estacionario globalmente estable y no negativo. Uno de los resusltados m´as importantes es el Teorema 2.2.3, el cual permite aplicar la Teor´ıa LCP al sistema Lotka-Volterra demostrando que si A∈Sw, entonces existe un ´unico punto estacionario factible para (L-V). Adem´as, la teor´ıa LCP nos permite estudiar la biodiversidad del ecosistema a trav´es de la geometr´ıa. La matriz Adivide RNen una serie de conos complementarios que forman una partici´on, y en funci´on de a qu´e cono pertenezca b, el sistema tendr´a m´as o menos biodiversidad. Esto nos proporciona un nuevo m´etodo para calcular la estabilidad (llamada estabilidad estructural), ya que cuanto mayor sea el cono mayor ser´a la perturbaci´on que puede soportar el vector bantes de que el sistema pierda biodiversidad. El an´alisis de los ecosistemas reales que hemos realizado ha demostrado que los ecosistemas con una cooperaci´on alta son m´as estables que los ecosistemas con menor cooperaci´on. Adem´as, hemos podido comprobar que cuanto mayor sea el n´umero de especies de un ecosistema, menor es su estabilidad estructural. Por ´ultimo, hemos planteado un nuevo modelo de Lotka-Volterra que incluye las relaciones indirectas entre las especies, y hemos analizado la estabilidad estructural de ´este con respecto al Lotka-Volterra original. Observamos que el papel que juega la cooperaci´on y la competici´on en el volumen del cono se invierte cuando incluimos los efectos indirectos, siendo creciente respecto a la competici´on y decreciente respecto a la cooperaci´on. A pesar de todo, considerar los efectos indirectos ha 71 Programas de MatLab 78 genPuntos 1function [M] = genPuntos(dim,npunt) 2% genPuntos(dim,npunt) 3% genera npunt puntos en la bola unidad en dimension dim con una ... distribucion 4% uniforme. Los puntos generados se guardan en forma de matriz de ... tamano 5% dim x npunt, donde cada columna representa las coordenadas de ... un punto. 6rng('shuffle') 7 8M = zeros(dim,npunt); 9 10 for i = 1:npunt 11 X = normrnd(0,1,dim,1); 12 Y = exprnd(1); 13 M(:,i) = X./sqrt(Y+norm(X,2)ˆ2); 14 end Programas de MatLab 79 volIntersecciones 1function [inter] = volIntersecciones(A,B,puntos) 2% volIntersecciones(A,B,puntos) calcula el el volumen de las ... intersecciones 3% de los conos complementarios de las matrices A y B. El tercer ... argumento 4% es el numero de puntos que se utiliza en el metodo de Montecarlo. 5if ¬isequal(size(A),size(B)) 6error('Las matrices no son del mismo tamano'); 7end 8[N,¬] = size(A); 9inter = zeros(N+1); 10 datosAnt = csvread('intersecciones.csv'); 11 M = genPuntos(N,puntos); 12 for i = 1:puntos 13 b = M(:,i); 14 [¬,za] = LCPSolve(-A,-b); 15 [¬,zb] = LCPSolve(-B,-b); 16 fila = sum(za>0)+1; 17 colu = sum(zb>0)+1; 18 inter(fila,colu) = inter(fila,colu) + 1; 19 end 20 inter = datosAnt+inter; 21 writematrix(inter,'intersecciones.csv'); Programas de MatLab 80 resetPart 1function [] = resetPart(dim) 2% resetPart(dim) borra los datos guardados del calculo del ... volumen de las 3% intersecciones, y en su lugar escribe una matriz nula de tamano ... dim x dim 4A = zeros(dim+1); 5writematrix(A,'intersecciones.csv'); Programas de MatLab 81 esNDD 1function [b] = esNDD(A) 2% esNDD(A) comprueba si la matriz A es diagonal dominante negativa. 3d = diag(A); 4neg = all(d<0); 5dom = all(sum(abs(A-diag(d)),2)<abs(d)); 6b = neg & dom; Programas de MatLab 82 VolConoRapido2 1function [V] = VolConoRapido2(A) 2% VolConoRapido(A) es un metodo alternativo al de Montecarlo que 3% calcula el volumen del cono de A mas rapido. 4% La matriz A debe ser cuadrada y diagonal dominante negativa. 5%% Inicializacion y parametros 6npunt = 100; % numero de puntos por ortante. 7[N,col] = size(A); % numero de especies. 8volumen = 0; 9muestra = round(Nˆ2/2); 10 11 if N6=col % comprobar que la matriz es cuadrada 12 error('La matriz no es cuadrada'); 13 end 14 15 16 rng('shuffle')% generador de numeros aleatorios 17 %% Algoritmo 18 if esNDD(A) 19 for k = N:-1:1 20 21 ncomb = nchoosek(N,k); % numero de ortantes de tamano k 22 if ncomb ≤muestra 23 comb = nchoosek(1:N,k); 24 nconos = ncomb; 25 aumento = 1; 26 for i = 1:nconos 27 v = indiceBinario(comb(i,:),N); % indice del ortante 28 for n = 1:npunt 29 X = normrnd(0,1,N,1); 30 Y = exprnd(1); 31 p = abs(X./sqrt(Y+norm(X,2)ˆ2)); 32 p=p.*((-1)*ones(N,1)).ˆ(v'+1); 33 [¬,bio] = LCPSolve(-A,-p); 34 if all(bio >0) 35 volumen = volumen+aumento; 36 end 37 end 38 end 39 else 40 nconos = muestra; 41 aumento = ncomb/nconos; 42 for i = 1:nconos 43 v = generaCombinacion(N,k); % indice del ortante 44 for n = 1:npunt 45 X = normrnd(0,1,N,1); Programas de MatLab 83 46 Y = exprnd(1); 47 p = abs(X./sqrt(Y+norm(X,2)ˆ2)); 48 p=p.*((-1)*ones(N,1)).ˆ(v'+1); 49 [¬,bio] = LCPSolve(-A,-p); 50 if all(bio >0) 51 volumen = volumen+aumento; 52 end 53 end 54 55 end 56 57 end 58 % disp(['Completado: ', num2str((1-k/N)*100,' %.0f'), ' %']); 59 end 60 61 % Volumen final 62 V = volumen/(npunt*2ˆN); 63 else 64 V = 0; 65 end 66 %% Funciones auxiliares 67 function [vv] = indiceBinario(uu,NN) 68 vv = zeros(1,NN); 69 l = length(uu); 70 for j = 1:l 71 vv(uu(j)) = 1; 72 end 73 end 74 function [vv] = generaCombinacion(NN,kk) 75 while 1>0 76 pp = binornd(1,kk/NN,1,NN); 77 if sum(pp) == kk 78 vv = pp; 79 break 80 end 81 end 82 83 end 84 end Programas de MatLab 84 Graficas2DInd 1function [] = Graficas2DInd(ady,vectorT,gmax,leer,ndiv,Bcomple,beta) 2% Graficas2DInd(ady,vectorT,gmax,leer,ndiv,Bcomple,beta) muestra ... la grafica 3% del volumen del cono en funcion de la cooperacion. 4% - ady: matriz de conexiones del ecosistema mutualista. 5% - vectorT: vector con los distintos valores de T para los ... cuales se 6% desea representar las graficas. P. ej vectorT = [0, ... 0.5, 0.8, 1.2] 7% - gmax: valor maximo de la cooperacion 8% - leer: es posible ejecutar este programa en varias sesiones. leer 9% indica el numero de veces que se ha ejecutado. ... Escribir leer = 0 10 % si no se desea ejecutar en varias tandas. 11 % - ndiv: numero de divisiones en el eje de abscisas (cooperacion) 12 % - Bcomple: valor booleano que determina como se construyen las ... matrices 13 % de competicion. si Bcomple = 0, solo estan en ... competicion 14 % aquellas especies que tengan un cooeperador en ... comun. Si 15 % Bcomple = 1, las submatrices B1 y B2 son completas ... (no tienen 16 % ceros) 17 % - beta: valor de competicion. 18 19 rng('shuffle') 20 21 ntes = length(vectorT); 22 gamma = linspace(0,gmax,ndiv); 23 [P,A] = size(ady); 24 N = P+A; 25 ∆ = 0; 26 27 if leer>0 28 ZZZ = csvread('2dcono.csv'); 29 [mz,nz] = size(ZZZ); 30 ZZ = NaN*ones(ntes,ndiv); 31 ZZ(1:mz,1:nz) = ZZZ; 32 end 33 Z = NaN*ones(ntes,ndiv); 34 35 for k = 1:ntes 36 T = vectorT(k); 37 for i = 1:ndiv Programas de MatLab 85 38 g = gamma(i); 39 M = generaMatriz(ady,Bcomple,g,g,1.5,beta,beta,∆); 40 Mind = (1-T)*M-2*T*eye(N)-T*Mˆ(-1); 41 if esNDD(Mind) 42 Z(k,i) = VolConoRapido2(Mind); 43 end 44 porcentaje = (i + (k-1)*ndiv)/(ntes*ndiv); 45 disp(['Completado: ', num2str(porcentaje*100,' %.1f'), ' %']); 46 end 47 48 if leer>0 49 Z = (ZZ*leer+Z)/(leer+1); 50 end 51 etiqueta = ['T = ', num2str(T,' %.2f')]; 52 plot(gamma,Z(k,:),'DisplayName',etiqueta); 53 hold on 54 xlabel('\gamma p = \gamma a'); 55 ylabel('Vˆ*(C A)'); 56 end 57 58 csvwrite('2dcono.csv',Z); 59 hold off 60 legend show 61 legend('location','northwest'), 62 sound(0.001*ones(100000,1),44100); 63 [y,fs] = audioread('sonidito.wav'); 64 soundsc(y,fs); Programas de MatLab 86 Graficas2D 1function [] = Graficas2D(ady,vectorB,gmax,leer,ndiv,Bcomple) 2% Graficas2D(ady,vectorB,gmax,leer,ndiv,Bcomple)) muestra la grafica 3% del volumen del cono en funcion de la cooperacion. 4% - ady: matriz de conexiones del ecosistema mutualista. 5% - vectorB: vector con los distintos valores de competicion ... para los 6% cuales se desea representar las graficas. P. ej ... vectorB = 7% [0,0.5, 0.8, 1.2] 8% - gmax: valor maximo de la cooperacion 9% - leer: es posible ejecutar este programa en varias sesiones. leer 10 % indica el numero de veces que se ha ejecutado. ... Escribir leer = 0 11 % si no se desea ejecutar en varias tandas. 12 % - ndiv: numero de divisiones en el eje de abscisas (cooperacion) 13 % - Bcomple: valor booleano que determina como se construyen las ... matrices 14 % de competicion. Si Bcomple = 0, solo estan en ... competicion 15 % aquellas especies que tengan un cooperador en comun. Si 16 % Bcomple = 1, Las submatrices B1 y B2 son completas ... (no tienen 17 % ceros) 18 19 rng('shuffle') 20 21 nbetas = length(vectorB); 22 gamma = linspace(0,gmax,ndiv); 23 24 if leer>0 25 ZZZ = csvread('2dcono.csv'); 26 [mz,nz] = size(ZZZ); 27 ZZ = NaN*ones(nbetas,ndiv); 28 ZZ(1:mz,1:nz) = ZZZ; 29 end 30 Z = NaN*ones(nbetas,ndiv); 31 32 for k = 1:nbetas 33 b = vectorB(k); 34 for i = 1:ndiv 35 g = gamma(i); 36 M = generaMatriz(ady,Bcomple,g,g,1.5,b,b,1.5); 37 if esNDD(M) 38 Z(k,i) = VolConoRapido2(M); 39 end Programas de MatLab 87 40 porcentaje = (i + (k-1)*ndiv)/(nbetas*ndiv); 41 disp(['Completado: ', num2str(porcentaje*100,' %.1f'), ' %']); 42 end 43 44 if leer>0 45 Z = (ZZ*leer+Z)/(leer+1); 46 end 47 etiqueta = ['\beta p = \beta a = ', num2str(b,' %.2f')]; 48 plot(gamma,Z(k,:),'DisplayName',etiqueta); 49 hold on 50 xlabel('\gamma p = \gamma a'); 51 ylabel('Vˆ*(C A)'); 52 end 53 54 csvwrite('2dcono.csv',Z); 55 hold off 56 legend show 57 legend('location','northwest'), 58 sound(0.001*ones(100000,1),44100); 59 [y,fs] = audioread('sonidito.wav'); 60 soundsc(y,fs);