Tema 7

Aplicaciones de los S.E.D.O. 7.1

Introducci´ on

Nota: APUNTES INCOMPLETOS Estudiaremos en este Tema algunos modelos de inter´es en las Ciencias Naturales que utilizan para su modelizaci´on sistemas de ecuaciones diferenciales ordinarias. Para su resoluci´on y/o an´alisis precisaremos las t´ecnicas aprendidas en el tema anterior, as´ı como una introducci´on al An´alisis Cualitativo de los sistemas (aut´onomos). Se incluir´an tambi´en brevemente algunos conceptos acerca del an´alisis num´erico de los sistemas.

7.2

An´ alisis del Plano de fases

Como ya se ha explicado en el tema anterior, una soluci´on de un sistema de dos ecuaciones diferenciales de primer orden: x = x(t), y = y(t), puede ser vista como las ecuaciones param´etricas de una curva en el plano determinado por las variables dependientes, x e y, llamado Plano de Fases del sistema de ecuaciones. Dicha curva se llama entonces ´orbita soluci´on. Las condiciones iniciales: x(t0 ) = x0 e y(t0 ) = y0 se reducen entonces a fijar un punto concreto de dicho plano que determina una u ´nica ´orbita de entre todas 1 las posibles . Para los sistemas aut´onomos es posible determinar directamente la ecuaci´on diferencial cuya integraci´on da como resultado las ´orbitas. Por otro lado, siempre es posible en este tipo de sistemas analizar en primer lugar las soluciones estacionarias, como veremos a continuaci´on. 1

Damos por supuesto que se verifican las hip´ otesis del Teorema de Picard para sistemas.

1

2

7.2.1

´ CIENCIAS AMBIENTALES / MODELOS MATEMATICOS / TEMA 7

Soluciones estacionarias

Consideremos un s.e.d.o. aut´onomo: dx dt dy dt

= f (x, y) = g(x, y)

)

Las soluciones estacionarias (o puntos singulares) del sistema son las soluciones del sistema algebraico: ) 0 = f (x, y) 0 = g(x, y) que supondremos puntos aislados. En el plano de fases, se trata evidentemente, de un punto soluci´on. Razonando de manera an´aloga a lo estudiado en los temas 3 y 4 acerca del an´alisis cualitativo de las soluciones es posible intuir el comportamiento de las mismas sin necesidad de integrar el sistema de ecuaciones, como veremos en los ejemplos siguientes. Ejemplo 1: Analicemos el sistema: dx dt dy dt

=y−2 =1−x

)

El sistema posee obviamente s´olo una soluci´on estacionaria, el punto P ≡ (1, 2). Dicho punto es adem´as la intersecci´on de las rectas: y = 2, donde se anula la derivada de x, y x = 1 donde lo hace dy ı dividir el plano de fases xy en cuatro regiones separadas por ambas rectas. dt . Podemos as´ y 3

IV

I

III

II

2

1

x 1

-1

2

3

-1

Figura 7.1: Regiones en las que se divide el plano de fases del sistema. dy , por tanto la direcci´on del vector tangente Para todos los puntos de la recta x = 1 se anula dx dy dx ~t = ( , ) a las ´orbitas soluci´on, en dichos puntos, ser´a horizontal. Evidentemente, apuntar´a dt dt en el sentido de las x positivas siempre y cuando y > 2, mientras que lo har´a hacia las x negativas si y < 2. Un razonamiento an´alogo para la recta y = 2 nos conduce al resultado mostrado en la Figura 7.2 (izquierda).

3

´ CIENCIAS AMBIENTALES / MODELOS MATEMATICOS / TEMA 7

En el interior de las regiones I, II, III y IV, tendremos que el vector tangente ~t “apuntar´ a” en las direcciones y sentidos indicadas en la Figura 7.2 (derecha). y

y

3

3

2

2

1

1

x 1

-1

2

3

x 1

-1

-1

2

3

-1

