Aproximaciones numéricas descentradas para sistemas de EDP con quimiotaxis

Universidad de Sevilla Facultad de Matemáticas Aproximaciones numéricas descentradas para sistemas de EDP con quimiotaxis Autor: Marcos Romero Madroñal Tutor: Francisco Guillén González Trabajo Fin de Grado Doble Grado en Física y Matemáticas Departamento de Ecuaciones Diferenciales y Análisis Numérico Sevilla, 2022 Abstract In this work, the phenomenon of chemotaxis is modeled, mainly by focusing on the Keller-Segel model, as well as other variants of it. Furthermore, we carry out approximations of finite differences in time and different schemes in space, such as: finite differences (FD), finite elements (FE) and finite volumes (FV). Subsequently, the schemes are particularized to a one-dimensional spatial domain and numerical simulations are performed. Analyzing the results, spurious oscillations are observed in the centered approximations and upwind schemes are proposed in each of the approximations in space: FD, FE and FV. It is proved that the three upwind schemes are equivalent in 1D and they correct the numerical oscillations of the centered schemes. In addition, the upwind scheme is proven to be solvable, and preserves the properties of conservation of total cell number, and non-negativity of cells and chemical substance, these properties are corroborated with numerical simulations. Finally, the discretization techniques developed in this work are applied to other models with chemotaxis. Resumen En este Trabajo de Fin de Grado se modeliza el fenómeno de la quimiotaxis principalmente centrándonos en el modelo de Keller-Segel, además de otras variantes del mismo. A su vez, realizamos aproximaciones numéricas del tipo diferencias finitas en tiempo y distintos esquemas en espacio, tales como: diferencias finitas (DF), elementos finitos (EF) y volúmenes finitos (VF). Posteriormente, se particularizan los esquemas a un dominio espacial unidimensional y se realizan simulaciones numéricas. Analizando los resultados, se observan oscilaciones espúreas en las aproximaciones centradas, y se proponen esquemas descentrados en cada una de las aproximaciones en espacio: DF, EF y VF. Se prueba que los tres esquemas descentrados son equivalentes en 1D y corrigen las oscilaciones numéricas de los esquemas centrados. Además, se prueba que el esquema descentrado tiene solución, y preserva las propiedades de conservación de la cantidad total de células y de no negatividad de células y química, estas propiedades se corroboran con las simulaciones numéricas. Finalmente se aplican las técnicas de discretización desarrolladas en este trabajo a otros modelos con quimiotaxis. Índice 1. El modelo de Keller-Segel 8 1.1. Desarrollo matemático del modelo . . . . . . . . . . . . . . . . . . . . 8 1.2. Propiedades................................ 12 2. La aproximación en tiempo 17 2.1. Presentación del esquema . . . . . . . . . . . . . . . . . . . . . . . . 17 2.2. Conservación de las propiedades . . . . . . . . . . . . . . . . . . . . . 18 3. Aproximación en espacio por elementos finitos y diferencias finitas 22 3.1. Diferencias finitas centradas . . . . . . . . . . . . . . . . . . . . . . . 22 3.2. Método de elementos finitos . . . . . . . . . . . . . . . . . . . . . . . 26 3.3. Simulación y análisis de los resultados. Oscilaciones espúreas . . . . . 33 4. Diferencias finitas descentradas 40 4.1. La técnica del descentrado . . . . . . . . . . . . . . . . . . . . . . . . 40 4.2. Análisis numérico del esquema descentrado (I) . . . . . . . . . . . . . 42 4.3. Simulaciones y análisis de resultados . . . . . . . . . . . . . . . . . . 44 5. Volúmenes finitos descentrados 50 6. Elementos finitos descentrados 54 6.1. Análisis numérico del esquema descentrado (II) . . . . . . . . . . . . 57 4ÍNDICE 7. Otros modelos con quimiotaxis 61 7.1. Keller-Segel con reacción logística . . . . . . . . . . . . . . . . . . . . 61 7.1.1. Simulación y análisis de los resultados . . . . . . . . . . . . . 62 7.2. Modelo de patrones de bacterias quimiotácticas en medio líquido. . . 65 7.2.1. Simulación y análisis de los resultados . . . . . . . . . . . . . 67 7.3. Modelo de patrones de bacterias quimiotácticas en medio semi-sólido. 69 7.3.1. Simulación y análisis de los resultados . . . . . . . . . . . . . 70 A. Mass lumping 75 B. Códigos de simulación 77 Bibliografía 92 Introducción Motivación Hace ya más de 300 años, en el año 1676, el neerlandés Anton van Leeuwenhoek observaba detenidamente a las pequeñas entidades que acababa de descubrir en un preparado de agua con pimienta, mientras trataba de encontrar el origen del picor que esta producía en la lengua humana. Fue entonces cuando, en completo asombro, observó que estas eran capaces de moverse por sí mismas y al intuir estaban dotadas de vida, las bautizó como “animáculos” (se cree que eran bacterias). Poco después, observó cómo su movimiento no era totalmente estocástico, sino que parecía estar dirigido por ciertas sustancias, su hipótesis quedó reforzada al verlos desplazarse hacia conchas abiertas de mejillones y devorarlos. Esta es la primera observación registrada de microorganismos moviéndose de acuerdo a la concentración de sustancias químicas en su medio ambiente, el fenómeno que contempla esta dinámica, cuando se habla de seres vivos, recibe el nombre de quimiotaxis. [3] La manifestación de la quimiotaxis no es para nada anecdótica en la naturaleza, existen multitud de ejemplos que atienden a este fenómeno, como pudieran ser: la migración de células del sistema inmune hacia tejidos dañados o inflamados, atraídas por proteínas como las quimiocinas; el movimiento de bacterias hacia sustancias nutritivas como la glucosa, o en contra de sustancias venenosas como el fenol; la atracción que generan los óvulos de los helechos, los cuales poseen ácido málico que motiva el movimiento de los espermatozoides; el proceso de formación de vasos sanguíneos (angiogénesis) incentivado por químicos producidos por células tumorales, etc. Debido a la gran cantidad de procesos mediados por el fenómeno de la quimiotaxis, resulta de gran interés matemático tratar de modelizar este suceso. El primer modelo matemático de la quimiotaxis data de los años 1950s presentado por el matemático Clifford S. Patlak. Posteriormente en los años 70s, este modelo fue mejorado por la física Evelyn Fox Keller y el matemático Lee Aaron Segel tratando de modelar el comportamiento de un tipo de ameba capaz de agregarse con otras de su misma especie, conformando un único ser conocido como moho del fango (Dictyostelium discoideum), esta atracción se lleva a cabo mediante una sustancia química atrayente, el adenosín monofostafo cíclico (AMFc). 6ÍNDICE Descripción del trabajo En el Capítulo 1, nos centraremos en modelizar el fenómeno de quimiotaxis deduciendo las ecuaciones que constituyen el modelo de Keller-Segel e interpretando físicamente los términos que las componen. Asimismo, se deducen varias propiedades del modelo, como son: la conservación de la cantidad total de células, la no negatividad de las densidades de células y química, y una ley de energía donde en particular el funcional de energía es decreciente. En el Capítulo 2, proponemos un esquema de discretización en tiempo mediante DF. Posteriormente, se estudia si nuestro esquema en tiempo conserva las propiedades del problema continuo, manteniéndose la conservación, la no negatividad y el decrecimiento de la energía. Además el esquema introduce un término de disipación extra, que podemos entender como disipación numérica. En el Capítulo 3, se lleva a cabo una aproximación centrada en espacio en un dominio unidimensional mediante DF y EF, y se demuestra que estos dos enfoques son equivalentes, se prueba además para el esquema centrado la propiedad de conservación. Más adelante, se realizan simulaciones numéricas donde se observan oscilaciones espúreas y se corrobora el hecho de que estas se producen por la discretización espacial centrada. Asimismo, se comprueba mediante las simulaciones si el nuevo esquema completamente discretizado mantiene las propiedades, teniéndose únicamente la propiedad de conservación. En el Capítulo 4, con vistas a eliminar las oscilaciones espúreas, proponemos un nuevo método de descentrado para la aproximación en espacio en el caso de las DF. A su vez, aprovechando esta formulación, se demuestra que el esquema descentrado impide la negatividad de las variables. Posteriormente, se realizan aproximaciones numéricas para corroborar estos nuevos resultados, observándose asimismo que el esquema parece conservar la cantidad total de células. Finalmente, se demuestra que el esquema descentrado introduce un tipo de disipación numérica que produce un amortiguamiento en la aproximación de máximos crecientes y mínimos decrecientes. En el Capítulo 5, realizamos una aproximación descentrada en espacio mediante un método de VF. Se concluye que, en el caso de dominios 1D, el esquema resultante de la aproximación centrada coincide con el esquema centrado de DF y EF, mientras que la aproximación descentrada da lugar al mismo esquema descentrado del caso de DF. En el Capítulo 6, se realiza una aproximación descentrada en espacio mediante un método de EF descentrado. Se concluye que dicho esquema descentrado es equivalente al caso de DF y VF descentrado en dominios 1D. Asimismo, se aprovecha esta nueva formulación para demostrar la propiedad de conservación y se prueba la existencia de solución del esquema descentrado completamente discretizado. En el Capítulo 7, se introducen nuevos modelos continuos con quimiotaxis, que se discretizan de forma descentrada teniendo en cuenta los resultados de los demás 7ÍNDICE capítulos. Además, se realizan simulaciones numéricas y se analiza el comportamiento de las aproximaciones de sus soluciones con respecto al caso de Keller-Segel, tratando de explicar las dinámicas observadas. Para concluir, cabe decir que el análisis y simulación numérica en modelos de quimiotaxis tiene una extensa literatura con muy diversos enfoques; sin embargo, en lo que se refiere al objetivo principal del trabajo, que es la comparación de los esquemas numéricos más conocidos en el caso unidimensional, no hemos encontrado referencias en la literatura, con lo que gran parte del desarrollo del trabajo es original. Capítulo 1 El modelo de Keller-Segel 1.1. Desarrollo matemático del modelo En primer lugar, será necesario especificar nuestro marco de trabajo. El medio material sobre el que se trabaja se modela de forma aproximada como un dominio abierto y conexo Ω⊂Rd, donde cada punto x∈Ω, al tomar un sistema de referencia, se corresponde con un punto espacial fijo. Asimismo, los valores de la variable temporal se consideran reales desde el instante presente t= 0 hasta un momento futuro t=Tcon T > 0. Por tanto t∈(0, T)⊂R. Sea u(x, t)la densidad de células y v(x, t)la concentración de sustancia química en un punto del espacio xy en el instante de tiempo t. Sea O ⊂ Ωun volumen de control. El balance de la cantidad total de células implica la siguiente ecuación d dt ZO u(x, t)dx =P(u)(O, t)−Z∂O (J(u)(x, t)·n(x))dS(x),∀O ⊂ Ω,(1.1) donde n(x)es el vector normal exterior a la frontera de Oen el punto xy (·) denota al producto escalar correspondiente. Primeramente, ROut(x, t)dx := d dt ROu(x, t)dx denota la variación temporal de la cantidad total de células contenidas en el volumen O. Estudiemos ahora cada uno de los dos sumandos de la derecha de la ecuación (1.1) por separado. En primer lugar, si anulamos el segundo sumando, obtenemos ROut(x, t)dx = P(u)(O, t), lo cual indica que la función P(u)(O, t), rige la variación de la cantidad total de uen el tiempo t. Por tanto, dependiendo de su signo, se encarga de crear (+) y destruir (−)células directamente. Al depender del volumen de control, escribiremos P(u)(O, t) = ZO s(u)(x, t)dx, (1.2) siendo s(u)(x, t)la densidad de variación de células en el instante t. Llamaremos a 15 CAPÍTULO 1. EL MODELO DE KELLER-SEGEL Demostración. En primer lugar sustituimos en la primera ecuación de (1.33) la función test φpor la parte negativa de u, esto es u−= m´ın{u, 0}. ZΩ utu−dx +ZΩ Du∇u· ∇u−dx −ZΩ λu+∇v· ∇u−dx = 0.(1.34) Desarrollando (1.34) y considerando las derivadas el sentido de las “derviadas generalizadas” tenemos ZΩ utu−dx +ZΩ Du∇u· ∇u−dx − :0 ZΩ λu+∇v· ∇u−dx = 0 ⇒1 2 d dt ZΩ (u−)2dx =−ZΩ Du|∇u−|2dx, (1.35) donde se ha tenido en cuenta que u+y∇u−tienen soportes disjuntos. Por tanto es evidente que se tiene d dt ZΩ (u−)2dx ≤0.(1.36) Sabiendo que u−(x, 0) = 0, por (1.36) tendremos que RΩ(u−)2= 0 para t > 0y por tanto, u−(x, t)=0en Ω×(0, T), lo que equivale a que u≥0. Consideramos ahora la segunda ecuación de (1.33) sustituyendo la función test por v−, se tendrá ZΩ vtv−dx +ZΩ Dv∇v· ∇v−dx =ZΩ αuv−dx −ZΩ βvv−dx. (1.37) De forma análoga a como hicimos para u, desarrollamos la ecuación (1.37). 1 2 d dt ZΩ (v−)2dx =ZΩ αuv−dx −ZΩ β(v−)2dx −ZΩ Dv|∇v−|2dx. Teniendo en cuenta que, como se dedujo antes u≥0, se tiene d dt ZΩ (v−)2dx ≤0.(1.38) Con un razonamiento equivalente al del caso de u−se deduce que v−(x, t)=0 en Ω×(0, T). Corolario 1.2.3. Suponiendo que (1.33) tiene solución, entonces existe una solución no negativa de Keller-Segel (formulación variacional (1.20)). 16 CAPÍTULO 1. EL MODELO DE KELLER-SEGEL Otra característica fundamental de este modelo es que presenta la siguiente ley de energía disipativa Proposición 1.2.4. Suponiendo que (u,v) es solución de la formulación variacional (1.20), se tendrá la siguiente relación: d dtE(u, v) = −D(u, v)≤0,(1.39) siendo E(u, v) = ZΩ (F(u) + λDv |∇v|2 2+βλv2 2−αλuv)dx, (1.40) un funcional al que denominaremos energía, con F(u) = αDu(ulog u−u).Y además D(u, v) = ZΩ (αu|∇(Dulog u−λv)|2+λv2 t)dx, (1.41) es un funcional no negativo, al que llamaremos disipación. Demostración. En primer lugar, se hace notar que la EDP (1.19a), se escribe de forma equivalente como ut−∇·(u∇(Dulog u−λv)) = 0,(1.42) donde hemos usado que u∇log u=∇u. Multiplicamos la ecuación (1.42) por α(Dulog u−λv)y (1.19b) por λvt, integremos en Ωy sumamos las ecuaciones. En las integrales siguientes se sobreentienden los dx y se omiten por no sobrecargar las ecuaciones. 0 = ZΩ αDulog u | {z } F′(u) ut−ZΩ αλutv−ZΩ ∇ · (u∇(Dulog u−λv))α(Dulog u−λv) +ZΩ λv2 t−ZΩ λDv∆vvt−ZΩ αλuvt+ZΩ βλvvt =d dt ZΩ F(u)−d dt ZΩ αλuv +d dt ZΩ βλv2 2+ZΩ λDv∇v· ∇vt−Z∂Ω  :0 (∇v·n)vt +ZΩ αu|∇(Dulog u−λv)|2−Z∂Ω u:0 (∇(Dulog u−λv)) ·n α(Dulog u−λv) +ZΩ λv2 t =d dt ZΩ (F(u) + λDv |∇v|2 2+βλv2 2−αλuv) + ZΩ (αu|∇(Dulog u−λv)|2+λv2 t), donde F(u) = αDu(ulog u−u)es una integral del logaritmo con sus constantes correspondientes. Finalmente, identificando los términos con sus definiciones se tiene el resultado. Capítulo 2 La aproximación en tiempo En este capítulo se trata de obtener una discretización temporal para nuestro modelo que preserve en la medida de lo posible las propiedades del problema continuo deducidas en la Sección 1.2. 2.1. Presentación del esquema En primer lugar se divide el intervalo espacial (0, T)en nodos equiespaciados {tn}Nt n=0 de paso dt =T/Nt, entonces tk=kdt y llamaremos uk≈u(·, tk). Definimos el siguiente esquema en tiempo:                  δtuk−Du∆uk+λ∇ · (uk∇vk) = 0.en Ω,∀k= 1, ..., Nt.(2.1a) δtvk−Dv∆vk−αuk−1+βvk= 0.en Ω,∀k= 1, ..., Nt.(2.1b) ∂uk ∂n = 0,∂vk ∂n = 0.en ∂Ω,∀k= 1, ..., Nt.(2.1c) u0=u0(x), v0=v0(x).en Ω.(2.1d) Se hace notar que los problemas (2.1a) y (2.1b) son lineales y desacoplados. Hemos utilizado el método de Euler implícito para la ecuación de u(1.19a) como se representa en (2.2). ut(tk)≈δtuk:= uk−uk−1 dt .(2.2) Sin embargo, para la ecuación de v(1.19b) se toma una aproximación semi-implícita donde se explicita el término de producción de vdebida a la variable u. Esta última consideración nos permitirá desacoplar las ecuaciones simplificando así nuestro esquema y además se comprueba en este capítulo que es una elección adecuada para tener conservación, no negatividad y una energía discreta decreciente. 18 CAPÍTULO 2. LA APROXIMACIÓN EN TIEMPO El algoritmo de un paso para la resolución del esquema (2.1a)-(2.1d) es el siguiente: 1: Begin 2: dados u0, v0de (2.1d) 3: for k= 1 to Ntdo 4: calcular vkde (2.1b) 5: calcular ukde (2.1a) 6: end for 7: End 2.2. Conservación de las propiedades En esta sección tratamos de comprobar si la aproximación resultante de la aplicación de nuestro esquema en tiempo conserva las propiedades de la solución del problema continuo. En primer lugar se mantiene la conservación de la cantidad total de u. Proposición 2.2.1. El esquema en tiempo conserva la cantidad total de u. Esto es, (uk,1) = (u0,1) para todo k= 1, ..., Nt. Demostración. Para comprobar esta propiedad, integramos la ecuación (2.1a) en Ω y tendremos 0 = ZΩ (uk−uk−1 dt −Du∆uk+λ∇ · (uk∇vk)) dx (1.25) =ZΩ uk−uk−1 dt dx +Z∂Ω (−Du∇uk+λuk∇vk)·n dS (2.1c) =ZΩ uk−uk−1 dt dx +Z∂Ω (−Du 70 δtuk ∂n +λuk  70 δtvk ∂n )dS = 0,∀k= 1, ..., Nt, de esto se sigue que ZΩ ukdx =ZΩ uk−1dx, ∀k= 1, ..., Nt.(2.3) De las ecuaciones (2.3) fácilmente se sigue el resultado. Corolario 2.2.1. Si uk≥0, se tendrá |uk∥L1(Ω) =∥u0∥L1(Ω).(2.4) Para la cantidad de v, que hemos llamado V(t) = RΩv(·, t)dx, tenemos el siguiente resultado: 19 CAPÍTULO 2. LA APROXIMACIÓN EN TIEMPO Proposición 2.2.2. El esquema en tiempo conduce al esquema de Euler implícito para la EDO que satisface Ven (2.1b). Esto es: δtVk+βV k−αU0= 0,∀k= 1, ..., Nt.(2.5) Demostración. Integramos la ecuación (2.1b) en Ωse tiene ZΩ (vk−vk−1 dt −Dv∆vk−αuk−1+βvk)dx (1.25),(2.1c) =Vk−Vk−1 dt −Z∂Ω Dv 70 δtvk ∂n dS −αUk−1+βV k∀k= 1, ..., Nt. Recordando la Proposición 2.2.1, se sigue que Uk−1=U0y se tiene el resultado. Veremos ahora si el esquema en tiempo impide la negatividad de las variables. En primer lugar análogamente al caso continuo, necesitamos la formulación variacional para nuestro esquema en tiempo (2.1a)-(2.1d). Esta se puede obtener directamente traduciendo (1.20) al esquema en tiempo. Asimismo truncaremos la velocidad quimiotáctica ukpor su parte positiva (uk)+.            ZΩ δtukφ dx +ZΩ (Du∇uk−λ(uk)+∇vk)· ∇φ dx = 0. ZΩ (δtvk−αuk−1+βvk)φ dx +ZΩ Dv∇vk· ∇φ dx = 0. u0=u0(x), v0=v0(x)en Ω. φ∈H1(Ω).(2.6) Proposición 2.2.3. Supongamos para el esquema truncado (2.6) que uk−1, vk−1≥0, entonces se tiene que uk, vk≥0. Demostración. Centrándonos en la primera ecuación de (2.6) tomamos como función test su parte negativa (uk)−, se tiene entonces ZΩ uk−uk−1 dt (uk)−dx +ZΩ Du∇uk· ∇(uk)−dx −ZΩ λ(uk)+∇vk· ∇(uk)−dx = 0 ⇒ZΩ ((uk)−)2dx =ZΩ uk−1(uk)−dx | {z } ≤0 −dt ZΩ Du|∇(uk)−|2dx | {z } ≤0 . Teniendo en cuenta la igualdad anterior, se tendrá que ZΩ ((uk)−)2dx ≤0 20 CAPÍTULO 2. LA APROXIMACIÓN EN TIEMPO ⇒(uk)−= 0.(2.7) Se tendrá finalmente que uk≥0en Ω. Veamos ahora la no negatividad de vk. Para ello procederemos de forma similar a la demostración de la propiedad del problema continuo y por tanto, tomando (vk)− como función test se sigue ZΩ vk−vk−1 dt (vk)−dx+ZΩ Dv∇vk·∇(vk)−dx =ZΩ αuk−1(vk)−dx −ZΩ βvk(vk)−dx, donde reordenando los términos tenemos ZΩ ((vk)−)2dx =ZΩ vk−1(vk)−dx | {z } ≤0 +dt(ZΩ αuk−1(vk)− | {z } ≤0 dx −ZΩ β((vk)−)2 | {z } ≤0 dx −ZΩ Dv|∇(vk)−|2dx |{z } ≤0 ). Se sigue que (vk)−= 0 y en consecuencia vk≥0en Ω. Corolario 2.2.2. Suponiendo que el esquema truncado (2.6) tiene solución, entonces el esquema sin truncar (2.1a)-(2.1d) posee solución no negativa. Finalmente, en lo referente a la ley de energía veremos que aparece un nuevo término; sin embargo, este término promueve la disipación de energía, lo que implica que la energía discreta en tiempo sigue siendo disipativa. Proposición 2.2.4. Si (uk, vk)son soluciones de (2.1a)-(2.1d), se tiene la siguiente relación: δtEk=−Dk−Dk num ≤0,(2.8) con Ek=ZΩ (F(uk) + λDv |∇vk|2 2+βλ(vk)2 2−αλukvk)dx, (2.9) y Dk=ZΩ (αuk|∇(Dulog uk−λvk)|2+λ(δtvk)2)dx, (2.10) las traducciones discretas de los funcionales de energía y disipación continuos. Pero además con otro nuevo término no negativo Dk num =dt 2ZΩ (F′′(ξ(x))(δtuk)2+λDv|δt∇vk|2+βλ(δtvk)2)dx, (2.11) donde ξ(x)∈(uk−1(x), uk(x)) para cada x∈Ω. 21 CAPÍTULO 2. LA APROXIMACIÓN EN TIEMPO Demostración. Multiplicamos (2.1a) por α(Dulog uk−λvk)y (2.1b) por λδtvk, integramos en Ωy sumamos las ecuaciones. En lo sucesivo omitimos los dx por no sobrecargar la notación. 0 = ZΩ αDulog uk | {z } F′(uk) δtuk−ZΩ αλδtukvk+ZΩ λ(δtvk)2−ZΩ λDv∆vkδtvk−ZΩ αλuk−1δtvk −ZΩ ∇ · (uk∇(Dulog uk−λvk))α(Dulog uk−λvk) + ZΩ βλvkδtvk =ZΩ F′(uk)δtuk−δtZΩ αλukvk+ZΩ αuk|∇(Dulog uk−λvk)|2+ZΩ λ(δtvk)2 +δtZΩ λDv |∇vk|2 2+dt ZΩ λDv |δt∇vk|2 2+δtZΩ βλ(vk)2 2+dt ZΩ βλ(δtvk)2 2, donde hemos usado las siguientes propiedades de la derivada discreta: (1) δtukvk+uk−1δtvk=δt(ukvk). (2) δtvkvk=1 2δt(vk)2+dt(δtvk)2. Es interesante señalar que para el uso de (1) ha sido imprescindible haber tomado el término de producción de vkcomo uk−1. Es aquí donde verdaderamente se justifica esta elección Finalmente, para relacionar F′(uk)δtukcon δtF(uk)realizamos una expansión en serie de Taylor hasta segundo orden de F(uk−1)con centro ukpara cada x∈Ω: F(uk−1) = F(uk)−F′(uk)(uk−uk−1) + F′′(ξ(x)) 2(uk−uk−1)2,(2.12) donde ξ(x)∈(uk−1(x), uk(x)) para cada x∈Ω. Reordenando los términos de (2.12), se obtiene fácilmente F′(uk)δtuk=δtF(uk) + dtF′′(ξ) 2(δtuk)2.(2.13) Con todo esto, identificando los sumandos se tiene fácilmente el resultado. Es importante notar que el término de disipación numérica Dk num que introduce el esquema en tiempo es de orden dt al igual que el orden de la aproximación de dicho esquema. Capítulo 3 Aproximación en espacio por elementos finitos y diferencias finitas En este capítulo se aborda el primer intento de discretizar en espacio guiados por la idea de obtener una discretización total de nuestro problema, en vistas de resolverlo numéricamente. Trabajaremos sobre un dominio espacial unidimensional y resolveremos el problema de forma numérica aplicando el método de elementos finitos. No se ha programado el método de diferencias finitas centradas como primer intento de discretización espacial ya que, como se demuestra en este capítulo, los dos métodos sobre nuestro modelo de Keller-Segel en 1D, son equivalentes. Para la realización de este capítulo nos ha sido útil [1]. 3.1. Diferencias finitas centradas Consideramos el problema de Keller Segel unidimensional con la discretización en tiempo. Para cada k= 1, ..., Ntse tiene (PCo)      δtuk−uk xx + ((uvx))k x= 0 δtvk−vk xx +vk−uk−1= 0 u0=u0(x), v0=v0(x). ,+ CC. Neumann. (3.1) Para la aproximación en diferencias finitas del esquema en tiempo (3.1), tomaremos una partición del intervalo de espacio Ωcon nodos {xn}N n=1 tal que xi= x1+ (i−1)hy se trata de definir ujcomo aproximación de u(xj,·). En primer lugar la aproximación de la derivada segunda se lleva a cabo de esta forma vk xx(xj)≈vk x(xj+1/2)−vk x(xj−1/2) h.(3.2) 23 CAPÍTULO 3. ELEMENTOS FINITOS Y DIFERENCIAS FINITAS Para evaluar (3.2), se precisa de las derivadas aproximadas en los puntos medios de nuestros nodos, estas se aproximan vk x(xj+1/2)≈vk j+1 −vk j h.(3.3) Con esta aproximación, (3.2) queda: vk xx(xj)≈vk j+1 −2vk j+vk j−1 h2.(3.4) Se tiene entonces la discretización completa de la ecuación para vkde (3.1) δtvk j−vk j+1 −2vk j+vk j−1 h2+vk j−uk−1 j= 0,(3.5) llevando los términos de conocidos al segundo término tenemos −1 h2vk j−1+ ( 1 dt +2 h2+ 1)vk j−1 h2vk j+1 =1 dtvk−1 j+uk−1 j.(3.6) Para los nodos frontera, consideramos las condiciones de contorno, (uk x(x1) = uk x(xN)=0, vk x(x1) = vk x(xN) = 0.(3.7) Para tenerlas en cuenta, utilizaremos las derivadas primeras centradas con puntos auxiliares x0=x1−hyxN+1 =xN+h. La aportación de estos nodos, más allá de permitirnos imponer las condiciones de contorno, será nula ya que los valores que toman las variables sobre ellos quedarán ligados a los valores de estas en los puntos interiores mediante las condiciones de contorno. Por tanto, el número de grados de libertad no aumentará y seguirá siendo N(ver Figura 3.1). ... x0x1x2xN−1xNxN+1 Figura 3.1: Intervalo espacial. Las derivadas primeras centradas serán por tanto vk x(xj)≈vk j+1 −vk j−1 2h.(3.8) Aplicando (3.8) a (3.7) junto con las condiciones de contorno, tendremos para el caso de j= 1, N, 0 = vk x(x1) = vk 2−vk 0 2h⇒vk 0=vk 2, 0 = vk x(xN) = vk N+1 −vk N−1 2h⇒vk N+1 =vk N−1. (3.9) 24 CAPÍTULO 3. ELEMENTOS FINITOS Y DIFERENCIAS FINITAS Por tanto, introduciendo (3.9) en (3.2), obtenemos vk xx(x1) = vk 2−2vk 1+vk 0 h2=2vk 2−2vk 1 h2, vk xx(xN) = vk N+1 −2vk N+vk N−1 h2=−2vk N+ 2vk N−1 h2. (3.10) Con lo que la ecuaciones de la discretización de los nodos frontera ya despejadas quedarán (1 dt +2 h2+ 1)vk 1−2 h2vk 2=uk−1 1+1 dtvk−1 1. −2 h2vk N−1+ ( 1 dt +2 h2+ 1)vk N=uk−1 N+ + 1 dtvk−1 N. (3.11) Escribiremos matricialmente el sistema de ecuaciones (3.6) ya que nos será útil posteriormente para comparar este esquema con el de elementos finitos. Avvk=Bv,(3.12) siendo vk:= (vk 1, ..., vk N)T,(3.13) el vector de los valores de las aproximaciones de vken los nodos. Además con Av= ( 1 dt + 1)I+Adiff ,(3.14) la matriz de coeficientes, compuesta por la identidad Iy Adiff =                     2 h2−2 h2 −1 h2 2 h2−1 h2 ... −1 h2 2 h2−1 h2 −2 h2 2 h2                     (3.15) la matriz que se corresponde con el término de autodifusión. Y finalmente Bv=1 dtvk−1+uk−1,(3.16) el segundo miembro. 31 CAPÍTULO 3. ELEMENTOS FINITOS Y DIFERENCIAS FINITAS de difusión se corresponderá con la matriz de rigidez K. Ahora bien, aparecerá el nuevo término de quimiotaxis cuya matriz resultante tendremos que estudiar como en casos anteriores. Por tanto el sistema para ukresulta Ak uuk=Bu,(3.50) siendo Ak u=1 dtM+K−Qk,(3.51) la matriz de coeficientes y KyMlas matrices de rigidez y masa descritas anteriormente. Qkes la matriz que tiene en cuenta el término de quimiotaxis, los elementos de esta matriz serán qk ij = ((φi)x, φj N X l=1 vk l(φl)x!) = X Ik N X l=1 vk lZIk (φi)xφj(φl)xdx. (3.52) Por otro lado, el segundo miembro será Bu=1 dtMuk−1.(3.53) Evaluemos entonces los elementos de la matriz Q. Es evidente que la matriz es tridiagonal, ya que si las funciones de base no tienen índices consecutivos o iguales, sus soportes son disjuntos. Esta reflexión nos lleva a que dado un i, los elementos no nulos se dan para j, l ∈ {i−1, i, i + 1}. Se calculan entonces los términos de la matriz para i= 1, N. qk ii = i X l=i−1 vk lZIi−1 (φi)xφi(φl)xdx + i+1 X l=i vk lZIi (φi)xφi(φl)xdx =1 h i X l=i−1 vk lZIi−1 φi(φl)xdx −1 h i+1 X l=i vk lZIi φi(φl)xdx =1 h−1 hvk i−1ZIi−1 φidx +1 hvk iZIi−1 φidx−1 h−1 hvk iZIi φidx +1 hvk i+1 ZIi φidx =1 h−1 hvk i−1 h 2+1 hvk i h 2−1 h−1 hvk i h 2+1 hvk i+1 h 2=−vk i+1 −2vk i+vk i−1 2h. qii−1= i X l=i−1 vk lZIi−1 (φi)xφi−1(φl)xdx =1 h i X l=i−1 vk lZIi−1 φi−1(φl)xdx =1 h−1 hvk i−1ZIi−1 φi−1dx +1 hvk iZIi−1 φi−1dx=1 h−1 hvk i−1 h 2+1 hvk i h 2 =vk i−vk i−1 2h. 32 CAPÍTULO 3. ELEMENTOS FINITOS Y DIFERENCIAS FINITAS qii+1 = i+1 X l=i vk lZIi (φi)xφi+1(φl)xdx =−1 h i+1 X l=i vk lZIi φi+1(φl)xdx =−1 h−1 hvk iZIi φi+1dx +1 hvk i+1 ZIi φi+1dx=−1 h−1 hvk i h 2+1 hvk i+1 h 2= =−vk i+1 −vk i 2h, donde se ha tenido en cuenta que las integrales RIiφidx corresponden al área del triángulo rectángulo de catetos hy1. Restan solo los elementos correspondientes a las fronteras del dominio. Las fórmulas anteriores son válidas para los valores de q12 yqNN−1. qk 12 =−vk 2−vk 1 2h. qk NN−1=vk N−vk N−1 2h. Los elementos diagonales son qk 11 =−vk 2−vk 1 2h. qk NN =vk N−vk N−1 2h. Con todo esto podemos presentar la matriz Qk. Qk=                           −vk 2−vk 1 2h−vk 2−vk 1 2h vk 2−vk 1 2h−vk 3−2vk 2+vk 1 2h−vk 3−vk 2 2h ... vk N−1−vk N−2 2h−vk N−2vk N−1+vk N−2 2h−vk N−vk N−1 2h vk N−vk N−1 2h vk N−vk N−1 2h                           (3.54) Tras conocer lo términos de esta última matriz podríamos desarrollar explícitamente el sistema de ecuaciones presente en (3.50). Proposición 3.2.2. Los esquemas de diferencias finitas centradas (3.27) de elementos finitos (3.36) con “Mass lumping” en 1D son equivalentes. 33 CAPÍTULO 3. ELEMENTOS FINITOS Y DIFERENCIAS FINITAS Demostración. Se ve fácilmente cumplen las siguientes relaciones: M=M I, K=M Adiff , Qk=−M Achem. (3.55) Las igualdades (3.55) indican que las ecuaciones resultantes de los dos métodos son equivalentes. A la vista de estas, las ecuaciones del método de elementos finitos son las de diferencias finitas multiplicadas por hpara los nodos interiores y por h/2 para los nodos frontera. Cabe destacar que el “Mass lumping” ha sido crucial para tener esta equivalencia. 3.3. Simulación y análisis de los resultados. Oscilaciones espúreas El código que se ha utilizado para conseguir las aproximaciones, se ha hecho con la metodología de elementos finitos, realizando un barrido por elementos para construir las diferentes matrices. El script de MATLAB en el que se ha implementado el método se encuentra en el Apéndice B Para el intervalo espacial tomamos Ω = (−1,1) y para el temporal (0,0,5). Tomamos los valores de los parámetros como Du=Dv=β= 1, λ = 10 yα= 50. Asimismo como valores iniciales tomaremos u0(x) = 10 (x+ 1)2(x−1)2yv0(x) = 10 sen2(2πx). Para un primer cálculo tomamos los valores de paso espacial h= 0,001 y paso temporal dt = 0,002. -1 -0.5 0 0.5 1 x 0 2 4 6 8 10 u (a) u0 -1 -0.5 0 0.5 1 x 0 2 4 6 8 10 v (b) v0 Figura 3.3: Condiciones iniciales. 34 CAPÍTULO 3. ELEMENTOS FINITOS Y DIFERENCIAS FINITAS Esta solución evolucionará en el tiempo sin errores apreciables a simple vista pero como ilustra la Figura 3.4, llega un punto en el que la aproximación empieza a tomar valores negativos. Aún así, estos son pocos y prácticamente inapreciables en la imagen. Además la escasez de estos permite al programa continuar aproximando la solución y volviendo incluso a eliminar la negatividad para siguientes pasos de tiempo. -1 -0.5 0 0.5 1 x -200 0 200 400 600 800 1000 1200 u (a) u174 -1 -0.5 0 0.5 1 x 40 60 80 100 120 140 v (b) v174 Figura 3.4: Aproximación para t174 ≈0,3. Posteriormente, en la Figura 3.5 se aprecia cómo el método empieza a aproximar la solución cometiendo grandes errores, los cuales se van acumulando y magnificando a media que se obtienen las aproximaciones para siguientes etapas de tiempo. Todo esto acaba dando lugar a aproximaciones erróneas, las cuales presentan oscilaciones espúreas e impropias de la solución exacta. -1 -0.5 0 0.5 1 x -2 -1 0 1 2 3 4 u 104 (a) u211 -1 -0.5 0 0.5 1 x 40 60 80 100 120 140 160 180 v (b) v211 Figura 3.5: Aproximación para t211 ≈0,4. 35 CAPÍTULO 3. ELEMENTOS FINITOS Y DIFERENCIAS FINITAS Además de las gráficas de las propias funciones uyv, es interesante estudiar cómo nuestro modelo aproxima las demás propiedades del problema continuo, los cuales nos servirán como indicadores adicionales de la bondad de la aproximación. Con este objetivo, se muestran las gráficas de la cantidad total de células y sustancia química, y las de la energía. Evidentemente, será necesario definir estos conceptos para nuestra aproximación. En primer lugar, como cada uk, vk∈Vh, denotaremos a la cantidad total de u(v) en el tiempo tkcomo Ck u=RΩukdx (Ck v=RΩvkdx). Como ukyvkson loclamente segmentos y globalmente continuas, la integral se puede calcular con la regla de los trapecios Ck u=h 2(u1+uN) + h N−1 X j=2 uj.(3.56) Representaremos entonces la cantidad total de uyven función del tiempo 0 0.1 0.2 0.3 0.4 0.5 t 10 10.5 11 11.5 Cu Aproximación Teórica (a) Cantidad total de u 0 0.1 0.2 0.3 0.4 0.5 t 0 50 100 150 200 250 Cv Aproximación Teórica (b) Cantidad total de v Figura 3.6: Cantidad total. Como se aprecia en la Figura 3.6, el esquema reproduce correctamente la conservación de la cantidad total de u, como se demostró en la Proposición 3.2.1 y a su vez, parece aproximar correctamente la dinámica que sigue la cantidad total de v, al menos con estos datos. Por otro lado, con una filosofía similar a la de la cantidad total, calcularemos la energía discreta Ekde la energía asociada al problema que se definió en (2.9). Para calcularla, encontramos un problema adicional, el valor de (vk x)2no está definido en los nodos. Para ello, se aproxima vk xen cada nodo de forma centrada como se 36 CAPÍTULO 3. ELEMENTOS FINITOS Y DIFERENCIAS FINITAS presentó en (3.8). Con esto, podemos aplicar la fórmula de los trapecios. Ek=h 2(F(u1) + F(uN) + βλv2 1+v2 N 2)−αλ(uk 1vk 1+uk Nvk N)) +h N−1 X j=2 (F(uj) + λDv 1 2vj+1 −vj−1 2h2 +βλv2 j 2−αλuk jvk j) (3.57) Asimismo aproximaremos la disipación de un modo similar: Dk=h 2λ(δtvk 1)2+ (δtvk N)2 +h N−1 X j=2 αuk j Du log(uk j+1)−log(uk j−1) 2h−λ(vk j+1)−(vk j−1) 2h!2 +λ(δtvk j)2. (3.58) Representaremos dos gráficas en una se representan en función del tiempo la aproximación de la energía Ek, y en la otra, una aproximación de la ley de energía, que llamaremos LEk=δtEk+Dkpara k > 0. Idealmente la primera deberá ser decreciente y la segunda no positiva, ya que como vimos, el esquema en tiempo introducía disipación numérica. Se toma un intervalo de tiempo donde la aproximación aún no toma valores negativos, ya que para estos, ni el término de energía ni el de disipación están definidos. 0 0.05 0.1 0.15 0.2 0.25 0.3 t -5 -4 -3 -2 -1 0 E 105 (a) Energía 0 0.05 0.1 0.15 0.2 0.25 0.3 t -2 0 2 4 6 8 10 LE 106 (b) Ley de energía Figura 3.7: Gráficas de la energía. En la gráfica 3.7a se observa cómo a priori sí se tiene que la energía es decreciente en las etapas anteriores a alcanzar valores negativos. Sin embargo, la gráfica 3.7b nos dice que la dinámica que representa la ley de energía discreta, está tomando valores positivos en algunos tiempos. Este hecho nos está indicando que el esquema en espacio está introduciendo fuente numérica respecto de la ley de energía. Si la fuente numérica que es capaz de introducir el esquema en espacio superase a nuestros dos 37 CAPÍTULO 3. ELEMENTOS FINITOS Y DIFERENCIAS FINITAS términos de disipación Dk, ocurriría que la energía aproximada no sería decreciente, aunque este efecto en nuestro ejemplo no se observa y tampoco se ha demostrado que pueda o no pueda ocurrir. En el Capítulo 2 vimos que la aproximación en tiempo impedía la negatividad, entonces es interesante comprobar numéricamente si el responsable de esta es, efectivamente, la aproximación en espacio. Para ello, observaremos la solución calculada cambiando el valor de h. -1 -0.5 0 0.5 1 x -6000 -4000 -2000 0 2000 4000 6000 u (a) uend -1 -0.5 0 0.5 1 x 50 100 150 200 v (b) vend Figura 3.8: Aproximación en tiempo final para h′= 5h. Al tomar un hmayor, como se aprecia en la Figura 3.8, se pone de manifiesto que la aproximación en espacio es la culpable de la negatividad, ya que esta presenta un comportamiento totalmente caótico. 0 0.1 0.2 0.3 0.4 0.5 t 10 10.5 11 11.5 Cu Aproximación Teórica (a) Cantidad total de u 0 0.1 0.2 0.3 0.4 0.5 t 0 50 100 150 200 250 Cv Aproximación Teórica (b) Cantidad total de v Figura 3.9: Cantidad total para h′= 5h. Es interesante señalar tal y como se aprecia en la Figura 3.9, que incluso dando 38 CAPÍTULO 3. ELEMENTOS FINITOS Y DIFERENCIAS FINITAS lugar a soluciones sin sentido, el esquema sigue reproduciendo la conservación de la cantidad total de u. Esto es, como se vio en Proposición 3.2.1 fruto de la propia naturaleza del esquema de elementos finitos. Asimismo, un comportamiento regular de la cantidad total de vtambién se reproduce. Para este caso, las gráficas de la energía no se representan ya que, para este valor de h, las variables son negativas en tiempos bastante cercanos al inicial. -1 -0.5 0 0.5 1 x 0 2000 4000 6000 8000 10000 u (a) uend -1 -0.5 0 0.5 1 x 60 80 100 120 140 160 180 200 v (b) vend Figura 3.10: Aproximación en tiempo final para h′=h/2. Tomando ahora un hmenor, como se ve en la Figura 3.10, ya no se observan valores negativos y se corrobora el hecho de que el causante de la negatividad es la aproximación en espacio. Las cantidades totales aproximadas siguen cumpliendo correctamente lo esperado para la solución del problema continuo, por lo que se omiten las gráficas ya que son similares al caso anterior. 0 0.1 0.2 0.3 0.4 0.5 t -7 -6 -5 -4 -3 -2 -1 0 E 105 (a) Energía 0 0.1 0.2 0.3 0.4 0.5 t -0.5 0 0.5 1 1.5 2 2.5 LE 107 (b) Ley de energía Figura 3.11: Gráficas de la energía. 39 CAPÍTULO 3. ELEMENTOS FINITOS Y DIFERENCIAS FINITAS Finalmente, como se aprecia en la Figura 3.11, se sigue teniendo que la energía es decreciente y que el esquema en espacio es capaz de introducir fuente numérica en algunos tiempos. Capítulo 4 Diferencias finitas descentradas Para lograr unas soluciones no negativas, salta a la vista por los resultados obtenidos en la Sección 3.3 que necesitamos modificar la aproximación en espacio de forma adecuada. 4.1. La técnica del descentrado Se propone una aproximación descentrada para el término de quimiotaxis. Este primer intento de descentrado se realiza en el caso de las diferencias finitas; esto es, sobre la misma partición del dominio introducida en la Sección 3.1 y considerando una partición uniforme con nodos auxiliares, ilustrada en la Figura 3.1. Asimismo, se considerarán las mismas ligaduras impuestas por las condiciones de contorno sobre las variables auxiliares, las cuales están recogidas en las ecuaciones (3.9) y (3.21). Por tanto, nuestra misión consiste en estimar de forma distinta el valor de (uvx)k x(xj). Comenzando del mismo modo que en el caso anterior, haremos uso de la ecuación (3.17). Ya con esto, el trabajo se reduce a conseguir una nueva aproximación descentrada del flujo (uvx)k j±1/2. Ya con todo esto se considera la siguiente aproximación descentrada (uvx)k j+1/2≈ vk j+1 −vk j h!+ uj+ vk j+1 −vk j h!− uk j+1.(4.1) Para el caso de ((uvx)k)j−1/2basta con considerar la ecuación anterior reemplazando jpor j−1. (uvx)k j−1/2≈ vk j−vk j−1 h!+ uk j−1+ vk j−vk j−1 h!− uk j.(4.2) 47 CAPÍTULO 4. DIFERENCIAS FINITAS DESCENTRADAS -1 -0.5 0 0.5 1 x 0 1000 2000 3000 4000 5000 6000 u (a) uend -1 -0.5 0 0.5 1 x 60 80 100 120 140 160 180 200 v (b) vend Figura 4.5: Aproximación descentrada en tiempo final para h′=h/2. En las Figuras 4.4 y 4.5, podemos observar que los máximos de las aproximaciones de uson claramente distintos. A continuación se desarrolla una idea de por qué el descentrado está ocasionando este efecto. Dado una etapa de tiempo k−1, supongamos que tenemos la configuración de la Figura 4.6 vk−1 uk−1 xi−1xixi+1 Figura 4.6: Configuración en tk−1. Supongamos, por simplicidad en la argumentación que aproximásemos la quimiotaxis de forma explícita y que además únicamente tenemos este término, ya que el término de autodifusión es el mismo en las dos aproximaciones. Se hace notar que en esta configuración se tiene que (vk−1 x)j−1/2>0y(vk−1 x)j+1/2<0. Con esto, para 48 CAPÍTULO 4. DIFERENCIAS FINITAS DESCENTRADAS la aproximación centrada, el valor de uk jse obtiene de (δtuk)cen =(vk−1 x)j−1/2 h uk−1 j+uk−1 j−1 2−(vk−1 x)j+1/2 h uk−1 j+1 +uk−1 j 2,(4.15) mientras que en el descentrado se tiene de (δtuk)upw =(vk−1 x)+ j−1/2 huk−1 j−1−(vk−1 x)− j+1/2 huk−1 j+1 .(4.16) Se ve claramente que: (δtuk)upw =(vk−1 x)j−1/2 huk−1 j−1−(vk−1 x)j+1/2 huk−1 j+1 <(vk−1 x)j−1/2 h uk−1 j+uk−1 j−1 2−(vk−1 x)j+1/2 h uk−1 j+1 +uk−1 j 2= (δtuk)cen. (4.17) Esto repercute en la aproximación de ukde una forma parecida a la que se observa en la Figura 4.7 uk cent uk upw xi−1xixi+1 Figura 4.7: Aproximación de ukdurante la agregación. Con esta idea, es fácil ver que en el caso de la formación de un valle ocurre lo que se representa en la Figura 4.8 49 CAPÍTULO 4. DIFERENCIAS FINITAS DESCENTRADAS uk upw uk cent xi−1xixi+1 Figura 4.8: Aproximación de uken un valle. A la vista de las Figuras 4.7 y 4.8, se ve que el descentrado produce en la solución una especie de “amortiguamiento” en la aproximación de máximos crecientes y mínimos decrecientes, esto es otra disipación numérica en un sentido distinto al caso de la ley de energía. Esta es la principal justificación de lo que se aprecia en las gráficas de uen las Figuras 4.4 y 4.5. Capítulo 5 Volúmenes finitos descentrados En este capítulo, se lleva a cabo un nuevo enfoque para la aproximación en espacio. Esta se basa en la adaptación a dimensión 1 del esquema de volúmenes finitos presentado en [9]. Consideramos junto con nuestra malla habitual {xi}N i=1 de equiespaciado h, una malla dual {yi}N+1 i=1 cuyos nodos interiores se corresponderán con los puntos medios de los nodos de la malla principal y en la frontera coincidirán las dos mallas, esto es    yi=xi+xi−1 2para i= 2, .., N. y1=x1, yN+1 =xN. (5.1) Además de esto, definimos nuestros “volúmenes” como Di= (yi, yi+1),(5.2) que se corresponden con los intervalos que unen a nuestros nuevos nodos yi(ver Figura 5.1). Esta es la reducción a dimensión 1de la idea de obtener una malla de volúmenes de Voronoi a partir de una malla triangular. D1D2... DN−1DN x1 y1y2x2y3yN−1xN−1yN yN+1 xN Figura 5.1: Dominio espacial y volúmenes. Como se puede observar fácilmente en la Figura 5.1, cada volumen Diestá asociado al nodo xi. Por este motivo aproximaremos las soluciones de nuestro problema en los nodos principales xjintegrando las ecuaciones de 3.1 en los diferentes volúmenes, en los que se considera que la aproximación es constante para cada uno de ellos. 51 CAPÍTULO 5. VOLÚMENES FINITOS DESCENTRADOS Comenzaremos como es habitual con la ecuación correspondiente a la variable vky la integraremos en el volumen Djcon j= 1, N. Aproximando por constantes vjen cada volumen Djtendremos ZDj (δtvk−vk xx +vk−uk) = 0 ⇒ h(δtvk j+vk j−uk j)−ZDj vk xx = 0.(5.3) Es por tanto necesario evaluar la integral RDjvk xx. Resulta conveniente, sabiendo que vk xx = (vk x)x, aplicar la regla de Barrow. ZDj vk xx =vk x(yj+1)−vk x(yj).(5.4) Esto plantea un nuevo problema, deberemos entonces aproximar el valor de vx en los puntos medios de nuestra malla principal. Una buena aproximación es tomar la derivada centrada vk x(yj)≈vk(xj)−vk(xj−1) h=vk j−vk j−1 h.(5.5) Introduciendo entonces 5.5 en 5.4 y esta a su vez en 5.3, tendremos h(δtvk j+vk j−uk j)− vk j+1 −vk j h−vk j−vk j−1 h!= 0 ⇒ h(δtvk j+vk j−uk j) + −vk j+1 + 2vk j−vk j−1 h= 0,(5.6) que es el esquema de volúmenes finitos centrados. Resta por tanto comprobar qué ocurre para los volúmenes frontera D1yDN, teniendo en cuenta que la longitud de estos es h/2. Comencemos con D1 ZD1 (δtvk−vk xx +vk−uk) = 0 ⇒h 2(δtvk 1+vk 1−uk 1)−ZD1 vk xx = 0 ⇒ ⇒h 2(δtvk 1+vk 1−uk 1) + vk x(y1)−vk x(y2) = 0. Teniendo en cuenta la condición de contorno de la frontera izquierda para v, se tendrá que vx(y1)=0. h 2(δtvk 1+vk 1−uk 1) + *0 vk x(y1)−vk x(y2) = 0 ⇒ 52 CAPÍTULO 5. VOLÚMENES FINITOS DESCENTRADOS ⇒h 2(δtvk 1+vk 1−uk 1)−vk 2−vk 1 h= 0.(5.7) Con la misma idea se tiene para el volumen DN h 2(δtvk N+vk N−uk N) + vk N−vk N−1 h= 0.(5.8) A la vista de los desarrollos, se llega al siguiente resultado: Proposición 5.0.1. El esquema de volúmenes finitos centrados para aproximar vk coincide con el de elementos finitos y con el de diferencias finitas en 1D. Procederemos ahora con la ecuación de u. En primer lugar, es evidente que la aproximación del término de autodifusión es análoga al caso de v. Por tanto, nuestra única labor será aproximar el término de quimiotaxis. Empecemos con j= 1, N. Con la misma idea que antes ZDj ((uvx))k x= (uvx)k(yj+1)−(uvx)k(yj).(5.9) Necesitamos entonces aproximar, (uvx)k(yj), para ello utilizaremos la técnica del descentrado “upwind” (uvx)(yj)≈vk x(yj)+uk j−1+vk x(yj)−uk j= vk j−vk j−1 h!+ uk j−1+ vk j−vk j−1 h!− uk j. (5.10) Introduciendo por tando (5.10) en (5.9), se tiene ZDj ((uvx))k x= (uvx)k(yj+1)−(uvx)k(yj) = vk j+1 −vk j h!+ uk j+ vk j+1 −vk j h!− uk j+1 − vk j−vk j−1 h!+ uk j−1− vk j−vk j−1 h!− uk j= − vk j−vk j−1 h!+ uk j−1+ vk j+1 −vk j h!+ − vk j−vk j−1 h!−!uK j+ vk j+1 −vk j h!− uk j+1. (5.11) Se comprueba ahora el caso para j= 1, donde se tiene en cuenta la condición de contorno en la frontera izquierda ya que y1=x1, y por tanto vx(y1) = vx(x1)=0. ZD1 ((uvx))k x= (uvx)k(y2)− :0 (uvx)k(y1) = vk 2−vk 1 h+ uk 1+vk 2−vk 1 h− uk 2.(5.12) 53 CAPÍTULO 5. VOLÚMENES FINITOS DESCENTRADOS Considerando ahora la condición de contorno de la frontera derecha, se tendrá de forma análoga ZDN ((uvx))k x=−vk N−vk N−1 h+ uk N−1−vk N−vk N−1 h− uk N.(5.13) Proposición 5.0.2. El esquema de volúmenes finitos descentrados para aproximar ukes equivalente al de diferencias finitas descentradas en 1D. Demostración. El sentido de esta equivalencia es el mismo que el de la Proposición 3.2.2, las ecuaciones de los nodos frontera están multiplicados por h/2y las de los nodos interiores por h. Capítulo 6 Elementos finitos descentrados En este capítulo se comenta y particulariza a dimensión 1en espacio el esquema presentado en [7]. La obtención de dicho esquema se lleva a cabo mediante un enfoque distinto a los casos anteriores. Comprobaremos si este nuevo esquema equivale al esquema anterior. En primer lugar partimos de la formulación débil de nuestro sistema presentada en (3.36). Y todos los elementos presentados para el método de elementos finitos (ver Figura 3.2), así como Vhel espacio generado por {φn}N n=1. Con esto en cuenta definimos una malla dual, la cual, casualmente, en contraposición al caso de dimensiones superiores, coincide con la de volúmenes finitos en 1D (ver Figura (5.1)). Consideramos ahora 1 Dila función característica de Di. Asimismo, introducimos el espacio de Hilbert Vh⊂L2(Ω) el espacio generado por { 1 Dn}N n=1; es decir, las funciones constantes a trozos en la malla dual. Con todo esto introducimos el operador Mh:Vh→Vhdefinido como Mhvh= N X i=1 vh(xi) 1 Di,(6.1) el cual es llamado operador de “lumping”. Y se define (uh, vh)h= (Mhuh, Mhvh),(6.2) un nuevo producto escalar sobre Vh. Por tanto, con estas definiciones aproximaremos la ecuación para vde (3.36) como (δtvk h,¯u)h+ ((vk h)x,¯ux)+(vk h,¯u)h−(uk−1 h,¯u)h= 0,∀¯u∈Vh.(6.3) A partir de aquí, tomando como funciones test las φi, obtenemos un sistema de ecuaciones igual al presentado en (3.39). 55 CAPÍTULO 6. ELEMENTOS FINITOS DESCENTRADOS En primer lugar matriz de rigidez Kes la misma que en FEM (3.44), ya que esta surge del término ((vh)x,¯ux), el cual es idéntico al que se tiene en el método de elementos finitos. Por otro lado, los elementos de la nueva matriz de masa Mprovienen de los términos de la forma (ah, φh)hy como Mhφj= 1 Dj, se tiene que mij = (Mhφj, Mhφi) = ( 1 Dj, 1 Di) = |Di|δij.(6.4) Como |D1|=|DN|=h/2y|Di|=hpara i= 1, N, se tendrá que la matriz de masa que obtenemos es la misma que la que se obtiene con la técnica de “Mass lumping”. Con todo esto, la discretización espacial de vkpor este método de elementos finitos es la misma que aplicando los demás métodos. En definitiva, este método para la solución vkpresenta una manera diferente de realizar la técnica del “Mass lumping”. Procedemos ahora para u, su ecuación en (3.36) aplicamos el “Mass lumping” en el sumando de la derivada temporal (δtuk h,¯u)h+ ((uk h)x,¯ux)−(uk h(vk h)x,¯ux) = 0,∀¯u∈Vh.(6.5) Cuando tomemos como funciones test las φi, se consiguen términos vistos en el sistema de ecuaciones (3.50). La matrices de masa My rigidez Kserán las mismas que para el caso anterior para vk, ya que sus elementos se calculan de la misma forma. Lo interesante es por tanto aproximar el término de quimiotaxis. −(uk h(vk h)x,¯ux) = ZΩ (uk h(vk h)x)x¯u dx (Integración por partes) ≈ N X j=1 ¯u(xj)ZDj (uk h(vk h)x)xdx (¯u≈Mh¯u) = N X j=1 ¯u(xj)(uk h(vk h)x)j+1/2−(uk h(vk h)x)j−1/2,(Regla de Barrow) (6.6) donde para los nodos frontera se tiene que (uk h(vk h)x)1/2= (uk h(vk h)x)1 (uk h(vk h)x)N+1/2= (uk h(vk h)x)N.(6.7) Ahora es cuando aplicamos nuestro descentrado “upwind” (uk h(vk h)x)upw j+1/2≈(uk h)j((vk h)x)+ j+1/2+ (uk h)j+1((vk h)x)− j+1/2.(6.8) 56 CAPÍTULO 6. ELEMENTOS FINITOS DESCENTRADOS Podemos evaluar los términos de (6.8) como elementos del espacio Vh. Por tanto, dado un j= 1, N se tendrá que: ((vk h)x)j+1/2= (vk h)xxj+xj+1 2 =vk j(φj)x+vk j+1(φj+1)x=vk j(−1 h) + vk j+1(1 h) =vk j+1 −vk j h, (6.9) en los nodos frontera, aplicando las condiciones de contorno (vk h)x= 0. Finalmente tenemos (uk h(vk h)x)upw j+1/2=uj vk j+1 −vk j h!+ +uk j+1 vk j+1 −vk j h!− .(6.10) En este momento es cuando sustituimos ¯upor φi∈Vhpara i= 1, ..., N en (6.6), y vamos obteniendo los elementos de la matriz Qupw de quimiotaxis. Tomando φi como función test e al introduciendo (6.10) en (6.6), se tiene N X j=1 δij (uk h(vk h)x)upw j+1/2−(uk h(vk h)x)upw j−1/2 =−vk i−vk i−1 h+ uk i−1+ vk i+1 −vk i h+ −vk i−vi−1 h−!uk i+vk i+1 −vk i h− uk i+1. (6.11) Esto se desarrolla a continuación. Introduciendo φ1como función test en (6.6) se tiene (uk h(vk h)x)upw 3/2−(uk h(vk h)x)upw 1CC = (uk h(vk h)x)upw 3/2 (6.9) =vk 2−vk 1 h+ uk 1+vk 2−vk 1 h− uk 2,(6.12) y para φNse tiene −vk N−vk N−1 h+ uk N−1−vk N−vk N−1 h− uk N.(6.13) Con todo esto, el esquema completamente discretizado de elementos finitos descentrado para la quimiotaxis queda      (δtuk h,¯u)h+ ((uk h)x,¯ux) + aupw h(vk h, uk h,¯u)=0, (δtvk h,¯u)h+ ((vk h)x,¯ux)+(vk h,¯u)h−(uk−1 h,¯u)h= 0, u0 h=PN j=1 u0(xj)φj(x), v0 h=PN j=1 v0(xj)φj(x). ∀¯u∈Vh.(6.14) 63 CAPÍTULO 7. OTROS MODELOS CON QUIMIOTAXIS -1 -0.5 0 0.5 1 x 0 2 4 6 8 10 12 u (a) u51 -1 -0.5 0 0.5 1 x 9 10 11 12 13 14 v (b) v51 Figura 7.1: Aproximación en t51 = 0,1. En la Figura 7.1, se observa cómo la logística provoca la creación de picos de altura similar en la aproximación de u. Aún así estos no son producto de la reacción logística exclusivamente, ya que se encuentran por encima de la capacidad δ= 5, por tanto este efecto es una combinación de la reacción logísitca junto con la quimiotaxis. -1 -0.5 0 0.5 1 x 0 5 10 15 20 25 u (a) u101 -1 -0.5 0 0.5 1 x 10 12 14 16 18 20 v (b) v101 Figura 7.2: Aproximación en t51 = 0,2. En la Figura 7.2, se observa cómo se produce una agregación de umucho menor que en el caso sin la logística, y lo mismo ocurre para v, la producción de ves mucho menor debido a que lo es u. Cabe decir que la reducción de la altura era un efecto esperable ya que la capacidad del medio es pequeña con respecto a los valores que alcanza la agregación. También cabe señalar que no se tiene a diferencia que en el caso centrado la negatividad. Por último, otra diferencia respecto del modelo de Keller-Segel es que en este caso, la agregación final se produce antes con el término logístico (t≈0,2) que sin él (t≈0,35). El momento en el que se produce la 64 CAPÍTULO 7. OTROS MODELOS CON QUIMIOTAXIS agregación final en Keller-Segel se puede observar en la gráfica de la energía de la Figura (4.3), cuando esta decrece bruscamente. -1 -0.5 0 0.5 1 x 0 2 4 6 8 u (a) uend -1 -0.5 0 0.5 1 x 11 12 13 14 15 16 v (b) vend Figura 7.3: Aproximación en tiempo final. En la última etapa de tiempo, como se aprecia en la Figura 7.3, se observa cómo el término logístico ha ido dominando desde que se produjo la agregación, disminuyendo el pico de u, además esto ha ocasionado que la degradación de vtambién domine. Aún así, se ha observado que la dinámica va disminuyendo, lo que parece indicar que estas aproximaciones están cerca de la solución estacionaria con estos parámetros. Esto nos lleva a pensar que podemos llegar a soluciones estacionarias no constantes, lo que se conoce como el fenómeno de formación de patrones. Asimismo, se ha comprobado que variando la tasa de crecimiento logístico r la solución estacionaria de uparece crear patrones distintos como se aprecia en la Figura 7.4, o incluso para valores de r≫λparece estar regida por la logística y corresponderse con la capacidad, este último efecto se aprecia en la Figura 7.5. 65 CAPÍTULO 7. OTROS MODELOS CON QUIMIOTAXIS -1 -0.5 0 0.5 1 x 0 5 10 15 20 25 u (a) uend -1 -0.5 0 0.5 1 x 32 34 36 38 40 42 v (b) vend Figura 7.4: Aproximación en tiempo final r= 100. -1 -0.5 0 0.5 1 x 4.5 5 5.5 6 u (a) uend -1 -0.5 0 0.5 1 x 100 100.5 101 101.5 v (b) vend Figura 7.5: Aproximación en tiempo final r= 1000. 7.2. Modelo de patrones de bacterias quimiotácticas en medio líquido. Este modelo se corresponde con el de Keller-Segel que hemos tratado habitualmente, modificándose la velocidad quimiotáctica para que esta se reduzca al aumentar v, eliminándose el término de autodegradación de vy sustituyéndose el término de producción de vpor otro que aunque sigue siendo creciente respecto de u, está 66 CAPÍTULO 7. OTROS MODELOS CON QUIMIOTAXIS acotado por 1 para valores grandes de u.                      ut−Du∆u+λ∇ · (u (1 + v)2∇v)=0,en Ω×(0, T). vt−Dv∆v−αu2 µ+u2= 0,en Ω×(0, T). ∂u ∂n = 0,∂v ∂n = 0,en ∂Ω×(0, T). u(x, 0) = u0(x), v(x, 0) = v0(x),en Ω. (7.6) Una propiedad fácil de ver de (7.6) es la siguiente: Proposición 7.2.1. El problema (7.6) conserva la cantidad total de u Demostración. Es bastante similar al caso de Keller-Segel, basta con integrar en Ω y aplicar las condiciones de contorno. Empecemos estudiando la ecuación de v. Aplicaremos una aproximación centrada en espacio para el laplaciano. Para conseguir positividad incondicional, tomaremos el término de difusión implícito y el término de producción explícito. Esto es −α(uk−1)2 µ+ (uk−1)2.(7.7) Esta última elección (7.7), incluso desacoplará las ecuaciones de ukde las ecuaciones de vk. La aproximación vk, por tanto, se reduce a resolver el sistema de ecuaciones vk j−vk−1 j dt −Dv∆hvk j−α(uk−1)2 µ+ (uk−1)2= 0,∀j= 1, .., N, (7.8) Por otro lado, para la ecuación de u, en lo relacionado a la discretización en espacio tomaremos el término de autodifusión centrado y el término de quimiotaxis descentrado. Para la aproximación en tiempo tomamos todos los términos implícitos. Por tanto el sistema de ecuaciones que resuelve ukqueda tal que así uk j−uk−1 j dt −Du∆uk j+λchemupw j= 0,∀j= 1, .., N, (7.9) donde el descentrado del término de quimiotaxis viene recogido en chemupw j=−((vx)k j−1/2)+ h uk j−1 (1 + vk j−1)2−((vx)k j−1/2)− h uk j (1 + vk j)2 +((vx)k j+1/2)+ h uk j (1 + vk j)2+((vx)k j+1/2)− h uk j+1 (1 + vk j+1)2. (7.10) El orden y la forma de resolver la ecuaciones será la misma que en la sección anterior. 67 CAPÍTULO 7. OTROS MODELOS CON QUIMIOTAXIS 7.2.1. Simulación y análisis de los resultados En esta sección utilizaremos los parámetros habituales junto con µ= 1 y compararemos las aproximaciones de este nuevo modelo con las aproximación descentrada de Keller-Segel para β= 0. Esta última se presenta en la Figura 7.6 -1 -0.5 0 0.5 1 x 0 1000 2000 3000 4000 5000 u (a) uend -1 -0.5 0 0.5 1 x 100 150 200 250 v (b) vend Figura 7.6: Aproximación de Keller-Segel en tiempo final β= 0. A continuación se muestran las aproximaciones del nuevo modelo -1 -0.5 0 0.5 1 x 0 2 4 6 8 10 12 14 u (a) u6 -1 -0.5 0 0.5 1 x 3.5 4 4.5 5 5.5 6 6.5 7 v (b) v6 Figura 7.7: Aproximación para t6= 0,01. En la Figura 7.7, se observa que la aproximación del udel nuevo modelo no genera picos tan pronunciados como en Keller-Segel o incluso en la logística con autodegradación de u. 68 CAPÍTULO 7. OTROS MODELOS CON QUIMIOTAXIS -1 -0.5 0 0.5 1 x 2 3 4 5 6 7 8 9 u (a) u26 -1 -0.5 0 0.5 1 x 6.8 7 7.2 7.4 v (b) v26 Figura 7.8: Aproximación para t26 = 0,05. Observando la Figura (7.8), se aprecia que la agregación final de use produce en un tiempo muy corto con respecto a Keller-Segel y esta es mucho menos prominente. -1 -0.5 0 0.5 1 x 4.5 5 5.5 6 u (a) uend -1 -0.5 0 0.5 1 x 28 28.5 29 29.5 v (b) vend Figura 7.9: Aproximación en tiempo final. La aproximación en tiempo final del nuevo modelo que observamos en la Figura 7.9 difiere bastante de la de Keller-Segel (ver Figura 7.6). Estas nuevas gráficas dan lugar a pensar que están dominando los términos de autodifusión y la solución de uparece tender a un valor fijo. Este será el valor medio de su condición inicial ˜u=1 2R1 −1u0(x)dx = 5,3, ya que por la Proposición 7.2.1 la cantidad total se conserva. Esto es diferente en el caso de vya que no conserva la cantidad total. Veamos por último las cantidades totales numéricamente 69 CAPÍTULO 7. OTROS MODELOS CON QUIMIOTAXIS 0 0.1 0.2 0.3 0.4 0.5 t 10 10.5 11 11.5 Cu (a) Cantidad total de u 0 0.1 0.2 0.3 0.4 0.5 t 10 20 30 40 50 60 Cv (b) Cantidad total de v Figura 7.10: Cantidad total. En la Figura (7.10), se corrobora que el esquema discretizado mantiene la conservación de la cantidad total de u. 7.3. Modelo de patrones de bacterias quimiotácticas en medio semi-sólido. Este es un nuevo modelo distinto al de Keller-Segel. Lo más notorio es que se introduce una variable adicional w, cuya interpretación física se corresponde con la concentración de nutrientes en el medio. La ecuación que gobierna la dinámica de esta nueva variable posee únicamente un término de autodifusión y otro de consumo por la variable uregulado al mismo tiempo por w. Por otro lado, con respecto al las ecuaciones de uyvdel esquema de la Sección 7.2 , se incorpora un término de reacción logístico a la variable u, cuya capacidad de carga dependerá de la variable w. Asimismo, se introduce para la ecuación de v un término de consumo por u, con dependencia lineal de vy se añade al término de 70 CAPÍTULO 7. OTROS MODELOS CON QUIMIOTAXIS producción una dependencia lineal de w.                                ut−Du∆u+λ∇ · (u (1 + v)2∇v)−ru δw2 1 + w2−u= 0,en Ω×(0, T). vt−Dv∆v−αwu2 µ+u2+βuv = 0,en Ω×(0, T). wt−Dw∆w+κu w2 1 + w2= 0,en Ω×(0, T). ∂u ∂n = 0,∂v ∂n = 0,∂w ∂n = 0,en ∂Ω×(0, T). u(x, 0) = u0(x), v(x, 0) = v0(x), w(x, 0) = w0(x),en Ω. (7.11) Comenzaremos con la ecuación de w. Como sabemos, tomamos la aproximación del laplaciano centrada en espacio e implícita en tiempo. El término de degradación, como hemos visto, es conveniente que tenga carácter semi-implícito en wpara conseguir positividad incondicional. Con esta idea, lo más simple es tomarlo de una forma que linealice la ecuación y la desacople. Por tanto la elección que tomamos dará lugar al siguiente sistema wk j−wk−1 j dt −Dw∆hwk j+κuk−1wk−1wk 1+(wk−1)2= 0,∀j= 1, .., N. (7.12) Con todas las ideas propuestas, y habiendo sido capaces de desacoplar el sistema que despeja wk, resulta directo para vkrealizar la siguiente elección vk j−vk−1 j dt −Dv∆hvk j−αwk(uk−1)2 µ+ (uk−1)2+βuk−1vk= 0,∀j= 1, .., N, (7.13) Finalmente, la ecuación de ues una combinación de las que se presentaron en las Secciones 7.1 y 7.2, donde ya wkyvkhan podido ser calculada independientemente. Por tanto para uktendremos uk j−uk−1 j dt −Du∆huk j+λchemupw j−ruk−1 j(δ(wk)2 1+(wk)2−uk j) = 0,∀j= 1, .., N, (7.14) Para obtener las aproximaciones, dados uk−1, vk−1ywk−1, obtendremos wkde (7.12), la cual introduciremos en (7.13) para calcular vky finalmente, sustituyendo wkyvken (7.14) obtendremos uk. 7.3.1. Simulación y análisis de los resultados Este modelo, al introducir una nueva variable es difícilmente comparable con los anteriores, por ello se estudia independientemente. 71 CAPÍTULO 7. OTROS MODELOS CON QUIMIOTAXIS Como dominio espacio-temporal tomamos el habitual Ω×T= (−1,1)×(0,0,5), y los pasos de espacio h= 0,001 y tiempo dt = 0,002. Por otro lado, tomamos los siguientes datos para los parámetros: Du=Dv=Dw=µ= 1,λ=α= 50,β= 10, κ= 20,r= 3 yδ= 5. -1 -0.5 0 0.5 1 x 0 2 4 6 8 10 u (a) u0 -1 -0.5 0 0.5 1 x 0 2 4 6 8 10 v (b) v0 -1 -0.5 0 0.5 1 x 0 2 4 6 8 10 w (c) w0 Figura 7.11: Condiciones iniciales. Las condiciones iniciales se corresponden con las gráficas de la Figura 7.11. estas son: u0(x) = 10 (x+ 1)2(x−1)2,v0(x) = 10 sen2(2πx)yw0(x) = 10 sen2(5π/2x). 72 CAPÍTULO 7. OTROS MODELOS CON QUIMIOTAXIS -1 -0.5 0 0.5 1 x 2 3 4 5 6 7 u (a) u26 -1 -0.5 0 0.5 1 x 2.2 2.4 2.6 2.8 3 3.2 v (b) v26 -1 -0.5 0 0.5 1 x 1 1.2 1.4 1.6 1.8 w (c) w26 Figura 7.12: Aproximación para t26 = 0,05. Estudiemos lo que ocurre en esta configuración. En primer lugar, las aproximaciones parecen acomodarse rápidamente a las funciones que se observan en la Figura (7.12) y con esta forma permanecen durante un tiempo. Pero la forma del sistema (7.11) sugiere que esta situación debe evolucionar a otra distinta, ya que vywson consumidas por uy las tres funciones poseen grandes valores en zonas cercanas en el espacio. 79 APÉNDICE B. CÓDIGOS DE SIMULACIÓN 75 76 for k=1: N1 77 78 Mq(k,k)=Mq(k,k)-Vx(k)/2; 79 Mq (k,k +1)=Mq(k ,k+1) - Vx(k)/2 ; 80 Mq(k+1,k)=Mq(k+1,k)+ Vx(k)/2 ; 81 Mq(k+1,k +1)=Mq(k+1,k+1)+Vx(k)/2; 82 83 end 84 85 % Resolucion de u 86 Au =Du*R +1/ dt *M- lambda * Mq; 87 Bu=M *(1/ dt*U0); 88 U=Au\Bu; 89 90 % Actualizacion 91 U0=U; 92 V0=V; 93 94 % Guardado de las soluciones 95 Mu (t+1 ,:)=U0; 96 Mv (t+1 ,:)=V0; 97 98 %% Graficas 99 subplot (1 ,2 ,1) 100 plot (X,Mu(t+1 ,:)) 101 xlabel('x'); 102 ylabel('u'); 103 title ('Densidad de celulas '); 104 shg 105 106 subplot (1 ,2 ,2) 107 plot (X,Mv(t+1 ,:)) 108 xlabel('x'); 109 ylabel('v'); 110 title ('Densidad de quimico '); 111 shg 112 113 end Cantidad total 1%% Cantidad total 2 80 APÉNDICE B. CÓDIGOS DE SIMULACIÓN 3Cu =zeros (1, Nt); 4Cv =zeros (1, Nt); 5 6for t=1: Nt 7 8Cu (t)= trapz (X,Mu(t ,:) ); 9Cv (t)= trapz (X,Mv(t ,:) ); 10 11 end 12 13 %% Graficas 14 15 figure 16 subplot (1 ,2 ,1) 17 plot (T,Cu ,T ,Cu (1) *ones (1 ,Nt)) 18 title ('Cantidad total de sustancia celulas ') 19 axis ([0 ,T(end ),Cu (1) -1,Cu (1) +1]) 20 xlabel('t') 21 ylabel('C_u') 22 legend('Aproximaci’on','Te’orica ','Location ','northwest') 23 shg 24 25 % V teorica 26 Vteo = @(x) Cv (1) *exp (- beta *x)+ alpha / beta *Cu (1) *(1 - exp (- beta *x)); 27 28 subplot (1 ,2 ,2) 29 plot (T,Cv ,T, Vteo (T)) 30 title ('Cantidad total de sustancia quimica ') 31 xlabel('t') 32 ylabel('C_v') 33 legend('Aproximacion','Teorica','Location ','northwest') 34 shg Energía 1%% Energia 2 3% E =\ integral (F(u) +1/2(( Vx )^2+ v ^2) -u*v), F(u)= ulogu -u-v 4F=@(u) alpha *Du*u.*( log(u) -1); 5E=zeros (1 ,Nt); 6D= zeros (1, Nt1 ); 7 8for t=1: Nt 81 APÉNDICE B. CÓDIGOS DE SIMULACIÓN 9 10 Vx =( Mv(t ,3: N)-Mv(t ,1:N -2)) /(2* h); 11 Vx =[0,Vx ,0]; %Condiciones de contorno 12 Fu =F(Mu(t ,:) ); 13 Fu (Mu (t ,:) ==0) =0; %Fu (0) =0 14 15 E(t)= trapz (X,Fu + lambda * Dv *1/2* Vx .^2+ lambda * beta *1/2* Mv (t ,:) .^2 - lambda * alpha *Mu(t ,:) .* Mv(t ,:) ); 16 17 %% Disipacion 18 19 if t ==1 20 Vt (1 ,:)= zeros (N ,1) ; 21 else 22 Vt =( Mv(t ,:) -Mv(t -1 ,:) )/dt; 23 end 24 25 logx =( log (Mu (t ,3: N))-log (Mu(t ,1:N -2))) /(2* h); 26 logx =[0 ,logx ,0]; %Condiciones de contorno 27 28 29 D(t)= trapz (X, alpha * Mu(t ,:) .*( Du *logx - lambda *Vx) .^2 + lambda *Vt .^2) ; 30 end 31 32 %% Ley de Energia 33 34 LE =(E(2: Nt)-E(1: Nt1))/dt+D(1: Nt1 ); 35 36 37 %% Graficas 38 39 figure 40 subplot (1 ,2 ,1) 41 plot (T,E) 42 title ('Energia') 43 xlabel('t') 44 ylabel('E') 45 shg 46 47 subplot (1 ,2 ,2) 48 plot (T(2: Nt),LE (2: Nt)) 49 title ('Ley de energia ') 50 xlabel('t') 51 ylabel('LE ') 52 shg 82 APÉNDICE B. CÓDIGOS DE SIMULACIÓN Diferencias finitas descentradas 1clear all 2clf; 3 4%% Parametros 5 6Du =1; Dv =1; lambda =10; beta =1; alpha =50; 7 8%% Datos iniciales 9 10 h =0.001; L =1; X=( -L:h:L); N= length (X); N1=N -1; 11 dt =0.002; T =(0: dt :0.5) ; Nt = length (T); Nt1 =Nt -1; 12 13 %% Condicion inicial 14 15 factorv0 =10; factoru0 =10; 16 U0 = factoru0 *((( X +1) .^2.*( X -1) .^2) ) '; 17 V0 = factorv0 *(( sin (4* pi /2* X)) .^2) '; 18 19 %% Graficas ( inicial ) 20 21 subplot (1 ,2 ,1) 22 plot (X,U0) 23 title ('Densidad de celulas '); 24 ylabel('u') 25 xlabel('x') 26 shg 27 28 subplot (1 ,2 ,2) 29 plot (X,V0) 30 title ('Densidad de quimico '); 31 ylabel('v') 32 xlabel('x') 33 % hold on 34 shg 35 36 %% Matrices para guardar las soluciones 37 38 Mu= zeros (Nt ,N); 39 Mv= zeros (Nt ,N); 40 Mu (1 ,:)=U0; 41 Mv (1 ,:)=V0; 42 43 %% Matrices de coeficientes 44 83 APÉNDICE B. CÓDIGOS DE SIMULACIÓN 45 % Identidad 46 I= spdiags ( ones (N ,1) ,0,N ,N); 47 48 % Difusion 49 v1 =( -1/h ^2) *ones (N ,1) ; 50 v2 =(2/ h^2) *ones (N ,1) ; 51 Adif = spdiags ([ v1 v2 v1],[ -1 0 1],N,N); 52 Adif (1 ,2)=- Adif (1 ,1); 53 Adif (N,N -1)=- Adif (N,N); 54 55 % Matriz de v 56 Av =(1/ dt+ beta )*I+Dv* Adif; 57 58 % Matris de u sin quimiotaxis 59 Au1 =1/ dt*I+ Du* Adif ; 60 61 cero = zeros (N ,1) ; 62 for t=1: Nt1 63 %% Resolucion de v 64 65 Bv =alpha *U0 +1/ dt* V0; 66 V=Av\Bv; 67 68 %% Quimiotaxis 69 Vx (2: N)=(V (2: N)-V (1: N1))/h; 70 71 Vxpos =Vx >0; 72 Vxneg =Vx <0; 73 74 Auchem = spdiags ([ cero cero cero ] ,[-1 0 1],N,N); 75 for k=2: N1 76 Auchem (k,k -1) = -Vxpos (k)*Vx(k)/h; 77 Auchem (k,k)= Vxpos (k +1) *Vx(k+1) /h - Vxneg (k)*Vx (k)/h; 78 Auchem (k,k +1) = Vxneg (k+1) *Vx (k +1) /h; 79 end 80 Auchem (1 ,1) =2* Vxpos (2) *Vx (2) /h; 81 Auchem (1 ,2) =2* Vxneg (2) *Vx (2) /h; 82 Auchem (N,N -1) = -2* Vxpos (N)*Vx(N)/h; 83 Auchem (N,N) = -2* Vxneg (N)* Vx(N)/h; 84 85 Au =Au1 + lambda * Auchem ; 86 Bu =1/ dt*U0 ; 87 88 %% Resolucion de u 89 U=Au\Bu; 90 %% Graficas 84 APÉNDICE B. CÓDIGOS DE SIMULACIÓN 91 92 subplot (1 ,2 ,1) 93 plot (X,U) 94 title ('Densidad de c’elulas '); 95 ylabel('u') 96 xlabel('x') 97 shg 98 99 subplot (1 ,2 ,2) 100 plot (X,V) 101 title ('Densidad de qu’imico '); 102 ylabel('v') 103 xlabel('x') 104 shg 105 106 107 %% Actualizacion 108 U0=U; 109 V0=V; 110 111 %% Guardado de las soluciones 112 Mu (t+1 ,:)=U0; 113 Mv (t+1 ,:)=V0; 114 end Keller-Segel con reacción logística 1clc; clear ; clf ; 2 3%% parametros 4Du =1; Dv =1; lambda =10; beta =1; alpha =50; r=1; delta =5; 5 6%% Datos iniciales 7h =0.001; L =1; X=( -L:h:L); N= length (X); N1=N -1; 8dt =0.002; T =(0: dt :5) ; Nt= length (T); Nt1 =Nt -1; 9 10 %% Condicion inicial 11 factoru0 =10; factorv0 =10; 12 U0 = factoru0 *((( X +1) .^2.*( X -1) .^2) ) '; 13 V0 = factorv0 *(( sin (4* pi /2* X)) .^2) '; 14 15 %% Graficas 16 subplot (1 ,2 ,1) 17 plot (X,U0) 85 APÉNDICE B. CÓDIGOS DE SIMULACIÓN 18 title ('Densidad de c’elulas '); 19 ylabel('u') 20 xlabel('x') 21 shg 22 23 subplot (1 ,2 ,2) 24 plot (X,V0) 25 title ('Densidad de qu’imico '); 26 ylabel('v') 27 xlabel('x') 28 29 shg 30 31 %% Matrices de rigiez y masa 32 cero = zeros (N ,1) ; 33 K=spdiags ([ cero cero cero ] ,[-1 0 1],N,N); 34 M= spdiags ( cero , 0,N ,N); 35 36 for k=1: N1 37 K(k,k)=K(k,k) +1/h; 38 K(k+1,k+1)=K(k+1,k+1) +1/ h; 39 K(k,k+1)= K(k,k+1) -1/h; 40 K(k+1,k)=K(k,k+1) ; 41 M(k,k)=M(k,k)+h/2; 42 M(k+1,k+1)=M(k+1,k+1)+h/2; 43 end 44 Av=Dv*K+(1/ dt+ beta)*M; 45 46 for t=1: Nt1 47 48 %% Resolucion de v 49 Bv =M*( alpha *U0 +1/ dt*V0); 50 51 V=Av\Bv; 52 53 %% Quimiotaxis 54 Vx (2: N)=(V (2: N)-V (1: N1))/h; 55 56 Vxpos =Vx >0; 57 Vxneg =Vx <0; 58 59 Auchem = spdiags ([ cero cero cero ] ,[-1 0 1],N,N); 60 for k=2: N1 61 Auchem (k,k -1) = -Vxpos (k)*Vx(k); 62 Auchem (k,k)= Vxpos (k +1) *Vx(k+1) - Vxneg (k)* Vx(k); 63 Auchem (k,k +1) = Vxneg (k+1) *Vx (k +1) ; 86 APÉNDICE B. CÓDIGOS DE SIMULACIÓN 64 end 65 Auchem (1 ,1)= Vxpos (2)*Vx (2) ; 66 Auchem (1 ,2)= Vxneg (2)*Vx (2) ; 67 Auchem (N,N -1) =- Vxpos (N)* Vx(N); 68 Auchem (N,N)=- Vxneg (N)*Vx(N); 69 70 %% Resolucion de u 71 Bu =(1/ dt+r* delta )*M*U0; %/r*delta -> produccion explicita ( logistica ) 72 Umod = spdiags (U0 ,0,N ,N); % degradacion semiimplicita ( logistica ) 73 Au =M *1/ dt +r*M* Umod +Du *K+ lambda * Auchem ; 74 U=Au\Bu; 75 76 %% Graficas 77 subplot (1 ,2 ,1) 78 plot (X,U) 79 title ('Densidad de c’elulas '); 80 ylabel('u') 81 xlabel('x') 82 shg 83 84 subplot (1 ,2 ,2) 85 plot (X,V) 86 title ('Densidad de qu’imico '); 87 ylabel('v') 88 xlabel('x') 89 shg 90 91 %% Actualizacion 92 U0=U; 93 V0=V; Modelo de patrones de bacterias quimiotácticas en medio líquido. 1clc; clear ; clf ; 2 3%% parametros 4Du =1; Dv =1; lambda =10; mu =1; alpha =50; 5 6%% Datos iniciales 7h =0.001; L =1; X=( -L:h:L); N= length (X); N1=N -1; 87 APÉNDICE B. CÓDIGOS DE SIMULACIÓN 8dt =0.002; T =(0: dt :0.3) ; Nt = length (T); Nt1 =Nt -1; 9 10 %% Condicion inicial 11 factoru0 =10; factorv0 =10; 12 U0 = factoru0 *((( X +1) .^2.*( X -1) .^2) ) '; 13 V0 = factorv0 *(( sin (4* pi /2* X)) .^2) '; 14 15 %% Graficas 16 subplot (1 ,2 ,1) 17 plot (X,U0) 18 title ('Densidad de c’elulas '); 19 ylabel('u') 20 xlabel('x') 21 shg 22 23 subplot (1 ,2 ,2) 24 plot (X,V0) 25 title ('Densidad de qu’imico '); 26 ylabel('v') 27 xlabel('x') 28 shg 29 30 31 %% Matrices de rigiez y masa 32 cero = zeros (N ,1) ; 33 K=spdiags ([ cero cero cero ] ,[-1 0 1],N,N); 34 M= spdiags ( cero , 0,N ,N); 35 36 for k=1: N1 37 K(k,k)=K(k,k) +1/h; 38 K(k+1,k+1)=K(k+1,k+1) +1/ h; 39 K(k,k+1)= K(k,k+1) -1/h; 40 K(k+1,k)=K(k,k+1) ; 41 M(k,k)=M(k,k)+h/2; 42 M(k+1,k+1)=M(k+1,k+1)+h/2; 43 end 44 Av =Dv *K+1/ dt*M; 45 46 for t=1: Nt1 47 48 %% Resolucion de v 49 Bv=M*( alpha *U0 .^2./( mu+U0 .^2) +1/ dt*V0); % alpha *... -> produccion explicita 50 51 V=Av\Bv; 52 88 APÉNDICE B. CÓDIGOS DE SIMULACIÓN 53 54 %% Quimiotaxis 55 Vx (2: N)=(V (2: N)-V (1: N1))/h; 56 57 Vxpos =Vx >0; 58 Vxneg =Vx <0; 59 60 Auchem = spdiags ([ cero cero cero ] ,[-1 0 1],N,N); 61 for k=2: N1 62 Auchem (k,k -1) =- Vxpos (k)* Vx(k); 63 Auchem (k,k)= Vxpos (k +1) *Vx(k+1) - Vxneg (k)* Vx(k); 64 Auchem (k,k +1) = Vxneg (k +1) *Vx(k +1) ; 65 end 66 Auchem (1 ,1)= Vxpos (2)*Vx (2) ; 67 Auchem (1 ,2)= Vxneg (2)*Vx (2) ; 68 Auchem (N,N -1) =- Vxpos (N)* Vx(N); 69 Auchem (N,N)=- Vxneg (N)*Vx(N); 70 71 %%Resolucion de u 72 Vmod = spdiags (1./(1+ V) .^2 ,0 ,N,N); %modificacion de quimiotaxis 73 Bu =1/ dt*M* U0; 74 Au =1/ dt *M+Du *K+ lambda * Auchem * Vmod ; 75 U=Au\Bu; 76 77 %% Graficas 78 subplot (1 ,2 ,1) 79 plot (X,U) 80 title ('Densidad de c’elulas '); 81 ylabel('u') 82 xlabel('x') 83 shg 84 85 subplot (1 ,2 ,2) 86 plot (X,V) 87 title ('Densidad de qu’imico '); 88 ylabel('v') 89 xlabel('x') 90 shg 91 92 93 %% Actualizacion 94 U0=U; 95 V0=V; 96 end