Unha introdución ás ecuacións en diferenzas

Traballo Fin de Grao Unha introdución ás ecuacións en diferenzas Martín Lema Pailos 2019/2020 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Unha introdución ás ecuacións en diferenzas Martín Lema Pailos 2019/2020 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Traballo proposto Área de Coñecemento: Análise matemática. Título: Unha introdución ás ecuacións en diferenzas. Breve descrición do contido As ecuacións en diferenzas resultan de gran interese, non só como ferramenta para a aproximación das solucións de ecuacións diferenciais, senón tamén porque aparecen de forma natural en diversos fenómenos económicos ou biolóxicos, entre outros. Neste traballo introducirase o concepto de ecuación en diferenzas, para estudar posteriormente os resultados que garanten a existencia de solución, así como os métodos explícitos de resolución en certos casos. Comprobarase deste xeito a clara analoxía existente entre este tipo de ecuacións e as ecuacións diferenciais ordinarias. Recomendacións Ter cursadas as materias de Introdución ás Ecuacións Diferenciais Ordinarias e Ecuacións Diferenciais Ordinarias. iii Índice xeral Resumo vii Introdución ix 1. Cálculo en diferenzas 1 1.1. Operadordiferenza................................ 1 1.2. Sumaindefinida.................................. 5 2. Ecuacións lineais en diferenzas 9 2.1. Resultados xerais de ecuacións en diferenzas lineais . . . . . . . . . . . . . . 10 2.2. Resolución de ecuacións lineais en diferenzas . . . . . . . . . . . . . . . . . . 15 2.2.1. Resolución de ecuacións lineais de primeira orde . . . . . . . . . . . . 15 2.2.2. Ecuacións lineais con coeficientes constantes . . . . . . . . . . . . . . 18 2.2.3. Ecuacións lineais con coeficientes non constantes . . . . . . . . . . . 24 3. Sistemas de ecuacións lineais en diferenzas 29 3.1. Sistemas lineais autónomos con coeficientes constantes . . . . . . . . . . . . 31 3.1.1. Resolución de sistemas lineais autónomos . . . . . . . . . . . . . . . 32 3.2. Sistemas lineais non autónomos . . . . . . . . . . . . . . . . . . . . . . . . . 36 4. Teoría da estabilidade 39 4.1. Métododaescaleira ............................... 41 4.2. FuncióndeLyapunov............................... 44 4.3. Omodeloloxístico ................................ 48 4.3.1. Crecemento estable (0< r ≤2) ..................... 52 4.3.2. Crecemento cíclico (2< r < 2,57).................... 54 4.3.3. Crecemento caótico (2,57 < r ≤3) ................... 56 Bibliografía 59 v Resumo Neste traballo estudaremos as ecuacións en diferenzas. Introduciremos o cálculo en diferenzas para así poder traballar coas ecuacións en diferenzas, vendo métodos explícitos para resolvelas en certos casos e resultados que garanten a existencia de solución. Tras isto, veremos os sistemas de ecuacións lineais en diferenzas de forma breve, comprobando a analoxía existente cos sistemas de ecuacións lineais diferenciais á hora da busca de solucións. Por último, pasaremos a estudar o comportamento asintótico das solucións dos sistemas e das ecuacións non lineais en diferenzas, que é o que se coñece como teoría da estabilidade. Abstract In this work we will study the difference equations. We will introduce the calculation in differences in order to work with equations in differences, seeing explicit methods to solve them in certain cases and results that guarantee the existence of a solution. After this, we will see the systems of linear equations in differences in a brief way, checking the existing analogy with the systems of differential linear equations when looking for solutions. Finally, we will study the asymptotic behavior of the solutions of the systems and of the nonlinear equations and differences, which is known as the stability theory. vii 2CAPÍTULO 1. CÁLCULO EN DIFERENZAS Pódense considerar diferenzas de orde maior. En particular, defínese a diferenza de orde 2 coma ∆2y(t) = ∆(∆y(t)), e, en xeral, a diferenza de orde n coma ∆ny(t) = ∆(∆n−1y(t)). Ademais do operador diferenza, tamén se acostuma a usar o operador shift. Definición 1.2. Dada unha función ycon valores reais, definimos o operador shift, E, como Ey(t) = y(t+ 1). Vemos que ∆ = E−I, onde I denota o operador identidade definido por I y(t) = y(t). Podemos expresar entón a diferenza de orde n do seguinte xeito ∆ny(t)=(E−I)ny(t) = n X k=0 n k(−I)kEn−ky(t) = n X k=0 n k(−1)ky(t+n−k). Procedemos agora á introdución dalgunhas propiedades importantes do operador ∆. Teorema 1.3. Sexan y e z dúas funcións. Tense que 1. ∆m(∆ny(t)) = ∆m+ny(t)para todo m, n ∈N. 2. ∆(y(t) + z(t)) = ∆y(t)+∆z(t). 3. ∆(c y(t)) = c∆y(t)para calquera constante c. 4. ∆(y(t)z(t)) = y(t) ∆z(t) + E z(t) ∆ y(t). 5. ∆y(t) z(t)=z(t) ∆ y(t)−y(t) ∆ z(t) z(t)E z(t). Demostración. 2 Por definición, ∆[y(t) + z(t)] = [y(t+ 1) + z(t+ 1)] −[y(t) + z(t)] = [y(t+ 1) −y(t)] + [z(t+ 1) −z(t)] = ∆y(t)+∆z(t). 1.1. OPERADOR DIFERENZA 3 3 Usando a definición, ∆[c y(t)] = c y(t+ 1) −c y(t) = c[y(t+ 1) −y(t)] = c∆y(t). 4 Aplicando a definición á función y(t)z(t), ∆[y(t)z(t)] = y(t+ 1) z(t+ 1) −y(t)z(t) =y(t+ 1) z(t+ 1) −y(t)z(t+ 1) + y(t)z(t+ 1) −y(t)z(t) = ∆ y(t)E z(t) + y(t) ∆ z(t). Podemos escribir agora, ∆2y(t) = ∆(∆ y(t)) = ∆(y(t+ 1) −y(t)) = ∆ y(t+ 1) −∆y(t) = (y(t+ 2) −y(t+ 1)) −(y(t+ 1) −y(t)) = y(t+ 2) −2y(t+ 1) + y(t). Do mesmo xeito, por indución, ∆ny(t) = y(t+n)−n y(t+n−1) + n(n−1) 2! y(t+n−2) + ···+ (−1)ny(t) = n X k=0 (−1)kn ky(t+n−k). Ademais destas propiedades, hai determinadas funcións que cumpren outras igualdades que poden resultar de utilidade. Antes de velas, imos definir unha función de gran importancia no estudo das ecuacións en diferenzas. Definición 1.4. Definiremos a función gamma,Γ, como Γ(z) = Z∞ 0 e−ttz−1dt, con z∈C, sempre que a parte real de zsexa positiva. Aplicando integración por partes, vemos que Γ(z+ 1) = Z∞ 0 e−ttzdt =−e−ttz∞ 0−Z∞ 0−e−tz tz−1dt =zZ∞ 0 e−ttz−1dt, polo que a función gamma satisfai a seguinte ecuación en diferenzas Γ(z+ 1) = zΓ(z), ou o que é o mesmo Γ(z) = Γ(z+ 1) z. 4CAPÍTULO 1. CÁLCULO EN DIFERENZAS Teorema 1.5. Sexa aunha constante. Entón, 1. ∆at= (a−1) at. 2. ∆ sen(a t) = 2 sen(a 2) cos(a(t+1 2)). 3. ∆ cos(a t) = −2 sen(a 2) sen(a(t+1 2)). 4. ∆ log(a t) = log(1 + 1 t). 5. ∆ log(Γ(t)) = log(t), sendo t>0. Demostración. 2. Aplicando a propiedade trigonométrica da suma de senos, temos que ∆ sen (a t) = sen(a(t+ 1)) −sen(a t) = 2 sen a(t+ 1) −a t 2cos a(t+ 1) + a t 2 = 2 sen a 2cos at+1 2. 3. Do mesmo xeito, usando a propiedade trigonométrica da suma de cosenos, temos que ∆ cos (a t) = cos (a(t+ 1)) −cos (a t) = −2 sen a(t+ 1) −t 2sen a(t+ 1) + a t 2 =−2 sen a 2sen at+1 2. Vexamos agora un exemplo onde combinamos o uso dos Teoremas 1.3 e 1.5. Exemplo 1.6. Calculemos ∆ sec(πt): ∆ sec(π t)=∆1 cos(π t)=(cos(π t)∆(1) −(1)(∆(cos(π t))) cos(π t) cos(π(t+ 1)) =2 sen(π 2) sen(π(t+1 2)) cos(π t) cos(π(t+ 1)) =2 (sen(π t) cos(π 2) + cos(π t) sen(π 2)) cos(π t)(cos(π t) cos(π)−sen(π t) sen(π)) =2 cos(π t) (cos(π t)(−cos(π t)) =−2 sec(π t). Unha das propiedades máis importantes empregadas no cálculo diferencial é d dttn=n tn−1. Desafortunadamente, isto non se cumpre no cálculo en diferenzas pois ∆ttn= (t+ 1)n−tn. Por este motivo, definimos a seguinte función, que xogará o papel das potencias para a derivada. 1.2. SUMA INDEFINIDA 5 Definición 1.7. Definimos, sempre que teña sentido, o factorial de orde r,tr, como tr=             t(t−1) ···(t−r+ 1), r ∈N, 1, r = 0, 1 (t+1)(t+2)···(t−r), r =−1,−2, ... Γ(t+1) Γ(t−r+1) , r /∈Z. O factorial de orde r cumpre as seguintes propiedades. Teorema 1.8. Tense que 1. ∆ttr=rtr−1. 2. ∆tt r=t r−1. 3. ∆tr+t t=r+t t+1. 1.2. Suma indefinida Introduciremos agora o inverso dereito do operador diferenza para poder realizar cálculos con el. Definición 1.9. Chamaremos suma indefinida (ou antidiferenza) de y, e denotarémolo por Py(t)ou ∆−1y(t), a calquera función que verifique que ∆Xy(t)=y(t) para todo t dentro do dominio de y. Vemos así que este operador terá o papel da integral no cálculo diferencial. É dicir, mentres que a derivada ten por inversa á integral, o operador diferenza terá por inversa a antidiferenza. Do mesmo xeito que acontece coa integral, a antidiferenza dunha función non é única, tal e como podemos ver no seguinte exemplo. Exemplo 1.10. Calculamos a suma indefinida P6t. Sabemos, pola propiedade 1 do Teorema 1.3, que ∆6t= 5 ·6t. Polo que X6t=6t 5. Ademais, se Cé unha función co mesmo dominio que 6te verifica que ∆C(t)=0, temos que 6CAPÍTULO 1. CÁLCULO EN DIFERENZAS ∆6t 5+C(t)= ∆ 6t 5= 6t. Daquela 6t 5+C(t)é unha suma indefinida de 6t. Diremos entón que X6t=6t 5+C(t), onde Cé unha función co mesmo dominio que 6te∆C(t)=0. Enunciamos así o seguinte teorema. Teorema 1.11. Se z é unha suma indefinida de y, entón toda suma indefinida de y virá dada por Xy(t) = z(t) + C(t), onde Cten o mesmo dominio que ye ademais ∆C(t)=0. No caso de que esteamos traballando en conxuntos discretos, temos o seguinte corolario. Corolario 1.12. Sexa yunha función definida nun conxunto da forma {a,a+1,a+2,... }, sendo aun número real, e sexa zunha suma indefinida de y. Entón toda suma indefinida de yvén dada por Xy(t) = z(t) + C, sendo C unha constante arbitraria. Empregando o Teorema 1.11 e o Corolario 1.12, podemos probar as propiedades que aparecen a continuación. Teorema 1.13. Sexa a unha constante. Entón, se C é unha función tal que ∆C(t) = 0, tense que 1. Pat=at a−1+C(t)se a6= 1. 2. Psen (a t) = −cos(a(t−1 2)) 2 sen(a 2)+C(t)se a 6= 2nπ. 3. Pcos (a t) = sen(a(t−1 2)) 2 sen(a 2)+C(t). 4. Plog(t) = log (Γ(t)) + C(t). 5. Pta=ta+1 a+1 +C(t). 6. Pt a=t a+1+C(t). 7. Pa+t t=a+t t−1+C(t). 1.2. SUMA INDEFINIDA 7 Teorema 1.14. Sexan yezfuncións con valores reais. Entón cúmprese que 1. P(y(t) + z(t)) = Py(t) + Pz(t). 2. PD y(t) = DPy(t), se D é unha constante. 3. P(y(t)∆z(t)) = y(t)z(t)−PE z(t)∆y(t). Demostración. 3. Pola propiedade 3 do Teorema 1.3 ∆(y(t)z(t)) = y(t)∆z(t) + E z(t) ∆y(t). Daquela X[y(t) ∆ z(t) + E z(t) ∆y(t)] = X[∆(y(t)z(t))] = y(t)z(t) + C(t). Para rematar, vexamos nas Táboas 1.1 e 1.2 as similitudes e diferenzas entre o cálculo diferencial e o cálculo en diferenzas, onde Drespresenta ao operador derivada. Cálculo en diferenzas Cálculo diferencial 1. ∆y(x) = y(x+h)−y(x). 2. ∆ny= ∆(∆n−1y), n = 1,2, ... 3. ∆(c y) = c∆y∀c∈R. 4. ∆(c1y1+c2y2) = c1∆y1+c2∆y2. 5. ∆xn=n xn−1. 6. ∆(u v) = (E u)∆v+v∆u. 7. ∆(u v) = v∆u−u∆v vEv . 8. Se yé un polinomio de grado n, entón ∆nyé unha constante e todas as diferenzas de orde maior son nulas. 1. Dy(x) = l´ımh→0∆y(x) h. 2. Dny= D(Dn−1y), n = 1,2, ... 3. D(c y) = cDy∀c∈R. 4. D(c1y1+c2y2) = c1Dy1+c2Dy2. 5. Dxn=n xn−1. 6. D(u v) = uDv+vDu. 7. D(u v) = vDu−uDv v2. 8. Se yé un polinomio de grado n, entón Dnyé unha constante e todas as derivadas de orde maior son nulas. Táboa 1.1: Comparación entre o operador diferenza e a derivada. 8CAPÍTULO 1. CÁLCULO EN DIFERENZAS Cálculo en diferenzas Cálculo diferencial 1. P(d y) = dPy∀d∈R. 2. P(d1y1+c2d2) = d1Py1+d2Py2. 3. Pxn=xn n−1+Csendo Cunha función tal que ∆C= 0. 4. P(u∆v)) = u v −PE v∆u. 1. R(d y) = dRy∀d∈R. 2. R(d1y1+c2d2) = d1Ry1+d2Ry2. 3. Rxndn=xn+1 n+1 +Csendo Ccalquera constante. 4. Rudv=u v −Rvdu. Táboa 1.2: Comparación entre o operador antidiferenza e a integral. Capítulo 2 Ecuacións lineais en diferenzas Neste capítulo falaremos de ecuacións en diferenzas, para iso seguiremos fundamentalmente o Capítulo 3 de [6] e as Seccións 2.1, 2.2 e 2.3 de [5]. Outras referencias onde se pode consultar son [3] e [7]. Antes de comezar co estudo e resolución de ecuacións en diferenzas, definamos os conceptos de ecuación en diferenzas e de solución da mesma, para despois poder buscar as solucións. Definición 2.1. Sexa yunha función con valores nun conxunto S. Chamaremos ecuación en diferenzas sobre Sa unha ecuación que conteña os valores dunha ou máis diferenzas ∆y, ∆2y, ... Normalemente considérase Scomo o conxunto dos número naturais, aínda que tamén se poden considerar intervalos de R. Por comodidade, no noso traballo empregaremos sempre un conxunto Sda forma {a, a + 1, a + 2, ... }, onde a∈R. Definición 2.2. Unha función yé unha solución dunha ecuación en diferenzas sobre un conxunto Sse os valores de yreducen a ecuación en diferenzas a unha identidade en S, é dicir, se os valores de yfan a ecuación en diferenzas unha afirmación certa en todo o conxunto S. Se unha función yé solución dunha ecuación en diferenzas, diremos que a función satisfai a ecuación en diferenzas. Exemplo 2.3. Consideremos a seguinte ecuación en diferenzas y(t+ 2) −2y(t+ 1) + y(t) = t(t+ 1), t ∈N, que tamén pode ser expresada como ∆2y(t) = t2. 9 10 CAPÍTULO 2. ECUACIÓNS LINEAIS EN DIFERENZAS Empregando o Teorema 1.8, sabemos que ∆2t4= ∆ 4 t3= 12 t2, polo que, y(t) = t4 12 é unha solución da nosa ecuación. 2.1. Resultados xerais de ecuacións en diferenzas lineais Definición 2.4. Unha ecuación en diferenzas sobre un conxunto Sélineal se se pode escribir como p0(t)y(t+n) + p1(t)y(t+n−1) + ···+pn−1(t)y(t+ 1) + pn(t)y(t) = r(t), t ∈S, (2.1) onde p0, p1, ... , pn−1, pnerson funcións con valores en S. Definición 2.5. Unha ecuación lineal nun conxunto Sdise de orde n se a podemos escribir da forma (2.1) con p0epnnon nulos en S. Definición 2.6. Diremos que unha ecuación en diferenzas lineal é homoxénea cando é da forma p0(t)y(t+n) + p1(t)y(t+n−1) + ···+pn−1(t)y(t+ 1) + pn(t)y(t)=0, t ∈S. (2.2) Notemos que a ecuación (2.5) tamén se pode escribir usando o operador shift pn(t)En+···+p0(t)E0y(t) = r(t). Ademais, como E= ∆+I, tamén podemos expresar a ecuación (2.1) en termos do operador diferenza, pero isto non é de utilidade pois con este operador non sabemos a orde da nosa ecuación, como podemos ver no seguinte exemplo. Exemplo 2.7. Consideremos a ecuación en diferenzas ∆3y(t) + 3 ∆2y(t)+∆y(t)−y(t) = r(t). Se substituimos ∆ = E−I, vemos que (E3−3E2+ 3 E−I)y(t) + 3(E2−2E+I)y(t)+(E−I)y(t)−y(t) = r(t). Operando e simplificando obtemos y(t+ 3) −2y(t+ 1) = r(t), que é de orde dous, aínda que na ecuación inicial puidésemos pensar que é de orde tres. 2.1. RESULTADOS XERAIS DE ECUACIÓNS EN DIFERENZAS LINEAIS 11 Definición 2.8. Chamaremos problema de valor inicial para unha ecuación da forma (2.1) á restricción de impoñer á nosa solución as seguintes condicións y(t+k) = yk para k= 0,1, ... , n −1. Exemplo 2.9. Consideremos o seguinte problema de valor inicial (P)(y(t+ 2) −4y(t+ 1) + 4 y(t) = 0, y(0) = 1 ey(1) = 6. A ecuación en diferenzas ten por solución xeral y(t)=2t(C1+C2t). Se lle impoñemos as condicións iniciais obtemos y(0) = C1ey(1) = 2 (C1+C2), de onde C1= 1 e2 (C1+C2)=6, é dicir, C1= 1 eC2= 2. Polo que a solución ao noso problema (P)será y(t) = 2t(1 + 2 t). O problema de valor inicial ten unha única solución, tal e como veremos no seguinte resultado. Teorema 2.10 (Teorema de existencia e unicidade de solución).Sexan p0, p1, ... , pn−1, pn erfuncións con valores en S={a, a + 1, a + 2, ... }con p0(t)6= 0 epn(t)6= 0 para todo t∈S. Entón, para todo t0∈Se calquera y0, ... , yn−1existe unha única solución yque satisfai a ecuación (2.1) para t∈Sey(t0+k) = ykpara k= 0, ... , n −1. Demostración. Por hipótese coñecemos os nprimeiros valores de yen t0, que serán y(t0+k) = yk,para k= 0, ... , n −1. Probaremos por indución que o valor de ypara todo t∈Sestá determinado de forma única. Se escribimos a ecuación (2.1) evaluada en t0temos que p0(t0)y(t0+n) + p1(t0)y(t0+n−1) + ···+pn−1(t0)y(t0+ 1) + pn(t0)y(t0) = r(t0). 18 CAPÍTULO 2. ECUACIÓNS LINEAIS EN DIFERENZAS 2.2.2. Ecuacións lineais con coeficientes constantes Trataremos agora de buscar nsolucións linealmente independentes da ecuación homoxénea (2.2) no caso de que os coeficientes sexan todos constantes. Como pn6= 0, podemos dividir ambos lados da ecuación por pne obtemos así u(t+n) + pn−1u(t+n−1) + ···+p0u(t)=0, t ∈S, (2.6) onde p0, ... , pn−1son constantes e p06= 0. Definición 2.22. 1. Chamaremos polinomio característico da ecuación (2.6) ao seguinte polinomio λn+pn−1λn−1+···+p0. 2. Chamaremos ecuación característica da ecuación (2.6) a λn+pn−1λn−1+···+p0= 0. 3. As solucións da ecuación característica, λ1, ... , λk,son os autovalores. Vexamos agora como obter un sistema de solucións linealmente independentes de (2.6). Para iso escribamos a ecuación (2.6) empregando o operador shift: En+pn−1En−1+···+p0u(t)=0, t ∈S. (2.7) Sexan λ1, ... , λkos autovalores asociados á súa ecuación característica e α1, ... , αkas multiplicidades correspondentes. Podemos escribir así a ecuación característica asociada a (2.7) como λn+pn−1λn−1+···+p0= (λ−λ1)α1···(λ−λk)αk= 0. Vemos así que (2.7) se pode expresar como (E−λ1)α1···(E−λk)αku(t) = 0, t ∈S, (2.8) onde α1+···+αk=ne a orde dos factores é arbitraria. Notemos que os autovalores son non nulos pois p06= 0. Se consideramos a ecuación (E−λ1)α1u(t)=0, t ∈S, (2.9) todas as solucións de (2.9) serán tamén solución de (2.8). Vexamos como calcular esta solución: Se α1= 1, a ecuación (2.9) redúcese a u(t+ 1) = λ1u(t), que ten como solución u(t) = λt 1,t∈S. 2.2. RESOLUCIÓN DE ECUACIÓNS LINEAIS EN DIFERENZAS 19 Se α1>1, consideremos u(t) = λt 1v(t)e substituamos en (2.9): (E−λ1)α1λt 1v(t) = α1 X i=0 α1 i(−λ1)α1−iEiλt 1v(t) = α1 X i=0 α1 i(−1)α1−iλα1−i 1Eiλt 1v(t) = α1 X i=0 α1 i(−1)α1−1λα1−i 1λt+i 1Eiv(t) =λα1+t 1 α1 X i=0 α1 i(−1)α1−iEiv(t) =λα1+t 1(E−1)α1v(t) = λα1+t 1∆α1v(t)=0, onde 1representa a función identidade. Polas propiedades do Teorema 1.5 sabemos que ∆α1v(t)=0se v(t)=1, t, t2, ... , tα1−1. Así, a ecuación (2.9) ten α1solucións que serán λt 1, t λt 1, t2λt 1, ... , tα1−1λt 1. Repetindo este proceso para todo αicon i∈ {1, ... , k}obtemos finalmente nsolucións da ecuación (2.6) que son linealmente independentes. Recollemos isto no seguinte teorema. Teorema 2.23. Sexan λ1, ... , λkcon multiplicidade α1, ... , αk, respectivamente, os autovalores da ecuación (2.6). Entón (2.6) ten nsolucións independentes que veñen dadas por λt 1, ... , tα1−1λt 1, λt 2, ... , tα2−1λt 2, ... , λt k, ... , tαk−1λt k. Observación 2.24.Se entre os autovalores temos un par complexo λ1=a+b i eλ2=a−b i, entón λt 1= (a+b i)teλt 2= (a−b i)t. Se expresamos os autovalores complexos en función das súas coordenadas polares temos que a=rcos(θ), b =rsen(θ), onde r=√a2+b2eθ= arctan b ase a6= 0,θ=π 2se a= 0 eb > 0eθ=−π 2se a= 0 e b < 0. Así c1λt 1+c2λt 2=c1(rcos(θ) + i r sen(θ))t+c2(rcos(θ)−i r sen(θ))t =rt[(c1+c2) cos(t θ) + i(c1−c2) sen(t θ)] =rt[a1cos(t θ) + a2sen(t θ)] . Polo tanto, en lugar de tomar como solucións independentes λt 1eλt 2, que serán funcións complexas, tomaremos rtcos(θ t)ertsen(θ t). No caso de que os autovalores complexos 20 CAPÍTULO 2. ECUACIÓNS LINEAIS EN DIFERENZAS λ=a±b i teñan multiplicidade α > 1temos que, de forma análoga aos autovalores reais, as solucións reais linealmente independentes virán dadas por rtcos(θ t), ... , tα−1rtcos(θ t), rtsen(θ t), ... , tα−1rtsen(θ t). Exemplo 2.25. Buscaremos todas as solucións da ecuación en diferenzas u(t+ 3) −7u(t+ 2) + 16 u(t+ 1) −12 u(t) = 0, t =a, a + 1, ... A ecuación característica asociada será λ3−7λ2+ 16 λ−12 = 0, ou o que é o mesmo (λ−2)2(λ−3) = 0. Así, aplicando o Teorema 2.23, sabemos que as tres solucións linealmente independentes da ecuación en diferenzas son u1(t)=2t, u2(t) = t2teu3(t) = 3t. Así, a solución xeral da ecuación vén dada por u(t) = C12t+C2t2t+C33t, con C1, C2, C3constantes arbitrarias. Agora imos ver como resolver unha ecuación non homoxénea con coeficientes constantes, y(t+n) + pn−1y(t+n−1) + ···+p0y(t) = r(t), t ∈S, (2.10) por un método coñecido como método do aniquilador, cando r(t)é a solución dunha ecuación homoxénea con coeficientes constantes. Recollemos a idea central deste método no seguinte teorema. Teorema 2.26. Supoñamos que yé solución da ecuación (2.10), é dicir, (En+pn−1En−1+···+p0)y(t) = r(t), t ∈S, e existen q0, ... , qm−1para os cales se satisfai (Em+qm−1Em−1+···+q0)r(t) = 0, t ∈S. Entón ysatisfai a ecuación (Em+···+q0)(En+···+p0)y(t)=0, t ∈S. 2.2. RESOLUCIÓN DE ECUACIÓNS LINEAIS EN DIFERENZAS 21 Demostración. Basta con aplicar o operador Em+···+q0a ambos lados da ecuación (2.10). Apliquemos este teorema a un caso particular. Exemplo 2.27. Consideremos a ecuación y(t+ 2) −7y(t+ 1) + 6 y(t) = t. (2.11) Comecemos escribindo esta ecuación empregado o operador shift (E2−7E+ 6)y(t) = t, ou equivalentemente (E−1)(E−6) y(t) = t. Podemos ver que a función tsatisfai a ecuación homoxénea (E−1)2t= ∆2t= 0, polo que (E−1)2será o noso aniquilador. Entón, aplicando o Teorema 2.26, a función y(t) é solución da seguinte ecuación (E−1)3(E−6) y(t)=0. Agora, empregando o Teorema 2.23, a solución xeral da ecuación homoxenénea anterior será da forma y(t) = C16t+C2+C3t+C4t2. Para determinar os coeficientes, imos substituir y(t)na ecuación orixinal. Como C16t+C2 satisfai a parte homoxénea da ecuación, abonda con substituír y(t) = C3t+C4t2, obtendo C3(t+ 2) + C4(t+ 2)2−7C3(t+ 1) −7C4(t+ 2)2+ 6 C3t+ 6 C4t2=t, ou o que é o mesmo t2[C4−7C4+ 6 C4] + t[4 C4+C3−14 C4−7C3+ 6 C3] + [4 C4+ 2 C3−7C4−7C3] = t. Se igualamos os coeficientes dos dous lados da igualdade, obtemos o seguite sistema de ecuacións −10 C4= 1 −5C3−3C4= 0 que ten por solución C4=−1 10 eC3=3 50. Así, a solución da ecuación (2.11) será y(t) = C16t+C2+3 50 t−1 10 t2. 22 CAPÍTULO 2. ECUACIÓNS LINEAIS EN DIFERENZAS Algunhas aplicacións das ecuacións lineais con coeficientes constantes Exemplo 2.28 (A sucesión de Fibonacci. O problema dos coellos).Este problema foi formulado por Leonardo de Pisa en [11] e xurdiu tratando de dar resposta á seguinte pregunta: Cantos pares de coellos haberá despois dun ano se se comeza cunha parella adulta de coellos e se cada parella ten un par de descendentes cada mes a partir da súa etapa adulta (o segundo mes)? O primeiro par de coellos reprodúcese ao final do primeiro mes, temos así dous pares. Ao final do segundo mes, só se pode reproducir a primeira parella e temos entón tres pares de coellos. Ao final do terceiro mes, a primeira e a segunda parella reprodúcense e temos cinco pares de coellos. Continuando con este procedemento, chegamos a que se F(n)é o número de pares de coellos ao final do mes n, entón a relación de recurrencia que representa este modelo vén dada pola ecuación linear en diferenzas de segunda orde seguinte F(n+ 2) = F(n+ 1) + F(n), F (0) = 1, F(1) = 2,0≤n≤10. Este exemplo é un caso particular da sucesión de Fibonacci, que vén dada por F(n+ 2) = F(n+ 1) + F(n), F (0) = 0, F(1) = 1, n ≥0.(2.12) A ecuación característica correspondente a (2.12) é λ2−λ−1 = 0. Os autovalores correspondentes son λ1=1+√5 2eλ2=1−√5 2. Así, a solución xeral de (2.12) é F(n) = a1 1 + √5 2!n +a2 1−√5 2!n , n ≥1. Empregrando os valores iniciais F(0) = 0 eF(1) = 1, obtemos que a1=1 √5, a2=−1 √5. Consecuentemente, F(n) = 1 √5" 1 + √5 2!n − 1−√5 2!n#, n ≥1. Exemplo 2.29 (Ley de Hardy-Weinberg, Sección 1.7 de [1]).Calquera característica dun individuo, como a cor de ollos, está determinado por un par de xens, un procedente do seu pai e outro procedente de súa nai. Cada xen ten dúas formas, a dominante denotada por Ae a recesiva denotada por a. Así, cada individuo pode ter tres combinacións diferentes de xens: dominantes AA, híbridos (Aa ou aA) ou recesivos aa. 2.2. RESOLUCIÓN DE ECUACIÓNS LINEAIS EN DIFERENZAS 23 Supoñamos que na xeración t-ésima, as proporcións de dominantes, híbridos e recesivos son p(t),q(t)er(t)respectivamente, de maneira que p(t) + q(t) + r(t)=1, p(t)≥0, q(t)≥0, r(t)≥0. Se consideramos que os apareamentos son aleatorios (é dicir, non depende de se os individuos son dominantes, híbridos ou recesivos), todos teñen a mesma supervivencia e igual fertilidade, vexamos como varían as proporcións de cada xenotipo ao longo das xeracións. Daquela, tendo en conta que p(t+ 1) é a probabilidade de que un individuo sexa dominante na xeración (t+ 1)-ésima, é dicir, a probabilidade de que recibise un xen dominante, A, do pai e outro da nai, podemos escribir que p(t+ 1) = p(t) + 1 2q(t)·p(t) + 1 2q(t)=p(t) + 1 2q(t)2 . Do mesmo xeito, q(t+ 1) = 2 p(t) + 1 2q(t)r(t) + 1 2q(t) e r(t+ 1) = t(t) + 1 2q(t)2 . Notemos que estas novas proporcións verifican que p(t+ 1) + q(t+ 1) + r(t+ 1) = p(t) + 1 2q(t) + 1 2q(t) + r(t)2 = 1. Se queremos ver as proporcións de dominantes da seguinte xeración p(t+ 2) = p(t+ 1) + 1 2q(t+ 1)2 = p(t) + 1 2q(t)2 +p(t) + 1 2q(t)r(t) + 1 2q(t)!2 =p(t) + 1 2q(t)2p(t) + 1 2q(t) + 1 2q(t) + r(t)2 =p(t) + 1 2q(t)2 =p(t+ 1). Obtemos así a seguinte ecuación en diferenzas p(t+ 2) = p(t+ 1). Esta ecuación tan simple demostranos que a proporción de individuos dominantes se matén constante ao longo das xeracións. Podemos realizar o mesmo proceso para a proporción de individuos híbridos e recesivos. Isto explica, por exemplo, porque aínda que o xen que determina a cor de ollos marrón é dominante, hoxe en día continua existindo xente con ollos de cor verde ou azul. 24 CAPÍTULO 2. ECUACIÓNS LINEAIS EN DIFERENZAS 2.2.3. Ecuacións lineais con coeficientes non constantes En xeral, as ecuacións de orde superior a un con coeficientes non constantes non poden ser resoltas na meirande parte dos casos. Por este motivo, veremos distintos métodos que servirán para obter solucións explícitas en determinadas situacións. Consideremos de novo a ecuación non homoxénea pn(t)y(t+n) + pn−1(t)y(t+n−1) + ···+p1(t)y(t+ 1) + p0(t)y(t) = r(t), t ∈S, (2.13) e a ecuación homoxénea asociada pn(t)y(t+n) + pn−1(t)y(t+n−1) + ···+p1(t)y(t+ 1) + p0(t)y(t)=0, t ∈S. (2.14) Método de variación de parámetros O método de variación de parámetros é o procedemento xeral para resolver (2.13). Se supoñemos que coñecemos nsolucións linealmente independentes de (2.14), a través deste método obtemos todas as solucións de (2.13) en termos de nsumas indefinidas. Veremos este método para o caso particular n= 2.Sexan u1eu2dúas solucións independentes de p2(t)y(t+ 2) + p1(t)y(t+ 1) + p0(t)y(t)=0, t ∈S. (2.15) Buscaremos unha solución da forma y(t) = a1(t)u1(t) + a2(t)u2(t), onde a1ea2teñen que ser determinados. Vemos que y(t+ 1) = a1(t+ 1) u1(t+ 1) + a2(t+ 1) u2(t+ 1) =a1(t)u1(t+ 1) + a2(t)u2(t+ 1) + ∆ a1(t)u1(t+ 1) + ∆ a2(t)u2(t+ 1). Se escollemos as funcións a1ea2de forma que ∆a1(t)u1(t+ 1) + ∆ a2(t)u2(t+ 1) = 0,(2.16) entón y(t+ 2) = a1(t+ 1) u1(t+ 2) + a2(t+ 1) u2(t+ 2) =a1(t)u1(t+ 2) + a2(t)u2(t+ 2) + ∆ a1(t)u1(t+ 2) + ∆ a2(t)u2(t+ 2). Se subsituímos agora as expresións de y(t),y(t+ 1) ey(t+ 2) na ecuación (2.15) temos que 2.2. RESOLUCIÓN DE ECUACIÓNS LINEAIS EN DIFERENZAS 25 p2(t)y(t+ 2) + p1(t)y(t+ 1) + p0(t)y(t) =a1(t)[p2(t)u1(t+ 2) + p1(t)u1(t+ 1) + p0(t)u1(t)] +a2(t)[p2(t)u2(t+ 2) + p1(t)u2(t+ 1) + p0(t)u2(t)] +p2(t)[u1(t+ 2) ∆ a1(t) + u2(t+ 2) ∆ a2(t)]. Como u1eu2satisfan (2.15), os dous primeiros termos da expresión son nulos e así, y(t)satisfai (2.13) se u1(t+ 2) ∆ a1(t) + u2(t+ 2) ∆ a2(t) = r(t) p2(t).(2.17) Podemos concluír así que y(t) = a1(t)u1(t) + a2(t)u2(t)é solución de (2.13) se ∆a1(t)e ∆a2(t)verifican as ecuacións (2.16) e (2.17). Este sistema de ecuacións lineais ten unha única solución se a matriz de coeficientes, W(t+ 1), ten determinante distinto de cero (en virtude do Teorema 2.15. Recollemos o caso para narbitrario no seguinte resultado. Teorema 2.30. Sexan u1, ... , unnsolucións linealmente independentes de (2.14). Entón y(t) = a1(t)u1(t) + ···+an(t)un(t) é unha solución de (2.13), sendo a1, ... , anfuncións que satisfan o seguinte sistema de ecuacións W(t+ 1)     ∆a1(t) . . . ∆an(t)     =    0 . . . r(t) pn(t)     . Método de factorización Escribamos agora a ecuación lineal de orde nen termos do operador shift (pn(t)En+···+p0(t)) y(t) = r(t), t ∈S. No caso de que o factor que involucra ao operador shift se poida factorizar en factores lineais, a solución á nosa ecuación en diferenzas pode obterse resolvendo unha serie de ecuacións de primeira orde. Ilustremos esta situación cun exemplo. Exemplo 2.31. Consideremos a ecuación en diferenzas E2−(t−1) E−ty(t)=0, t ∈S, 26 CAPÍTULO 2. ECUACIÓNS LINEAIS EN DIFERENZAS ou equivalentemente (E+ 1) (E−t)y(t)=0. Se consideramos entón a ecuación de primeira orde dada por (E+ 1) v(t)=0, que ten como solución v(t)=(−1)tC, para obter a solución da ecuación orixinal teremos que resolver a ecuación dada por (E−t)y(t)=(−1)tC. A ecuación homoxénea asociada será (E−t)y(t)=0. Desenvolvendo isto, temos que y(t+ 1) −t y(t) = 0, ou o que é o mesmo, y(t+ 1) = t y(t). Vemos así que esta ecuación en diferenzas satisfaina a función gamma, polo que a solución xeral da ecuación homoxénea é da forma DΓ(t+1) t, onde Dé unha constante arbitraria. Empregando o Teorema 2.21 temos que y(t) = DΓ(t+ 1) t+CΓ(t+ 1) tX(−1)tt Γ(t+ 2). Se consideramos que ttoma os valores discretos 0,1,2, ... entón y(t) = D(t−1)! + C(t−1)! t−1 X k=0 (−1)k(k−1) (k+ 1)! , onde CeDson constantes arbitrarias. Método de redución de orde Outro caso particular no que podemos calcular a solución dunha ecuación de orde superior a un é cando somos quen de atopar unha solución non nula da ecuación homoxénea asociada. Neste caso, poderemos reducir nunha unidade a orde da nosa ecuación. Así, se traballamos cunha ecuación de segunda orde, poderemos encontrar unha solución independente da primeira e conseguir así unha solución xeral da ecuación. Vexamos primeiro a ecuación de primeira orde que satisfai o Casoratiano. 2.2. RESOLUCIÓN DE ECUACIÓNS LINEAIS EN DIFERENZAS 27 Lema 2.32. Sexan u1, ... , unsolucións da ecuación pn(t)u(t+n) + ···+p0(t)u(t)=0, t ∈S, e sexa w(t)o Casoratiano correspondente. Entón w(t)satisfai w(t+ 1) = (−1)np0(t) pn(t)w(t), t ∈S. (2.18) Demostración. Denotemos por Fia fila i-ésima da matriz de Casorati. O valor de w(t+ 1) non varía se cambiamos a derradeira fila por Fn+p1(t) pn(t)F1+···+pn−1(t) pn(t)Fn−1. Así, esta nova fila será −p0(t) pn(t)u1(t),··· ,−p0(t) pn(t)un(t). Temos así que w(t+ 1) = det        u1(t+ 1) ··· un(t+ 1) . . ..... . . u1(t+n−1) ··· un(t+n−1) −p0(t) pn(t)u1(t)··· −p0(t) pn(t)un(t)        = (−1)np0(t) pn(t)w(t). Supoñamos agora que u1é unha solución non nula de p2(t)y(t+ 2) + p1(t)y(t+ 1) + p0(t)y(t) = 0, t ∈S, (2.19) e sexa u2outra solución. Recordando que ∆u2(t) u1(t)=u1(t) ∆ u2(t)−u2(t) ∆ u1(t) u1(t)u1(t+ 1) =w(t) u1(t)u1(t+ 1), obtemos que u2(t) u1(t)=Xw(t) u1(t)u1(t+ 1), ou o que é o mesmo u2(t) = u1(t)Xw(t) u1(t)u1(t+ 1). Obtemos así o seguinte teorema, que se coñece como método de redución de orde para ecuacións en diferenzas de orde dous. 34 CAPÍTULO 3. SISTEMAS DE ECUACIÓNS LINEAIS EN DIFERENZAS lineal de M0, ... , Mipara i= 0, ... , n −1,e empregando a Observación 3.5, deducimos que isto é certo para calquera potencia de A. Podemos escribir entón At= n−1 X i=0 ci+1(t)Mi,para t∈Z, onde ci+1(t)está por determinar. Tendo en conta que At+1 =A At, n−1 X i=0 ci+1(t+ 1)Mi=At+1 =A n−1 X i=0 ci+1(t)Mi= n−1 X i=0 ci+1(t) [Mi+1 +λi+1 Mi] = n−1 X i=1 ci(t)Mi+ n−1 X i=0 ci+1(t)λi+1 Mi, (3.8) aplicando a igualdade (3.7) no penúltimo paso. Esta ecuación satisfaise se ci(t)verifican o sistema     c1(t+ 1) . . . cn(t+ 1)     =           λ10 0 ··· 0 1λ20··· 0 0 1 λ3··· 0 . . ........... . . 0. . .0 1 λn               c1(t) . . . cn(t)     .(3.9) Ademais, como A0=I=c1(0)I+···+cn(0)Mn−1, temos que        c1(0) c2(0) . . . cn(0)        =       1 0 . . . 0        ,(3.10) é dicir, c1(t+ 1) = λ1c1(t), c1(0) = 1, ci(t+ 1) = λici(t) + ci−1(t), ci(0) = 0, i = 2, ... , n. Empregando o Teorema 3.7 sabemos que o problema de valor inicial dado por (3.9) e (3.10) ten unha única solución. Chegamos así ao seguinte resultado. Teorema 3.9 (O algoritmo de Putzer).A solución da ecuación (3.4) con vector de valores iniciais u0é u(t) = n−1 X i=0 ci+1(t)Miu0=Atu0, onde as Miveñen dadas por (3.7) e os ci(t)están determinados de forma única polas ecuacións (3.9) e(3.10). 3.1. SISTEMAS LINEAIS AUTÓNOMOS CON COEFICIENTES CONSTANTES 35 Exemplo 3.10. Resolveremos o seguinte sistema de ecuacións empregando o algoritmo de Putzer u(t+ 1) = "1 1 −1 3 #u(t), u(0) = "α β#. A ecuación característica correspondente é λ2−4λ+ 4 = 0, polo que a nosa matriz ten por autovalor λ= 2 con multiplicidade dous. Aplicando a ecuación (3.7), temos que M0=I, M1=A−2I="−1 1 −1 1 #. Polas ecuacións (3.9) e (3.10) temos que c1(t+ 1) = 2 c1(t), c1(0) = 1, de onde c1(t)=2t. Empregando de novo (3.9) e (3.10), c2(t+ 1) = 2 c2(t)+2t, c2(0) = 0, de onde c2(t) = t2t−1. Así, empregando o Teorema 3.9, temos que u(t) = (c1(t)I+c2(t)M1)"α β#= 2t"1 0 0 1 #+t2t−1"−1 1 −1 1 #!" α β# = 2t"1−t 2 t 2 −t 21 + t 2#" α β#. Sistemas non homoxéneos: variación de parámetros Finalmente, consideremos agora o sistema non homoxéneo u(t+ 1) = A u(t) + f(t).(3.11) O seguinte teorema establece unha fórmula de variación de parámetros para resolver este tipo de sistemas. Teorema 3.11. A única solución do sistema (3.11) satisfacendo a condición inicial u(0) = u0 é u(t) = Atu0+ t−1 X s=0 At−s−1f(s).(3.12) 36 CAPÍTULO 3. SISTEMAS DE ECUACIÓNS LINEAIS EN DIFERENZAS Demostración. Polo Teorema 3.7, abonda con comprobar que (3.12) satisfai o problema de valor inicial. Temos primeiro que, para t= 0 −1 X s=0 A−s−1f(s) = 0 por convención, daquela u(0) = u0. Para t≥1, u(t+ 1) = At+1 u0+ t X s=0 At−sf(s) = At+1 u0+ t−1 X s=0 At−sf(s) + f(t) =A"Atu0+ t−1 X s=0 At−s−1f(s)#+f(t) = A u(t) + f(t). 3.2. Sistemas lineais non autónomos Aínda que non imos entrar en detalle (para máis información véxase a Sección 3.2 de [3] e a Sección 4.4 de [6]) estudaremos agora sistemas con matriz de coeficientes non constantes. Recordemos que os sistemas non autónomos eran da forma u(t+ 1) = A(t)u(t) + f(t), t ∈S, e o seu sistema homoxéneo asociado era u(t+ 1) = A(t)u(t), t ∈S. Comecemos cunhas definicións básicas que nos permitan despois introducir o concepto de matriz fundamental, imprescindible no estudo de sistemas non autónomos. Definición 3.12. Diremos que as solucións u1, ... , unde (3.2) son linealmente independentes se c1u1(t) + ···+cnun(t)=0para todo t∈S⇒ci= 0,1≤i≤n. Se consideramos Φ(t)∈ Mn×ntal que as súas columnas son solución de (3.2), é dicir, Φ(t)=[u1(t), ... , un(t)] , entón, Φ(t+ 1) = [u1(t+ 1), ... , un(t+ 1)] = [A(t)u1(t), ... , A(t)un(t)] = A(t) [u1(t), ... , un(t)] =A(t) Φ(t). Vemos así que Φ(t)satisfai a ecuación Φ(t+ 1) = A(t) Φ(t).(3.13) 3.2. SISTEMAS LINEAIS NON AUTÓNOMOS 37 Definición 3.13. Chamaremos matriz fundamental do sistema (3.2) a unha matriz Φ(t)non singular para todo t∈Ssatisfacendo a ecuación (3.13). Podemos caracterizar as matrices fundamentais mediante o seguinte teorema. Teorema 3.14. Se Φ(t)é unha matriz fundamental de (3.2), entón Ψ(t)é outra matriz fundamental se e só se existe unha matriz non singular Ctal que Ψ(t) = Φ(t)Cpara todo t∈S. Demostración. ⇒Sexa Ψ(t) = Φ(t)C, onde Φ(t)é unha matriz fundamental de (3.2) e C é non singular. Entón Ψ(t)é non singular para todo t∈Se Ψ(t+ 1) = Φ(t+ 1) C=A(t) Φ(t)C=A(t) Ψ(t). Polo tanto, Ψ(t)é unha matriz fundamental de (3.2). ⇐Supoñamos que Ψ(t)eΦ(t)son matrices fundamentais de (3.2). Para algún t0∈S, sexa C= Φ−1(t0)Ψ(t0). Entón Ψ(t)eΦ(t)Cson solucións de (3.13), satisfacendo a mesma condición inicial e por unicidade temos que Ψ(t) = Φ(t)Cpara todo t∈S. Teorema 3.15. Se Φ(t)é unha solución de (3.13), entón det Φ(t)6= 0 para todo t∈Sou det Φ(t) = 0 para todo t∈S. Demostración. Como Φ(t)é solución de (3.13) para todo t∈S, Φ(t+ 1) = A(t) Φ(t). Polo tanto, det Φ(t+ 1) = det A(t) det Φ(t)para todo t∈S. Como por definición, det A(t)6= 0 para todo t∈Sentón a única posibilidade de que se cumpra a igualdade é que det Φ(t)6= 0 para todo t∈Sou que det Φ(t)=0para todo t∈S. En efecto, se det Φ(t0)=0para t0∈S, temos que det(t0+ 1) = det A(t) det Φ(t0)=0, det(t0) = det A(t) det Φ(t0−1) = 0, polo que det Φ(t)=0para todo t∈S. 38 CAPÍTULO 3. SISTEMAS DE ECUACIÓNS LINEAIS EN DIFERENZAS Teorema 3.16. Se Φ(t)é unha matriz fundamental de (3.2), entón a solución xeral de (3.2) estará dada por u(t) = Φ(t)c, onde cé un vector constante arbitrario. Os sistemas completos tamén se poden resolver empregando matrices fundamentais, tal e como se recolle no seguinte resultado. Teorema 3.17 (Método de variación de parámetros).Se Φ(t)é unha matriz fundamental de (3.13), a única solución de (3.5) que satisfai a condición inicial u(t0) = u0vén dada por u(t) = Φ(t) Φ−1(t0)u0+ Φ(t) t−1 X s=t0 Φ−1(s+ 1) f(s)para t≥t0. Vemos así a situación é análoga ao que acontecía no caso dos sistemas lineais de ecuacións diferenciais ordinarias. Capítulo 4 Teoría da estabilidade Ata agora estudamos ecuacións e sistemas lineais, pasemos a estudar as ecuacións e sistemas non lineais. Como normalmente estas ecuacións e sistemas non se poden resolver de forma explícita, imos estudar como se comportan as súas solucións, a medida que aumenta t, que é o que se coñece como comportamento asintótico e da lugar á teoría da estabilidade. Para iso seguiremos a Sección 1.3 de [3] e a Sección 4.5 de [6]. Tamén se pode seguir na Sección 4.1 de [5] ou en [7]. Centrarémonos nas ecuacións de orde un, pois podemos ver que calquera ecuación de orde maior, n, se pode transformar nun sistema de necuacións de orde un, como se ve a continuación. Exemplo 4.1. Consideremos a seguinte ecuación non linear de orde tres y(t+ 3) = 4 y4(t+ 2) y(t+ 1) −2 cos y(t). Podemos reescribir esta ecuación como o seguinte sistema de ecuacións en diferenzas de primeira orde, considerando u1(t) = y(t),u2(t) = y(t+ 1) eu3(t) = y(t+ 2), temos     u1 u2 u3     (t+ 1) =     u2(t) u3(t) 4u4 3(t)u2(t)−2 cos u1(t)     . Por este motivo, todos os resultados enunciados para sistemas de necuacións non lineais son tamén válidos para ecuacións non lineais de orde n. Estudaremos a estabilidade para sistemas autónomos non lineais da forma u(t+ 1) = f(u(t)), t ∈S={a, a + 1, a + 2, ... },(4.1) onde u(t)∈Rnpara t∈Sefé unha función continua de Rncon valores en Rn.Para iso, veremos algún resultado teórico e un método gráfico, coñecido como método da escaleira, pero primeiro introduciremos algúns conceptos básicos sobre estabilidade. 39 40 CAPÍTULO 4. TEORÍA DA ESTABILIDADE Definición 4.2. Un vector u∈Rné un punto fixo de fse f(u) = u. Un vector v∈Rn é un punto periódico de fse existe un enteiro positivo ktal que fk(v) = v(entendendo fkcomo evaluar kveces f) e ké o período de v. Así, podemos ver que os puntos fixos son solucións constantes do sistema (4.1) e os puntos periódicos proporcionan solucións periódicas. O período knon é único pois calquera múltiplo enteiro de ké tamén un período, por isto, chamaremos primeiro período ao menor período positivo. Definición 4.3. Chamaremos p-ciclo a unha solución periódica de primeiro periodo igual a p. Vexamos agora algunhas definicións sobre conceptos básicos da teoría da estabilidade. Definición 4.4. 1. Sexa u∈Rner > 0. Chamaremos bóla aberta centrada en u con radio rao conxunto B(u, r) = {v∈Rn:kv−uk< r}. 2. Sexa vun punto fixo de f. Diremos que véestable se dada calquera bóla B(v, ε), existe unha bóla B(v, δ)tal que se u∈B(v, δ), entón ft(u)∈B(v, ε)para t≥0. Diremos que véinestable se non é estable. 3. Diremos que véasintoticamente estable se é estable e existe unha bóla B(v, r) tal que ft(u)→vcando t→ ∞ para todo u∈B(v, r). 4. Sexa wun punto periódico de fcon período k. Diremos que wéestable (asintoticamente estable) se w, f(w), ... , fk−1(w)son estables (asintoticamente estables) como puntos fixos de fk. Observación 4.5.Así, un punto fixo vé estable se os puntos próximos a vnon se afastan de vao ser evaluados sucesivamente por f. Se ademais é asintoticamente estable, entón toda solución de (4.1) que comeza cerca de vconverxe a v. Vexamos un resultado para saber cando un punto fixo é asintoticamente estable para o caso escalar. Teorema 4.6. Sexa f:R→Runha función, cun punto fixo ve a súa primeira derivada continua nun intervalo aberto contendo a v. Entón: 1. Se |f0(v)|<1,vé asintoticamente estable. 2. Se |f0(v)|>1,vé inestable. 4.1. MÉTODO DA ESCALEIRA 41 Demostración. 1. Supoñamos que |f0(v)|<1. Por ser f0continua existe αtal que |f0(u)| ≤ α < 1en algún intervalor da forma I= (v−δ, v +δ)con δ > 0.Polo Teorema do Valor Medio temos que se u, w ∈I, entón |f(u)−f(w)|=|f0(c)||u−w| ≤ α|u−w|. Ademais, para cada u∈I, |f(u)−v| ≤ α|u−v|< δ, polo que f(u)∈Ie podemos concluir que vé estable. Por outra parte, |ft+1(u)−v| ≤ α|ft(u)−v| para cada t≥0eu∈I. Entón, por inducción, |ft(u)−v| ≤ αt|u−v|para t≥0eu∈I. Como l´ım t→∞ αt= 0, cada solución da ecuación (4.1) que comeza en Iconverxe a v cando t→ ∞ polo que vé asintoticamente estable. 2. Supoñamos que |f0(v)|>1. Entón, como f0é continua, existen λ > 1eε > 0tales que |f0(c)| ≥ λ > 1para todo λ > 1ec∈I= (v−ε, v +ε). Entón, polo Teorema do Valor Medio, |f(u)−f(w)|=|f0(c)||u−w| ≥ λ|u−w|,para todo u, w ∈I. Por inducción, |ft(u)−v| ≥ λt|u−w|, mentres ft(u)∈I. Como λ > 1, todas as solucións da ecuación (4.1) que comezan en I, excepto a solución constante u(t) = v, abandonan o intervalo Ipara to suficientemente grande. Concluimos así que vé inestable. 4.1. Método da escaleira Trátase dun método gráfico para estudar a estabilidade de ecuacións non lineais de primeira orde. Consideremos unha ecuación non lineal da seguinte forma u(t+ 1) = f(u(t)), t ∈S. 42 CAPÍTULO 4. TEORÍA DA ESTABILIDADE Este método consiste en debuxar as gráficas das funcións y=f(u)ey=unun mesmo plano de coordenadas. Así, os puntos fixos de fserán as interseccións destas dúas gráficas. Escollemos un valor inicial u(0) no eixo de abscisas e subimos verticalmene ata o punto (u(0), f(u(0)) = (u(0), u(1)) na gráfica de f. Agora, movémonos horizontalmente ao punto (u(1), u(1)) na liña y=u. Tras isto, consideramos u(2) = f(u(1)) e movémonos verticalmente a (u(1), u(2)) na gráfica de f. Así, alternando sucesivamente movementos verticais cara a gráfica de fe movementos horizontais cara a liña y=uxeramos unha sucesión {u(t)}, o cal nos permite observar o comportamento asintótico da solución. Ilustremos este método co seguinte exemplo. Exemplo 4.7. Consideremos a seguinte ecuación non lineal en diferenzas u(t+ 1) = 1 2u(t) (4 −u(t)). A función f(u) = 1 2u(4 −u)ten dous puntos fixos, o 0e o 2. En efecto, f(u) = u⇔1 2u(4 −u) = u⇔       u= 0 ou 1 2(4 −u)=1        ⇔       u= 0 ou u= 2        . Vexamos o comportamento das solucións mediante o método da escaleira en función de onde está a condición inicial u0que lle impoñemos a nosa ecuación. Comecemos considerando u0= 0,5e apliquemos o método da escaleira como se ve na seguinte gráfica. u y y=1 2u(4 −u) y=u u0u1u2u3u420 Figura 4.1: Método da escaleira para a función f(u) = 1 2u(4 −u)tomando como punto inicial u0∈(0,2). Podemos observar así na Figura 4.1 que se tomamos u0∈(0,2), a recta está por debaixo 4.1. MÉTODO DA ESCALEIRA 43 da parábola. Así, realizando este método iremos avanzando cara a dereita aproximándonos cada vez máis a 2. Consideremos agora u0= 3,5e apliquemos de novo o método da escaleira. u y y=1 2u(4 −u) y=u u0 u1u2u3 u4 20 Figura 4.2: Método da escaleira para a función f(u) = 1 2u(4 −u)tomando como punto inicial u0∈(2,4). Vemos así que se consideramos u0∈(2,4), como na Figura 4.2 a recta y=uestá enrriba da parábola, polo que na primeira iteración, cando avancemos horizontalmente, iremos cara a esquerda e despois estaremos na mesma situación que o caso anterior. Vexamos o que acontece cando tomamos como punto inicial u0/∈(0,4). Para iso, comecemos considerando a condición inicial u0=−0,5e apliquemos o método da escaleira. u y y=1 2u(4 −u) y=u u0 u1 u2 20 Figura 4.3: Método da escaleira para a función f(u) = 1 2u(4 −u)tomando como punto inicial u0∈(−∞,0). 50 CAPÍTULO 4. TEORÍA DA ESTABILIDADE Este modelo depende de r=α−β, que é o que se coñece como taxa de crecemento. Reescribimos así a ecuación (4.6) como Nt+1 =Nt+r Nt. Se r < 0, o número de nacementos en cada tempada é menor que o número de defuncións e polo tanto a poboación diminúe ata a súa extinción. Pola contra, se r > 0, o número de nacementos é maior que o de defuncións e a poboación incrementa infinitamente. Vexamos isto plasmado na Figura 4.6. t Nt 5 10 15 20 100 200 300 400 r= 0,1 r=−0,1 Figura 4.6: Crecemento poboación loxístico lineal para r= 0,1en vermello e para r=−0,1 en azul. Esta ecuación lineal non se pode xeneralizar para moitas poboacións pois, como podemos ver na natureza, o aumento das poboacións non é ilimitado, xa que o ecosistema no que viven os individuos ten o que se coñece como capacidade de carga, que é o tamaño máximo de poboación dunha especie biolóxica que se pode soster nese ecosistema específico, dado o alimento, o hábitat, a auga e outros recursos dispoñibles. Así, realizando estudos nos laboratorios puidose observar como a medida que aumenta a poboación tamén o fai o número de falecementos, mentres que se reduce o número de nacementos. Isto débese ao sobrepoboamento e a competencia polo alimento. Por este motivo, a partir da ecuación lineal de crecemento poboacional Nt+1 =Nt+r Nt, chegamos a unha ecuación non lineal ao incorporarlle o sobrepoboamento, obtendo Nt+1 =Nt+R(Nt)Nt,(4.7) 4.3. O MODELO LOXÍSTICO 51 onde R(Nt)é a ratio de crecemento, unha función que depende do tamaño da poboación Nt. Temos que impoñer certas condicións a esta función de crecemento da poboación: Debido ao sobrepoboamento, o número de mortes aumenta e o de nacementos diminúe, polo que R(Nt)ten que decrecer cando Ntaumenta. Se Nt=K, sendo Ka capacidade de carga, o crecemento é nulo, é dicir, R(K)=0. Cando Nt→0os efectos de sobrepoboación diminúen e o ratio de crecemento tende a unha constante, r, que representa o ratio de crecemento sen restriccións. Así, R(0) = r. Existen numerosas funcións que satisfan estas condicións, pero nós escolleremos a seguinte R(Nt) = −r KNt+r. (4.8) Substituíndo (4.8) en (4.7) temos que Nt+1 =Nt+r Nt1−Nt K,(4.9) que é o que se coñece como ecuación loxística discreta. Esta ecuación non lineal, como moitas outras, aínda non foi resolta de forma xeral, pero si se realizaron importantes estudos mediante iteracións para diferentes valores dos parámetros r(ratio de crecemento sen restrición), K(capacidade de carga) e N0(poboación inicial), como veremos a continuación. Comecemos buscando os puntos fixos da ecuación loxística, para iso, substituimos en (4.9) pola seguinte incognita Nt+1 =Nt=s, e obtemos s=s+r s 1−s K, ou o que é o mesmo, s1−s K= 0. Esta ecuación ten dúas solucións, s= 0 es=K. A primeira delas danos un resultado que parece obvio, se a poboación inicial é nula, continuará sendo nula na seguinte tempada de cría, pois non haberá ningún individuo que poida ter descendencia. A segunda solución constante, s=K, correspóndese coa capacidade de carga, é dicir, se a poboación inicial é N0=K, entón esta mantense estable no mesmo valor. Se escribimos agora a ecuación (4.9) como Nt+1 −Nt=r Nt1−Nt K(4.10) 52 CAPÍTULO 4. TEORÍA DA ESTABILIDADE podemos ver como é a taxa de variación da poboación entre dúas tempadas de cría sucesivas. Vemos así que se Nt< K a poboación total aumentará, pola contra, se Nt> K a poboación total diminuirá. A maior limitación deste modelo é que pode producir valores negativos de poboación se r > 3. Vexamos agora un caso particular dun ecosistema que para unha determinada especie ten unha capacidade de carga K= 100 e iremos variando o ratio de crecemento sen restricións, r, entre 0e3. Estes distintos tipos de comportamento pódense clasificar en tres grupos: crecemento estable, crecemento cíclico e crecemento caótico. 4.3.1. Crecemento estable (0< r ≤2) Comecemos primeiro considerando r= 0,8e empreguemos o método da escaleira para dúas poboacións iniciais diferentes, unha inferior á capacidade de carga e outra superior á capacidade de carga. Nt 10 140 K 0 Figura 4.7: Método da escaleira aplicado a ecuación Nt+1 =Nt+0,8Nt1−Nt 100tomando como condicións iniciais 10, en cor vermella, e 140, en cor azul. Vemos así na Figura 4.7 que cando a poboación inicial é menor que a capacidade de carga, esta vai aumentando aproximándose cada vez máis á capacidade de carga. Se a poboación inicial é maior que a capacidade de carga, xa na primeira iteración diminúe a poboación ata un valor menor que a capacidade de carga e estamos na mesma situación que o caso anterior. Representemos agora a poboación en cada época de cría nunha gráfica (Figura 4.8) para ver como se vai acercando a 100 e se estabiliza nese valor. 4.3. O MODELO LOXÍSTICO 53 t Nt K 100 123456789 Figura 4.8: Representación do número de individuos en cada época de cría para a ecuación loxística Nt+1 =Nt+ 0,8Nt1−Nt 100. Vexamos agora como se consideramos un ratio de crecemento sen restriccións máis próximo a 2, a poboación aumenta de forma máis rápida e tarda menos épocas de cría en achegarse á capacidade de carga, K= 100. Empregamos así o método da escaleira con r= 1,6considerando de novo as cantidades iniciais, N0,10 e140. Nt 10 140 K 0 Figura 4.9: Método da escaleira aplicado a ecuación Nt+1 =Nt+1,6Nt1−Nt 100tomando como condicións iniciais 10, en cor vermella, e 140, en cor azul. Vemos así na Figura 4.9 que acontece o mesmo que no caso anterior, se tomamos como poboación inicial N0= 10 a poboación vai aumentando ata estabilizarse entorno á capacidade de carga, K= 100. Se tomamos como poboación inicial N0= 140, que é maior 54 CAPÍTULO 4. TEORÍA DA ESTABILIDADE que capacidade de carga, a poboación diminúe na primeira iteración ata un valor inferior á capacidade de carga e despois estamos na mesma situación que o caso anterior. t Nt 100 K 123456789 Figura 4.10: Representación do número de individuos en cada época de cría para a ecuación loxística Nt+1 =Nt+ 1,6Nt1−Nt 100. A única diferenza resaltable entre r= 0,8er= 1,6é que no primeiro caso a poboación nunca aumenta ata superar a capacidade de carga pois a recta corta a parábola antes de que esta acade o seu máximo, entón polo método da escaleira nunca imos obter un Nt>100. Pola contra no segundo caso a recta corta á parábola despois do seu máximo, polo que cando nos imos aproximando á capacidade de carga, vemos como o método da escaleira nos leva a unha espiral entorno ao valor 100, no que nos imos aproximando á capacidade de carga por defecto e por exceso alternadamente (como podemos observar na Figura 4.10), polo que si se chega a ter Nt>100 en varias etapas. Diferenciamos así dous tipos de comportamento estable, mentres que para 0< r ≤1 a poboación crece ata estabilizarse entorno a 100 sen chegar a sobrepasalo en ningún momento, para 1< r ≤2a poboación tamén se estabiliza en 100, pero primeiro fluctúa arredor deste valor. 4.3.2. Crecemento cíclico (2< r < 2,57) Vexamos agora un novo comportamento para este rango de valores, para iso tomamos a ecuación loxística con r= 2,4e aplicamos o método da escaleira (Figura 4.11). 4.3. O MODELO LOXÍSTICO 55 Nt 10 K 0 Figura 4.11: Método da escaleira aplicado a ecuación Nt+1 =Nt+2,4Nt1−Nt 100tomando como condición inicial N0= 10. Vemos así que a poboación oscila entre valores superiores á capacidade de carga e valores inferiores a esta. Así, calculando as sucesivas iteracións vemos que a nosa órbita tende a un 2-ciclo. Para calcular os dous valores entre os que oscila este ciclo, resolvemos a ecuación Nt+2 =Nt.(4.11) Como Nt+2 =Nt+1 + 2,4Nt+1 1−Nt+1 K, podemos escribir Nt+1 + 2,4Nt+1 1−Nt+1 K=Nt, e substituindo Nt+1 =Nt+ 2,4Nt1−Nt 100obtemos Nt+ 2,4Nt1−Nt 100+2,4Nt+ 2,4Nt1−Nt 100 1−Nt+ 2,4Nt1−Nt 100  K!=Nt, que ten por solución 119,305 e64,0281, que serán os valores do 2-ciclo ao que tenderá a nosa órbita. Vemos na Figura 4.12 como os valores oscilan cada dúas épocas de cría e se van aproximando as solucións da ecuación (4.11). 56 CAPÍTULO 4. TEORÍA DA ESTABILIDADE t Nt 100 K 119.305 64.0281 5 10 15 20 25 30 35 40 45 50 55 60 Figura 4.12: Representación do número de individuos en cada época de cría para a ecuación loxística Nt+1 =Nt+ 2,4Nt1−Nt 100. Así, N6= 120,40665, N8= 118,2973818, ... , N24 = 119,173124, ... , N60 = 119,3042752 eN7= 61,43633727, N9= 66,39251878, ... , N23 = 63,63511996, ... , N59 = 64,02527594, ... e vemos como os valores pares sen van aproximando paulatinamente a 119,305 e os valores impares a 64,0281. Se tomasemos r≥2,5poderiamos ver como o comportamento do número de individuos seguiría ciclos de orde cada vez maior a medida que nos aproximamos a2,57. 4.3.3. Crecemento caótico (2,57 < r ≤3) Dentro deste rango de valores o comportamento deixa de seguir un patrón fixo como nos casos anteriores, é o que se coñece como caos. Apliquemos o método da escaleira á nosa ecuación para r= 2,7e vexamos o que acontece. Nt 10 K 0 Figura 4.13: Método da escaleira aplicado a ecuación Nt+1 =Nt+2,7Nt1−Nt 100tomando como condición inicial N0= 10. 4.3. O MODELO LOXÍSTICO 57 Vemos así na Figura 4.13 que este tipo de crecemento poboacional non segue ningún patrón. Un dos primerios en descubrir que un modelo tan sinxelo podía propiciar un comportamento tan complexo foi Rober May en diversos artículos que publicou na segunda parte do século XX (ver [8] e [9]). Porén, temos que destacar que para determinados valores dentro deste intervalo o crecemento poboacional sí presenta algún tipo de patrón recoñecible, como pode ser un 3-ciclo. Outra das peculiaridades deste comportamento é a gran variación ao longo das xeracións que supón cambiar un pouco o valor inicial da poboación. Por isto falamos de caos, ademais de non estabilizarse a poboación entorno a un comportamento concreto, para condicións iniciais moi próximas teñense comportamentos moi distintos. Comprobémolo na Figura 4.14 tomando como valores iniciais 10 individuos, en cor vermella, e 13 individuos, en cor azul. t Nt 100 K 1 2 3 4 5 6 7 8 9 10 11 12 13 14 Figura 4.14: Representación do número de individuos en cada época de cría para a ecuación loxística Nt+1 =Nt+2,7Nt1−Nt 100tomando como valor inicial N0= 10, en cor vermella, eN0= 13, en cor azul. Por este motivo, neste tipo de modelos é complexo saber cal vai ser o comportamento da poboación para un período longo de tempo. 58 CAPÍTULO 4. TEORÍA DA ESTABILIDADE Bibliografía [1] R. P. Agarwal, Difference Equations and Inequalities. Theory, Methods, and Applications. Marcel Dekker, Inc., 2000. [2] L. Edelstein-Keshet, Mathematical Models in Biology. SIAM, 1988. [3] S. N. Elaydi, An Introduction to Difference Equations. Springer, 1996. [4] G. Fulford, P. Forrester and A. Jones, Modelling with Differential and Difference Equations. Cambridge University Press, 1997. [5] S. Goldberg, Introduction to Difference Equations with Illustrative Examples from Economics, Psychology, and Sociology. Dover Publications, Inc., New York, 1986. [6] W. G. Kelley and A. C. Peterson, Difference Equations: An Introduction with Applications. Academic Press, 1991. [7] V. Lakshmikantham and D. Trigiante, Theory of Difference Equations: Numerical Methods and Applications. CRC Press, 1988. [8] R. M. May, Biological Populations obeying Difference Equations: stable points, stable cycles, and chaos. Journal of Theoretical Biology, 1975. [9] R. M. May, Simple mathematial models with very complicated dynamics. Nature, 1976. [10] J. Maynard Smith, Mathematical Ideas in Biology. Cambridge University Press, 1968. [11] Leonardo de Pisa, Liber abaci. 1202. 59