Figura 7.2: Direcciones de las tangentes en cada uno de los sectores (la inclinaci´on y magnitud del vector tangente no se ha determinado con precisi´ on, se trata por tanto de gr´ aficas aproximadas, que s´ olo indican aspectos cualitativos). En este ejemplo concreto, es f´acil resolver tanto el sistema directamente (es un sistema lineal con coeficientes constantes, no homog´eneo), como la ecuaci´on diferencial de las ´orbitas del sistema: dx − dy = ⇒ (x − 1)2 + (y − 2)2 = 2C y−2 x−1 es decir, una familia de circunferencias centrada en el punto estacionario (1, 2)

Ejemplo 2: Analicemos ahora el plano de fases del sistema: dx dt dy dt

= y − x2 = x − y2

)

restringi´endonos al primer cuadrante, x > 0 e y > 0. ............

En el siguiente apartado estudiaremos tambi´en c´omo se puede aproximar el sistema en los alrededores de una soluci´on estacionaria por medio de un sistema lineal con coeficientes constantes asociado, ello permitir´a, adem´as, clasificar la estabilidad de las soluciones estacionarias.

7.2.2

Estabilidad Lineal. Clasificaci´ on de los puntos estacionarios

Dada una soluci´on estacionaria de un sistema, veremos a continuaci´ on c´omo estudiar el comportamiento de las soluciones cercanas a la misma. El criterio con el que decidiremos si un punto estacionario es estable o no (estabilidad lineal) consistir´a simplemente en determinar lo que ocurre cuando t → ∞ en el sistema linealizado en el punto. Si las variables tiende a infinito tendremos un punto de equilibrio inestable, si tienden al punto de equilibrio, o al menos se mantienen acotadas en sus alrededores, tendremos un punto de equilibrio estable.

4

´ CIENCIAS AMBIENTALES / MODELOS MATEMATICOS / TEMA 7

En las cercan´ıas de un punto estacionario (xs , ys ), podemos aproximar las funciones f (x, y) y g(x, y) por sus correspondientes polinomios de Taylor de primer orden en el punto (xs , ys ), de manera que: ) ∂f ∂f dx ' f (x , y ) + (x , y ) (x − x ) + (x , y ) (y − y ) s s s s s s s s dt ∂x ∂y dy ∂g ∂g ' g(x , y ) + (x , y ) (x − x ) + s s s s s dt ∂x ∂y (xs , ys ) (y − ys ) que teniendo en cuenta que (xs , ys ) es singular, y redefiniendo las variables: x1 = x − xs , x2 = y − ys , puede escribirse como: ) dx1 = ax + bx 1 2 dt dx2 = cx + dx 1 2 dt es decir, un sistema lineal homog´eneo con coeficientes constantes, que sabemos resolver. Concluimos por tanto que en las cercan´ıas de un punto estacionario siempre podemos aproximar el sistema que estamos estudiando por uno lineal con coeficientes constantes homog´eneo. En el siguiente ejemplo se muestra esta situaci´on. Ejemplo : Estudiaremos a continuaci´on el siguiente sistema aut´onomo en el que estamos interesados en la regi´on x > 0 e y > 0: dx dt dy dt

= x(1 − y) = y 2 (x2 − 1)

)

Calcularemos en primer lugar los puntos estacionarios del sistema. ) f (x, y) = x(1 − y) = 0 ⇒ P1 = (0, 0), P2 = (1, 1), P3 = (−1, 1) g(x, y) = y 2 (x2 − 1) = 0 Nos interesa por tanto el punto P2 = (1, 1). Podemos tambi´en calcular la ecuaci´on de las ´orbitas soluci´on: dx dy 1 1 = 2 2 ⇒ x2 − ln |x| + + ln |y| = C x(1 − y) y (x − 1) 2 y En la figura 7.3 (izq.) pueden observarse algunas de estas ´orbitas. Calcularemos a continuaci´on la linealizaci´on del sistema en los alrededores de P2 . Para ello calcularemos las derivadas parciales: ∂f (x, y) =1−y ∂x

,

∂f (x, y) = −x ∂y

∂f (x, y) (1, 1) = 0 ∂x

