Aproximación numérica de variedades invariantes en ... · vamos a dibujar el diagrama de fases,...

Julen Sánchez Oroz

Víctor Lanchares Barrasa

Facultad de Ciencias, Estudios Agroalimentarios e Informática

Grado en Matemáticas

2013-2014

Título

Director/es

Facultad

Titulación

Departamento

TRABAJO FIN DE GRADO

Curso Académico

Aproximación numérica de variedades invariantes en sistemas dinámicos

Autor/es

© El autor© Universidad de La Rioja, Servicio de Publicaciones, 2014

publicaciones.unirioja.esE-mail: [email protected]

Aproximación numérica de variedades invariantes en sistemas dinámicos,trabajo fin de grado

de Julen Sánchez Oroz, dirigido por Víctor Lanchares Barrasa (publicado por la Universidad de La Rioja), se difunde bajo una Licencia

Creative Commons Reconocimiento-NoComercial-SinObraDerivada 3.0 Unported. Permisos que vayan más allá de lo cubierto por esta licencia pueden solicitarse a los

titulares del copyright.

Aproximación numérica devariedades invariantes en

sistemas dinámicos

Autor: Julen Sánchez Oroz

Tutor: Víctor Lanchares Barrasa

Curso académico: 2013/2014

Facultad de Ciencias, Estudios Agroalimentariose Informática

Grado en Matemáticas

Universidad de La Rioja

Resumen

En esta memoria nos hemos interesado, en primera instancia, en dar unasnociones generales de lo que son los sistemas dinámicos. En particular, noshemos centrado en la existencia de variedades invariantes asociadas a los puntoscríticos, aunque también hablamos en menor medida de otro tipo de variedadesinvariantes que aparecen en los sistemas dinámicos.

A continuación, nos hemos preocupado en el modo en el que se calculanlas variedades invariantes, primero en sistemas lineales, luego en sistemas nolineales de manera analítica, y finalmente, ante la imposibilidad de extender losmétodos analíticos para calcular las variedades en su totalidad, nos centramosen los métodos numéricos, que constituyen la segunda parte de la memoria.

En dicha parte, analizamos punto a punto toda la problemática existente a lahora de calcular estas variedades. Así, mediante un ejemplo sencillo, ponemos demanifiesto los diferentes problemas y como deberíamos solventarlos. Finalmente,se da un algoritmo existente en la literatura que recoge toda esta problemática.

SummaryIn this report we were interested, first of all, in giving general ideas of whatthe dynamic systems are. Specifically, we have been focused on the existence ofinvariant manifolds associated to the critical points, although we also talk to alesser extent about other kinds of invariant manifolds appearing in dynamicalsystems.

After that we have been concerned about how the invariant manifolds arecalculated, first in linear systems, then in analytically non-linear systems, andfinally, facing the impossibility of extending the analytical methods to calcu-late the manifolds as a whole, we have focused in numerical methods, whichconstitute the second part of the report.

In this part we analyze case by case all the existing problems on calculatingthese manifolds. Thus, using a simple example we highlight the different pro-blems and how we should solve them. Finally an existing algorithm that includesall these problems is given.

iii

Índice general

1. Nociones básicas 11.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2. Espacio de fases . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.3. Conjunto límite y trayectorias . . . . . . . . . . . . . . . . . . . . 6

2. Ecuaciones diferenciales lineales 92.1. Variedades invariantes . . . . . . . . . . . . . . . . . . . . . . . . 92.2. Linealización de sistemas no lineales . . . . . . . . . . . . . . . . 122.3. Resultados de existencia . . . . . . . . . . . . . . . . . . . . . . . 142.4. Persistencia de puntos hiperbólicos . . . . . . . . . . . . . . . . . 18

3. Cálculo numérico de variedades 213.1. Sistemas 2−dimensionales . . . . . . . . . . . . . . . . . . . . . . 213.2. Sistemas 3−dimensionales . . . . . . . . . . . . . . . . . . . . . . 233.3. Aproximación por curvas geodésicas . . . . . . . . . . . . . . . . 30

4. Apendice 35

Bibliografía 39

v

Capítulo 1

Nociones básicas desistemas dinámicos

1.1. IntroducciónDe sobra se sabe que las ecuaciones diferenciales se usan en la ciencia para

modelizar muchos procesos dinámicos. Sin embargo, para multitud de ecua-ciones diferenciales no lineales no somos capaces de encontrar su solución. Envez de esto, lo que intentaremos será obtener información cualitativa sobre elcomportamiento de sus soluciones; si son periódicas, eventualmente periódicas,atractivas etc...

Esto ha dado lugar a lo que se conoce como teoría cualitativa de ecuacionesdiferenciales o sistemas dinámicos. Comencemos dando la definición de lo quese entiende por sistema dinámico.

Definición 1 Un sistema dinámico en E, abierto de Rn, es una función de

clase C1:

ϕ : R × E → E (1.1)

tal que si ϕt(x) = ϕ(t, x) entonces se verifica:

(i) ϕ0(x) = x ∀x ∈ E,

(ii) ϕt(x) ◦ ϕs(x) = ϕt+s(x) ∀s, t ∈ R , ∀x ∈ E.

Un sistema dinámico puede interpretarse como la evolución en el tiempo de unaserie de variables que, a veces, se denominan de estado, según sea el problema.

Ejemplo 1 Si A es una matriz n×n entonces ϕ(t, x) = eAtx define un sistemadinámico en R

n, ya que, como puede verse, se verifican de inmediato (i) y (ii)dados en la Definición 1. Además, para cada x0 ∈ R

n, ϕ(t, x0) es solución delproblema de valor inicial:

1

2 CAPÍTULO 1. NOCIONES BÁSICAS

