Simulación numérica de modelos de ecuaciones diferenciales en Neurociencia.

FACULTAD DE MATEMÁTICAS DEPARTAMENTO DE ECUACIONES DIFERENCIALES Y ANÁLISIS NUMÉRICO Trabajo Fin de Grado Simulación numérica de modelos de ecuaciones diferenciales en Neurociencia. Elena Pérez Castilla ___________________________ Dirigido por: Soledad Fernández García María Macarena Gómez Mármol 2020 2 Índice general Resumen I Introducción III 1. Modelos matemáticos clásicos en neurociencia 1 1.1. Conceptosprevios............................. 1 1.2. Modelos .................................. 2 1.2.1. El modelo de Hodgkin-Huxley . . . . . . . . . . . . . . . . . . 2 1.2.2. El sistema de FitzHugh-Nagumo . . . . . . . . . . . . . . . . 4 1.2.3. Sistemas tipo Bursting ...................... 4 2. Análisis Matemático 7 2.1. Resultadosteóricos............................ 7 2.1.1. Funciones lipschitzianas . . . . . . . . . . . . . . . . . . . . . 7 2.1.2. Existencia y unicidad de solución . . . . . . . . . . . . . . . . 8 2.1.3. Órbitas .............................. 10 2.1.4. Estabilidad ............................ 11 2.1.5. Bifurcaciones ........................... 13 2.2. Análisis del sistema de FitzHugh-Nagumo . . . . . . . . . . . . . . . 16 2.3. Análisis del sistema Hindmarsh-Rose . . . . . . . . . . . . . . . . . . 19 2.4. Análisis del sistema de FitzHugh-Rinzel . . . . . . . . . . . . . . . . . 21 3. Métodos numéricos 25 3.1. Resolución numérica de problemas ligados a EDO’s . . . . . . . . . . 25 3.2. Métodos de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.3. Método de Gauss de dos etapas . . . . . . . . . . . . . . . . . . . . . 31 3.4. MétodoTR-BDF2 ............................ 33 3.5. Método Rosenbrock ROS2 . . . . . . . . . . . . . . . . . . . . . . . . 34 3.6. Aplicación de los métodos . . . . . . . . . . . . . . . . . . . . . . . . 36 3.6.1. El sistema de FiztHugh-Nagumo . . . . . . . . . . . . . . . . 36 3.6.2. El sistema de Hindmarsh-Rose . . . . . . . . . . . . . . . . . . 38 3.6.3. El sistema de FitzHugh-Rinzel . . . . . . . . . . . . . . . . . . 39 4. Simulaciones 41 4.1. Testacadémico .............................. 42 4.2. Testaplicado ............................... 45 4.2.1. Sistema de FitzHugh-Nagumo . . . . . . . . . . . . . . . . . . 45 4.2.2. Sistema de Hindmarsh-Rose . . . . . . . . . . . . . . . . . . . 49 4.2.3. Sistema de FitzHugh-Rinzel . . . . . . . . . . . . . . . . . . . 56 3 4ÍNDICE GENERAL Códigos MATLAB 65 4.3. MétododeGauss............................. 65 4.4. MétodoTR-BDF2 ............................ 66 4.5. Método Rosenbrock ROS2 . . . . . . . . . . . . . . . . . . . . . . . . 66 Bibliografía 68 Resumen En este trabajo, consideraremos los sistemas de FitzHugh-Nagumo, Hindmarsh- Rose y FitzHugh-Rinzel, que modelan el comportamiento de las neuronas al recibir un impulso eléctrico. Los estudiaremos desde el punto de vista de los sistemas dinámicos y la teoría de bifurcaciones. Para el sistema de FitzHugh-Nagumo nos centraremos en la bifurcación de Hopf, mientras que para los sistemas de Hindmarsh-Rose y FitzHugh-Rinzel veremos las bifurcaciones que provocan la dinámica bursting de los mismos. Además, realizaremos simulaciones numéricas que nos permitirán comprobar los resultados teóricos obtenidos para estos sistemas. Para ello, utilizaremos distintos métodos numéricos para problemas stiff y los implementaremos en el programa MATLAB. Abstract In this work, we will consider the FitzHugh-Nagumo, Hindmarsh-Rose and FitzHugh-Rinzel systems, which model the behavior of neurons when they experience a change of electric charge. We will study them from the point of view of dynamical systems and bifurcation theory. For the FitzHugh-Nagumo system we will focus on the Hopf bifurcation, while for the Hindmarsh-Rose and FitzHugh- Rinzel systems we will focus on the bifurcations that cause bursting. Moreover, we perform numerical simulations that will allow us to check the theoretical results obtained for these systems. Thus, we will study different numerical methods for stiff problems and implement them in MATLAB program. i ii RESUMEN Introducción La biomatemática es la rama de la matemática aplicada que se encarga de crear y estudiar modelos que reproduzcan algún fenómeno biológico con el fin de comprender los mecanismos subyacentes. En este trabajo, trataremos la neurociencia desde el punto de vista de la matemática aplicada. La neurociencia es la ciencia que se ocupa de cada uno de los diversos aspectos y funciones especializadas del sistema nervioso. La compresión del funcionamiento del sistema nervioso es uno de los grandes retos de la comunidad científica actual. Su estudio se divide en distintos niveles de complejidad: neuronas aisladas (escala microscópica), conjuntos de neuronas (escala mesoscópica) y áreas completas del cerebro y la médula espinal (escala macroscópica). A escala microscópica, uno de los aspectos más relevantes es la sinapsis neuronal. Existen numerosos modelos matemáticos que tratan de explicar este proceso, dentro de ellos el más relevante es el modelo de Hodgkin-Huxley, realizado por Alan Lloyd Hodgkin y Andrew Fielding Huxley en 1952 [1], a los que se les concedió el premio Nobel de medicina en 1963. Debido a la complejidad de este sistema, a lo largo de los años han surgido varias simplificaciones del mismo. Dentro de estos sistemas destaca el sistema de FitzHugh-Nagumo, que mantiene las características del sistema original. Asímismo, encontramos sistemas tipo bursting, como el sistema de Hindmarsh-Rose y el sistema de FitzHugh-Rinzel. Estos sistemas son lento-rápidos y es difícil encontrar soluciones exactas de los mismos. Sin embargo, mediante la resolución numérica se pueden hallar aproximaciones de la mismas. Para comenzar el trabajo debemos introducir algunos conceptos básicos sobre neurociencia, esto lo haremos en el Capítulo 1, en el que también presentaremos cada uno de los sistemas mencionados con más detenimiento. En el Capítulo 2, enunciaremos los resultados teóricos necesarios para llevar a cabo el estudio de estos sistemas, que enfocaremos desde el punto de vista de los sistemas dinámicos y la teoría de bifurcaciones. Para el sistema de FitzHugh-Nagumo nos centraremos en la bifurcación de Hopf, y para los sistemas tipo bursting en las bifurcaciones responsables de este fenómeno. En el Capítulo 3, se estudiarán distintos métodos numéricos para problemas stiff. Comenzaremos estudiando los métodos de Runge-Kutta, centrándonos después en los métodos de Gauss, TR-BDF2 y Rosenbrock ROS2. En el Capítulo 4 realizaremos distintos tipos de simulaciones, en primer lugar, iii iv INTRODUCCIÓN realizaremos un test académico que nos permitirá comprobar los resultados numéricos de los métodos estudiados en el Capítulo 3. Una vez comprobada la validez de los mismos, realizaremos simulaciones sobre los sistemas objeto de estudio, las cuales nos permitirán observar la dinámica de estos. Por último, en el apéndice se añaden los códigos de MATLAB desarrollados en el trabajo. Capítulo 1 Modelos matemáticos clásicos en neurociencia En este capítulo se presentan algunos modelos matemáticos del campo de la neurociencia, que analizaremos y simularemos numéricamente en los siguientes capítulos. Previamente, veremos algunos conceptos imprescindibles para la comprensión de los mismos. 1.1. Conceptos previos Los modelos matemáticos con los que trabajaremos modelan la sinapsis neuronal. Para una mejor compresión de los mismos conviene describir qué es una neurona, cómo funciona la sinapsis neuronal y en qué parte de este proceso nos centraremos. Comenzamos introduciendo el concepto de neurona. Una neurona es una célula nerviosa que consta de cuatro partes: el núcleo, el pericarión o soma, las dendritas y el axón. Estas células son de suma importancia, pues son únicas en el sentido de que solo ellas pueden transmitir señales eléctricas a largas distancias. La sinapsis neuronal es el proceso mediante el cual una neurona recibe, a través de contactos entre sus dendritas, señales eléctricas transmitidas por otras neuronas, ver Figura 1.1. Al recibir dichas señales eléctricas se producen cambios en las concentraciones iónicas de la membrana de la neurona, es decir, ocurre un cambio del potencial de membrana. Estos cambios se llaman potenciales postsinápticos (PSP). Pequeñas señales producen pequeños PSP y grandes señales producen importantes PSP, que pueden amplificarse por los canales sensibles al voltaje de la membrana neuronal y conducir a la generación de lo que se conoce como spike, (oscilación rápida simple, pulso en español), un cambio abrupto y transitorio del voltaje de la membrana que se propaga a otras neuronas a través del axón; cuando esto ocurre decimos que la neurona se dispara. Los spike son los principales medios de comunicación entre neuronas. 1 8CAPÍTULO 2. ANÁLISIS MATEMÁTICO usuales de suma de funciones y producto de un número real por una función, es un espacio vectorial sobre R. En particular, todas las funciones constantes pertenecen aLip(x, Ω). Definición 2.2. Supongamos que Ωes abierto. Se dice que f: Ω →Rnes una función localmente lipschitziana respecto de la variable xen Ω, si para cada punto (t0, x0)∈Ωexiste una bola abierta de centro dicho punto, contenida en Ω, y tal que la restricción de f a dicha bola es globalmente lipschitziana respecto de la variable x. Es decir, f es localmente lipschitziana respecto de la variable x en Ω, si para cada punto (t0, x0)∈Ωexisten ε(t0, x0), y una constante L(t0, x0)>0, tales que: B((t0, x0), ε(t0, x0)) ⊂Ω, y para todo (t, x1),(t, x2)∈B((t0, x0), ε(t0, x0)), |f(t, x1)−f(t, x2)| ≤ L(t0, x0)|x1−x2|. Si Ωes abierto, se denota Liploc(x, Ω) al conjunto de todas las funciones f: Ω → Rnque son localmente lipschitzianas respecto de la variable xen Ω. Se comprueba de manera inmediata que Liploc(x, Ω), con las operaciones usuales de suma de funciones y producto de un número real por una función, es un espacio vectorial sobre R, y Lip(x, Ω) ⊂Liploc(x, Ω) de forma estricta. 2.1.2. Existencia y unicidad de solución Para Ω⊆Rn+1 enunciaremos un resultado de existencia y unicidad de solución global para el problema de Cauchy: (PC)(x0(t) = f(t, x(t)),(t, x(t)) ∈Ω, x(t0) = x0 donde la prima denota la derivada con respecto de la variable t. Antes, debemos enunciar el Teorema de Picard de existencia y unicidad de solución para el (PC) en un intervalo suficientemente pequeño: Teorema 2.3 (de Picard).Sean Ω⊆Rn+1 un conjunto abierto no vacío, y f: Ω→Rntal que: f∈C(Ω;Rn)∩Liploc(x, Ω). Con estas condiciones, para cada (x0, t0)∈Ω, existe un δ > 0, tal que si denotamos Iδ= [t0−δ, t0+δ], existe una y solo una solución del (PC) en Iδ. Bajo las mismas condiciones, se tienen resultados de existencia y unicidad de solución maximal, es decir, definida en un intervalo lo más grande posible. 2.1. RESULTADOS TEÓRICOS 9 Definición 2.4. Sean Ω⊆Rn+1 un conjunto abierto no vacío, y f: Ω →Rntal que f∈C(Ω;Rn)∩Liploc(x, Ω). Para cada (t0, x0)∈Ω, denotemos por (I(t0, x0), ϕ(·;t0, x0)) a la solución maximal del problema (PC). Se define el conjunto Θ = {(t, t0, x0)∈Rn+2;(t0, x0)∈Ω, t ∈I(t0, x0)}, y la función ϕ: (t, t0, x0)∈Θ7→ ϕ(t, t0, x0)∈Rn. A la función ϕasí definida se la denomina la solución (maximal) del (PC) expresada en función de los datos iniciales. Teorema 2.5 (Lema de Gronwall).Se tienen los siguientes resultados: 1. Supongamos dados −∞ < x0< x1<+∞, dos funciones u, k ∈C([x0, x1]), y una constante h∈R, tales que k(x)≥0, y u(x)≤h+Zx x0 k(s)u(s)ds, ∀x∈[x0, x1]. En tal caso, también se satisface u(x)≤heRx x0k(s)ds,∀x∈[x0, x1]. 2. Supongamos dados −∞ < x1< x0<+∞, dos funciones u, k ∈C([x0, x1]), y una constante h∈R, tales que k(x)≥0, y u(x)≤h+Zx0 xk(s)u(s)ds, ∀x∈[x1, x0]. En tal caso, también se satisface u(x)≤heRx0 xk(s)ds,∀x∈[x1, x0]. Con una aplicación sencilla del Lema de Gronwall, obtenemos el siguiente resultado de unicidad global para (PC): Teorema 2.6 (de unicidad global).Sean Ω⊆Rn+1 un conjunto abierto no vacío, yf: Ω →Rntal que f∈C(Ω;Rn)∩Liploc(x, Ω), y consideremos fijado un punto (t0, x0)∈Ω. En estas condiciones, si (I1, ϕ1)e(I2, ϕ2)son dos soluciones del (PC), entonces, ϕ1(t) = ϕ2(t),∀t∈I1∩I2. Teorema 2.7 (de existencia y unicidad de solución global).Sean Ω⊂Rn+1 un abierto no vacío y f: Ω →Rntal que f∈C(Ω;Rn)∩Liploc(x, Ω). En estas condiciones, para cada (t0, x0)∈Ωdado, existe una y sólo una solución maximal o global del problema de Cauchy (PC) que denotaremos (I(t0, x0), ϕ(·;t0, x0)). Además, el intervalo I(t0, x0)de definición de la solución global es abierto. Los modelos que vamos a estudiar son sistemas autónomos, por ello nos centramos en sistemas de este tipo a partir de ahora. Comenzamos introduciendo el concepto de órbita en la siguiente sección. 10 CAPÍTULO 2. ANÁLISIS MATEMÁTICO 2.1.3. Órbitas Sean Ω⊆Rnun abierto no vacío, y f: Ω →Rntal que f∈C0(Ω;Rn)∩ Liploc(x, Ω). Consideremos el sistema autónomo: x0=f(x) =       f1(x) f2(x) . . . fn(x)       .(2.1) Analizaremos las propiedades cualitativas de las soluciones del sistema autónomo (2.1). Comenzaremos definiendo el concepto de órbita que pasa por un punto x0∈Ω: Definición 2.8. Sea x0∈Ωy denotemos I(x0) = I(0, x0). Se llama órbita del sistema autónomo (2.1) asociada (o, que pasa por) x0al conjunto γ(x0)dado por γ(x0) = {ϕ(t;0, x0) : t∈I(x0)} ⊂ Ω. De la definición se deduce que x0∈γ(x0), y que en general γ(x0)es una curva en Rn. Definición 2.9. Sea x0∈Ω. Se dice que x0es un punto crítico o punto de equilibrio para el sistema (2.1) si se tiene f(x0)=0. Definición 2.10. El conjunto Γ = {x∈Ω : fi(x) = 0}se denomina nulclina de la variable xi. Observación 2.1. Un punto es crítico (de equilibrio) si y sólo si pertenece a la instersección de todas las nulclinas. Definición 2.11. Sea x0∈Ω, se dice que la órbita γ(x0)de (2.1) es cíclica, periódica o cerrada si I(x0) = Ry existe T > 0tal que ϕ(t+T;0, x0) = ϕ(t;0, x0),∀t∈R. En este caso se dice que T es un periodo de la órbita. Además, si una órbita periódica se reduce a un punto, se llamará degenerada. En un sistema dinámico puede haber dos o más puntos de equilibrio y, por tanto, también es posible que existan órbitas conectando dichos equilibrios. Definición 2.12. A las órbitas que conectan puntos de equilibrio se las llama separatrices, y pueden darse dos casos: Si una órbita conecta dos puntos de equilibrio diferentes, se la denomina órbita heteroclina. Si una órbita conecta un punto de equilibrio consigo mismo, se la denomina órbita homoclina. Los equilibrios y órbitas periódicas son elementos organizadores de la dinámica del sistema. Así, una vez identificados, es importante conocer su estabilidad. Intuitivamente, el concepto de estabilidad nos dice si el resto de órbitas del sistema se acercan o alejan de ellos, cuando el tiempo evoluciona. Veamos la definición de este concepto en la siguiente sección. 2.1. RESULTADOS TEÓRICOS 11 2.1.4. Estabilidad Definición 2.13. Un punto de equilibrio del sistema (2.1) se dice estable si ∀ε > 0, existe δ > 0tal que, para cualquier x con kx−x0k< δ, la solución ϕ(t, x)del sistema (1.1) que pasa por x en t= 0 safistace la desigualdad kϕ(t, x)−x0k< ε, ∀t≥0. El punto de equilibrio x0se dirá inestable si no es estable. Definición 2.14. Un punto de equilibrio x0se dice que es asintóticamente estable si es estable y, además, existe r > 0tal que kϕ(t, x)−x0k → 0cuando t→+∞, para todo x satisfaciendo kx−x0k< r. A continuación, nos centraremos en la estabilidad de los puntos de equilibrio de los sistemas lineales. Consideremos el sistema lineal x0=Ax, (2.2) donde A=      a11 a12 ··· a1n a21 a22 ··· a2n . . .. . ..... . . an1an2··· ann       con aij ∈R. Teorema 2.15. Si todos los autovalores de la matriz de coeficientes Adel sistema lineal x0=Ax tienen parte real negativa, entonces el punto de equilibrio x0= 0 es asintóticamente estable. Más aún, existen constantes positivas K y αtal que keAtx0k ≤ Ke−αtkx0k,∀t≤0, x0∈Rn. Si algún autovalor de la matriz Atiene parte real positiva, entonces el punto de equilibrio x0= 0 es inestable. El tipo de estabilidad de un punto de equilibrio es una propiedad local. Por consiguiente, es de esperar que bajo ciertas condiciones el tipo de estabilidad de x0 se pueda aproximar por el campo vectorial fcon sus derivadas, el cual es un campo vectorial lineal. Para ello, supongamos f= (f1, f2, ..., fn)∈C1(Ω) y sea Df(x) =        ∂f1 ∂x1(x)∂f1 ∂x2(x)··· ∂f1 ∂xn(x) ∂f2 ∂x1(x)∂f2 ∂x2(x)··· ∂f2 ∂xn(x) . . .. . ..... . . ∂fn ∂x1(x)∂fn ∂x2(x)··· ∂fn ∂xn(x)        la matriz Jacobiana de fen el punto x. Definición 2.16. Si x0es un punto de equilibrio de x0=f(x), entonces la ecuación diferencial lineal x0=Df(x0)x, se denomina ecuación lineal variacional o linealización del campo vectorial f en el punto de equilibrio x0. 12 CAPÍTULO 2. ANÁLISIS MATEMÁTICO Teorema 2.17. Sea f∈C1(Ω). Si todos los autovalores de la matriz Jacobiana Df(x0)tiene parte real negativa, entonces el punto de equilibrio x0de la ecuación diferencial (2.1) es asintóticamente estable. Teorema 2.18. Sea f∈C1(Ω). Si al menos uno de los autovalores de la matriz Jacobiana Df(x0)tiene parte real positiva, entonces el punto de equilibrio x0de la ecuación diferencial (2.1) es inestable. Definición 2.19. Un punto de equilibrio x0del sistema(2.1) se dice hiperbólico si todos los autovalores de la matriz Jacobiana Df(x0)tienen parte real no nula. El siguiente Teorema es un resultado central de la teoría de Sistemas Dinámicos, pues nos permite conocer el comportamiento del sistema en un entorno de los puntos de equilibrio a partir de su linealización, siempre que el equilibrio sea hiperbólico. Teorema 2.20 (Grobman- Hartman).Si x0es un punto de equilibrio hiperbólico de x0=f(x), entonces existe un entorno de x0en el cual f es topológicamente equivalente al campo vectorial x0=Df(x0)x. Como el sistema de FN es plano, y en los sitemas tipo Bursting es imprescindible conocer la dinámica del subsistema rápido, que también es plano, estudiaremos a continuación el caso de los sistemas lineales planos. Consideremos el sistema lineal x0=f(x), en el que f:R2→R2. En concreto, consideremos los sistemas de la forma: (x0 1=a11x1+a12x2, x0 2=a21x1+a22x2, o equivalententemente x0=Ax, donde A= a11 a12 a21 a22 ! con aij ∈R. Daremos un resultado mediante el cual se pueden clasificar topológicamente los puntos de equilibrio de este sistema, pero antes debemos recordar algunos conceptos: Definición 2.21. Un número (real o complejo) λse denomina autovalor de la matriz Asi existe un vector no nulo (real o complejo) vtal que Av =λv. El vector vse denomina autovector de Aasociado al autovalor λ. En esta definición es necesario que el sistema lineal (A−λI) = 0 sea homogéneo con solución no nula. Esto es equivalente, por el teorema de Rouché Frobenius, a que det(A−λI) = 0. Desarrollando esta expresión nos queda el polinomio cuadrático en λ λ2−τλ + ∆ = 0,(2.3) 2.1. RESULTADOS TEÓRICOS 13 donde τ=traza(A) y ∆ =det(A). Por tanto, los autovalores de Ason las raíces de este polinomio (2.3), que se denomina polinomio característico de A. Los autovalores son λ1=τ+√τ2−4∆ 2, λ2=τ−√τ2−4∆ 2. La naturaleza y la estabilidad de los puntos de equilibrio de un sistema autónomo lineal se pueden describir a partir de la traza (τ) y el determinante (∆) de la matriz A: Lema 2.1. Dados ∆yτde la matriz Adel sistema (2.2), tenemos que: Si ∆<0, los autovalores son reales y tienen signos opuestos; los puntos de equilibrio se denominan puntos de silla. Si ∆>0y: •Si τ2−4∆ >0, los autovalores son reales y tienen el mismo signo; en este caso los puntos de equilibrio se denominan nodos y son: ◦Estables si τ < 0. ◦Inestables si τ > 0. •Si τ2−4∆ <0, los autovalores son complejos conjugados; los puntos de equilibrio se denominan focos o espirales y son: ◦Estables si τ < 0. ◦Inestables si τ > 0. ◦Centros si τ= 0. •Si τ2−4∆ <0y los autovalores son: ◦Iguales y negativos, el punto de equilibrio es un nodo impropio asintóticamente estable. ◦Iguales y positivos, el punto de equilibrio es un nodo impropio inestable. Si ∆ = 0, al menos un autovalor es cero. Entonces el punto de equilibrio no es un punto de equilibrio aislado. Hay una recta de puntos fijos, o un plano de puntos fijos, si A= 0. 2.1.5. Bifurcaciones En referencia a los sistemas dinámicos, una bifurcación se da cuando una pequeña variación en los valores de los parámetros causa un cambio brusco cualitativo o topológico en su comportamiento. 14 CAPÍTULO 2. ANÁLISIS MATEMÁTICO Para el posterior análisis de los sistemas que hemos planteado, es imprescindible el estudio de la bifurcación de Andronov-Hopf, la bifurcación silla-nodo, la bifurcación homoclina y la bifurcación silla-nodo de órbitas periódicas. Para empezar, nos centramos en la bifurcación de Andronov-Hopf, que se da cuando un equilibrio tipo foco cambia de estabilidad, dando lugar al nacimiento de una órbita periódica. Definición 2.22. Una bifurcación de un equilibrio del sistema x0=F(λ, x), para el cual sus autovalores son puramente imaginarios, es decir, λ1,2=±iβ con β > 0, se denomina bifurcación Poincaré-Andronov-Hopf, o bifurcación del nacimiento de un ciclo límite. El siguiente teorema establece algunas condiciones de existencia de órbitas periódicas de la bifurcación de Poincaré-Andronov-Hopf. Teorema 2.23. (Poincaré-Andronov-Hopf) Sea x0=A(λ)x+F(λ, x)de clase Ckcon k≥3un campo vectorial plano dependiendo del parámetro λtal que F(λ, 0) = 0 yDxF(λ, 0) = 0 para todo |λ|suficientemente pequeño. Supongamos que la parte lineal A(λ)tiene, en el origen, los autovalores α(λ)±iβ(λ)con α(0) = 0 y β(0) 6= 0. Además, supongamos que estos autovalores cruzan el eje imaginario con velocidad distinta de cero, es decir, dα dλ(0) 6= 0. Entonces, en cualquier entorno del origen U⊂R2y dado cualquier λ0>0, existe λcon |λ|< λ0tal que la ecuación diferencial x0=A(λ)x+F(λ, x)tiene una órbita periódica no trivial en U. Usando un cambio lineal de variables, podemos transformar la matriz linealizada del campo vectorial en el origen en la matriz A(λ) = α(λ)β(λ) −β(λ)α(λ)!. A partir de la suposición β(0) 6= 0, podemos cambiar la variable de tiempo de modo que β(λ) = 1 para |λ|pequeño. Además, dado que dα dλ (0) 6= 0 , en conjunción con el Teorema de la Función Inversa, implica que hay una correspondencia uno a uno entre α(λ)yλ. Esto nos permite usar α(λ)como el parámetro en lugar de λ. Como resultado de todas estas transformaciones, podemos suponer que la parte lineal del campo vectorial en el punto de equilibrio es de la forma A(λ) = λ1 −1λ!. Sea Funa función Ck,k≥3, definida de la forma 2.1. RESULTADOS TEÓRICOS 15 F:R×R2→R2 (λ, x)7→ F(λ, x), satisfaciendo F(λ, 0) = 0,DxF(λ, 0) = 0. Bajo estas condiciones, consideramos el sistema de ecuaciones (x0 1=λx1+x2+F1(λ, x1, x2), x0 2=−x1+λx2+F2(λ, x1, x2).(2.4) Para el sistema (2.4), el Teorema de Poincaré-Andronov-Hopf se puede reformular de la siguiente manera, donde la variable adebemos verla como la amplitud aproximada de la solución correspondiente a la órbita periódica: Teorema 2.24. Para el sistema (2.4), existen constantes a0>0,λ0>0,δ0>0, y funciones reales C1,λ?(a)yT?(a)de una variable real a, y una función vectorial T?(a)-periódica x?(t, a)con las siguientes propiedades: para 0≤a<a0, λ?(0) = 0,T?(0) = 2π,kx?(0, a)k=a. La función x?(t, a)es una solucón del sistema (2.4) con el valor del parámetro λ=λ?(a)y sus componentes vienen dadas por x? 1(t, a) = acost+o(|a|), x? 2(t, a) = −asint+o(|a|)cuando a→0. Para |λ|< λ0y|T−2π|< δ0toda solución T-periódica x(t)del sistema (2.6) satisfaciendo kx(0)k=aykx(t)k< a0deben venir dadas por la función x?(t, a), excepto por una posible traslación en fase. Veamos ahora la estabilidad de las órbitas del Teorema 2.18, que se puede deducir de la derivada de la función λ?(a)siempre que esta no sea nula. En concreto, se tiene el siguiente resultado: Teorema 2.25. Sea λ?(a)la función dada en el Teorema 2.18 y sea Γala órbita periódica correspondiente al sistema (2.4). Entonces, para a=alo suficientemente pequeña, la órbita periódica Γaes orbitalmente asintóticamente estable si dλ?(a) da >0, e inestable si dλ?(a) da <0. Definición 2.26. Si el ciclo límite que bifurca es estable, la bifurcación de Hopf se llama supercrítica. En otro caso, se dirá subcrítica. Ahora que hemos visto lo que es una bifurcación de Hopf, pasamos a estudiar lo que es una bifurcación silla-nodo, en principio, de equilibrios. El siguiente teorema se centra en la existencia de la bifurcación silla-nodo. 16 CAPÍTULO 2. ANÁLISIS MATEMÁTICO Teorema 2.27. Sea Funa función Ck,k≥1, dependiendo de un parámetro λ, F:R×R2→R2 (λ, x)7→ F(λ, x) satisfaciendo F(0, x) = f(x),f(0) = 0 yDf(0) = 0. Consideramos el sistema (˙x1=F1(λ, x1, x2) ˙x2=−x2+F1(λ, x1, x2), Si ∂F1 ∂λ (0,0,0) 6= 0 y∂2F1 ∂x12(0,0,0) 6= 0, entonces existe una bifurcación silla-nodo para λ= 0, esto es, cuando λ∂F1 ∂λ ∂2F1 ∂x12<0, existen dos equilibrios hiperbólicos, uno es de tipo silla y otro es un nodo asintóticamente estable, y si λ∂F1 ∂λ ∂2F1 ∂x12>0, entonces no existen puntos de equilibrio. Como acabamos de ver, una bifurcación de silla-nodo de equilibrios se da cuando dos equilibrios de estabilidades cambiadas colapsan y desaparecen. Asímismo, una bifurcación de silla-nodo de órbitas periódicas se da cuando dos ciclos de estabilidades cambiadas chocan y desaparecen. Dicho fenómeno, que es global, se puede estudiar a través de los puntos fijos de la llamada aplicación de Poincaré. Finalmente, la bifurcación Homoclina es también un fenomeno global, que se da cuando una órbita periódica se convierte en homoclina y desaparece. Estos fenómenos más complejos, se visualizarán mediante simulaciones numéricas de los sistemas. Una vez introducidos estos resultados procedemos con el análisis de los sistemas antes citados. 2.2. Análisis del sistema de FitzHugh-Nagumo En esta sección analizaremos el sistema de FitzHugh-Nagumo. Para la elaboración de esta sección se han consultado las referencias [4],[5] y [7]. Recordamos las ecuaciones del sistema:      ˙v=v−v3 3+w+I, ˙w=ε(µ−v−νw), (2.5) donde la variable vrepresenta el potencial de membrana de la neurona, wes la variable iónica de recuperación, µ,I∈R, siendo Ila corriente aplicada, ν > 0y 0< ε 1. Sea v= (v, w)∈R2, como      f(v, w) = v−v3 3+w+I, g(v, w) = ε(µ−v−νw), 2.2. ANÁLISIS DEL SISTEMA DE FITZHUGH-NAGUMO 17 son C∞, también son localmente lipschitzianas, f, g ∈Liploc(v,R2), esto implica que f, g ∈C(R2;R)∩Liploc(v,R2), por lo que el Teorema 2.7 (de existencia y unicidad local de solución) se aplica para todo (v0, w0)∈R2. Sin embargo, no podemos asegurar que la función sea globalmente Lipschitz, por tanto, podrían producirse fenómenos de explosión en tiempo finito, cuyo estudio queda fuera de los objetivos de este trabajo. Una vez sabemos que existe solución local y es única, igualando fyga cero y despejando wtenemos las nulclinas:      w=v3 3−v−I(nulclina v), w=µ−v ν(nulclina w). (2.6) Podemos observar que la nulclina v, asociada a la variable rápida, es cúbica (ver Figura 2.1), tiene un máximo y un mínimo que nos permiten diferenciar tres partes de la curva: a la izquierda del máximo, la parte comprendida entre los puntos extremos y a la derecha del mínimo. Mientras que la segunda es una recta y está asociada a la variable lenta. Figura 2.1: Representación gráfica nulclina v Estas curvas pueden intersecarse en uno, dos o tres puntos, que son los equilibrios del sistema. Dado que la pendiente de la nulclina wdepende del parámetro ν, podemos elegirlo de manera que tengamos uno, dos o tres equilibrios, ver Figuras 2.2 - 2.4. De igual manera, variando el parámetro Ien la nulclina vdesplazamos hacia arriba o hacia abajo la curva. Figura 2.2: Un punto de equilibrio Figura 2.3: Dos puntos de equilibrio Figura 2.4: Tres puntos de equilibrio 24 CAPÍTULO 2. ANÁLISIS MATEMÁTICO Capítulo 3 Métodos numéricos En este capítulo se presentan los métodos numéricos que se usarán para realizar las simulaciones de los sistemas. En la primera sección, se explica de manera breve como pasar del problema continuo al problema discreto y se introducen los llamados métodos de un paso, centrándonos en los métodos de Runge-Kutta. Para el desarrollo de esta primera sección se han seguido las referencias [11] y [12], donde se pueden encontrar las demostraciones de los resultados. 3.1. Resolución numérica de problemas ligados a EDO’s Nota 3.1. En lo que sigue nos centramos en EDO’s de primer orden, ya que de tener una EDO de orden mayor podemos reducirla a un sistema de EDO’s de primer orden. Para resolver una EDO numéricamente necesitamos fijar un problema que tenga una única solución. Para lograr esto, completamos la EDO con una condición. Consideremos el problema de valor inicial (PVI): (PV I)(y0(t) = f(t, y(t)), t ∈I0⊆R(I0abierto), y(t0) = y0.(3.1) El siguiente teorema nos garantiza existencia y unicidad de solución. Teorema 3.1 (Existencia y unicidad de solución).Supongamos f:I0×Rm→ Rmy∃l∈L1(I0)tal que (f(t, y)−f(t, z), y −z)Rm≤l(t)ky−zk2 Rm,∀y, z ∈ Rm,∀t∈I0. Entonces: Existe un único y∈C1(I0,Rm)solución de (PVI). Teorema 3.2 (Dependecia continua respecto a los datos).Supongamos que el (PVI) posee una única solución. Sea g∈C(I0,Rn)y sea y0g∈Rm. Dado el 25 26 CAPÍTULO 3. MÉTODOS NUMÉRICOS problema (y0 g(t) = f(t, yg(t)) + g(t), t ∈I0⊆R(I0abierto), yg(t0) = y0g.(3.2) Entonces: |y(t)−yg(t)|≤|y0−y0g|eL(t)+Zt t0 eL(t)−L(s)|g(s)|ds, ∀t∈I0(3.3) siendo L(t) = Rt t0l(s)ds. Con este resultado tenemos que al añadir una perturbación pequeña (g(t)) a la solución, la diferencia de las soluciones con y sin perturbación también es pequeña. Para discretizar el problema, en primer lugar realizamos una partición de I0en N+ 1 puntos: I0= [t0, T] = {t0< t1< ... < tN=T}.(3.4) Siendo hn=tn+1 −tn,∀n= 0,1, ..., N −1yh=max 0<n<N−1hnes el diámetro de la partición. Para que la partición sea buena es necesario que lim N→∞h= 0. En adelante consideraremos particiones uniformes para simplificar la notación, es decir, h=tn+1 −tn,∀n= 0,1, ..., N −1. Una vez realizada la partición, nuestro objetivo es calcular una aproximación de la solución en dichos puntos. A estas aproximaciones las denotaremos por {y(tk)}N k=0 ∈ RN+1. Entonces, partimos del problema (3.1), en cada subintervalo [tk, tk+1],∀k= 0, ..., N −1. (y0(t) = f(t, y(t)),∀t∈[tk, tk+1],∀k= 0, ..., N −1, y(t0) = y0,(3.5) intengrando y aplicando la Regla de Barrow llegamos al problema:      y(tk+1) = y(tk) + Ztk+1 tk f(t, y(t))dt y(t0) = y0 (3.6) A partir de la ecuación (3.6) hay distintas formas de obtener las aproximaciones, una de ellas consiste en aproximar la integral mediante una fórmula de integración numérica. También es muy usual utilizar la fórmula de Taylor, entre otras. En nuestro caso, nos centramos en los que llamamos métodos de un paso (M1P), en los que para calcular yk+1 solo se necesita yk. Estos vienen dados por: (M1P)(yk+1 =ytk+hφ(t, y, h), y0dado, siendo φ(t, y, h)la función que define el método. 3.1. RESOLUCIÓN NUMÉRICA DE PROBLEMAS LIGADOS A EDO’S 27 Todos los métodos que se pueden escribir de esta manera son (M1P) y pueden ser explícitos, implícitos o semi-implícitos. Las propiedades de estos métodos varían según la función φ(t, y, h)que los defina. A continuación se dan algunos resultados relacionados con las propiedades de los métodos de un paso, tales como el error, la estabilidad o la consistencia. Consistencia Definición 3.3. Sea ϕ(tk)la solución exacta del problema continuo. Se define el error de consistencia asociado a un (M1P) como el vector {εk}N−1 k=0 =εN, donde cada εkviene dado por: εk=ϕ(tk+1)−ϕ(tk) + hφ(tk, ϕ(tk), h). Definición 3.4. Se dice que un (M1P) es consistente si lim h→0kεNk=lim N→∞kεNk= 0. El siguiente teorema nos facilita ver la consistencia en los (M1P). Teorema 3.5. El (M1P) es consistente ⇔φ(t, y, 0) = f(t, y),∀y∈Rm,∀t∈I0. Estabilidad Definición 3.6. Un método de un paso es estable si ∃C > 0tal que kyk−zkk ≤ C(ky0−z0k+ N−1 X k=0 |ηk|), donde {yk},{zk}son las soluciones de (yk+1 =ytk+hφ(t, y, h), y0dado, y(zk+1 =ztk+hφ(t, z, h) + ηk, z0dado, respectivamente. En la práctica, para ver si un método es estable se usa la función φ. Se tiene el siguiente resultado sobre la estabilidad: Teorema 3.7. Si φ(t, y, h)es lipschitziana respecto al segundo argumento ⇒el (M1P) es estable. Convergencia y orden Definición 3.8. Se dice que el (M1P) es convergente si lim h→0ky−yhk= 0, donde y es la solución del problema continuo e yhes la aproximación obtenida mediante el método. Teorema 3.9. Si (M1P) es consistente y estable ⇒es convergente. 28 CAPÍTULO 3. MÉTODOS NUMÉRICOS Definición 3.10. Se define el error de discretizacion, que denotamos EDD, como: EDD =ky−yhk. La convergencia del método nos garantiza que estamos aproximando realmente la solución del problema continuo, lo que no nos proporciona es cómo de rápido o lento se acerca la aproximada a la exacta. Para conocer la velocidad de convergencia estudiamos el orden del método. Definición 3.11. Se dice que el (M1P) es de orden (de convergencia) p > 0si ∃C > 0tal que kεNk ≤ Chp. En general es difícil calcular el orden de convergencia de manera teórica, tenemos el siguiente resultado: Teorema 3.12. Supongamos que f∈Cp([t0, T]×Rm)y que φ, ∂φ ∂h ,∂2φ ∂h2, ..., ∂pφ ∂hpexisten y son continuas. Entonces, el (M1P) es de orden p si y sólo si se verifica: ∂kφ(t, y, 0) ∂hk=1 k+ 1f(k)(t, y) ∀t∈[t0, T],∀y∈Rm,∀k= 0,1, ..., p −1. Notación. f(k)(t, y) = ∂f(k−1)(t,y) ∂t +∂f(k−1)(t,y) ∂y f(t, y),∀k≥1. Nota 3.2. Si el (M1P) es consistente, entonces es de al menos orden uno. 3.2. Métodos de Runge-Kutta Dentro de los (M1P) encontramos los métodos de Runge-Kutta (M.R-K), nos centraremos en este tipo de métodos ya que los problemas que tratamos durante el trabajo tienen escalas temporales distintas, es decir, al ser lento-rápido, son problemas stiff o rígidos, y estos métodos están diseñados para este tipo de problemas. Los métodos de Runge-Kutta vienen dados por:      yk+1 =ytk+h q X i=1 bif(tk,i, yk,i),∀k= 0, ..., N −1, y0dado,para ciertos bi, (3.7) siendo tk,i =tk+cih∈[tk, tk+1]con ci∈[0,1],∀i= 1, ..., q, yk,i =yk+h q X j=1 aijf(tk,j, yk,j)∀i= 1, ..., q, y cumpliéndose la siguiente relación q X j=1 aij =ci∀i= 1, ..., q. (3.8) 3.2. MÉTODOS DE RUNGE-KUTTA 29 Las constantes aij ,bi,cise recogen en una tabla de coeficientes que se conoce como Tablero de Butcher. Usualmente los (M.R-K) vienen dados por su Tablero de Butcher. c1a11 a12 ... a1q c2a21 a22 ... a2q ... ... ... ... ... cqaq1aq2... aqq b1b2... bq Cuadro 3.1: Tablero de Butcher de un método Runge-Kutta. Definiendo los vectores b= [b1, b2, ..., bq]T∈Rq,c= [c1, c2, ..., cq]T∈Rqy la matriz A= [aij]∈Rq×Rq, tenemos: c A bT Cuadro 3.2: Representación matricial del tablero de Butcher. Nota 3.3. Todo método que se pueda escribir como tablero de Butcher es un método de tipo Runge-Kutta. Nota 3.4. Los métodos de Runge-Kutta pueden ser explícitos, implícitos o semiimplícitos. Se tiene la siguiente clasificación según la matriz Adel tablero de Butcher del método: Si Aes estrictamente triangular inferior, el (M.R-K) es explícito. Si Aes triangular inferior, el (M.R-K) es semi-implícito. Si Aes no es triangular inferior, el (M.R-K) es implícito. En los (M.R-K) implícitos debemos resolver un sistema no lineal en yk,i. A continuación, estudiaremos las propiedades de estos métodos. Al ser métodos de un paso se les pueden aplicar todos los resultados vistos anteriormente. Se tienen, además, los siguientes resultados. Nota 3.5. En adelante, kAkFdenota la norma fila de la matriz A, y Ldenota la constante de Lipschitz de la función φ, que en los métodos de Runge-Kutta es: φ(t, y, h) = q X j=1 bjf(t+cjh, T(t, y, h)), donde T(t, y, h)es la única solución de y=τ(y), y τes la aplicación dada por: τ:Rq→Rq (x1, ..., xq)t7→ (yk, ..., yk)t+hA(f(tk,1, x1), ..., f(tk,q, xq))t. 30 CAPÍTULO 3. MÉTODOS NUMÉRICOS Comenzamos por ver cuándo el método está bien definido. En caso de que el método sea explícito siempre está bien definido, para los otros casos el siguiente resultado nos proporciona bajo que condiciones el método está bien definido: Lema 3.1. Si hL kAkF<1, entonces el (M.R-K) está bien definido. Una vez vista cual es la condición para que el método esté bien definido, veamos bajo que condiciones es estable y consistente, y por tanto, convergente: Lema 3.2. Si hL kAkF<1, entonces el (M.R-K) es estable. Nota 3.6. Si hL kAkF<1, entonces la función φde la Nota 3.5 es lipschitziana y por el Teorema 3.7 el (M.R-K) es estable. Nota 3.7. Si el método está bien definido ⇒es estable. Lema 3.3. Si bTe= 1, entonces el (M.R-K) es consistente. Nota 3.8. Si bTe= 1, entonces φ(t, y, 0) = f(t, y),∀y∈Rm,∀t∈I0, y por el Teorema 3.5 el método es consistente. Lema 3.4. Si hL kAkF<1ybTe= 1, entonces el (M.R-K) es convergente. Nota 3.9. Por el Teorema 3.9, si el método es estable y consistente entonces es convergente. Respecto al orden de convergencia de los (M.R-K), tenemos el siguiente resultado: Lema 3.5. El método de Runge-Kutta es de orden: Al menos uno si bTe= 1. Al menos dos si bTe= 1 ybTAe =bTc=1 2. Para que el método sea de orden de convergencia al menos tres hay dos condiciones además de las tres ya dadas, es decir, en total hay cinco condiciones. Conforme aumenta el orden aumenta el número de condiciones. La siguiente tabla recoge las condiciones necesarias y suficientes para que el método sea de orden al menos 4. 3.3. MÉTODO DE GAUSS DE DOS ETAPAS 31 Orden Condiciones Nºde condiciones Nºtotal de condiciones 1 q X i=1 bi= 1 1 1 2 q X i=1 bici=1 22 3 q X i,j=1 biaij =1 2 3 q X i=1 bic2 i=1 32 5 q X i,j=1 biaijcj=1 6 4 q X i=1 bic3 i=1 44 9 q X i,j=1 biciaijcj=1 8 q X i,j=1 biaijc2 j=1 12 q X i,j,k=1 biaijajkck=1 24 Cuadro 3.3: Condiciones para alcanzar orden entre 1 y 4. 3.3. Método de Gauss de dos etapas En esta sección se presenta el método de Gauss de dos etapas. Este es un método de Runge-Kutta implícito de cuarto orden que nos ofrece gran precisión a la hora de realizar simulaciones numéricas. Para elaborar esta sección se han consultado las referencias [13] y [20]. El Método de Gauss de dos etapas viene dado por el tablero de Butcher 3−√3 6 1 4 3−2√3 12 3+√3 6 3+2√3 12 1 4 1 2 1 2 Cuadro 3.4: Tablero de Butcher del Método de Gauss. Para construir este método se imponen las condiciones de orden recogidas en el Cuadro 3.3. Fijándonos primero en las condiciones en las que solo están implicados los vectores byc, podemos hallar los valores b1, b2, c1yc2resolviendo el sistema de ecuaciones no lineales 32 CAPÍTULO 3. MÉTODOS NUMÉRICOS                        b1+b2= 1, b1c1+b2c2=1 2, b1c2 1+b2c2 2=1 3, b1c3 1+b2c3 2=1 4, ⇐⇒                        b1= 1 −b2, b2(c2−c1) = 1 2−c1, b2c2(c2−c1) = 1 3−c1 2, b2c2 2(c2−c1) = 1 4−c1 3. Dividiendo la tercera ecuación entre la segunda, y la cuarta entre la tercera, se obtienen dos expresiones para c2, e igualándolas y despejando c1tenemos c1=3±√3 6. Eligiendo c1=3−√3 6(para que c1< c2), nos queda                            b1+b2= 1, b2(c2−3−√3 6) = √3 6, b2c2(c2−3−√3 6) = 1 + √3 12 , b2c2 2(c2−3−√3 6) = 3+2√3 36 . Sumando a la tercera ecuación la segunda multiplicada por −c2y despejándo, obtenemos c2=3 + √3 6. Sustituyendo c2en la segunda ecuación, obtenemos b2=1 2, y sustituyendo ahora b2en la primera se tiene b1=1 2. Una vez conocidos los valores b1, b2, c1yc2, establecemos el siguiente sistema de ecuaciones lineales en aij                    a11 +a12 =c1, a21 +a22 =c2, b1(a11c1+a12c2) + b2(a21c1+a22c2) = 1 6, b1c1(a11c1+a12c2) + b2c2(a21c1+a22c2) = 1 8, ⇐⇒                      a11 +a12 =c1, a11c1+a12c2=4c2−3 12(c2−c1), a21 +a22 =c2, a21c1+a22c2=4c1−3 12(c2−c1). Resolviendo este sistema se obtiene a11 =a22 =1 4, a12 =3−2√3 12 , a21 =3+2√3 12 . 3.4. MÉTODO TR-BDF2 33 Se comprueba con un sencillo cálculo que los valores hallados satisfacen las dos condiciones de orden que no hemos utilizado. Además, aplicando los lemas anteriores tenemos que el método está bien definido y es estable si h < 3−√3 L. 3.4. Método TR-BDF2 Para elaborar esta sección se han consultado las referencias [18] y [19]. El método TR-BDF2 viene dado por: un+γ−γh 2fn+γ=un+γh 2fn, un+1 −1−γ 2−γhfn+1 =1 γ(2 −γ)un+γ−(1 −γ)2 γ(2 −γ)un. Este método está compuesto por dos etapas, en la primera se aplica el método trapezoidal para avanzar del paso tnal paso tn+γ=tn+γh, y en la segunda, se usa el método BDF2 (Backward Differentiation Formulae) para pasar de tn+γa tn+1 =tn+h, siendo γ < 1. Además, este método puede ser escrito como un método de Runge-Kutta semiimplícito, cuyo tablero de Butcher asociado es 0 0 0 0 γγ 2 γ 20 11 2(2−γ) 1 2(2−γ) 1−γ 2−γ 1 2(2−γ) 1 2(2−γ) 1−γ 2−γ Cuadro 3.5: Tablero de Butcher del Método TR-BDF2. Para ver las propiedades de este método aplicaremos los resultados dados anteriormente. Se tiene que si h < 1 Lel método está bien definido y es estable. Además, bte= 1, luego el método es consistente. Por tanto, el método converge con orden al menos uno siempre que h < 1 L. Se comprueba fácilmente que se cumplen las condiciones: bTAe =bTc=1 2. Luego el método es de al menos orden dos. Veamos si puede ser de orden tres, para esto deben cumplirse una serie de condiciones, siendo una de ellas: 3 X i=1 bic2 i=1 3 40 CAPÍTULO 3. MÉTODOS NUMÉRICOS Para el método TR-BDF2 se tiene: xn+1 =xn+√2 4h(f(xk,1) + f(xk,2)) + 2−√2 2hf(xk,3), xk,1=xk, xk,2=xk+2−√2 2h(f(xk,1) + f(xk,2)), xk,3=xk+√2 4h(f(xk,1) + f(xk,2)) + 2−√2 2hf(xk,3). (3.29) Y finalmente, aplicando el método Rosenbrock ROS2 tenemos: xn+1 =xn+3 2k1+1 2k2, (I−γhJ)k1=hf(xn), (I−γhJ)k2=hf(xn+k1)−2k1, (3.30) -1.5 -1.4 -1.3 -1.2 -1.1 -1 -0.9 -0.8 -0.7 -0.5 0 0.5 1 1.5 Figura 3.3: Determinante para ε= 0,08. donde (I−γhJ) =    1−γh(1 −x2 1n)γh −γh −δγh 1+0,8δγh 0 εγh 0 1 + εγh    y se tiene det(I−γhJ) = (1 −γh(1 −x2 1n))(1 + 0,8δγh)(εγh)+ +γ2h2ε(1 + 0,8δγh) + δγ2h2(1 + εγh)6= 0, ∀x1n∈[−3/2,−7/10]. En la Figura 3.3 se representa el determinante anterior en el intervalo [-3/2, -7/10] para ε= 0,0001 y δ= 0,08. Capítulo 4 Simulaciones En este capítulo compararemos los resultados que nos ofrecen los métodos presentados en el Capítulo 3 aplicados a los sistemas de FitzHugh-Nagumo, FitzHugh- Rinzel y Hindmarsh-Rose. Nos centraremos en cada de uno de ellos por separado y veremos gráficamente las soluciones de los mismos, así como los resultados obtenidos en el Capítulo 2. Para obtener los resultados implementaremos los métodos en el programa MATLAB. Dado que apenas se conocen soluciones exactas de los sistemas que nos ocupan, para realizar la comparación entre los métodos debemos fijar qué solución consideraremos exacta. Aunque el programa que utilizamos para realizar las simulaciones cuenta con una serie de solvers específicos para problemas stiff, usaremos los métodos estudiados en el trabajo, que nos permiten fijar un paso de tiempo concreto. Consideraremos como solución exacta la proporcionada por el método de Gauss construido en la Sección 3.3 que, como hemos probado, es un método de tipo Runge-Kutta de cuarto orden. Previamente, con el fin de validar nuestros métodos realizaremos un test académico sobre ellos. Es decir, compararemos las soluciones que nos ofrecen cada uno de ellos con la solución exacta de un problema stiff, para ello calcularemos el error que cometen cada uno de ellos y comprobaremos gráficamente que la soluciones aproximadas son buenas. Además, comprobaremos que los métodos son del orden esperado. Nota 4.1. Los códigos desarrollados pueden encontrarse en el apéndice de este trabajo. A lo largo del capítulo calcularemos el error cometido por distintos métodos frente al método de Gauss, para ello usaremos las normas l2yl∞que nos ofrece MATLAB, estas son: kvk2= N X k=1 |vk|2!1 2 ykvk∞= m´ax k(|vk|). Cuyos comandos en MATLAB son norm(v,2) ynorm(v, inf), respectivamente. 41 42 CAPÍTULO 4. SIMULACIONES 4.1. Test académico En esta sección vamos a realizar un test académico con el fin de validar nuestros métodos. Dado que los sistemas que nos ocupan son stiff, realizaremos el test sobre un problema stiff del que sí conocemos solución. Sea (PS): (PS)     u0(t) = −2000u(t) + 999,75v(t) + 1000,25, v0(t) = u(t)−v(t), u(0) = 0, v(0) = −2, (4.1) cuya solución exacta es: (u(t)=1−1,499875e−0,5t+ 0,499875e−2000,5t, v(t)=1−2,99975e−0,5t−0,00025e−2000,5t.(4.2) Para que se cumplan las condiciones de estabilidad es necesario implementar el método con paso de tiempo h < 0,00138. Hemos calculado los errores en norma l2que cometen el método de Gauss, el método ROS2 y el método TR-BDF2 respecto a la solución exacta del sistema (PS) al implementar los métodos con distintos pasos de tiempo en el intervalo de tiempo [0,1.e−2]. Dado que estamos trabajando con un problema stiff es importante conocer la talla de la solución exacta. Además, estamos calculando las soluciones en un intervalo de tiempo muy pequeño, por lo que calculamos también los errores relativos, que nos servirán para determinar el orden de los métodos. Para calcular el error relativo en norma l2hacemos ERROR RELATIV O =Error en norma l2 Norma l2de la solución exacta hM. GAUSS M.ROS2 M.TR-BDF2 Error E.Rel Error E.Rel Error E.Rel 1.e−3 3.9123e−3 5.7788e−4 1.0281e−1 1.5186e−2 3.4264e−2 5.0610e−3 1.e−3/2 3.5926e−4 3.8346e−5 5.1676e−2 5.5156e−3 1.1400e−2 1.2168e−3 1.e−3/4 3.1109e−5 2.3746e−6 2.1668e−2 1.6540e−3 3.7816e−3 2.8866e−4 1.e−3/8 2.7239e−6 1.4788e−7 8.2627e−3 4.4858e−4 1.2968e−3 7.0402e−5 Cuadro 4.1: Error y error relativo cometido por los métodos respecto a la solución exacta de (PS) en norma l2. 4.1. TEST ACADÉMICO 43 2 3 4 5 6 7 8 9 10 h10-4 10-6 10-5 10-4 10-3 10-2 10-1 100 Error Figura 4.1: Error en norma l2. 2 3 4 5 6 7 8 9 10 h10-4 10-7 10-6 10-5 10-4 10-3 10-2 10-1 Error Figura 4.2: Error relativo en norma l2. Tanto en el Cuadro 4.1 como en las Figuras 4.1-4.2 se observa como el error cometido por el método de Gauss es menor que el cometido por los otros dos métodos, siendo el método TR-BDF2 más preciso que el método ROS2. Este hecho se observa de manera clara al representar las soluciones. En las Figuras 4.3-4.5 se representa la primera componente de la solución proporcionada por cada uno de los métodos junto a la solución exacta, podemos observar como para h=1.e−3 el método ROS2 tiene un pico que se aleja de la solución para t≈0,001, debido a que el paso de tiempo es grande, ver Figura 4.3. Sin embargo, al disminuir el paso de tiempo es casi imperceptible a simple vista la diferencia entre las soluciones, ver Figuras 4.4-4.5. De igual forma, en las Figuras 4.6-4.8 se representa la segunda componente de la solución proporcionada por cada método junto con la solución exacta, en este caso podemos observar como la solución del método ROS2 está bastante alejada de la solución que proporcionan los otros dos métodos. Por último, comprobamos el orden de convergencia de cada método, haciendo ORDEN =log2 Error relativo con paso h Error relativo con paso h/2! obtenemos los siguientes datos: hM. GAUSS M.ROS2 M.TR-BDF2 E.Relativo Orden E.Relativo Orden E.Relativo Orden 1.e−3 5.7788e−4 - 1.5186e−2 - 5.0610e−3 - 1.e−3/2 3.8346e−5 4.0001 5.5156e−3 1.9450 1.2168e−3 2.0153 1.e−3/4 2.3746e−6 4.0005 1.6540e−3 1.9735 2.8866e−4 2.0068 1.e−3/8 1.4788e−7 3.9915 4.4858e−4 1.9870 7.0402e−5 2.0032 Cuadro 4.2: Estimación del orden de los métodos. Para el método de Gauss obtenemos orden 4, y para los métodos ROS2 y TRBDF2 obtenemos orden 2, por tanto nuestros métodos convergen con el orden esperado. 44 CAPÍTULO 4. SIMULACIONES 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 Sol. exacta Sol. M.GAUSS Sol. M.ROS2 Sol. M.TR-BDF2 Figura 4.3: u(t)para h=1.e-3. 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 -0.5 -0.45 -0.4 -0.35 -0.3 -0.25 -0.2 -0.15 -0.1 -0.05 0 Sol. exacta Sol. M.GAUSS Sol. M.ROS2 Sol. M.TR-BDF2 Figura 4.4: u(t)para h=1.e-3/2. 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 -0.5 -0.45 -0.4 -0.35 -0.3 -0.25 -0.2 -0.15 -0.1 -0.05 0 Sol. exacta Sol. M.GAUSS Sol. M.ROS2 Sol. M.TR-BDF2 Figura 4.5: u(t)para h=1.e-3/8. 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 -2 -1.995 -1.99 -1.985 -1.98 -1.975 -1.97 -1.965 -1.96 -1.955 -1.95 Sol. exacta Sol. M.GAUSS Sol. M.ROS2 Sol. M.TR-BDF2 Figura 4.6: v(t)para h=1.e-3. 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 -2 -1.995 -1.99 -1.985 -1.98 -1.975 -1.97 Sol. exacta Sol. M.GAUSS Sol. M.ROS2 Sol. M.TR-BDF2 Figura 4.7: v(t)para h=1.e-3/2. 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 -2 -1.998 -1.996 -1.994 -1.992 -1.99 -1.988 -1.986 -1.984 -1.982 Sol. exacta Sol. M.GAUSS Sol. M.ROS2 Sol. M.TR-BDF2 Figura 4.8: v(t)para h=1.e-3/8. 4.2. TEST APLICADO 45 4.2. Test aplicado Una vez visto que los métodos funcionan de manera adecuada, procedemos a implementarlos para los sistemas que nos ocupan, considerando como exacta la solución del método de Gauss. Comenzamos por el sistema de FitzHugh-Nagumo. 4.2.1. Sistema de FitzHugh-Nagumo Para comenzar, presentamos la siguiente tabla en la que se recoge el error del método Rosenbrock ROS2 y el método TR-BDF2 respecto al método de Gauss, cuando resolvemos el problema con condición inicial v0= (0,0) en el intervalo de tiempo [0,1]: FN ERROR M.ROS2 ERROR M.TR-BDF2 hl∞l2l∞l2 1.e−1 8.9451e−4 1.1720e−3 3.5169e−4 4.2946e−4 1.e−2 9.3173e−6 3.3610e−5 3.5674e−6 1.1418e−5 1.e−3 9.3636e−8 1.0527e−6 3.5724e−8 3.5422e−7 1.e−4 9.3681e−10 3.3256e−8 3.5731e−10 1.1180e−8 Cuadro 4.3: Error en norma l∞yl2al resolver el sistema de FN con paso de tiempo h. 10-4 10-3 10-2 10-1 h 10-10 10-8 10-6 10-4 10-2 Error Figura 4.9: Error sistema FN en norma l∞yl2. 46 CAPÍTULO 4. SIMULACIONES Se observa que el error cometido en el método TR-BDF2 es ligeramente menor que en el método ROS2 al resolver el sistema con el mismo paso de tiempo, ver Figura 4.9. Dado que los errores son muy pequeños para cada paso de tiempo, en lo que sigue fijamos el paso de tiempo h= 0,01. A continuación, veremos gráficamente la bifurcación de Hopf, que como estudiamos en el Capítulo 2 ocurre para I=−0,7. En las Figuras 4.10-4.12, representamos las nulclinas del sistema, (nulclina ven rojo y nuclina wen verde), junto con algunas órbitas del plano de fase (en negro) para distintos valores del parámetro I, y en las Figuras 4.13-4.15, las soluciones del sistema para dichos valores de I. En la Figura 4.10, para I=−0,6el sistema tiene un equilibrio tipo foco estable, que corresponde al estado de reposo de la neurona, como se puede observar en la Figura 4.13. Disminuyendo el valor del parámetro, para I=−0,7, ocurre la bifurcación de Hopf, ver Figura 4.11, el sistema comienza a oscilar levemente, ver Figura 4.14. Por último, tras producirse la bifurcación nace un ciclo límite estable, ver Figura 4.12, que corresponde al estado de actividad neuronal o spike, ver Figura 4.15. -3 -2 -1 0 1 2 3 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Figura 4.10: I=−0,6. 4.2. TEST APLICADO 47 -3 -2 -1 0 1 2 3 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Figura 4.11: I=−0,7. -3 -2 -1 0 1 2 3 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Figura 4.12: I=−0,8. 48 CAPÍTULO 4. SIMULACIONES 0 50 100 150 200 250 300 350 400 450 500 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 Potencial de membrana Recuperacion ionica Figura 4.13: I=−0,6. 0 50 100 150 200 250 300 350 400 450 500 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 Potencial de membrana Recuperacion ionica Figura 4.14: I=−0,7. 0 50 100 150 200 250 300 350 400 450 500 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 Potencial de membrana Recuperacion ionica Figura 4.15: I=−0,8. 4.2. TEST APLICADO 49 4.2.2. Sistema de Hindmarsh-Rose En la siguiente tabla se recoge el error del método Rosenbrock ROS2 y del método TR-BDF2 respecto al método de Gauss, cuando resolvemos el problema con condición inicial x0= (0,0,0) en el intervalo de tiempo [0,1]: HR ERROR M.ROS2 ERROR M.TR-BDF2 hl∞l2l∞l2 1.e−1 1.6875e−1 2.5248e−1 7.6499e−2 1.1972e−1 1.e−2 3.0020e−3 1.3804e−2 7.7610e−4 3.7238e−3 1.e−3 3.2096e−5 4.6627e−4 7.7687e−6 1.1741e−4 1.e−4 3.2315e−7 1.4844e−5 7.7692e−8 3.7117e−6 Cuadro 4.4: Error en norma l∞yl2al resolver el sistema de HR con paso de tiempo h. 10-4 10-3 10-2 10-1 h 10-8 10-6 10-4 10-2 100 Error Figura 4.16: Error Sistema HR en norma l∞yl2. Como ocurría para el sistema de FN el error cometido por el método TR-BDF2 es menor que el cometido por el método ROS2, ver Figura 4.16. Sin embargo, en este caso los errores, aun siendo muy pequeños para un paso de tiempo adecuado, son mayores que los cometidos para FN. En las Figuras 4.17-4.19, se representa la solución dada por cada uno de los métodos, al ser el error entre los métodos tan pequeño no es posible distinguir a simple vista la solución proporcionada por cada uno de ellos. 56 CAPÍTULO 4. SIMULACIONES 4.2.3. Sistema de FitzHugh-Rinzel Presentamos la siguiente tabla en la que se recoge el error del método Rosenbrock ROS2 y del método TR-BDF2 respecto al método de Gauss, cuando resolvemos el problema con condición inicial x0= (0,0,0) en el intervalo de tiempo [0,1]: FR ERROR M.ROS2 ERROR M.TR-BDF2 hl∞l2l∞l2 1.e−1 1.1456e−4 2.5342e−4 1.4118e−4 2.5967e−4 1.e−2 1.0757e−6 7.3863e−6 1.4125e−6 7.7808e−6 1.e−3 1.0690e−8 2.3169e−7 1.4126e−8 2.4467e−7 1.e−4 1.0685e−10 7.3216e−9 1.4124e−10 7.7319e−9 Cuadro 4.5: Error en norma l∞yl2al resolver el sistema de FR con paso de tiempo h. 10-4 10-3 10-2 10-1 h 10-12 10-10 10-8 10-6 10-4 10-2 Error Figura 4.26: Error Sistema FR en norma l∞yl2. A diferencia de los errores obtenidos para los sistema de FN y HR, en este caso observamos que ambos métodos cometen un error similar en norma infinito, siendo además algo mayor el cometido por el método TR-BDF2 que el cometido por el método ROS2 para la norma l2, ver Figura 4.26. En las Figuras 4.27-4.29 se representa la solución obtenida con cada uno de los métodos, como ocurría en el sistema de HR, al ser los errores tan pequeños a simple 4.2. TEST APLICADO 57 vista no se diferencian las soluciones proporcionadas por cada uno de los métodos. Las siguientes gráficas se han elaborado con condición inicial x0= (−0,9,0,0,009) y paso de tiempo h= 0,01 en el intervalo [0,500]. Se representa individualmente cada una de las componentes de la solución debido a que sus rangos son distintos y al representarlas unidas se pierde información sobre las mismas. 0 50 100 150 200 250 300 350 400 450 500 -1.5 -1.4 -1.3 -1.2 -1.1 -1 -0.9 -0.8 -0.7 -0.6 x1 ROS2 x1 TRBDF2 x1 GAUSS Figura 4.27: Soluciones para el sistema de FR. 58 CAPÍTULO 4. SIMULACIONES 0 50 100 150 200 250 300 350 400 450 500 -0.4 -0.35 -0.3 -0.25 -0.2 -0.15 -0.1 -0.05 0 x2 ROS2 x2 TRBDF2 x2 GAUSS Figura 4.28: Soluciones para el sistema de FR. 0 50 100 150 200 250 300 350 400 450 500 0.009 0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 ROS2 TRBDF2 GAUSS Figura 4.29: Soluciones para el sistema de FR. 4.2. TEST APLICADO 59 A continuación veremos la dinámica bursting del sistema, que recordamos es del tipo elíptico, ver figura 4.30. En la Figura 4.31 se tiene una representación tridimesional de las órbitas tipo bursting del sistema de FR, podemos observar como para µ≈0,018 sucede la bifurcación de Hopf subcrítica que da lugar al comienzo de las oscilaciones y para µ≈0,01 sucede la bifurcación silla-nodo de óbitas periódicas, provocando el cese de las oscilaciones. Al aumentar de nuevo µse completa el ciclo de histéresis. La Figura 4.32 nos ofrece otro punto de vista de las órbitas tipo bursting, además la línea roja representa el único equilibrio del sistema como función de µ, que ocurre en x1≈ −0,9674, cuando x1<−0,9674 el equilibrio es un foco estable y cuando x1>−0,9674 el equilibrio es un foco inestable. 0 500 1000 1500 2000 2500 3000 3500 4000 -3 -2 -1 0 1 2 3 Figura 4.30: Bursting elíptico en el sistema de FR. 60 CAPÍTULO 4. SIMULACIONES Figura 4.31: Representación tridimensional de las órbitas tipo bursting del sistema de FR. Para finalizar, en las Figuras 4.33-4.35 se representan las nulclinas del subsistema rápido (nulcina x1en rojo y nulclina x2en verde), junto con algunas órbitas del plano de fases. Estas figuras nos permiten ver la dinámica bursting de forma más concreta. En principio, el sistema se encuentra en estado de reposo o resting, el equilibrio se comporta como un foco atractivo, ver Figuras 4.33. Conforme µaumenta el sistema sigue al equilibrio estable, hasta que ocurre la bifurcación de Hopf subcrítica, ver Figura 4.34, desaparece el ciclo límite inestable y el equilibrio se vuelve inestable, lo cual provoca que el sistema siga al ciclo límite estable exterior y, por tanto, las oscilaciones. A continuación, µdecrece y el sistema sigue al ciclo límite estable, hasta que para µ≈0,01 sucede la bifurcación silla-nodo de órbitas periódicas, ver Figura 4.35, provocando el cese de las oscilaciones. El sistema recupera el estado de reposo y sigue al equilibrio estable. Finalmente, µaumenta y se completa el ciclo de histéresis. 4.2. TEST APLICADO 61 0 0.005 0.01 0.015 0.02 0.025 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 Figura 4.32: Esquema global de bursting en el sistema de FR -4 -3 -2 -1 0 1 2 -3 -2 -1 0 1 2 3 Figura 4.33: µ= 0 62 CAPÍTULO 4. SIMULACIONES -4 -3 -2 -1 0 1 2 -3 -2 -1 0 1 2 3 Figura 4.34: µ= 0,01878 -4 -3 -2 -1 0 1 2 -3 -2 -1 0 1 2 3 Figura 4.35: µ= 0,01 4.2. TEST APLICADO 63 A lo largo del capítulo observamos que los tres métodos nos proporcionan buenas aproximaciones de las soluciones exactas, sirviendo cualquiera de ellos para hacer las simulaciones y permitiéndonos ver los resultados de manera clara. Sin embargo, si hemos de elegir uno de ellos debemos tener en cuenta ciertos aspectos a nivel de trabajo informático. Por un lado, tenemos el método de Gauss de orden 4, el cual es totalmente implícito, lo que implica mayor trabajo en cada iteración al tener que recurrir a un método auxiliar que nos permita resolver el problema, en nuestro caso, el método de aproximaciones sucesivas. Por otro lado, tenemos el método TR-BDF2 de orden 2 que, pese a ser semi-implícito, a nivel informático funciona de la misma manera que el método de Gauss, por tanto, en este caso a la desventaja que supone el coste informático debemos añadirle que al ser de un orden inferior a el método de Gauss necesitariamos más iteraciones para conseguir el mismo resultado. Por último, tenemos el método ROS2 de orden 2, el cual para aproximar las soluciones de los sistemas resuelve un sistema de ecuaciones lineales, por tanto, encontramos en este método una ventaja a nivel de trabajo informático. En definitiva, la elección del método depende de los sistemas a resolver y las capacidades tecnológicas de las que dispongamos. 64 CAPÍTULO 4. SIMULACIONES Códigos MATLAB 4.3. Método de Gauss function [tt,v]=GAUSS(h,T,f,v0) %Método de Gauss de dos etapas hmed=h/2; tt=T(1):h:T(2); Nt= length(tt); c1= 1/4; c2=(3-2*sqrt(3))./12; c3=(3+2*sqrt(3))./12; c4=(3-sqrt(3))./6; c5=(3+sqrt(3))./6; v=[]; v(1,:)=v0; for k=1:(Nt-1) tk1=tt(k)+c4*h; tk2=tt(k)+c5*h; g1=@(x) v(k,:)+c1*h*f(tk1, x(:,1))+c2*h*f(tk2,x(:,2)); g2=@(x) v(k,:)+c3*h*f(tk1, x(:,1))+c1*h*f(tk2,x(:,2)); g=@(x) [g1(x); g2(x)]’; x0=[v(k,:);v(k,:)]’; vk= mas(g,x0); $Método de aproximaciones sucesivas v(k+1,:)= v(k,:)+hmed*(f(tk1,vk(:,1))+f(tk2,vk(:,2))); end end Método de aproximaciones sucesivas (mas) function x= mas(g,x0) %Método de aproximaciones sucesivas epsi= 1.e-8; Nmax= 1000; x = x0; for k = 1:Nmax 65