,

∂f (x, y) (1, 1) = −1 ∂y

∂g(x, y) = 2xy 2 ∂x

, ,

∂g(x, y) (1, 1) = 2 ∂x

Tendremos as´ı, en P2 = (1, 1), el sistema lineal: dx dt dy dt

,

= 0(x − 1) − 1(y − 1) = 2(x − 1) + 0(y − 1)

)

∂g(x, y) = 2y(x2 − 1) ∂y ,

∂g(x, y) (1, 1) = 0 ∂y

5

´ CIENCIAS AMBIENTALES / MODELOS MATEMATICOS / TEMA 7

y 3

1.4 2.5

1.2

2

1

1.5

0.8 0.6

1

0.4 0.5

0.2 0.5

1

1.5

2

2.5

3

1

0.5

1.5

2

x

Figura 7.3: a) Gr´aficas de algunas ´orbitas del sistema. b) Gr´afica de la ´orbita exacta con condiciones iniciales x(0) = y(0) = 0.8 y de la ´ orbita del sistema lineal asociado (elipse) con las mismas condiciones iniciales. Si definimos variables nuevas x ¯ = (x − 1) e y¯ = (y − 1), tendremos: d¯ x dt d¯ y dt

= =

−¯ y 2¯ x

)

cuya resoluci´on es f´acil: x ¯ = C1 cos x = 1 + C1 cos

√

√ 2t + C2 sen 2t

,

y¯ =

√

√

√ 2t + C2 sen 2t

,

y =1+

2C1 sen √

√

2t −

√

√ 2C2 cos 2t

√ √ √ 2C1 sen 2t − 2C2 cos 2t

Consideramos ahora las condiciones iniciales (cercanas al punto P2 ): x(0) = y(0) = 0.8. 0.8 = 1 + C1

x = 1 − 0.2 cos

√

,

0.8 = 1 −

√

2C2 ⇒ C1 = −0.2

√ 0.2 2t + √ sen 2t , 2

,

0.2 C2 = √ 2

√ √ √ y = 1 − 0.2 2 sen 2t − 0.2 cos 2t

Es f´acil deducir que esta soluci´on tiene como ´orbita a una elipse centrada en el punto (1, 1) √ √ y de semi-ejes 0.06 y 0.12. Este tipo de puntos estacionarios se denomina centros. dx −dy 1 = ⇒ (x − 1)2 + (y − 1)2 = C (y − 1) 2(x − 1) 2 que teniendo en cuenta que debe de pasar por el punto (0.8, 0.8): 1 (x − 1)2 (y − 1)2 0.22 + 0.22 = 0.06 ⇒ + =1 2 0.06 0.12 En la figura 7.3 (derecha) est´a representada la ´orbita (exacta) correspondiente a las condiciones iniciales dadas y la elipse que acabamos de calcular.

6

´ CIENCIAS AMBIENTALES / MODELOS MATEMATICOS / TEMA 7

Clasificaci´ on de los Puntos Estacionarios Casos cr´ıticos: la clasificaci´on afecta s´olo a la aproximaci´ on lineal, la estabilidad o inestabilidad de ´esta no queda garantizada para el sistema completo. Casos principales o gen´ericos: No est´an afectados por peque˜ nos errores de medida o de redondeo, son todos ellos casos no cr´ıticos (Hay tambi´en casos no principales y no estrictamente cr´ıticos). • Dos ra´ıces reales distintas λ1 y λ2 : – Del mismo signo: a) ambas positivas: Equilibrio inestable (principal). Nodo. b) ambas negativas: Equilibrio estable (principal). Nodo. – De diferente signo: Equilibrio inestable (principal). Punto de Silla. – Una ra´ız nula: a) La otra positiva: equilibrio inestable. b) La otra negativa: equilibrio estable (cr´ıtico). • Una ra´ız real doble: – Diferente de cero: a) Positiva: equilibrio inestable. Nodo tipo 2. b) Negativa: equilibrio estable. Nodo tipo 2. – Ra´ız nula: equilibrio inestable (cr´ıtico). • Dos ra´ıces complejas – Parte real no nula: a) positiva: equilibrio inestable (principal). Espiral. b) negativa: equilibrio estable (principal). Espiral. – Imaginarias puras: equilibrio estable (cr´ıtico). Centro.