{x = Ax,x(0) = x0.

(1.2)

En este ejemplo vemos que podemos asociar a un sistema dinámico un sis-tema de ecuaciones diferenciales. Esto es así en general ya que, dado ϕ(t, x),sistema dinámico en E, para cada x0 ∈ E, ϕ(t, x0) es solución del problema devalor inicial:

⎧⎨⎩

x = f(x), f(x) =d

dtϕ(t, x)|t=0,

x(0) = x0.

(1.3)

Recíprocamente, dado un sistema de ecuaciones diferenciales x = f(x), conf ∈ C1(E), la solución de éste, ϕ(t, x0), define (junto con una condición inicialx(0) = x0) un sistema dinámico si y solo si ∀x ∈ E el intervalo maximal deexistencia de ϕ(t, x0) es R.

En este sentido, podemos identificar sistemas dinámicos y ecuaciónes dife-renciales. La función ϕ se conoce a menudo como flujo y representa la evolucióncon el tiempo de una cierta condición inicial.

Ejemplo 2 Un ejemplo típico de sistéma dinámico es el que describe el ángulode un péndulo simple, θ, respecto al tiempo. Este satisface la ecuación diferencial

d2θ

dt2 + sen θ = 0. (1.4)

Este es un problema en el que no puede encontrarse explícitamente la funciónϕ(t, θ0, θ

′0) (con θ0 una condición inicial) en términos de funciones elementales.

Por eso, este sistema es en ocasiones aproximado por otro lineal

d2θ

dt2 + θ = 0, (1.5)

cuando se consideran oscilaciones pequeñas, es decir cuando sen θ ≈ θ. En estecaso sí podemos determinar el flujo, que vendrá definido en E = [0, 2π) × R

como ϕ(t, θ0, θ′0) = (θ0 cos t + θ

′0 sin t, −θ0 sin t + θ

′0 cos t).

1.2. Espacio de fasesNótese que un sistema dinámico no solo es el flujo, también juega un papel

importante el conjunto E donde está definido. Este conjunto se denomina espa-cio de fases y es en dicho espacio donde puede representarse la evolución de lascondiciones iniciales bajo la acción del flujo.

En este punto nos centraremos en un tipo concreto de ecuaciones diferencia-les.

1.2. ESPACIO DE FASES 3

Definición 2 Llamamos ecuaciones diferenciales autónomas a aquellas que sonde la forma

x = f(x), x ∈ Rn, (1.6)

donde la variable independiente no aparece explícitamente en el segundo miem-bro de la ecuación.

La solución de estas ecuaciones, cuando disponemos de una condición inicialx(0) = x0, es una curva en R

n parametrizada en el tiempo, a la que llamaremoscurva integral, órbita o trayectoria a través de x0.

Definición 3 La curva (x1(t), . . . , xn(t)) ∈ Rn es una curva integral de la ecua-

ción (1.6) si (x1(t), . . . , xn(t)) = f(x1(t), . . . , xn(t)) ∀t ∈ R.

Es evidente que una curva integral es solución de la ecuación diferencial.

Si tenemos un sistema de n ecuaciones diferenciales autónomas, x = f(x),con f : Rn −→ R

n, dicho sistema puede escribirse componente a componente:

xi = fi(x1, . . . , xn), i = 1, . . . , n con x = (x1, . . . , xn). (1.7)

Así, una curva integral se puede encontrar haciendo un cambio de varia-ble independiente. Si tomamos x1 como nueva variable, el sistema (1.7) seríaequivalente a:

dxk

dx1=

fk(x1, . . . , xn)f1(x1, . . . , xn)

k = 2, 3, . . . , n (1.8)

Nótese que es posible obtener la ecuación de las curvas integrales siemprey cuando ∃k ∈ {1, .., n} tal que fk(x1, ...., xn) �= 0. Si esto no sucede, entoncestenemos una solución de equilibrio.

Definición 4 x ∈ Rn es un punto estacionario o de equilibrio para un flujo

ϕ ⇐⇒ ϕ(x, t) = x, ∀t ∈ R. Luego las soluciones de f(x) = 0 dan lugar a lospuntos estacionarios.

Los conceptos introducidos en las definiciones anteriores quedan reflejadosen el siguiente ejemplo.

Ejemplo 3 Sea el sistema dinámico dado por la ecuación diferencial x = ax +bx3. Para los siguientes valores de los parámetros,

(a, b) = (1, 1), (−1, 1), (1, −1), (−1, −1),

vamos a dibujar el diagrama de fases, haciendo hincapié en los puntos de equi-librio.

En este caso tenemos que el espacio de fases es E = R2 ≡ {(x, x)|x, x ∈ R},

y podemos obtener las curvas integrales introduciendo una nueva variable y = x.Haciendo esto la ecuación diferencial se transforma en el sistema


�1.0 �0.5 0.0 0.5 1.0

�1.0

�0.5

0.0

0.5

1.0

�1.0 �0.5 0.0 0.5 1.0

�1.0

�0.5

0.0

0.5

1.0

Figura 1.1: Diagrama de fases para el sistema dinámico del Ejemplo 3, x =ax+bx3. A la izquierda el caso en que (a, b) = (1, 1), mientras que en la derecha(a, b) = (−1, −1)

.

x = y, y = ax + bx3, (1.9)

y las óribtas quedan determinadas por la siguiente ecuación

dx

dy=

y

ax + bx3 ⇒ (ax + bx3)dx = ydy ⇒ ax2

2+

bx4

4− y2

2= C (1.10)

Es evidente que diferentes valores de a y b van a dar lugar a diferentesfamilias de curvas integrales.

Para el caso (a, b) = (1, 1) obtenemos las curvas integrales:

x2

2+

x4

4− y2

2= C, C ∈ R. (1.11)

Es claro que, para este caso, el único punto estacionario es x = 0, y = 0,ya que

x + x3 = 0 ⇒ x(1 + x2) = 0 ⇒ x = 0. (1.12)

Del mismo modo, obtenemos que también para el caso (a, b) = (−1, −1) elúnico punto estacionario es x = 0, y = 0, pues

− x − x3 = 0 ⇒ x(1 + x2) = 0 ⇒ x = 0. (1.13)

Para estos dos casos el diagrama de fases queda representado en la Figura1.1. Nótese que, aunque (0, 0) es el único punto crítico, las órbitas se comportande forma distinta en cada caso.

Pasamos a calcular las curvas integrales para los valores (a, b) = (−1, 1), queresultan ser

1.2. ESPACIO DE FASES 5

�1.0 �0.5 0.0 0.5 1.0

�1.0

�0.5

0.0

0.5

1.0

�1.0 �0.5 0.0 0.5 1.0

�1.0

�0.5

0.0

0.5

1.0

Figura 1.2: Diagrama de fases para el sistema dinámico del Ejemplo 3, x =ax + bx3. A la izquierda el caso en que (a, b) = (−1, 1), mientras que en laderecha (a, b) = (1, −1)

.

− x2

2+

x4

4− y2

2= C. (1.14)

En este caso los puntos estacionarios del sistema dinámico son tres, x =0, ±1, y = 0, lo que se deduce de resolver

− x + x3 = 0 ⇒ x(−1 + x2) = 0 ⇒ x = 0, ±1. (1.15)

Del mismo modo, en el caso (a, b) = (1, −1), obtenemos los mismos puntosestacionarios, siendo ahora las curvas integrales

x2

2− x4

4− y2

2= C. (1.16)

Para estos dos casos el diagrama de fases queda representado en la Figura1.2. Nótese que, aunque de nuevo los dos sistemas dinámicos dispongan de losmismos puntos críticos, las órbitas se comportan de nuevo de forma distinta encada caso.

Como se aprecia en los ejemplos, los puntos de equilibrio condicionan, engran medida, el resto de las trayectorias. Estos puntos gozan de una importantepropiedad, como es la de ser invariantes respecto al flujo.

No obstante, no son los únicos subconjuntos del espacio de fases con estapropiedad. En este sentido, podemos dar la siguiente definición.

Definición 5 Un conjunto M ⊆ E es invariante ⇐⇒ ∀x ∈ M , ϕ(x, t) ∈M, ∀t ∈ R. Un conjunto es positivamente (negativamente) invariante si ∀x ∈ M ,ϕ(x, t) ∈ M, ∀t > 0 (∀t < 0).

Obviamente, las órbitas son conjuntos invariantes. No obstante existen otrosconjuntos invariantes no triviales.


Figura 1.3: Ejemplo de conjunto invariante para el sistema dinámico definidopor r = r(1 − r), θ = 1. En rojo la circunferencia de centro (0, 0) y radio 1, yen azul una corona circular y dos trayectorias que parten de la corona circular,y vemos como se mantienen en su interior.

Ejemplo 4 Sea el sistema de ecuaciones diferenciales, en coordenadas polares,dado por

r = r(1 − r), θ = 1. (1.17)Además de las órbitas, como la circunferencia de centro (0, 0) y radio 1 o el

origen, que corresponden a una órbita y a un punto crítico, respectivamente, elsistema (1.17) tiene otros conjuntos invariantes. Por ejemplo, cualquier coronacircular de radios a y b con a > 1 y b < 1 es positivamente invariante, ya que sir > 0 ⇒ r < 0 y si r < 1 ⇒ r > 0. Por tanto, cualquier solución que empiece enla corona se mantiene en la corona si t > 0. Este hecho se muestra en la Figura1.3 para un caso particular de corona circular con dichas características.

Parece tarea complicada encontrar conjuntos invariantes de un sistema diná-mico. El siguiente lema define algunas propiedades de éstos, que pueden servirpara obtener conjuntos invariantes a partir de otros ya conocidos.

Lema 1 : Sea M ⊂ Rn

(i) M es invariante ⇐⇒ Rn \ M es invariante.

(ii) Sea Mi una colección contable de subconjuntos invariantes de Rn. Entonces

∪iMi y ∩iMi son también invariantes.

1.3. Conjunto límite y trayectoriasPodemos ver que cualquier órbita del Ejemplo 4 se acerca a la circunferencia

r = 1 si t → ∞ y además dicha circunferencia es invariante. Vamos a generalizar

1.3. CONJUNTO LÍMITE Y TRAYECTORIAS 7

este hecho.

Supongamos que disponemos de la ecuación diferencial (1.6) que define elflujo ϕ(x, t), que asumimos que está definido ∀x ∈ R

n y ∀t ∈ R. En este casopodemos dar una definición rigurosa de algunos conceptos.

Definición 6 La trayectoria a través de x es el conjunto:

γ(x) =⋃t∈R

ϕ(x, t). (1.18)

Análogamente, la semi-trayectoria positiva γ+(x) y la semi-trayectoria ne-gativa γ−(x) vienen dadas por

γ+(x) =⋃t≥0

ϕ(x, t) , γ−(x) =⋃t≤0

ϕ(x, t) (1.19)

La trayectoria a través de x es precisamente la curva integral que pasa porx.

Como hemos visto las trayectorias pueden tener un límite, lo que motiva lasiguiente definición.

Definición 7 Dada una trayectoria γ(x) =⋃

t∈R

ϕ(x, t), el conjunto ω-límite de

x, Λ(x), y el conjunto α-límite, A(x), de la trayectoria que pasa por x, son losconjuntos:

Λ(x) = {y ∈ Rn|∃(tn), tn → ∞, y ϕ(x, tn) → y n → ∞}

A(x) = {y ∈ Rn|∃(sn), sn → −∞, y ϕ(x, sn) → y n → ∞}

(1.20)

En el Ejemplo 4 es fácil identificar los conjuntos límite. Así, la circunferenciade radio 1 actúa como conjunto ω-límite de todas las trayectorias excepto elorigen, mientras que el origen es el conjunto α-límite de todas las trayectoriasinteriores a la circunferencia unidad. Para las exteriores no hay conjunto α-límite. Como vemos, los conjuntos límite en este caso son invariantes. Esto esasí en general como consecuencia del siguiente teorema.

Teorema 1 [2] El conjunto Λ(x) es invariante, y si γ+(x) está acotado enton-ces Λ(x) es no vacío y compacto.

De manera análoga se cumple el Teorema 1 para el conjunto A(x).

Capítulo 2

Ecuaciones diferencialeslineales

Como ya hemos comentado en el Ejemplo 1, un sistema de n ecuacionesdiferenciales lineales x = Ax define un sistema dinámico en R

n. Notar que elorigen x = 0 es un punto estacionario del flujo, y además es único si | A |�= 0.

Se sabe, por la teoría de ecuaciones diferenciales, que realizando un cambiode coordenadas podemos transformar el sistema inicial en y = Jy con J matrizdel sistema en forma canónica de Jordan, lo que permite resolver el sistema demanera sencilla, a través de la exponencial de J .

No entraremos más en detalle en el cálculo de soluciones de este tipo deproblemas, ya que corresponde a la teoría de ecuaciones diferenciales. En elsiguiente apartado tan solo veremos las propiedades de los subespacios funda-mentales asociados a los valores propios de la matriz J .

2.1. Variedades invariantesRecordemos que si para una matriz A, λk es un valor propio de multiplicidad

nk, entonces el subespacio fundamental asociado a λk se define como:

Ek = {x ∈ Rn | (A − λkI)nk x = 0} (2.1)

Teorema 2 [1] Si los valores propios de una matriz A n × n de coeficientesreales son distintos, entonces R

n se descompone en suma directa de subespaciosuno y dos-dimensionales. Cada uno de estos subespacios es invariante bajo elflujo que define x = Ax.En el caso en que A tenga algún valor propio de multiplicidad mayor que uno, elsubespacio fundamental asociado a éste también será invariante, con dimensiónigual a su multiplicidad.

9

10 CAPÍTULO 2. ECUACIONES DIFERENCIALES LINEALES

A continuación distinguiremos 3 tipos de subespacios fundamentales diná-micamente importantes. Denotemos por (ui), 1 ≤ i ≤ ni, a los vectores propiosde A asociados a los valores propios con parte real positiva, (cj), 1 ≤ j ≤ nc, alos vectores propios de A asociados a los valores propios con parte real cero, y(sk), 1 ≤ k ≤ ns, a los vectores propios de A asociados a los valores propios conparte real negativa (ni + nc + ns = n).

Definición 8 El espacio inestable del punto estacionario x = 0, Eu(0), es elespacio invariante generado por los vectores propios ui.

El espacio centro del punto estacionario x = 0, Ec(0), es el espacio inva-riante generado por los vectores propios cj.

El espacio estable del punto estacionario x = 0, Es(0), es el espacio inva-riante generado por los vectores propios sk.

Nótese que los espacios Eu(0) y Es(0) pueden caracterizarse a partir delflujo.

Teorema 3 Dado un sistema lineal de ecuaciones diferenciales x = Ax, dondeA es una matriz real n × n, entonces:

Eu(0) = {x ∈ Rn | etAx → 0 cuando t → −∞}. (2.2)

Es(0) = {x ∈ Rn | etAx → 0 cuando t → ∞}. (2.3)

Ejemplo 5 Consideremos el sistema x = Ax donde la matriz A tiene dos va-lores propios complejos conjugados con parte real positiva (ρ ± iω), y un valorpropio real negativo, λ < 0. Entonces, en forma canónica, podemos escribir elsistema de la siguiente manera

⎡⎣x

yz

⎤⎦ =

⎡⎣ρ −ω 0

ω ρ 00 0 λ

⎤⎦

⎡⎣x

yz

⎤⎦ . (2.4)

Podemos calcular Es mediante la Definición 8

⎡⎣ρ − λ −ω 0

ω ρ − λ 00 0 0

⎤⎦

⎡⎣x

yz

⎤⎦ =⇒

{(ρ − λ)x − ωy = 0,ωx − (ρ − λ)y = 0 =⇒ x = y = 0. (2.5)

por lo tanto Es es el eje z. Del mismo modo calculamos el espacio inestable, Eu,con la diferencia de que en este caso estará generado por los vectores propiosasociados a los valores propios con parte real positiva, que, en este caso comohemos dicho se trata de dos valores complejos conjugados

2.1. VARIEDADES INVARIANTES 11

⎡⎣ρ − ρ − iω −ω 0

ω ρ − ρ − iω 00 0 λ − ρ − iω

⎤⎦

⎡⎣x

yz

⎤⎦ =⇒

⎧⎨⎩

−iωx − ωy = 0,ωx − iωy = 0,(λ − ρ − iω)z = 0,

=⇒ z = 0.

(2.6)

Concluimos que el espacio inestable Eu queda determinado por el plano z =0. También podemos caracterizar los espacios estable e inestable mediante elTeorema 3. Sabemos que la solución del sistema es de la forma

⎧⎨⎩

x(t) = x0eρt cos(ωt) − y0eρt sen(ωt),y(t) = x0eρt sen(ωt) + y0eρt cos(ωt),z(t) = z0eλt.

(2.7)

Así, a partir del Teorema 3, obtenemos que Eu(0) viene dado por

Eu(0) = {(x0, y0, z0) ∈ R3 | lım

t→−∞(x(t), y(t), z(t)) = (0, 0, 0)} ={(x0, y0, z0) ∈ R

3 | z0 = 0}.(2.8)

Del mismo modo obtenemos el espacio estable

Es(0) = {(x0, y0, z0) ∈ R3 | lım

t→∞(x(t), y(t), z(t)) = (0, 0, 0)} ={(x0, y0, z0) ∈ R

3 | x0 = y0 = 0}.(2.9)

La importancia de estos espacios radica en que cualquier trayectoria se puedever como suma de una órbita del espacio estable y otra del inestable como semuestra en la Figura 2.1. Así en el espacio estable las órbitas son espiralesque se acercan al origen en sentido horario, mientras que en el inestable sontrayectorias rectilíneas que se alejan del origen. La suma de los dos movimientoses una espiral que se va enroscando en el eje z al mismo tiempo que se aleja delorigen.

Así vemos que estos espacios tienen una gran importancia, ya que caracte-rizan todo el movimiento del sistema (consecuencia de que como ya sabemos deteoría de ecuaciones diferenciales, R2 = Es ⊕ Eu). Además separan el espaciode fases en diferentes zonas, de manera que no se puede pasar de una a otra.En este ejemplo los conjuntos

A = {(x, y, z) ∈ R3|z > 0}, B = {(x, y, z) ∈ R

3|z < 0}, (2.10)

quedan separados por Es(0).


Figura 2.1: Representación de una órbita para el Ejemplo 5. La segunda imagencorresponde a la órbita proyectada en el plano xy

.

2.2. Linealización de sistemas de ecuaciones di-ferenciales no lineales

Consideremos ahora una ecuación diferencial no lineal x = f(x), x ∈ Rn, con

f suficientemente diferenciable, que tiene un punto estacionario x0. Mediante unsimple cambio de coordenadas podemos conseguir trasladar el punto estacionarioal origen, por lo que asumimos que f(0) = 0. Desarrollando f(x) en serie deTaylor en un entorno del origen se tiene

x = Df(0)x + O(|x|2), (2.11)donde Df(0) es la matriz jacobiana de f evaluada en el origen; Df(0)i,j =∂fi

∂xj(0). Así, si ignoramos términos de orden |x|2 y mayores, obtenemos la ecua-

ción diferencial lineal:

x = Df(0)x, (2.12)que representa una aproximación de la ecuación diferencial original en un en-torno del origen. Podemos pensar que, en dicho entorno, (2.12) es una buenaaproximación de las soluciones de la ecuación diferencial. En este sentido seríabueno poder encontrar una equivalencia entre el comportamiento del sistemalineal y del sistema no lineal, equivalencia en el sentido de encontrar un cambiode coordenadas que transforme un sistema en otro. Esta cuestión depende delos valores propios de A.

2.2. LINEALIZACIÓN DE SISTEMAS NO LINEALES 13

Definición 9 Sean λ1, . . . , λn los valores propios de Df(0). Decimos que Df(0)es resonante si existen enteros no negativos m1, . . . , mn cumpliendo

∑k

mk ≥ 2

tales que:

(m, λ) =n∑

k=1mkλk = λs (2.13)

para algún s ∈ {1, 2, ...., n}. El valor |m| =n∑

k=1mk se llama orden de la reso-

nancia.

Esta definición de resonancia aparece en la tesis de Henri Poincaré, en laque también encontramos el siguiente resultado acerca de la linealización desistemas de ecuaciones diferenciales no lineales.

Teorema 4 (Linealización de Poincaré) [2] Supongamos que x = f(x), con funa función analítica en un entorno del origen, f(0) = 0 y Df(0) no resonante.Supongamos también que Reλ > 0, ∀λ, o Reλ < 0, ∀λ, con λ valor propio deDf(0). Entonces existe un cambio de coordenadas y = x + · · · que transformala ecuación diferencial en y = Df(0)y.

La demostración se reduce a probar que existe una sucesión de cambios decoordenadas, en los que se van eliminando, con cada cambio, los términos deorden j, con j = 2, 3, . . . respectivamente, y finalmente comprobar que la seriede potencias definida en el cambio de coordenadas converge.

Sin embargo, se puede dar un resultado de linealización más general queabarque el caso en el que Df(0) tenga valores propios con parte real negativay valores propios con parte real positiva. Cumplirá un papel determinante lasiguiente condición, dada por Siegel.

Definición 10 La n-tupla λ = (λ1, . . . , λn) satisface la condición de Siegel siexisten C > 0 y ν > 0 tales que

|λi − (m, λ)| ≥ C

|m|ν , ∀i = 1, . . . , n (2.14)

∀m = (m1, . . . , mn) donde mi son enteros no negativos con |m| =n∑

i=1mi ≥ 2

Uniendo los resultados de Poincaré y Siegel obtenemos el siguiente resultado:

Teorema 5 (Siegel) [2] Supongamos que x = f(x), con f una función analítica,f(0) = 0. Si todos los valores propios de Df(0) cumplen las condiciones delTeorema 4, o si todos los valores propios de Df(0) cumplen la condición deSiegel, entonces existe un cambio de coordenadas y = x + · · · que transforma laecuación diferencial en y = Df(0)y, y además la serie de potencias definida enel cambio de coordenadas converge en un entorno del origen.


0.0002 0.0004 0.0006 0.0008 0.0010

�6.�10�6

�5.�10�6

�4.�10�6

�3.�10�6

�2.�10�6

�1.�10�6

1.�10�6

Figura 2.2: Representación de las funciones y = x2 e y = x2(Log|x| + k) en unentorno del origen.

Ejemplo 6 Para aplicar tanto el Teorema 4, como el Teorema 5, necesitamosanaliticidad del sistema de ecuaciones diferenciales y también de su solución. Esdecir, la solución admite un desarrollo en serie de potencias en un entorno delpunto de equilibrio. La presencia de resonancias puede hacer que la condiciónde analiticidad de la solución se pierda. Veámoslo en el siguiente ejemplo. Seael sistema de ecuaciones diferenciales

x = x , y = 2y + x2. (2.15)

Tenemos que la linealización en el origen viene dada por x = x, y = 2y. Losvalores propios de la linealización son (λ1, λ2) = (1, 2), que son resonantes deorden dos. La solución de la ecuación lineal es:

dy

dx=

2y

x=⇒ y = kx2 (2.16)

Sin embargo, es fácil ver que la ecuación de las órbitas, para el sistema(2.15),viene dada por

y = x2(log|x| + k). (2.17)

que no es analítica en un entorno del origen.En la Figura 2.2 se muestra que aún acercándonos al punto estacionario, la

solución del sistema original y la solución del sistema lineal difieren considera-blemente.

Todos estos teoremas son muy útiles, debido a que caracterizar los sistemascerca de los puntos de equilibrio puede ayudar a comprender el comportamientoglobal del sistema.

2.3. Resultados sobre la existencia de varieda-des invariantes

Si únicamente nos interesan propiedades topológicas del flujo, podemos sus-tituir las condiciones del teorema anterior por otras más débiles.

2.3. RESULTADOS DE EXISTENCIA 15

Definición 11 Un punto estacionario x es hiperbólico si Df(x) no tiene valorespropios con parte real 0.

Definición 12 Sea U un entorno del punto estacionario x. Análogamente a ladefinición de espacio estable (inestable) para sistemas lineales, podemos definirla variedad estable local de x, W s

loc(x), y la variedad inestable local de x, W uloc(x),

como

W sloc(x) = {y ∈ U |ϕ(x, t) → x, t → ∞, ϕ(y, t) ∈ U, ∀t ≥ 0} (2.18)

W uloc(x) = {y ∈ U |ϕ(x, t) → x, t → −∞, ϕ(y, t) ∈ U, ∀t ≤ 0} (2.19)

Definición 13 Supongamos que x0 es un punto estacionario hiperbólico de x =f(x), y Df(x0) denota la matriz Jacobiana de f evaluada en x0. Entonces x0 esun sumidero si todo valor propio de Df(x0) tiene parte real negativa (W u

loc(x0) =∅); y fuente si todo valor propio de Df(x0) tiene parte real positiva(W s

loc(x0) =∅). En otro caso se trata de una silla.

La existencia de estas variedades se puede asegurar gracias al siguiente teo-rema:

Teorema 6 (Variedad estable) [2] Supongamos que el origen es un punto esta-cionario hiperbólico de x = f(x) y Es y Eu son los espacios estable e inestabledel sistema lineal x = Df(0)x. Entonces existen variedades locales estable einestable, W s

loc(0) y W uloc(0), de la misma dimensión que Es y Eu respectiva-

mente. Éstas además son (respectivamente) tangentes a Es y Eu en el origen ytan suaves como la función original f .

Notar que para puntos estacionarios hiperbólicos, la variedad central es vacíay que el comportamiento local de las soluciones se mantiene, lo que no se puedeasegurar en otro caso.

La demostración del Teorema 6 proporciona una manera de obtener las va-riedades estable e inestable mediante desarrollos en serie de potencias al estilode Poincaré. En este sentido, consideramos una ecuación diferencial x = f(x)tal que x = 0 es un punto estacionario hiperbólico.Si elegimos el sistema de coordenadas para que la parte lineal de la ecuacióndiferencial en el origen esté en forma canónica, siempre podemos organizarlapara escribirla de la forma:

x = Ax + g1(x, y) y = −By + g2(x, y) (2.20)

donde x ∈ Rnu e y ∈ Rns (con nu la dimensión de la variedad inestable localy ns de la estable; nu + ns = n) y las matrices A y B tienen valores propioscon parte real positiva. Las funciones gi(x, y), i = 1, 2, contienen las partes nolineales de la ecuación.


Despreciando, inicialmente, la contribución de g1 y g2 tenemos, asociado alorigen, los espacios estable e inestable, dados por

Es(0, 0) = {(x, y)|x = 0} Eu(0, 0) = {(x, y)|y = 0} (2.21)

Como Eu, Es son tangentes a W uloc y W s

loc, podemos caracterizar W sloc a

partir de

xi = Si(y),∂Si

∂yj(0) = 0, 1 ≤ i ≤ nu, 1 ≤ j ≤ ns, (2.22)

donde está claro que la variedad es tangente a Es en el origen. Del mismo modopodemos describir la variedad inestable como:

yj = Uj(x),∂Uj

∂xj(0) = 0, 1 ≤ i ≤ nu, 1 ≤ j ≤ ns. (2.23)

Esto nos permite aproximar variedades estables e inestables desarrollandolas funciones Si y Uj en serie de potencias al estilo de Poincaré en el Teorema4. Tomando, por ejemplo, el desarrollo de Uj , tenemos

Uj(x) =∑r≥2

∑m∈Mr

Umj xm, (2.24)

donde Mr es simplemente un vector de índices, de tal modo que la suma detodos ellos es igual a r, y m ∈ Z

n, xm = xm11 · · · x

mnsns

Si B es diagonal con valores propios (λi), i = 1, . . . , ns entonces:

yj = −λjyj + g2j(x, y) (2.25)

y en la variedad inestable:

yj = −λjUj(x) + g2j(x, U(x)). (2.26)

Por otro lado tenemos también que:

yj =d

dtUj(x) =

nu∑k=1

xk∂

∂xkUj(x). (2.27)

Igualando las expresiones en (2.26) y (2.27) tenemos

− λjUj(x) + g2j(x, U(x)) =nu∑

k=1xk

∂

∂xkUj(x). (2.28)

Finalmente, si sustituimos (2.24) en (2.28) e igualamos coeficientes obtene-mos un conjunto de ecuaciones lineales que determinan los coeficientes Umi

.

2.3. RESULTADOS DE EXISTENCIA 17

Ejemplo 7 Sea el sistema de ecuaciones diferenciales

x = x , y = −y + x2. (2.29)

El único punto estacionario de este sistema es el origen (x, y) = (0, 0).

El sistema linealizado es x = x, y = −y, que tiene un punto de silla en elorigen con subespacios invariantes:

Es(0, 0) = {(x, y)|x = 0}, Eu(0, 0) = {(x, y)|y = 0}. (2.30)

Por el teorema de la variedad estable, sabemos que el sistema no lineal tieneuna variedad de la forma:

y = U(x),∂U

∂x(0) = 0. (2.31)

Buscamos un desarrollo en serie de potencias para U, U(x) =∑

k≥2Uix

i. Así:

y = −y + x2 = −∑k≥2

Uixi + x2 (2.32)

y en la variedad inestable:

y = x∂U

∂x(x) =

∑k≥2

kUkxk. (2.33)

Igualando los coeficientes de las diferentes potencias de x, se tienen las ecua-ciones

− u2 + 1 = 2u2, −uk = kuk, k ≥ 3. (2.34)

Así u2 = 13 , uk = 0, ∀k ≥ 3 y por tanto:

W uloc(0, 0) = {(x, y)|y =

13

x2}. (2.35)

Del mismo modo se obtiene que W sloc(0, 0) = Es(0, 0).

El próximo paso es extender las variedades locales para obtener variedadesestables e inestables globales, definidas por:

W u(0) =⋃t≥0

ϕ(W uloc(0), t), W s(0) =

⋃t≤0

ϕ(W sloc(0), t). (2.36)

Nótese que en el Ejemplo 7 se tiene que W u = W uloc y W s = W s

loc. Paraverlo basta con probar que las ecuaciones que hemos encontrado correspondena variedades invariantes. Así, por ejemplo, en el caso de la variedad inestableresulta


y =13

x2 =⇒ y =23

xx (y = −y + x2, x = x) =⇒

−y + x2 =23

x2 =⇒ −13

x2 + x2 =23

x2,

(2.37)

lo que demuestra la invariancia de la curva. También podemos ver que la va-riedad inestable local es la variedad inestable global mediante la resolución delsistema. En efecto,

{x = x =⇒ x = Aet,y = −y + x2 =⇒ y = −y + A2e2t.

(2.38)

Usando ahora el método de coeficientes indeterminados, obtenemos la solu-ción general del sistema de ecuaciones diferenciales

⎧⎨⎩

x = Aet

y = Be−t +A2

3e2t.

(2.39)

Imponiendo condiciones iniciales cuando t = 0, x = x0 e y = y0, la solucióndel problema de valores iniciales será

⎧⎨⎩

x = x0et

y = (y0 − x20

3)e−t +

x20

3e2t.

(2.40)

Finalmente aplicando la Definición 12, obtenemos la ecuación de la variedadinestable

lımt→−∞(x(t), y(t)) = (0, 0) =⇒ y0 − x2

03

= 0. (2.41)

Y del mismo modo obtenemos la variedad estable

lımt→∞(x(t), y(t)) = (0, 0) =⇒ x0 = 0. (2.42)

2.4. Persistencia de puntos estacionarios hiper-bólicos

Otra característica importante sobre los puntos estacionarios hiperbólicoses el hecho de que persisten bajo pequeñas perturbaciones de la ecuación di-ferencial. Así, si el origen es un punto estacionario hiperbólico de x = f(x) yv es un campo vectorial suave en R

n, y ε suficientemente pequeño, entoncesx = f(x) + εv(x) tiene un punto estacionario hiperbólico, cerca del origen, delmismo tipo que el de la ecuación no perturbada.

Ahora bien, tendremos que ver cuándo podemos asegurar que las perturba-ciones son suficientemente pequeñas para que los flujos de la ecuación original

2.4. PERSISTENCIA DE PUNTOS HIPERBÓLICOS 19

y la perturbada sean equivalentes. Dependerá de la noción de equivalencia y deltipo de perturbaciones que se permitan. Estas cuestiones se alejan de nuestroestudio y solo enunciaremos el siguiente teorema:

Teorema 7 (Grobman-Hartman): [3] Sea f : Rn −→ Rn un campo vectorial

C1, y sea x0 un punto de equilibrio. Supongamos que A = Df(x0) es un operadorhiperbólico (i.e. todos sus valores propios tienen parte real no nula). Entoncesexisten V, U entornos abiertos de x0 y de 0 respectivamente tales que los camposf|V y A|U son topológicamente conjugados (i.e existe un homeomorfismo ϕ :U −→ V tal que ϕ ◦ etA = Φt ◦ ϕ, donde Φ es el flujo definido por f).

Capítulo 3

Cálculo numérico de lavariedad estable e inestable

El Ejemplo 7 es un ejemplo académico, de manera que la variedad establee inestable pueden calcularse de forma exacta. Pero, en general, los métodosanalíticos, como el dado para el cálculo de W s

loc y W uloc son limitados y no

permiten siempre calcular W s y W u. Por este motivo es necesario recurrir a losmétodos numéricos.

3.1. Sistemas de ecuaciones diferenciales de di-mensión dos

Empezamos considerando el caso más sencillo, como es el caso de la deter-minación de las variedades estable e inestable para un sistema de ecuacionesdiferenciales 2-dimensional. De este modo introducimos el procedimiento que sesigue a la hora de realizar los cálculos.

Sea el sistema del Ejemplo 7:{

x = x,y = −y + x2.

(3.1)

En la Figura 3.1 se muestra el flujo definido por este sistema de ecuacionesdiferenciales, donde se ve que el origen es un punto de silla, que tiene asociadaslas variedades estable e inestable.

El cálculo de estas variedades empieza por obtener los correspondientes Eu

y Es. Calculamos los valores y vectores propios del sistema lineal asociado alpunto (0, 0). Así, se tienen los valores propios −1, 1, y como vectores propiosasociados (0, 1) y (1, 0).Recordemos que la variedad inestable es tangente en el origen al espacio inestable

21

22 CAPÍTULO 3. CÁLCULO NUMÉRICO DE VARIEDADES

Figura 3.1: Flujo generado por el sistema de ecuaciones diferenciales del Ejemplo7.

Eu(0) (Teorema 6), que está generada por los vectores propios asociados a losvalores propios con parte real positiva. Así la variedad inestable será tangenteal eje x. Por el mismo razonamiento, la variedad estable será tangente al eje y.

Procedemos, en primer lugar, a calcular numéricamente la variedad inestable.La idea es tomar dos condiciones iniciales sobre el espacio inestable (una encada sentido de la dirección definida por el vector (0, 1), (x0, y0) = ε(1, 0) y(x0, y0) = −ε(1, 0), ε > 0), e integrar numéricamente con tiempo positivo ambascondiciones iniciales, con lo que se espera obtener una buena aproximación dela variedad inestable, tanto mejor cuanto menor sea el valor de ε.

Esto nos da una curva que, en la Figura 3.2, comparamos con la variedadinestable exacta del problema, calculada en el Ejemplo 7. Como vemos, la aproxi-mación numérica de la variedad inestable se ajusta realmente bien a la variedadinestable real. Esto es así como consecuencia de (2.36).

Actuando de modo similar, podemos calcular la variedad estable. Nos si-tuamos sobre la dirección del vector (0, 1), muy cerca del origen, e integramosnuméricamente hacia atrás dos condiciones iniciales (una en cada sentido). Eneste caso, como se aprecia en la Figura 3.3, la variedad estable coincide con elespacio estable. Visto el caso de sistemas de dimensión 2, pasemos a analizarlos de dimensión mayor. El paso natural es ir a dimensión 3. A priori, la ideadesarrollada para el cálculo de variedades de dimensión uno parece que deberíapoder generalizarse a dimensión mayor, sin mayor inconveniente. No obstante,

3.2. SISTEMAS 3−DIMENSIONALES 23

�5 5

5

10

15

20

�4 �2 2 4 6 8

5

10

15

20

Figura 3.2: En la izquierda tenemos el cálculo numérico de la variedad inestabledel Ejemplo 7. En la derecha tenemos la comparación de ésta con la variedadinestable exacta y = x2

3 (en azul).

veremos que las cosas no son tan sencillas y tendremos que tener en cuentadiferentes aspectos. Para ello formulamos una generalización del problema delEjemplo 7 a dimensión 3.

3.2. Sistemas de ecuaciones diferenciales de di-mensión 3

Consideremos el sistema dinámico dado por las ecuaciones diferenciales⎧⎨⎩

x = x,y = αy,z = −z + 3x2 + (1 + 2α)y2,

(3.2)

donde α es un parámetro tal que α > 0. El único punto crítico de este sistemaes el origen de coordenadas.

Calculando los valores y vectores propios del sistema lineal asociados al punto(0, 0, 0), obtenemos como valores propios −1, 1 y α, y como vectores propiosasociados (0, 0, 1), (1, 0, 0) y (0, 1, 0). Por tanto tenemos que Es ≡ x = y = 0(eje z), y Eu ≡ z = 0 (el plano xy).

Podemos determinar la variedad inestable local mediante desarrollos en seriede potencias, aplicando el Teorema de la variedad estable. De este modo comola variedad inestable es tangente al espacio Eu, podemos escribirla como:


Figura 3.3: Flujo generado por el sistema de ecuaciones diferenciales del Ejemplo7 junto con la variedad inestable en el origen (en rojo), y la variedad estable (enverde).

z = a20x2 + a11xy + a02y2 + a30x3 + a21x2y + a12xy2 + a03y3 + . . . (3.3)

En el primer apartado del Apéndice se muestran el resto de cálculos necesa-rios para obtener la variedad inestable local. Vemos que ésta es

z = x2 + y2, (3.4)

independientemente del parámetro α.

Pasemos a calcular numéricamente la variedad inestable del sistema. Estaserá 2 dimensional, mientras que la variedad estable será 1-dimensional, por loque no mostraremos su cálculo numérico debido a que es semejante a lo hechoen la sección anterior. Calculemos la variedad inestable para el caso α = 1. Elprocedimiento a seguir es una generalización del caso 1-dimensional:

1. Tomamos una circunferencia sobre Eu de centro el punto crítico y conradio infinitesimal

{ε cos θ, ε sin θ, 0}, (3.5)

y, sobre ella, un conjunto de N condiciones iniciales.

{ε cos2kπ

5N, ε sin

2kπ

N, 0} k = 0, . . . , N (3.6)


Mi

Mi+1

ri,j r

i,j+1

ri+1,j r

i+1,j+1

Figura 3.4: Esquema de la construcción de la malla que aproxima a la variedadinestable.

(el que el primer y último punto sean iguales es por simplificar la progra-mación).

2. Integramos numéricamente las ecuaciones diferenciales para cada una delas ecuaciones iniciales escogidas sobre la circunferencia un tiempo Δtprefijado. Así, obtenemos un nuevo conjunto de puntos situados sobre unanueva circunferencia. Uniendo los puntos de la circunferencia inicial y losde la nueva vamos dando textura a la superficie solución.Esto se hace de la siguiente forma. Supongamos que disponemos de unacircunferencia Mi y en ella dos puntos consecutivos ri,j , ri,j+1. Por otrolado sea Mi+1 la circunferencia obtenida al integrar los puntos de Mi, yri+1,j , ri+1,j+1 los puntos de Mi+1 obtenidos mediante integración de lospuntos mencionados anteriormente de la circunferencia Mi. Así, medianteestos 4 puntos creamos un polígono como se muestra la Figura 3.4.

3. El punto 2 se repite, partiendo de los nuevos puntos obtenidos hasta al-canzar un tiempo t prefijado. De este modo se va creando una malla queaproxima a la variedad inestable. Esto puede verse en la Figura 3.5, dondese representa la primera etapa del proceso.

La implementación de este algoritmo se ha hecho con Mathematica y elcódigo del mismo puede verse en el segundo apartado del Apéndice. Para elcaso del sistema (3.2), con α = 1, hemos tomado ε = 0,01 y N = 50 y unaseparación entre circunferencias correspondiente a t = 0,1. Para ver la bondadde la aproximación obtenida, en la Figura 3.6 hemos representado a la izquierdala variedad inestable, dada por la ecuación (3.4) y en la parte derecha esta mismasuperficie junto con la aproximación numérica de la misma. Puede verse que elajuste es bastante bueno. De este modo, podemos pensar que este algoritmo,que es la generalización natural del dado para variedades de dimensión uno, vaa funcionar para aproximar variedades de dimensión mayor.

Sin embargo vamos a ver a continuación que este algoritmo no es adecuadopara calcular numéricamente las variedades estables e inestables de cualquiersistema de ecuaciones diferenciales.


Figura 3.5: Primer paso en el algoritmo de obtención de la variedad inestable delsistema (3.2), en el que se aprecian la circunferencia inicial y la circunferenciaformada por la integración de las condiciones iniciales , y los rectángulos queforman la malla.

Figura 3.6: A la izquierda se muestra la variedad inestable del sistema (3.2). Ala derecha, la variedad inestable junto con la aproximación numérica.


Figura 3.7: En la izquierda tenemos la variedad inestable aproximada numéri-camente para el sistema (3.2), con α = 1, 1. A la derecha, ésta junto con lavariedad inestable.

Calculemos ahora la variedad inestable para el caso α = 1, 1. Como hemosvisto, la variedad inestable no depende del valor del parámetro α, luego será lamisma que hemos obtenido antes z = x2 + y2. Sin embargo, al intentar obteneruna aproximación numérica utilizando el mismo método vemos, en la Figura 3.7,cómo, aunque el ajuste entre la aproximación numérica y la variedad inestablees buena, la construcción numérica se deforma. Esto es debido a que la variedadno crece uniformemente en todas las direcciones, ya que los valores propiosasociados al sistema de ecuaciones diferenciales no son iguales, 1 y 1,1. Cuantomayor sea la diferencia entre los valores propios, mayor será la diferencia decrecimiento en las distintas direcciones.

Podemos dar una solución a este problema forzando a la variedad a creceruniformemente en todas las direcciones. Conseguimos esto pasando a resolverun problema de contorno. La idea es que la separación entre dos circunferenciasno venga dada por un tiempo fijo, sino que la separación entre ambas sea unadistancia prefijada. Así, integrando a partir del punto (xi, yi, zi) exigiremos queel punto (xi+1, yi+1, zi+1), cumpla

√x2

i+1 + y2i+1 + z2

i+1 = id, (3.7)

siendo d una distancia prefijada.De este modo la superficie crecerá uniformemente en todas las direcciones. El


Figura 3.8: En la izquierda tenemos la variedad inestable aproximada numérica-mente mediante el segundo algoritmo, cuando α = 1,1. A la derecha, ésta juntocon la variedad inestable.

código de este segundo algoritmo aparece en el tercer apartado del Apéndice. Esuna sencilla adaptación del primer algoritmo en el que se han usado las opcionesde la función NDSolve para transformar el problema de valores iniciales en unode contorno. Para ello se ha usado la opción EventLocator, gracias a la cualse asegura que la distancia del nuevo punto al origen se corresponde con ladistancia prefijada.

En la Figura 3.8 vemos que, para el caso en el que α = 1, 1, el problemaqueda solventado.

Esta solución vemos que es buena en este caso, pero no necesariamente loserá siempre. En la Figura 3.9 mostramos la aproximación numérica de la va-riedad inestable del sistema (3.2) cuando α = 2. El cociente entre los valorespropios es mayor, y lo que sucede es que, aunque los puntos estén uniforme-mente distribuidos en la circunferencia inicial, las trayectorias se concentran enla dirección del vector propio asociado al valor propio de mayor módulo, porlo que los polígonos que forman la malla de la variedad inestable no son uni-formes, y así partes de la malla están excesivamente densas, y por el contrariootras apenas disponen de puntos.

Este hecho se ve plasmado en la Figura 3.10, en donde se muestra el flujode un sistema de ecuaciones diferenciales proyectado en el plano xy [4]. Dichosistema tiene valores propios 1, 10 y -1, y debido a esto vemos cómo al integrarciertas condiciones iniciales, éstas se propagan más en una dirección, que es laque corresponde al vector propio asociado al valor propio de mayor magnitud.

La forma de solucionar este problema es proporcionando un algoritmo deinserción y eliminación de puntos, de modo que, a partir de una circunferencia


Figura 3.9: Variedad inestable aproximada numéricamente mediante el segundoalgoritmo, con α = 2.

con puntos igualmente espaciados, obtengamos otra circunferencia en la que lospuntos estén situados de modo uniforme, con el fin de obtener una malla losuficientemente buena.

De este modo, se propone un tercer algoritmo donde se define una horquillade distancias entre dos puntos de una circunferencia dada por dos valores δ yΔ, de manera que la distancia entre cada dos puntos verifique

δ < dist(ri,j , ri,j+1) < Δ (3.8)

Esta distancia no es exactamente la distancia euclídea, debe ser una distanciamedida sobre la circunferencia y normalizada, para que no dependa del radio otamaño de la misma. Por ejemplo, aquí hemos usado

dist(ri,j , ri,j+1) =

√√√√√⎛⎝ xi,j√

x2i,j + y2

i,j

− xi,j+1√x2

i,j+1 + y2i,j+1

⎞⎠

2

+

⎛⎝ yi,j√

x2i,j + y2

i,j

− yi,j+1√x2

i,j+1 + y2i,j+1

⎞⎠

2

.

Para la implementación de este tercer algoritmo, se van calculando los puntosde la nueva circunferencia Mi+1 uno a uno a partir de los de la circunferenciaMi. Si dos puntos de la misma dan lugar a puntos a distancia inferior a δ sesuprime el segundo punto y si la distancia es mayor que Δ se inserta un nuevopunto, que se calcula a partir de dos punto de la circunferencia Mi, mediante


Figura 3.10: Flujo de un sistema de ecuaciones diferenciales proyectada en elplano xy. Los valores propios del sistema difieren considerablemente, y es poresto por lo que vemos que en una dirección se propagan muchas más condicionesiniciales.

un algoritmo de bisección. El código de este algoritmo podemos encontrarlo enel Apéndice, y cabe destacar que presenta mayor coste computacional que losanteriores.

Este método tiene un claro inconveniente, tenemos que ir moviendo el puntoescogido sobre la circunferencia Mi de modo que al integrarlo, este sea el óptimo,al tiempo que hay que llevar una ordenación paralela para saber cómo asociarlos puntos de Mi con los de Mi+1.

En la Figura 3.11 se muestra el cálculo de la variedad inestable del sistema(3.2), para el valor α = 2, utilizando el tercer algoritmo.

Este último algoritmo tiene en cuenta la mayoría de la problemática quepodemos encontrarnos para calcular las variedades estable e inestable.

3.3. Aproximación numérica mediante curvas geo-désicas de nivel

A pesar de todo hay un problema añadido, y es que no siempre, en unasuperficie, la distancia radial y la distancia mínima entre dos puntos, medidasobre la superficie, son iguales. En este sentido es conveniente introducir ladistancia geodésica entre dos puntos x e y de una superficie, que se representapor dg(x, y), como la longitud de arco del camino más corto en la superficieconectando x e y.

El tercer algoritmo dado en la sección anterior debería modificarse para teneren cuenta la distancia geodésica, en lugar de una distancia de tipo euclídeo. Estealgoritmo, denominado algoritmo de aproximación numérica mediante curvasgeodésicas de nivel, es el que encontramos en el trabajo de Krauskopf y Osinga[5] que es el que exponemos a continuación.

3.3. APROXIMACIÓN POR CURVAS GEODÉSICAS 31

Figura 3.11: En la izquierda tenemos la variedad inestable aproximada numéri-camente mediante el tercer algoritmo, para el valor α = 2. A la derecha, éstajunto con la variedad inestable.

Presentaremos el caso de una variedad inestable 2-dimensional de un puntohiperbólico, en un espacio 3-dimensional. Sin embargo el método se puede utili-zar para calcular variedades (in)estables k−dimensional en campos vectorialesde R

n.

Consideremos un sistema de ecuaciones diferenciales 3-dimensional, dondex0 ∈ R

3 es un punto hiperbólico. El método de las curvas geodésicas de nivelaproxima una variedad inestable global como una secuencia de círculos geodé-sicos (con círculo denotamos a cualquier curva suave y cerrada). Así la parame-trización geodésica de W u(x0) vendrá dada por

W u(x0) = {Sη}η>0, Sη = {x ∈ W u(x0)|dg(x, x0) = η}. (3.9)

El dato inicial es de nuevo un mallado uniforme M0 de un círculo geodésicopequeño Sη0 = Sδ ⊂ Eu(x0). El método consiste en un conjunto de pasos, talque para cada paso i se obtiene una nueva lista circular de puntos, Mi+1, queaproximará un nuevo conjunto de nivel Sηi+1 , de tal modo que se irá definien-do un mallado que consistirá en la variedad inestable. Los nuevos puntos delmallado se calculan resolviendo un problema de contorno apropiado.

Consideremos la tarea de encontrar Mi+1 a partir de la lista circular cono-cida Mi, que es una aproximación de Sηi

. Para ello previamente se da valor alparámetro Δi, que indicará la distancia que separa las distintas curvas geodé-sicas. La lista circular Mi+1 se construye punto a punto. Sea xi,j ∈ Mi; paraobtener el correspondiente punto en Mi+1 necesitaremos también los puntosxi,j−1, xi,j+1 ∈ Mi. De este modo, definimos, Fxi,j

, el plano que contiene al pun-to xi,j y cuyo vector normal viene dado por la media de los vectores −−−−−−→xi,j−1xi,j


Figura 3.12: Esquema de obtención de una variedad (in)estable de un sistemade ecuaciones diferenciales, mediante el método de las curvas geodésicas denivel. En la Figura r = xi,j , br(tr) denota el nuevo punto xi+1,j ∈ Mi+1, yqr(tr) denota el punto que hay que escoger sobre el circulo Mi, de modo que alintegrarlo obtengamos xi+1,j .

y −−−−−−→xi,jxi,j+1, que queda determinado de manera única. Este plano va a cortara la variedad inestable en una cierta curva. Ahora, para llegar a la siguientecurva geodésica, lo que hacemos es movernos a partir de xi,j sobre la circunfe-rencia Mi, buscando un punto sobre dicha curva geodésica de tal modo que latrayectoria que pase por él, corte al plano en un punto xi+1,j , de manera que ladistancia de este punto al punto xi,j sea la distancia prefijada Δi. Así, el puntoxi+1,j será un nuevo punto de Mi+1. Repitiendo esto con todos los puntos deMi obtendremos el nuevo círculo geodésico Mi+1. En la Figura 3.12 se muestraesta situación.

Podemos expresar la obtención de los nuevos puntos de una curva geodésicamediante el siguiente problema de contorno [6]:

φ0(qr(t)) = qr(t) ∈ Mi,

φt(qr(t)) = br(t) ∈ Fr,(3.10)

donde el tiempo de integración, t, es un parámetro libre, y φ denota el flujodefinido por el sistema de ecuaciones diferenciales.

Una vez que hemos encontrado todos los puntos candidatos de Mi+1 paraΔi, entonces hay que decidir si el tamaño de paso Δi es apropiado. Para ello secomprueba que la curvatura geodésica a través de todos los puntos de Mi no es

3.3. APROXIMACIÓN POR CURVAS GEODÉSICAS 33

demasiado grande. Esto se comprueba del siguiente modo; denotemos por αxi,j

el ángulo formado por la recta que pasa por xi,j y xi+1,j y la recta que pasa porxi,j y xi−1,j ∈ Mi−1. Así, diremos que Δi es aceptable si

αxi,j< αmax,

Δiαxi,j< (Δα)max,

(3.11)

para todo xi,j ∈ Mi. En este caso Mi+1 se acepta y el paso i finaliza. Si por elcontrario existe algún xi,j ∈ Mi de modo que no se cumple una de las desigual-dades de (3.11), entonces reducimos Δi a la mitad y repetimos el paso i con elnuevo valor de Δi. Del mismo modo Δi puede duplicarse si para todo xi,j ∈ Mi

se cumplen las dos condiciones dadas en (3.11). Nótese que tanto αmax como(Δα)max serán constantes prefijadas.

También es importante asegurar que Mi+1 es una buena aproximación deSηi+1 . En otras palabras, puntos vecinos de Mi+1 pueden estar demasiado cercao demasiado lejos entre ellos. Por tanto, del mismo modo que con la distanciaradial, necesitaremos algoritmos de inserción y eliminación de puntos del malla-do. Para ello se introduce una nueva constante ΔF . Cuando dos puntos vecinosde Mi proporcionan dos puntos vecinos en Mi+1 que están a distancia mayor deΔF , entonces se añade un nuevo punto entre ellos. La forma de añadir puntoses la misma que usábamos con la distancia radial, se añade un punto en Mi demodo que al integrarlo nos de un punto en Mi+1 que se encuentre a mitad dedistancia entre los puntos que estaban a una distancia mayor que ΔF . Con elfin de garantizar una relación de orden adecuada entre los puntos vecinos deMi+1, se elimina un punto si dos puntos vecinos en Mi+1 se encuentran a unadistancia menor que δF , que también es una distancia prefijada.

Este algoritmo ha sido implementado de manera satisfactoria para calcularvariedades invariantes, incluso en casos complejos, como es el caso del sistemade Lorenz [5]. En la Figura 3.13 se muestra la variedad estable asociada al origendel sistema de Lorenz.

Cabe mencionar que existen otros algoritmos diferentes para aproximar va-riedades invariantes, basados en otras ideas. Aquí nos hemos limitado a exponerun método, desarrollado a partir de una idea simple, como es el de construir lavariedad invariante a partir de la evolución de condiciones iniciales apropiadaspor medio del flujo. Para una revisión de los diferentes métodos y algoritmospuede consultarse la referencia [5].


Figura 3.13: Variedad estable del sistema de Lorenz.

Capítulo 4

Apendice

A continuación vamos a presentar algunos de los algoritmos que hemos uti-lizado a lo largo de la memoria para realizar los diferentes cálculos.

Algoritmo 1: Algoritmo general para calcular la variedad inestable localhasta orden n, siempre que la aproximación lineal sea z = 0.

Funciones auxiliares para introducir los coeficientes y el desarrollo en serie dela variable z.

In[1]:=SetAttributes[a, Constant]

coeficientes[n_] :=Flatten[Table[a[j, i - j], {i, 2, n}, {j, i, 0, -1}]]

serie2[n_] :=Apply[Plus,Flatten[Table[a[j, i - j] x^j y^(i - j), {i, 2, n}, {j, i, 0, -1}]]]

Función para calcular el desarrollo en serie de potencias de una función f(x, y)en torno a (0, 0) hasta un orden prefijado n.

In[2]:=desarrollo[f_, n_] :=Normal[Series[f /. {x -> e x, y -> e y}, {e, 0, n}]] /. e -> 1

Función para calcular el desarrollo en serie de la variedad inestable local hastaorden n.

In[3]:=inestable[eq1_, eq2_, eq3_, n_] :=serie2[n] /.Flatten[Solve[

LogicalExpand[

35

36 CAPÍTULO 4. APENDICE

Series[desarrollo[eq3 - Dt[serie2[n], t] /. {Dt[x, t] -> eq1,

Dt[y, t] -> eq2} /. z -> serie2[n], n], {x, 0, n}, {y, 0,n}] == 0], coeficientes[n]]]

Algoritmo 2: Código del primer algoritmo para obtener la variedadinestable local del sistema de ecuaciones diferenciales (48) numéri-camente.

Algoritmo basado en la resolución de problemas de valor inicial,mediante la función NDsolve. En cada paso se parte de una cir-cunferencia Mi y se construye otra Mi+1 = φt(Mi), donde t = 0,1en este caso. Es decir la circunferencia Mi+1 es la que se obtiene deMi cuando ha transcurrido un tiempo t = 0,1.

In[4]:=poligonos = {};puntos = Table[{.01 Cos[2 k Pi/50], .01 Sin[2 k Pi/50], 0}, {k, 0,

50}];For[i = 1, i <= 55, i = i + 1,puntos2 =Table[{x[t], y[t], z[t]} /.

Flatten[NDSolve[{x’[t] == x[t], y’[t] == y[t],z’[t] == -z[t] + 3 x[t]^2 + 3 y[t]^2, x[0] == puntos[[k, 1]],y[0] == puntos[[k, 2]], z[0] == puntos[[k, 3]]}, {x[t], y[t],z[t]}, {t, 0, .5}]] /. t -> .1, {k, 1, 51}];

AppendTo[poligonos,Table[Polygon[{puntos[[k]], puntos2[[k]], puntos2[[k + 1]],

puntos[[k + 1]]}], {k, 1, 50}]]; puntos = puntos2]

Algoritmo 3: Código Mathematica donde se solventa el problemadel crecimiento no uniforme de la superficie.

Algoritmo muy similar al anterior. En este caso se resuelve un pro-blema de contorno gracias a la instrucción EventLocator, un locali-zador de eventos.In[5]:=

poligonos3 = {};

37

h = 0.2;puntos = Table[{.01 Cos[2 k Pi/50], .01 Sin[2 k Pi/50], 0}, {k, 0,

50}];For[i = 1, i <= 30, i = i + 1, puntos2 = {};For[k = 1, k <= 51, k = k + 1,rr = Flatten[

Reap[NDSolve[{x’[t] == x[t], y’[t] == 1.1 y[t],z’[t] == -z[t] + 3 x[t]^2 + 3.2 y[t]^2, x[0] == puntos[[k, 1]],y[0] == puntos[[k, 2]], z[0] == puntos[[k, 3]]}, {x[t], y[t],

z[t]}, {t, 0, 3},Method -> {"EventLocator",

"Event" -> x[t]^2 + y[t]^2 + z[t]^2 - (i h)^2,"EventAction" :> Sow[t]}]]];

AppendTo[puntos2, ({x[t], y[t], z[t]} /. Rest[Reverse[rr]]) /.t -> Last[rr]]];

AppendTo[poligonos3,Table[Polygon[{puntos[[k]], puntos2[[k]], puntos2[[k + 1]],


Algoritmo 4: Código Mathematica donde se solventa el problemade la no uniformidad del mallado.

Como ya hemos visto anteriormente, hemos aprovechado en granparte las herramientas que dispone Mathematica, como por ejemplotodas las herramientas numéricas, como es el caso de NDsolve pararesolver problemas de valores iniciales y problemas de contorno. Aparte de esto, nosotros nos hemos definido algunas funciones au-xiliares con el fin de que la programación fuese más sencilla. Deeste modo, la siguiente función calcula a partir de un punto de lacircunferencia Mi, un punto de la circunferencia Mi+1.In[6]:=

iteracion[{x0_, y0_, z0_}, d_] :=Block[{pp},pp = Flatten[

Reap[NDSolve[{x’[t] == x[t], y’[t] == 2 y[t],z’[t] == -z[t] + 3 x[t]^2 + 5 y[t]^2, x[0] == x0, y[0] == y0,z[0] == z0}, {x[t], y[t], z[t]}, {t, 0, 3},

Method -> {"EventLocator","Event" -> Norm[{x[t], y[t], z[t]}] - d,"EventAction" :> Sow[t]}]]]; ({x[t], y[t], z[t]} /.

Rest[Reverse[pp]]) /. t -> Last[pp]]

Definimos la siguiente función auxiliar para calcular la distancianormalizada entre dos puntos.


In[7]:=distancia[a_, b_] :=Norm[Rest[Reverse[a]]/Norm[Rest[Reverse[a]]] -

Rest[Reverse[b]]/Norm[Rest[Reverse[b]]]]

Introducimos la colección inicial de puntos.In[8]:=

puntos = Table[{.01 Cos[2 k Pi/50], .01 Sin[2 k Pi/50], 0}, {k, 0,50}];

Difinimos una horquilla para controlar la distancia entre puntos.In[9]:=

Delta = distancia[puntos[[2]], puntos[[1]]] (1 + 10^(-3))delta = distancia[puntos[[2]], puntos[[1]]] (1 - 10^(-3))

Bloque de cálculo de la malla. Partiendo de las condiciones inicialesse va añadiendo punto a punto, los siguientes puntos de la malla. Secalcula un candidato y es aceptado si la distancia al punto anteriorestá en la horquilla. Si no, se inserta un punto a distancia adecuadao si no se pasa a un nuevo candidato.In[10]:=

poligonos4 = {};h = 0.2;puntos = Table[{.01 Sin[2 k Pi/50], .01 Cos[2 k Pi/50], 0}, {k, 0,

50}];For[i = 1, i <= 30, i = i + 1,aux = puntos;puntos2 = {};AppendTo[puntos2, iteracion[puntos[[1]], i h]];For[k = 2, k <= Length[aux], k = k + 1,ite = iteracion[aux[[k]], i h];If[distancia[ite, Last[puntos2]] > Delta, izda = aux[[k - 1]];dcha = aux[[k]];medio = (izda + dcha)/2;candidato =iteracion[medio, i h];

While[(distancia[candidato, Last[puntos2]] >Delta) || (distancia[candidato, Last[puntos2]] < delta),

If[distancia[candidato, Last[puntos2]] > Delta, dcha = medio,If[distancia[candidato, Last[puntos2]] < delta, izda = medio]];

medio = (izda + dcha)/2;candidato = iteracion[medio, i h]];

aux = Insert[aux, medio, k];AppendTo[puntos2, candidato],If[distancia[ite, Last[puntos2]] > delta, AppendTo[puntos2, ite]]]

]; If[Length[puntos2] < 51,For[l = 1, l <= 51 - Length[puntos2],AppendTo[puntos2, iteracion[Last[puntos], i h]]]];

AppendTo[poligonos4,Table[Polygon[{puntos[[k]], puntos2[[k]], puntos2[[k + 1]],


Conclusiones

Para la realización de esta Memoria he revisado los conceptos bá-sicos de los sistemas dinámicos continuos, con especial hincapié enlas variedades invariantes asociadas a los puntos críticos, y el modode calcularlas. Dicha memoria me ha aportado distintos aspectosrespecto a los estudios de la titulación. Ha supuesto familiarizarmecon algunos conceptos de la teoría cualitativa de ecuaciones diferen-ciales y de la dinámica no lineal, como por ejemplo el hecho de quealgunas estructuras dentro de las ecuaciones diferenciales puedenhacer comprender el comportamiento de las soluciones sin llegar acalcularlas. Además he podido ver la utilidad de los métodos numé-ricos aprendidos en la titulación.

En la parte de la memoria en la que nos hemos centrado en laultilización de métodos numéricos, hemos trabajado en un problemamuy concreto, poniendo de relieve que algunas ideas elementalespara problemas de una dimensión, cuando se intentan generalizara más dimensiones dejan de ser triviales. Así, hemos sido capacesde poner de manifiesto prácticamente todas las dificultades que nosencontramos a la hora de calcular las variedades invariantes.

39

Bibliografía

[1] G.F. Simmons, Ecuaciones Diferenciales con Aplicaciones y no-tas Históricas. McGraw-Hill, Madrid, 1988.

[2] P. Glendinning, Stability, Instability and Chaos. Cambridge Uni-versity Press, Cambridge, 1994.

[3] F. Dumortier, J. Llibre,J.C. Artés Qualitative Theory of PlanarDifferential Systems. Springer, New York, 2006.

[4] J. Guckenheimer y A. Vladimirsky, A fast method for approxi-mating invariant manifolds, SIAM Journal on Applied DynamicalSystems, 3, 232–260, 2004.

[5] B. Krauskopf, H.M. Osinga, E.J. Doedel, M.E. Henderson, J.Guckenheimer, A. Vladimirsky, M. Dellnitz y O. Jungle, A sur-vey of methods for computing (un)stable manifold of vectorfields, International Journal of Bifurcation and Chaos, 15, 763–791, 2005.

[6] B. Krauskopf y H.M. Osinga, Computing geodesic level sets onglobal (un)stable manifolds of vector fields, SIAM Journal onApplied Dynamical Systems, 2, 546–569, 2003.

41

Aproximación numérica de variedades invariantes en ... · vamos a dibujar el diagrama de fases,...

Documents

Transcript of Aproximación numérica de variedades invariantes en ... · vamos a dibujar el diagrama de fases,...