7.3

Modelos de Din´ amica de poblaciones con m´ as de una especie

Analizaremos en esta secci´on diferentes modelos de Din´amica de poblaciones con m´as de una especie. Existen tres casos particularmente relevantes, se trata de modelos tipo depredador-presa, donde una de las dos especies se beneficia del contacto con la otra,

´ CIENCIAS AMBIENTALES / MODELOS MATEMATICOS / TEMA 7

7

que resulta perjudicada, los del tipo competencia mutua, en el que ambas resultan perjudicadas, y finalmente los de tipo simbiosis, donde las dos especies resultan favorecidas por la interacci´on entre ellas. El modelo quiz´as m´as conocido y pionero hist´oricamente es el de Lotka-Volterra, con el que comenzaremos este estudio.

7.3.1

Modelo de Lotka-Volterra

Breve historia del modelo: . . . . . . . . . x= poblaci´on presa. y= poblaci´on depredadora. dx dt dy dt

= ax − bxy = −cy + dxy

)

Explicaci´on: en ausencia de depredador: presa malthusiana creciente; en ausencia de presas: depredador malthusiano decreciente. xy da cuenta del n´ umero de contactos entre los individuos de las dos especies. Soluciones estacionarias: claramente dos, una trivial: x = y = 0 y la otra: x = dc , y = ab .

Figura 7.4: Al ser un sistema aut´onomo, es posible analizar directamente la ecuaci´on de las ´orbitas soluci´on: dx dy = x(a − by) y(dx − c) que adem´as es de variables separables, su integraci´ on nos lleva a: d x + b y − ln(xc y a ) = C Mientras que la soluci´on estacionaria x = y = 0 es claramente una soluci´on singular de la ecuaci´on de las ´orbitas, la otra soluci´on estacionaria puede verse como soluci´on particular (extendiendo l´ogicamente el dominio de los par´ametros) para el valor: C = c + a − c ln dc − a ln ab . Analicemos la estabilidad (lineal) de las soluciones estacionarias:

8

´ CIENCIAS AMBIENTALES / MODELOS MATEMATICOS / TEMA 7

1. Soluci´ on: x = y = 0 El sistema se linealiza a: dx dt dy dt

= ax = −cy

)

es decir, las ecuaciones malthusianas correspondientes, y as´ı las soluciones ser´an: x = C1 eat

y = C2 e−ct

;

Se trata claramente de un punto de silla inestable y en consecuencia, el comportamiento de las soluciones en las cercan´ıas del punto (0, 0) es el mostrado en las figuras 7.5 y 7.6.

x

y

t

t

Figura 7.5: Gr´aficas de algunas soluciones en las cercan´ıas de x = 0, y = 0.

y

x

Figura 7.6: Gr´aficas de algunas ´orbitas en las cercan´ıas de (0, 0). 2. Soluci´ on: x = dc , y =

a b

El sistema linealizado resulta ahora: d¯ x dt d¯ y dt

donde x ¯ = x − dc , y¯ = y − ab , α = Las soluciones son :

bc d

= −α¯ y = βx ¯

>0yβ=

√ √ x ¯ = C1 cos act + C2 sen act;

ad b

)

> 0.

√ √ √ ¢ d a¡ y¯ = √ C1 sen act − C2 cos act b c

9

´ CIENCIAS AMBIENTALES / MODELOS MATEMATICOS / TEMA 7

y el punto resulta ser entonces un “centro”. (Ver figuras 7.7 y 7.8 izq.). Aunque no es dif´ıcil observar directamente en las expresiones anteriores que las ´orbitas de las soluciones son elipses, procederemos a calcularlas directamente a partir del sistema linealizado: d¯ x d¯ y β 2 α 2 = ⇒ βx ¯ d¯ x + α¯ y d¯ y=0 ⇒ x ¯ + y¯ = C −α¯ y βx ¯ 2 2 Deshaciendo los cmabios y definiendo los semiejes horizontal y vertical de la elipse, A y B, de la forma: r r 2bC 2dC A= , B= ad bc tendremos finalmente la ecuaci´on: (x − dc )2 (y − ab )2 + =1 A2 B2 Si se conocen unas condiciones iniciales x(0) = x0 , e y(0) = y0 , entonces la constante de integraci´on C quedar´a determinada: C=

(x0 − dc )2 2b ad

x

+

(y0 − dc )2 2d bc

y

t

t

Figura 7.7: Gr´aficas de algunas soluciones del sistema linealizado en las cercan´ıas de x = dc ,

y = ab .

Finalmente, la resoluci´on num´erica del sistema para diferentes puntos iniciales puede verse representada en la Figura 7.8 (derecha). Estimaci´ on de los par´ ametros En general no es f´acil deducir de los datos experimentales los valores de los diferentes par´ametros que intervienen en el Modelo de Lotka-Volterra. En el caso en el que nos encontremos cerca del punto de equilibrio antes deducido, y dando por buena la aproximaci´ on que la elipse nos proporciona, podemos razonar del siguiente modo:

10

´ CIENCIAS AMBIENTALES / MODELOS MATEMATICOS / TEMA 7

y

y

x

x

Figura 7.8: a) Gr´aficas de algunas ´orbitas (elipses) en las cercan´ıas de ( dc , ad ). b) Gr´aficas de algunas ´ orbitas obtenidas por resoluci´ on num´erica del sistema completo.

El periodo del ciclo de las poblaciones T es T = √2πac . Las poblaciones medias las tomaremos iguales a las de los valores de equilibrio: xe = dc , ye = ab . Las poblaciones m´aximas y m´ınimas del depredador y de la presa nos vendr´ an dadas l´ogicamente por los semi-ejes de las elipses. xM = xe + A, xm = xe − A, yM = ye + B, ym = ye − B. De esta manera, a partir de las constantes del modelo a, b, c, d, y de los valores iniciales (x0 , y0 ), podemos determinar (en la aproximaci´ on lineal) el periodo y los valores m´aximos y m´ınimos previstos para las dos especies: a 2π c , ye = T =√ , xe = d b ac r r r r c c a a 2bC 2bC 2dC 2dC xM = + , xm = − , yM = + , ym = − d ad d ad b bc b bc Siendo C la constante de la ´orbita: (x0 − dc )2 (y0 − dc )2 C= + 2b 2d ad

bc

Es posible determinar tambi´en la soluci´on del problema inverso, es decir: conocidos de manera aproximada los valores del periodo, xM , xm , yM e ym , calcular las constantes originales del modelo: a, b, c y d: En primer lugar, los valores de xe e ye ser´an , obviamente: xM + xm yM + ym xe ' , ye ' 2 2 mientras que los semiejes de la elipse y la raz´on entre ambos, R, ser´an: yM − ym A xM − xm xM − xm , B= , R= = A= 2 2 B yM − ym Con estos ingredientes es posible despejar de las identidades anteriores para obtener: 2πR 2πxe 2π 2πRye , b= , c= , d= a= xe T xe T ye RT ye RT

,

´ CIENCIAS AMBIENTALES / MODELOS MATEMATICOS / TEMA 7

Mejoras del modelo L-V Modelo de Lotka-Volterra con comportamiento log´ıstico en la especie “presa”. .........

7.3.2

Modelos de competencia dx dt dy dt

= ax − bxy = cy − dxy

)

Figura 7.9: dx dt dy dt

= ax − bxy − γx2 = cy − dxy − δy 2

)

Casos: . . . . . . . . .

7.3.3

Modelos de simbiosis dx dt dy dt

= ax + bxy − γx2 = cy + dxy − δy 2

Casos: . . . . . . . . .

7.4

Otros Modelos

.........

7.5

M´ etodo de Euler para S.E.D.O.

.........

)

11