Análisis teórico y numérico de problemas diferenciales con quimiotaxis repulsiva.
Full text
Programa de Doctorado en Matem´aticas Tesis Doctoral An´alisis te´orico y num´erico de problemas diferenciales con quimiotaxis repulsiva Autor Diego Armando Rueda G´omez Directores Francisco Guill´en Gonz´alez Mar´ıa ´ Angeles Rodr´ıguez Bellido Dpto. Ecuaciones Diferenciales y An´alisis Num´erico Facultad de Matem´aticas Universidad de Sevilla Julio de 2018
. ii
A mis padres. A mi esposa y mi hijo. iii
. iv
Agradecimientos En primer lugar, agradezco a Dios por conducir mi vida y por permitirme alcanzar este gran logro profesional. Agradezco a mis padres: Jes´us Rueda y Sandra G´omez, por su apoyo incondicional en cada una de las experiencias vividas hasta hoy. Tambi´en quiero agradecer a mi esposa, Mayra Alejandra, y a mi hijo, Diego Alejandro, por ser mi mayor soporte durante estos a˜nos de estudio y la principal motivaci´on en todo lo que hago. Agradezco profundamente a mis directores, el profesor Francisco Guill´en Gonz´alez y la profesora Mar´ıa ´ Angeles Rodr´ıguez Bellido, no solo por su empe˜no y dedicaci´on en la direcci´on de este trabajo, sino por su gran calidad humana y simpat´ıa, y por la hospitalidad que han tenido conmigo durante mi estad´ıa aqu´ı en Sevilla. Quiero agradecer tambi´en al profesor ´ Elder Jes´us Villamizar, a quien tuve el privilegio de tener como maestro durante mi M´aster, y ahora es colega y amigo en la Universidad Industrial de Santander. Gracias por todo su apoyo y ayuda desinteresada, y por todos los consejos que me ha brindado cuando lo he necesitado. Agradezco a todos los profesores del Departamento EDAN por acogerme durante estos a˜nos. A todos, muchas gracias. v
. vi
Resumen Esta tesis se enmarca en el ´ambito del an´alisis te´orico y num´erico de Ecuaciones en Derivadas Parciales, con aplicaciones a otras ciencias. Concretamente, aborda el estudio de algunos problemas diferenciales de quimiotaxis de tipo repulsiva-productiva. Los primeros tres cap´ıtulos est´an dedicados al estudio de un modelo de quimiotaxis repulsiva con t´ermino de producci´on cuadr´atico, y los restantes dos cap´ıtulos se enfocan en modelos con t´erminos de producci´on lineal y potencial (con potencia superlineal y subcuadr´atica). En el Cap´ıtulo 1, se presentan dos esquemas num´ericos discretos solamente en tiempo, energ´eticamente estables, para un modelo de quimiotaxis repulsiva con t´ermino de producci´on cuadr´atico, y se estudian algunas propiedades adicionales de estos esquemas tales como la conservaci´on de la cantidad total de masa, positividad, resolubilidad, convergencia hacia soluciones d´ebiles y estimaciones de error. En el Cap´ıtulo 2, se estudia un esquema num´erico completamente discreto con elementos finitos, energ´eticamente estable, asociado al modelo estudiado en el Cap´ıtulo 1, basado en la introducci´on de una variable auxiliar. Nuevamente, se estudian algunas propiedades como resolubilidad, conservaci´on de masa, convergencia hacia soluciones d´ebiles, estimaciones de error, y estimaciones d´ebiles y fuertes del esquema. Adicionalmente, como el esquema bajo estudio es no lineal, se proponen dos m´etodos iterativos para aproximar las soluciones y se prueba la resolubilidad y la convergencia de ambos esquemas hacia el esquema no lineal. En el Cap´ıtulo 3, se estudia el comportamiento asint´otico de las soluciones del modelo estudiado en los Cap´ıtulos 1and 2. En la primera parte, se analiza el comportamiento en tiempo infinito de soluciones d´ebiles del problema continuo y se prueba convergencia exponencial hacia un estado constante. En la segunda parte, se estudia este mismo comportamiento para dos esquemas num´ericos completamente discretos asociados a este modelo. Finalmente, los Cap´ıtulos 4y5se centran en el estudio de modelos de quimiotaxis repulsiva con t´erminos de producci´on lineal y potencial, respectivamente. Aqu´ı, usando una t´ecnica de regularizaci´on, se proponen algunos esquemas num´ericos completamente discretos con elementos finitos, energ´eticamente estables, asociados a estos modelos, y se prueban algunas propiedades adicionales tales como la resolubilidad, conservaci´on de la cantidad total de masa, y positividad aproximada de las soluciones. vii
Abstract This PhD thesis falls within the scopes of Theoretical and Numerical analysis of Partial Differential Equations, with applications to other sciences. Specifically, it addresses the study of some differential problems of repulsive-productive chemotaxis. The first three chapters are devoted to study a chemo-repulsion model with quadratic production, and other two chapters are focused on models with linear and potential (with a superlinear and subquadratic power) production. In Chapter 1, we present two unconditionally mass-conservative and energy-stable timediscrete numerical schemes for a chemo-repulsion model with quadratic production, and study some additional properties of the schemes such as positivity, solvability, convergence towards weak solutions and error estimates of these schemes. In Chapter 2, we study an unconditionally mass-conservative and energy-stable fully discrete FE scheme associated to the problem studied in Chapter 1, in which an auxiliary variable is introduced. Again, we study some properties like solvability, convergence towards weak solutions, error estimates, and weak, strong and more regular a priori estimates of the scheme. Additionally, as the scheme is nonlinear, we propose two different linear iterative methods to approach the solutions and we prove solvability and convergence of both methods to the nonlinear scheme. In Chapter 3, we focus on the study of the asymptotic behaviour of the solutions of the model studied in Chapters 1and 2. In the first part, we analyze the large-time behavior of the global weak-strong solutions and we prove the exponential convergence to a constant state as time goes to infinity; and in the second part, we study this same behaviour for two fully discrete FE numerical schemes associated to this model. Finally, in Chapters 4and 5we focus on the study of chemo-repulsion models with linear and potential (superlinear and subquadratic) production, respectively. Here, by using a regularization technique, we propose some unconditionally energy-stable and massconservative fully discrete FE schemes associated to these models, and we prove some additional properties such as solvability and approximated positivity of the solutions. viii
. ix
.
Introduction Chemotaxis is a biological phenomenon which describes the oriented movement of living organisms in response to a chemical stimulus. The chemotaxis is called attractive when the organisms move towards regions with higher chemical concentration, while if the motion is towards lower concentrations, the chemotaxis is called repulsive. At the same time, the presence of living organisms can produce or consume chemical substance. This process allows the bacteria to find food, moving towards the highest concentration of food molecules, or move away from poisonous substances. The typical example for chemotaxis is the amoebae Dictyostelium, which is a species of soil-living amoeba belonging to the phylum Mycetozoa. When they are moving towards the nutrients, they produce a chemical substance, cyclic Adenosine Monophosphate, attracting other amoebae. The main purpose of this thesis is the theoretical and numerical study of repulsiveproductive chemotaxis models, given by the following parabolic PDE’s system: ∂tu−∆u=∇·(u∇v) in Ω, t > 0, ∂tv−∆v+v=f(u) in Ω, t > 0,(1) where Ω ⊂Rn,n= 2,3, is a bounded domain with boundary ∂Ω, and u(x, t)≥0 and v(x, t)≥0 denote the cell density and the chemical concentration, respectively. Moreover, f(u)≥0 is the production term. Model (1) possesses some properties among which we focus on: (i) The blow-up phenomenon is not expected to take place here. (ii) Mass-conservation: This problem is conservative in u, because the total mass RΩu(·, t) remains constant in time, that is, d dt ZΩ u(·, t)= 0,i.e. ZΩ u(·, t) = ZΩ u0:= m0,∀t > 0.(2) (iii) In some cases (for instance, when the production is given by a power of u, that is, 1
f(u) = up), this problem satisfies an energy inequality in the form: E(u(t1), v(t1)) −E(u(t0), v(t0)) + Zt1 t0 D(u(s), v(s)) ds ≤0,(3) for a.e. t0, t1:t1≥t0≥0, where E(u, v) and D(u, v) denote, respectively, the free energy and the physical dissipation terms of the problem. In particular, (3) implies that this problem has a decreasing in time energy. The motivation of this work was initially based on the study of Finite Element (FE) numerical approximations of problem (1), with production term f(u) = u, conserving at the discrete level the main properties of the continuous model, such as mass-conservation, energy-stability and positivity of the variables. Moreover, as part of the numerical analysis to be developed, we set out the study of the well-posedness, convergence and error estimates of the schemes, among others. However, we find that this is not an easy task, and the main difficulty lies in how to approximate at the discrete level the energy inequality (3) for this case, taking into account that, in the continuous problem, this energy inequality is obtained by testing the u-equation by the nonlinear function F(u) = ln u. Once these difficulties have been pointed out, we decided to focus first in the study of the problem (1) with quadratic production term, that is, f(u) = u2, mainly for the following two reasons: (a) the energy inequality in the continuous problem is obtained by testing the u-equation by uand the v-equation by −∆v, which it is not so difficult to reproduce for fully discrete FE approximations; and (b) the quadratic production term allows to control an energy in L2(Ω)-norm for u, which is very useful for performing numerical analysis. In this part, we studied some mass-conservative and energy-stable (in the sense that a discrete energy decreases in time) schemes associated to this problem, for which we also analyzed well-posedness, positivity, convergence towards weak solutions, error estimates and convergence at infinite time. Later, with the experience obtained in the case of quadratic production, we focused in the case of linear production, that is, f(u) = u. In this case, by using a regularization technique used in previous works, we construct some unconditionally energy-stable fully discrete FE approximations, for which we proved well-posedness and some additional properties such as mass-conservation and approximated positivity of the variables (in the sense that the negative part of the variable tends to 0 as the regularization parameter tends to 0). However, since we could not obtain uniform estimates independent of the discrete and regularization parameters that allowed us to take limits on the discrete problem, the convergence towards weak solutions was not proved. 2
Finally, in the last part of this work, we have studied the intermediate case, in which f(u) = upwith p∈(1,2). Here, we have adapted the ideas used in the case of linear production, in order to obtain mass-conservative and energy-stable fully discrete FE schemes. This PhD thesis is organized in five chapters, which we expect to correspond to five different papers: Chapter 1focuses on the study of numerical approximations of model (1) in the case of f(u) = u2. We present two unconditionally mass-conservative and energy-stable first order time schemes: the (nonlinear) Backward Euler scheme and a linearized coupled version. We analyze positivity, solvability, convergence towards weak solutions and error estimates of these schemes. In particular, uniqueness of the nonlinear scheme is proved assuming small time step with respect to a strong norm of the scheme. This hypothesis is simplified in 2D domains where a global in time strong estimate is proved. Finally, some numerical simulations are made in order to compare the behavior of the schemes. Chapter 2is devoted to study a fully discrete FE scheme associated to problem (1) in the case of f(u) = u2. By following the ideas presented in Chapter 1, we introduce σ=∇vas an auxiliary variable, and then the corresponding FE backward Euler scheme is unconditionally mass-conservative and energy-stable. For this nonlinear scheme, we study some properties like solvability, convergence towards weak solutions, error estimates, and weak, strong and more regular a priori estimates of the scheme. Additionally, we propose two different linear iterative methods to approach the nonlinear scheme: an energy-stable Picard’s method and the Newton’s method. We prove solvability and convergence of both methods to the nonlinear scheme. Finally, we provide some numerical results in agreement with our theoretical analysis about the error estimates. Chapter 3is focused on the study of the asymptotic behaviour of the problem (1) in the case of f(u) = u2. In the first part, we analyze the large-time behavior of the global weak-strong solutions and we prove the exponential convergence to a constant state as time goes to infinity. In the second part, we study this same behaviour for two fully discrete numerical schemes associated to this model: the FE backward Euler associated to (1) (with f(u) = u2) and the nonlinear scheme defined in Chapter 2. On the way, in order to analyze the asymptotic behaviour for the backward Euler scheme, we prove its solvability and unconditional energy-stability. Finally, we compare the numerical schemes throughout several numerical simulations. Chapter 4is devoted to study unconditionally energy-stable and mass-conservative numerical schemes for problem (1) in the case of f(u) = u. By using a regularization 3
technique (in which, some regularized functions approximating F(s) = u(ln u−1) and its first and second derivatives are introduced), we propose three fully discrete FE approximations. The first one is a nonlinear approximation in the variables (u, v); the second one is another nonlinear approximation obtained introducing σ=∇vas an auxiliary variable; and the third one is a linear approximation constructed by mixing the regularization procedure with the so called Energy Quadratization strategy, in which the energy of the system is transformed into a quadratic form by introducing new auxiliary variables. In addition, we prove the well-posedness of the numerical schemes. In fact, unconditional existence of solution, but conditional uniqueness (for the nonlinear schemes) are proved. Finally, we compare the behavior of the schemes throughout several numerical simulations. At last, in Chapter 5we focus on the study of problem (1) in the case of f(u) = up, with p∈(1,2). In the first part, by using a regularization technique, we prove the existence of solutions of the model. In the second part, we propose three fully discrete FE nonlinear approximations, where the first one is defined in the variables (u, v), and the second and third ones by introducing σ=∇vas an auxiliary variable. We prove some unconditional properties such as mass-conservation, energy-stability and solvability of the schemes. Finally, we compare the behavior of the schemes throughout several numerical simulations. 4
Chapter 1 On a chemo-repulsion model with quadratic production: The continuous problem and time-discrete numerical schemes 1.1 Introduction Chemotaxis is understood as the biological process of the movement of living organisms in response to a chemical stimulus which can be given towards a higher (attractive) or lower (repulsive) concentration of a chemical substance. At the same time, the presence of living organisms can produce or consume chemical substance. A repulsive-productive chemotaxis model can be given by the following parabolic PDE’s system: ∂tu−∆u=∇·(u∇v) in Ω, t > 0, ∂tv−∆v+v=f(u) in Ω, t > 0,(1.1) where Ω ⊂Rn,n= 2,3, is an open bounded domain with boundary ∂Ω. The unknowns for this model are u(x, t)≥0, the cell density, and v(x, t)≥0, the chemical concentration. Moreover, f(u) is a function which is nonnegative when u≥0. In this paper, we consider the particular case in which the production term is quadratic, that is f(u) = u2, and then we focus on the following initial-boundary problem: ∂tu−∆u=∇·(u∇v) in Ω, t > 0, ∂tv−∆v+v=u2in Ω, t > 0, ∂u ∂n=∂v ∂n= 0 on ∂Ω, t > 0, u(x,0) = u0(x)≥0, v(x,0) = v0(x)≥0 in Ω. (1.2) 5
The quadratic production term allows to control an energy in L2(Ω)-norm for u(see (1.20)- (1.21)), which is very useful for performing numerical analysis. Other production terms will be studied in Chapters 4and 5. In the case of linear production term, in [2] the authors proved that model (1.1) with f(u) = uis well-posed in the following sense: there exist global in time weak solutions (based on an energy inequality) and, for 2Ddomains, there exists a unique global in time strong solution. However, as far as we know, there are not works studying problem (1.2) with quadratic production. In addition, some papers on numerical analysis for chemotaxis models with linear production are [3,6,10,12,15]. In order to develop our analysis, we reformulate (1.2) introducing the new variable σ= ∇v. Then, we rewrite the model (1.2) as follows: ∂tu−∇·(∇u) = ∇·(uσ) in Ω, t > 0, ∂tσ−∇(∇·σ) + rot(rot σ) + σ=∇(u2) in Ω, t > 0, ∂u ∂n= 0 on ∂Ω, t > 0, σ·n= 0,[rot σ×n]tang = 0 on ∂Ω, t > 0, u(x,0) = u0(x)≥0,σ(x,0) = ∇v0(x) in Ω, (1.3) where (1.3)2was obtained applying the gradient to equation (1.2)2and adding the term rot(rot σ) using the fact that rot σ= rot(∇v) = 0. Once solved (1.3), we can recover v from u2solving ∂tv−∆v+v=u2in Ω, t > 0, ∂v ∂n= 0 on ∂Ω, t > 0, v(x,0) = v0(x)≥0 in Ω. (1.4) We will use the variable σin order to simplify the notation throughout the chapter. Moreover, for fully discrete schemes by using the Finite Elements Method (which will be analyzed in Chapter 2), it will be very convenient to use the variables (u, σ) in order to obtain an unconditionally energy-stable scheme. This chapter is organized as follows: In Section 1.2, we give the notation and some preliminary results that will be used along this paper. In Section 1.3, we analyze the continuous problem (1.2), obtaining global in time weak regularity for both two and three dimensions, and global in time strong regularity of the model assuming the regularity criteria (1.32), which is satisfied in 2D domains. In Section 1.4, we analyze the Backward Euler scheme corresponding to problem (1.3)-(1.4), including mass-conservation, unconditional energystability, solvability, positivity and error estimates of the scheme. In particular, uniqueness of solution of the scheme is proved under a hypothesis that relates the time step and a strong norm of the scheme (the discrete version of (1.32)), which can be simplified in the case of 2D 6
domains owing to the strong estimates already obtained for the scheme. Moreover, we prove the existence of weak solutions of model (1.2) throughout the convergence of this scheme when the time step goes to 0. In Section 1.5, we propose a linearized coupled scheme for model (1.3)-(1.4), and again we analyze some properties of this linear scheme as in Section 1.4, comparing to the previous nonlinear scheme. Finally, in Section 1.6, we show some numerical simulations using Finite Elements spatial approximations associated to both time schemes, in order to verify numerically the theoretical results obtained in terms of positivity and unconditional energy-stability. 1.2 Notations and preliminary results We recall some functional spaces which will be used throughout this paper. We will consider the usual Sobolev spaces Hm(Ω) and Lebesgue spaces Lp(Ω),1≤p≤ ∞,with the usual notations for norms k·kmand k·kLp, respectively. In particular, the L2(Ω)-norm will be represented by k·k0. Corresponding Sobolev spaces of vector valued functions will be denoted by H1(Ω),L2(Ω),and so on; and we denote by H1 σ(Ω) := {σ∈H1(Ω) : σ·n= 0 on ∂Ω} and H1 ∗(Ω) := {h∈H1(Ω) : RΩh= 0}. From now on, we will use the following equivalent norms in H1(Ω) and H1 σ(Ω), respectively (see [11] and [1, Corollary 3.5], respectively): kuk2 1=k∇uk2 0+ZΩ u2 ,∀u∈H1(Ω),(1.5) kσk2 1=kσk2 0+krot σk2 0+k∇·σk2 0,∀σ∈H1 σ(Ω).(1.6) In particular, (1.6) implies that, for all σ=∇v∈H1 σ(Ω), k∇vk2 1=k∇vk2 0+k∆vk2 0.(1.7) If Zis a Banach space, then Z0will denote its topological dual. Moreover, the letters C, K will denote different positive constants always independent of the time step. We define the linear elliptic operators Av =g⇔ −∆v+v=gin Ω, ∂v ∂n= 0 on ∂Ω,(1.8) and Bσ=f⇔(−∇(∇·σ) + rot (rot σ) + σ=fin Ω, σ·n= 0,[rot σ×n]tang =0on ∂Ω.(1.9) 7
The corresponding variational forms are given by A:H1(Ω) →H1(Ω)0and B:H1 σ(Ω) → H1 σ(Ω)0such that hAv, ¯vi= (∇v, ∇¯v)+(v, ¯v),∀v, ¯v∈H1(Ω), hBσ,¯ σi= (σ,¯ σ)+(∇·σ,∇· ¯ σ) + (rot σ,rot ¯ σ),∀σ,¯ σ∈H1 σ(Ω). We assume the H2and H3-regularity of problem (1.8) (see for instance [4]). Consequently, we assume the existence of some constants C > 0 such that kvk2≤CkAvk0∀v∈H2(Ω); kvk3≤CkAvk1∀v∈H3(Ω).(1.10) Then, if the right hand side of problem (1.9) is given by f=∇hwith h∈H1(Ω), then taking σ=∇v, we will prove the H2-regularity of problem (1.9) as follows: Lemma 1.2.1 If f=∇hwith h∈H1(Ω), then the solution σof problem (1.9) belongs to H2(Ω). Moreover, kσk2≤Ck∇hk0.(1.11) Proof. First, we assume that h∈H1 ∗(Ω), hence khk1≤Ck∇hk0. Then, from H3- regularity of problem (1.8) taking g=h, we have that v∈H3(Ω) with −∆v+v=hand kvk3≤Ckhk1≤Ck∇hk0. Then, taking σ=∇v, we have that σ∈H2(Ω) solves (1.9), and (1.11) holds. Finally, in the general case of h∈H1(Ω), we consider g=h−1 |Ω|RΩhin (1.8), deducing again that v∈H3(Ω) and kvk3≤Ckh−1 |Ω|RΩhk1≤Ck∇(h−1 |Ω|RΩh)k0= Ck∇hk0. Then, taking σ=∇v, we have that σ∈H2(Ω) solves (1.9) with f=∇(h− 1 |Ω|RΩh) = ∇h, and (1.11) holds. Along this paper, we will use repeatedly the classical interpolations inequalities kukL4≤Ckuk1/2 0kuk1/2 1∀u∈H1(Ω) (in 2D domains),(1.12) kukL3≤Ckuk1/2 0kuk1/2 1∀u∈H1(Ω) (in 3D domains).(1.13) Finally, in order to obtain uniform in time strong estimates for the continuous problem and the numerical schemes, we will use the following results (see [14] and [13], respectively): Lemma 1.2.2 (Uniform Gronwall Lemma) Let g=g(t),h=h(t),z=z(t)be three positive locally integrable functions defined in (0,+∞)with z0also locally integrable in (0,+∞), such that z0(t)≤g(t)z(t) + h(t)a.e. t≥0. 8
If for any r > 0there exist a1(r),a2(r)and a3(r)such that Zt+r t g(s)ds ≤a1(r),Zt+r t h(s)ds ≤a2(r),Zt+r t z(s)ds ≤a3(r),∀t≥0 then, z(t+r)≤a2(r) + a3(r) rexp(a1(r)),∀t≥0. Lemma 1.2.3 (Uniform discrete Gronwall lemma) Let k > 0and dn, gn, hn≥0such that dn+1 −dn k≤gndn+hn,∀n≥0.(1.14) If for any r∈N, there exist a1(tr),a2(tr)and a3(tr)depending on tr=kr, such that k n0+r−1 X n=n0 gn≤a1(tr), k n0+r−1 X n=n0 hn≤a2(tr), k n0+r−1 X n=n0 dn≤a3(tr),∀n0≥0, then dn≤a2(tr) + a3(tr) trexp {a1(tr)},∀n≥r. As consequence of Lemma 1.2.3 (with an estimate for dnfor any n≥r) and the classical Discrete Gronwall Lemma (estimating dnfor n= 0, . . . , r −1), we will prove the following result: Corollary 1.2.4 Assume hypothesis of Lemma 1.2.3. Let k0>0be fixed, then the following estimate holds for all k≤k0 dn≤C(d0, k0)∀n≥0.(1.15) Proof. We fix T= 2k0and choose r0∈Nsuch that k(r0−1) < T ≤kr0:= tr0. Then, from Lemma 1.2.3 we have dn≤a2(tr0) + a3(tr0) tr0exp {a1(tr0)} ≤a2(tr0) + a3(tr0) Texp {a1(tr0)}:= C1(k0),∀n≥r0.(1.16) On the other hand, applying the Discrete Gronwall Lemma to (1.14), one has dn≤a2(tr0) + d0exp {a1(tr0)}:= C2(d0, k0),∀n < r0.(1.17) Therefore, from (1.16)-(1.17) we deduce (1.15). 9
Therefore, since k∇vk4 1and k∇uk4 0∈L1(0,+∞) (owing to (1.30)2and (1.32)), Lemma 1.2.2 and (1.40) imply (∂tu, ∂tv)∈L∞(0,+∞;L2(Ω) ×H1(Ω)).(1.43) Moreover, integrating in time (1.42) and using (1.40) and (1.43), we deduce (∂tu, ∂tv)∈L2(0,+∞;H1(Ω) ×H2(Ω)).(1.44) Finally, applying time derivative to equations (1.2)1and (1.2)2, and taking into account (1.32)-(1.35), (1.38) and (1.43)-(1.44), we can deduce the following regularity for ∂ttuand ∂ttv: (∂ttu, ∂ttv)∈L2(0,+∞;H1(Ω)0×L2(Ω)).(1.45) By following with a bootstrap argument, it is possible to obtain more regularity for (u, v). However, the regularity obtained so far is sufficient to guarantee the hypothesis required later to prove error estimates (see Theorem 1.4.21 and Theorem 1.5.6). 1.3.4 Proof of (1.32) in 2D domains In order to prove (1.32) in 2D domains, we make (∇(1.2)1,∇u)+ 1 2(∆(1.2)2,∆v), integrating by parts, and arguing as in (1.33), but in this case using the 2Dinterpolation inequality (1.12), one has 1 2 d dt k∇uk2 0+1 2k∆vk2 0+k∆uk2 0+1 2k∇vk2 0+1 2k∇(∆v)k2 0 ≤Ck∇uk2 L4k∇2vk0≤Ck∇uk0k∇uk1k∇vk1≤εk∇uk2 1+Cεk∇uk2 0k∇vk2 1. (1.46) Therefore, if we add (1.22) and (1.46), use (1.7) and take εsmall enough, we have d dt kuk2 1+1 2k∇vk2 1+k∇uk2 1+k∇vk2 2≤Ck∇uk2 0k∇vk2 1 and thus, since k∇uk2 0∈L1(0,+∞), Lemma 1.2.2 implies (1.32). 1.4 Euler time discretization In this section, we study the Euler time discretization for the problem (1.3)-(1.4), and we analyze the unconditional stability (in weak norms, see Definition 1.4.5 below) and solvability 16
of the scheme, as well as its convergence towards weak solutions. We also study some additional properties as positivity of the cell and chemical variables, mass-conservation of cells, and error estimates. Additionally, we prove uniqueness of solution of the scheme under a hypothesis that relates the time step and strong norms of the scheme, which is simplified in the case of 2D domains owing to a strong estimate obtained for the scheme (the discrete version of (1.32)). Let us consider a fixed partition of the time interval [0,+∞) given by tn=nk, where k > 0 denotes the time step (that we take constant for simplicity). First, taking into account the (u, v)-problem (1.2), we consider the following first order, nonlinear and coupled scheme (Backward Euler): •Scheme UV : Initialization: We fix u0=u(0) and v0=v(0). Time step n: Given (un−1, vn−1)∈H1(Ω) ×H2(Ω) with un−1≥0 and vn−1≥0, compute (un, vn)∈H1(Ω) ×H2(Ω) with un≥0 and vn≥0 and solving ((δtun,¯u)+(∇un,∇¯u)+(un∇vn,∇¯u)=0,∀¯u∈H1(Ω), δtvn+Avn−u2 n= 0,a.e. x∈Ω,(1.47) where, in general, we denote δtan=an−an−1 k. On the other hand, we can consider the following Backward Euler scheme related to the reformulation in the (u, σ)-problem (1.3). Then, one also has the following first order, nonlinear and coupled scheme: •Scheme US: Initialization: We fix v0=v(0) and (u0,σ0) = (u(0),σ(0)), with σ0=∇v0. Time step n: Given (un−1,σn−1)∈H1(Ω) ×H1 σ(Ω) with un−1≥0, compute (un,σn)∈H1(Ω) ×H1 σ(Ω) with un≥0 and solving ((δtun,¯u)+(∇un,∇¯u)+(unσn,∇¯u) = 0,∀¯u∈H1(Ω), (δtσn,¯ σ) + hBσn,¯ σi−2(un∇un,¯ σ) = 0,∀¯ σ∈H1 σ(Ω).(1.48) Once the Scheme US is solved , given vn−1∈H2(Ω) with vn−1≥0, we can recover vn=vn(u2 n)∈H2(Ω) (with vn≥0) solving: δtvn+Avn=u2 n,a.e. x∈Ω.(1.49) 17
Remark 1.4.1 (Positivity and regularity of vn)It is not difficult to prove that, given un∈H1(Ω) and vn−1∈H2(Ω), there exists a unique vn∈H2(Ω) solution of (1.49). Even more, using the H3-regularity of problem (1.8), we can prove that vn∈H3(Ω). Moreover, if vn−1≥0then vn≥0. Indeed, testing by (vn)−= min{vn,0} ≤ 0in (1.49), and taking into account that (vn)−= 0 if (vn)≥0, as well as (vn)−∈H1(Ω) with ∇((vn)−) = ∇(vn)if (vn)≤0, and ∇((vn)−) = 0 if (vn)>0, we obtain 1 kk(vn)−k2 0−1 kZΩ vn−1(vn)−+k∇((vn)−)k2 0+k((vn)−)k2 0=ZΩ u2 n((vn)−)≤0.(1.50) Then, as vn−1≥0, from (1.50) we conclude k(vn)−k2 1≤0, which implies that (vn)−≡0a.e. x∈Ω, and thus, vn≥0in Ω. Lemma 1.4.2 If σn−1=∇vn−1, the schemes UV and US are equivalents in the following sense: If (un, vn)is a solution of scheme UV then (un,σn)with σn=∇vnsolves scheme US; and reciprocally, if (un,σn)is a solution of the scheme US and vn=vn(u2 n)is the unique solution of (1.49), then σn=∇vn, and therefore (un, vn)is a solution of the scheme UV. Proof. Suppose that (un, vn) is a solution of the scheme UV, then testing (1.47)2by −∇· ¯ w, for any ¯ w∈H1 σ(Ω), and taking into account that rot(∇vn) = 0, we obtain (δt∇vn,¯ w) + hB∇vn,¯ wi= 2(un∇un,¯ w),∀¯ w∈H1 σ(Ω).(1.51) Then, defining σn=∇vnand assuming the hypothesis σn−1=∇vn−1, from (1.51) we conclude that (un,σn) is solution of the scheme US. On the other hand, if (un,σn) is a solution of the scheme US and vnsatisfies (1.49), reasoning as above, we conclude that ∇vn satisfies (1.51). Therefore, from (1.48)2and (1.51), we obtain (δt(σn−∇vn),¯ σ) + hB(σn−∇vn),¯ σi= 0,∀¯ σ∈H1 σ(Ω).(1.52) Then, taking ¯ σ=σn−∇vnin (1.52) and using the formula a(a−b) = 1 2(a2−b2)+ 1 2(a−b)2, we deduce δt1 2kσn−∇vnk2 0+k 2kδt(σn−∇vn)k2 0+kσn−∇vnk2 1= 0, which implies that σn=∇vnusing that σn−1=∇vn−1. Thus, we conclude that (un, vn) is solution of the scheme UV. 18
Although, as it was said at the beginning, both time schemes UV and US are equivalent, we will study the scheme US with the variable σnin order to facilitate the notation throughout the paper. Moreover, both schemes furnish different fully discrete schemes considering for instance the spatial approximation by Finite Elements, which will be analyzed in Chapter 2. In fact, it will be necessary to use the variable σnin order to obtain a fully discrete unconditional energy-stable scheme. 1.4.1 Solvability, Energy-Stability and Convergence Taking ¯u= 1 in (1.48)1we see that the scheme US is conservative, that is: ZΩ un=ZΩ un−1=···=ZΩ u0.(1.53) Moreover, integrating (1.49) in Ω, we deduce the following discrete in time equation for ZΩ vn: δtZΩ vn+ZΩ vn=ZΩ u2 n.(1.54) Theorem 1.4.3 (Unconditional existence and conditional uniqueness) There exists (un,σn)∈H1(Ω) ×H1 σ(Ω) solution of the scheme US, such that un≥0. Moreover, if kk(un,σn)k4 1is small enough (1.55) then the solution of the scheme US is unique. Proof. Let (un−1,σn−1)∈H1(Ω) ×H1 σ(Ω) be given, with un−1≥0. The proof is divided into two parts. Part 1: In order to prove the existence of (un,σn)∈H1(Ω) ×H1 σ(Ω) solution of the scheme US, such that un≥0, we consider the following auxiliary problem: ((δtun,¯u)+(∇un,∇¯u) + ((un)+σn,∇¯u) = 0,∀¯u∈H1(Ω), (δtσn,¯ σ) + hBσn,¯ σi−2(un∇un,¯ σ)=0,∀¯ σ∈H1 σ(Ω),(1.56) where (un)+= max{un,0}. In fact, it is the same scheme US but changing unby (un)+in the chemotaxis term. 19
A. Positivity of un:First, we will see that if (un,σn) is a solution of (1.56), then un≥0. Testing by ¯u= (un)−= min{un,0} ≤ 0 in (1.56), and using that ((un)+σn,∇((un)−)) = 0, we obtain 1 kk(un)−k2 0−1 kZΩ un−1(un)−+k∇((un)−)k2 0= 0. Then, using the fact that un−1≥0, one has that k(un)−k2 1≤0, and thus (un)−≡0 a.e. x∈Ω. Therefore, un≥0 in Ω. B. Existence of solution of (1.56): It can be proven by using the Leray-Schauder fixed-point theorem. (See Appendix A). Then, from parts A and B, we conclude that there exists (un,σn) solution of (1.56) with un≥0. In particular, taking into account that un= (un)+, we conclude that (un,σn) is also a solution of the scheme US, with un≥0. Part 2: In order to prove the uniqueness of solution (un,σn) of the scheme US, we suppose that there exist (u1 n,σ1 n),(u2 n,σ2 n)∈H1(Ω) ×H1 σ(Ω) two possible solutions of (1.48). Then, defining un=u1 n−u2 nand σn=σ1 n−σ2 n, we have that (un,σn)∈H1(Ω) ×H1 σ(Ω) satisfies 1 k(un,¯u)+(∇un,∇¯u)+(u1 nσn,∇¯u)+(unσ2 n,∇¯u)=0,∀¯u∈H1(Ω),(1.57) 1 k(σn,¯ σ) + hBσn,¯ σi−2(u1 n∇un,¯ σ)−2(un∇u2 n,¯ σ) = 0,∀¯ σ∈H1 σ(Ω).(1.58) Taking ¯u=un,¯ σ=1 2σnin (1.57)-(1.58) and adding the resulting expressions, the terms (u1 nσn,∇un) cancel, and using the fact that ZΩ un= 0, we obtain 1 2kk(un,σn)k2 0+1 2k(un,σn)k2 1 ≤ kunkL3kσ2 nkL6k∇unk0+kunkL3k∇u2 nk0kσnkL6 ≤Ckunk1/2 0kσ2 nkL6kunk3/2 1+Ckunk1/2 0kunk1/2 1k∇u2 nk0kσnk1 ≤1 4k(un,σn)k2 1+Ckunk2 0kσ2 nk4 1+Ckunk2 0k∇u2 nk4 0, which implies that 1 2k(un,σn)k2 0+k 4k(un,σn)k2 1≤C k k∇u2 nk4 0+kσ2 nk4 1k(un,σn)k2 0. Therefore, since k(k∇u2 nk4 0+kσ2 nk4 1) is small enough (from hypothesis (1.55)), we conclude that k(un,σn)k1= 0, thus u1 n=u2 nand σ1 n=σ2 n. 20
Remark 1.4.4 In the case of 2Ddomains, using estimate (1.75) (see Theorem 1.4.20 below), the uniqueness restriction (1.55) can be relaxed to k K2 0small enough, where K0is the constant appearing in (1.75) which depends on data (Ω, u0,σ0), but is independent of n. Definition 1.4.5 A numerical scheme with solution (un,σn)is called energy-stable with respect to the energy E(u, σ)given in (1.26) if this energy is time decreasing, that is E(un,σn)≤E(un−1,σn−1),∀n. (1.59) In the next Lemma, we obtain unconditional energy-stability for the scheme US. Lemma 1.4.6 (Unconditional stability) The scheme US is unconditionally energy-stable with respect to E(u, σ). In fact, for any (un,σn)solution of scheme US, the following discrete energy law holds δtE(un,σn)+k 2kδtunk2 0+k 4kδtσnk2 0+k∇unk2 0+1 2kσnk2 1= 0.(1.60) Proof. Taking ¯u=unin (1.48)1and ¯ σ=1 2σnin (1.48)2and adding the resulting expressions, the chemotaxis and production terms cancel, obtaining (1.60). Remark 1.4.7 Comparing the energy law (1.22) of the continuous problem, and the discrete version (1.60), we can say that the scheme US introduces the following two first order “numerical dissipation terms”: k 2kδtunk2 0and k 4kδtσnk2 0. From the (local in time) discrete energy law (1.60), we deduce the following global in time estimates for any (un,σn) solution of the scheme US: Theorem 1.4.8 (Uniform Weak estimates) Let (un,σn)be a solution of the scheme US. Then, the following estimates hold k(un,σn)k2 0+k2 n X m=1k(δtum, δtσm)k2 0+k n X m=1k(∇um,σm)k2 L2×H1≤C0,∀n≥1,(1.61) k n0+n X m=n0+1k(um,σm)k2 1≤C0+C1(nk),∀n≥1,(1.62) where n0≥0is any integer and C0, C1are positive constants depending on the data (u0,σ0) and (Ω, u0)respectively, but independent of n0,kand n. 21
Proof. Observe that from the discrete energy law (1.60), we have 1 4k(un,σn)k2 0+k2 4 n X m=1k(δtum, δtσm)k2 0+k 2 n X m=1k(∇um,σm)k2 L2×H1≤1 2k(u0,σ0)k2 0, which implies (1.61). Moreover, starting again from (1.60), but now summing for mfrom n0+ 1 to n+n0, using (1.61) and the Poincar´e inequality for the zero-mean value function um−m0, where m0=1 |Ω|ZΩ u0=1 |Ω|ZΩ um, we have k n0+n X m=n0+1k(um−m0,σm)k2 1≤C0, and thus, we deduce (1.62). Remark 1.4.9 The proofs of solvability (without positivity) and unconditional energy-stability of the scheme US (see Theorems 1.4.3 and 1.4.8, and Lemma 1.4.6) can be followed almost line by line if we consider a fully discrete scheme corresponding to a Finite Element approximation of US, that is, if we take any finite-dimensional subspaces Uh⊂H1(Ω) and Σh⊂H1 σ(Ω) instead of H1(Ω) and H1 σ(Ω) respectively. Corollary 1.4.10 (Estimates for vn)If vn=vn(u2 n)is the solution of (1.49), it holds kvnkL1≤K0,∀n≥0,(1.63) where K0>0depends on the data (u0,σ0, v0), but independent of kand n. Moreover, the following estimates hold kvnk2 1≤K0,and k n0+n X m=n0+1kvmk2 2≤K0+K1(nk),∀n≥1,(1.64) with K1>0depending on the data u0,σ0, v0,Ω, but independent of n0,kand n. Proof. From (1.54) and (1.61) we have (1 + k)kvikL1−kvi−1kL1=kkuik2 0≤kC0.(1.65) Then, multiplying (1.65) by (1 + k)i−1and summing for i= 1,···, n, we obtain (1 + k)nkvnkL1≤ kv0kL1+k C0 n−1 X i=0 (1 + k)i 22
and taking into account that n−1 X i=0 (1 + k)i=1−(1 + k)n 1−(1 + k)=1 k((1 + k)n−1) ≤1 k(1 + k)n, we conclude kvnkL1≤(1 + k)−nkv0kL1+C0≤ kv0kL1+C0, which implies (1.63). Finally, taking into account the relation σn=∇vn, from (1.61)-(1.63), we can deduce (1.64). Starting from the previous stability estimates, following the ideas of [9] we can prove the convergence towards weak solutions. For this, let us to introduce the functions: •(euk,e σk) are continuous functions on [0,+∞), linear on each interval (tn−1, tn) and equal to (un,σn) at t=tn,n≥0; •(ur k,σr k) as the piecewise constant functions taking values (un,σn) on (tn−1, tn], n≥1. Theorem 1.4.11 (Convergence) There exists a subsequence (k0)of (k), with k0↓0, and a weak solution (u, σ)of (1.3) in (0,+∞), such that (euk0,e σk0)and (ur k0,σr k0)converge to (u, σ) weakly-?in L∞(0,+∞;L2(Ω) ×L2(Ω)), weakly in L2(0, T;H1(Ω) ×H1(Ω)) and strongly in L2(0, T;L2(Ω) ×L2(Ω)), for any T > 0. Proof. Observe that (1.48) can be rewritten as: d dteuk(t),¯u+ (∇ur k(t),∇¯u)+(ur k(t)σr k(t),∇¯u)=0,∀¯u= ¯u(t),for t∈[0,+∞)\{tn}, d dt e σk(t),¯ σ+hBσr k(t),¯ σi−2(ur k(t)∇ur k(t),¯ σ)=0,∀¯ σ=¯ σ(t),for t∈[0,+∞)\{tn}. (1.66) where ¯u(t)|In= ¯un∈H1(Ω) and ¯ σ(t)|In=¯ σn∈H1 σ(Ω), with In:= [tn−1, tn]. From Theorem 1.4.8 we have that (euk,e σk) and (ur k,σr k) are bounded in L∞(0,+∞;L2(Ω)×L2(Ω))∩ L2(0, T;H1(Ω) ×H1(Ω)). Moreover, using (1.61), it is not difficult to prove that euk−ur k and e σk−σr kconverge to 0 in L2(0, T;L2(Ω)) as k→0, for any T > 0. More precisely, we have keuk−ur k,e σk−σr kkL2(0,T;L2(Ω)) ≤(C0k/3)1/2. Therefore, there exists a subsequence (k0) of (k) and limit functions uand σverifying the following convergence as k0→0: (euk0,e σk0)→u, (ur k0,σr k0)→uin (L∞(0,+∞;L2(Ω) ×L2(Ω))-weak* L2(0, T;H1(Ω) ×H1(Ω))-weak. 23
Moreover, we can deduce d dt (euk0,e σk0) is bounded in L4/3(0, T;H1(Ω)0×H1(Ω)0). Therefore, a compactness result of Aubin-Lions type implies that the sequence is compact in L2(0, T;L2(Ω) ×L2(Ω)). This implies the strong convergence of both subsequences (euk,e σk) and (ur k,σr k) in L2(0, T;L2(Ω) ×L2(Ω)); and passing to the limit in (1.66), we obtain that (u, σ) satisfies (1.23)-(1.24). Now, in order to obtain that (u, σ) satisfies energy inequality (1.25), we test (1.66)1by ¯u=ur k(t) and (1.66)2by ¯ σ=1 2σr k(t), and taking into account that euk|In=un+tn−t k(un−1− un) (and e σkis defined in the same way), we obtain d dt 1 2keuk(t)k2 0+1 4ke σk(t)k2 0+ (tn−t)k(δtun, δtσn)k2 0+k∇ur k(t)k2 0+1 2kσr k(t)k2 1= 0, for any t∈(tn−1, tn), which implies that d dtE(euk(t),e σk(t)) + k∇ur k(t)k2 0+1 2kσr k(t)k2 1≤0,for t∈[0,+∞)\{tn}.(1.67) Then, integrating (1.67) in time from t0to t1, with t0, t1∈[0,+∞), and taking into account that Zt1 t0 d dtE(euk(t),e σk(t)) = E(euk(t1),e σk(t1)) −E(euk(t0),e σk(t0)) ∀t0< t1 since E(euk(t),e σk(t)) is continuous in time, we deduce E(euk(t1),e σk(t1))−E(euk(t0),e σk(t0))+Zt1 t0 (k∇ur k(t)k2 0+1 2kσr k(t)k2 1)dt ≤0,∀t0< t1.(1.68) Finally, we will prove that E(euk0(t),e σk0(t)) →E(u(t),σ(t)),a.e. t∈[0,+∞).(1.69) Indeed, for any T > 0, kE(euk0(t),e σk0(t)) −E(u(t),σ(t))kL1(0,T)=ZT 0|E(euk0(t),e σk0(t)) −E(u(t),σ(t))|dt =ZT 0 1 2keuk0(t)k2 0−ku(t)k2 0+1 4ke σk0(t)k2 0−kσ(t)k2 0dt ≤1 2keuk0−ukL2(0,T;L2)(keuk0kL2(0,T;L2)+kukL2(0,T;L2)) +1 4ke σk0−σkL2(0,T;L2)(ke σk0kL2(0,T;L2)+kσkL2(0,T;L2)),(1.70) and taking into account that (euk0,e σk0)→(u, σ) strongly in L2(0, T;L2(Ω)) for any T > 0, from (1.70) we conclude that E(euk0(t),e σk0(t)) →E(u(t),σ(t)) strongly in L1(0, T) for all 24
T > 0, which implies (1.69). Then, taking into account that (ur k0,σr k0)→(u, σ) weakly in L2(0, T;H1(Ω) ×H1(Ω)), we deduce lim inf k0→0Zt1 t0 (k∇ur k0(t)k2 0+1 2kσr k0(t)k2 1)dt ≥Zt1 t0 (k∇u(t)k2 0+1 2kσ(t)k2 1)dt ∀t1≥t0≥0 and, owing to (1.69), lim inf k0→0hE(euk0(t1),e σk0(t1)) −E(euk0(t0),e σk0(t0))i=E(u(t1),σ(t1)) −E(u(t0),σ(t0)), for a.e. t1, t0:t1≥t0≥0. Thus, taking lim inf as k0→0 in the inequality (1.68), we deduce the energy inequality (1.25) for a.e. t0, t1:t1≥t0≥0. Analogously, if we introduce the functions: •evkare continuous functions on [0,+∞), linear on each interval (tn−1, tn) and equal to vn,at t=tn,n≥0; •vr kas the piecewise constant functions taking values vnon (tn−1, tn], n≥1, the following result can be proved: Lemma 1.4.12 There exists a subsequence (k0)of (k), with k0↓0, and a strong solution vof (1.4) in (0,+∞), such that evk0and vr k0converge to vweakly-?in L∞(0,+∞;H1(Ω)), weakly in L2(0, T;H2(Ω)) and strongly in L2(0, T;H1(Ω)) ∩C([0, T]; Lp(Ω)), for 1≤p < 6 and any T > 0. Proof. Observe that (1.49) can be rewritten as: d dtevk(t) + Avr k(t)=(ur k(t))2,for t∈[0,+∞)\{tn},(1.71) From estimate (1.64) we have that evkand vr kare bounded in L∞(0,+∞;H1(Ω))∩L2(0, T;H2(Ω)). Moreover, it is not difficult to prove that evk−vr kconverge to 0 in L2(0, T;H1(Ω)) as k→0, for any T > 0. Therefore, there exists a subsequence (k0) of (k) and a limit function v verifying the following convergence as k0→0: evk0→v, vr k0→vin (L∞(0,+∞;H1(Ω))-weak* L2(0, T;H2(Ω))-weak. Moreover, we can deduce that d dt evk0is bounded in L4/3(0, T;L2(Ω)). Therefore, a compactness result of Aubin-Lions type implies that the sequence is compact in L2(0, T;H1(Ω)) ∩ C([0, T]; Lp(Ω)), for 1 ≤p < 6. This implies the strong convergence of both subsequences evkand vr kin L2(0, T;H1(Ω)) ∩C([0, T]; Lp(Ω)); and passing to the limit in (1.71), we obtain that vsatisfies (1.19). 25
Theorem 1.4.21 Let (un,σn)be a solution of the scheme US and assume the following regularity for the exact solution (u, σ)of (1.3): (u, σ)∈L∞(0, T;H1(Ω) ×H1(Ω)) and (utt,σtt)∈L2(0, T;H1 ∗(Ω)0×H1 σ(Ω)0).(1.99) Assuming that kk(∇u, ∇·σ)k4 L∞(L2)is small enough, (1.100) then the a priori error estimate k(en u, en σ)kl∞L2∩l2H1≤C(T)k(1.101) holds, where C(T) = K1exp(K2T), with K1, K2>0independent of k. Proof. Since u0=u(t0), then RΩen u=RΩe0 u= 0. Moreover, taking ¯u=en uin (1.97), ¯ σ=1 2en σin (1.98), and adding the resulting expressions, the terms (unen σ,∇en u) cancel, and using the H¨older and Young inequalities and (1.13), we obtain δt1 2ken uk2 0+1 4ken σk2 0+k 4k(δten u, δten σ)k2 0+1 2k(en u, en σ)k2 1 ≤1 8k(en u, en σ)k2 1+Ck(ξn 1, ξn 2)k2 (H1 ∗)0×(H1 σ)0+1 8k(en u, en σ)k2 1+Ck(∇u(tn),∇·σ(tn))k4 0ken uk2 0, and therefore δt1 2ken uk2 0+1 4ken σk2 0+1 4k(en u, en σ)k2 1 ≤Ck(ξn 1, ξn 2)k2 (H1 ∗)0×(H1 σ)0+Ck(∇u(tn),∇·σ(tn))k4 0ken uk2 0.(1.102) Now, taking into account that (ξn 1, ξn 2) = −1 kZtn tn−1 (t−tn−1)(utt(t),σtt(t))dt, using the H¨older inequality, we can deduce Ck(ξn 1, ξn 2)k2 (H1 ∗)0×(H1 σ)0≤Ck Ztn tn−1k(utt(t),σtt(t))k2 (H1 ∗)0×(H1 σ)0dt. (1.103) Therefore, from (1.102) and (1.103) we deduce δt1 2ken uk2 0+1 4ken σk2 0+1 4k(en u, en σ)k2 1 ≤Ck Ztn tn−1k(utt(t),σtt(t))k2 (H1 ∗)0×(H1 σ)0dt +Ck(∇u(tn),∇·σ(tn))k4 0ken uk2 0.(1.104) 32
Then, multiplying (1.104) by kand adding from n= 1 to n=r, we obtain (recall that e0 u=e0 σ= 0): 1 4−Ckk(∇u, ∇·σ)k4 L∞L2k(er u, er σ)k2 0+k 4 r X n=1k(en u, en σ)k2 1 ≤Ck2Ztr 0k(utt(t),σtt(t))k2 (H1 ∗)0×(H1 σ)0dt +Ck(∇u, ∇·σ)k4 L∞L2k r−1 X n=0k(en u, en σ)k2 0.(1.105) Therefore, if hypothesis (1.100) is satisfied, using the Discrete Gronwall Lemma (see [7], p. 369) in (1.105), and taking into account (1.99), we conclude (1.101). Remark 1.4.22 From (1.101), we deduce k(un,σn)k1≤C(T), for all n= 1, ..., N. In particular, this implies that in 3D domains, for finite time, the hypothesis (1.55) assuring the uniqueness of solution (un,σn)can be relaxed to k C(T)4small enough. Error estimates for vnsolution of (1.49) We will obtain error estimates for vnsolution of (1.49) with respect to a sufficiently regular solution vof (1.4). For this, we introduce the following notation for the error in t=tn: en v=v(tn)−vnand for the discrete in time derivative of this error: δten v=en v−en−1 v k. Then, subtracting (1.4) at t=tnand (1.49), we obtain that en vsatisfies δten v+Aen v= (u(tn) + un)en u+ξn 3,a.e. x∈Ω,(1.106) where ξn 3is the consistency error associated to (1.49), that is, ξn 3=δt(v(tn)) −vt(tn) = −1 kZtn tn−1 (t−tn−1)vtt(t)dt. (1.107) Theorem 1.4.23 Under hypothesis of Theorem 1.4.21. Let vnbe the solution of (1.49) and assume the following regularity for the exact solution vof (1.4): vtt ∈L1(0, T;H1(Ω)0).(1.108) Then the a priori error estimate holds ken vkl∞H1∩l2H2≤C(T)k(1.109) where C(T) = K1exp(K2T), with K1, K2>0independent of k. 33
Proof. From the relation en σ=∇en v, taking into account (1.101), we only need to prove the following estimate ZΩ en v≤C(T)k. (1.110) With this aim, if we integrate (1.106) in Ω, δtZΩ en v+ZΩ en v=ZΩ (u(tn) + un)en u+ZΩ ξn 3,(1.111) and therefore, multiplying (1.111) by kand using (1.107), we have (1 + k)ZΩ en v−ZΩ en−1 v≤kZΩ (u(tn) + un)en u+Ztn tn−1ZΩ (t−tn−1)vtt(x, t)dxdt ≤kku(tn) + unk0ken uk0+k|Ω|1/2Ztn tn−1kvtt(t)k(H1)0.(1.112) Then, adding from n= 1 to n=rin (1.112) and taking into account that u(tn) + unis bounded in l∞L2, we obtain (recall that e0 v= 0) ZΩ er v+k r X n=1 ZΩ en v≤Ckkvtt(t)kL1(H1)0+kC r X n=1ken uk0.(1.113) Thus, using (1.108) and (1.101) in (1.113), we deduce (1.110). 1.5 A linear scheme In this section, we propose the following first order in time, linear coupled scheme for model (1.3): •Scheme LC : Initialization: We fix v0=v(0) and (u0,σ0) = u(0),σ(0)), with σ0=∇v0. Time step n: Given (un−1,σn−1)∈H1(Ω) ×H1 σ(Ω), compute (un,σn)∈H1(Ω) × H1 σ(Ω) solving ((δtun,¯u)+(∇un,∇¯u) = −(un−1σn,∇¯u),∀¯u∈H1(Ω), (δtσn,¯ σ) + hBσn,¯ σi= 2(un−1∇un,¯ σ),∀¯ σ∈H1 σ(Ω).(1.114) Again, once solved the scheme LC, given vn−1∈H2(Ω) with vn−1≥0, we can recover vn=vn(u2 n)∈H2(Ω) solving (1.49). 34
1.5.1 Unconditional energy-stability and Unique Solvability Observe that scheme LC is also conservative in u(satisfying (1.53)) and also has the behavior for RΩvngiven in (1.54). Theorem 1.5.1 (Unconditional Unique Solvability) There exists a unique (un,σn) solution of the scheme LC. Proof. Let (un−1,σn−1)∈H:= H1(Ω) ×H1 σ(Ω) be given, and consider the following bilinear form a:H×H→R, and the linear operator l:H→Rgiven by a((un,σn),(¯u, ¯ σ)) = 2 k(un,¯u) + 1 k(σn,¯ σ) + 2(∇un,∇¯u) + hBσn,¯ σi +2(un−1σn,∇¯u)−2(un−1∇un,¯ σ), l((¯u, ¯ σ)) = 2 k(un−1,¯u) + 1 k(σn−1,¯ σ), for all (un,σn),(¯u, ¯ σ)∈H. Then, using the H¨older inequality and Sobolev embeddings, we can verify that a(·,·) is continuous and coercive on H, and l∈H0. Thus, from Lax-Milgram theorem, there exists a unique (un,σn)∈Hsuch that a((un,σn),(¯u, ¯ σ)) = l((¯u, ¯ σ)),∀(¯u, ¯ σ)∈H. Finally, taking first ¯ σ= 0 and then ¯u= 0, implies that (un,σn)∈H1(Ω) ×H1 σ(Ω) is the unique solution of (1.114). Moreover, following the proof of Lemma 1.4.6, we can prove unconditional energy-stability of the scheme LC. Lemma 1.5.2 (Unconditional energy-stability) The scheme LC is unconditionally energystable for E(u, σ). In fact, the same discrete energy law (1.60) holds. Remark 1.5.3 If we consider the fully discrete scheme corresponding to LC via a spatial approximation by using the Finite Elements method, i.e. taking Uh⊂H1(Ω) and Σh⊂ H1 σ(Ω) instead of H1(Ω) and H1 σ(Ω) respectively, then the proofs of solvability and unconditional energy-stability of this fully discrete scheme can be followed line by line from Theorem 1.5.1 and Lemma 1.5.2. Remark 1.5.4 We can prove weak estimates for the solution (un,σn)of the scheme LC analogously to Theorem 1.4.8. Moreover, assuming the H2-regularity for problem (1.9) in the case that the right hand side is not the gradient of a function, we can deduce strong and more regular estimates for this solution (un,σn)as in Subsection 1.4.2. 35
Remark 1.5.5 Unlike the scheme US, in the scheme LC it is not clear how to prove the nonnegativity of un. In fact, in some numerical simulations, very negative cell densities are obtained when h→0, where his the spatial parameter (see Subsection 1.6.1). 1.5.2 Error estimates in weak norms Theorem 1.5.6 (Error estimates for the scheme LC) Let (un,σn)be a solution of the scheme LC, and assume the regularity (1.99). Then, the a priori error estimate (1.101) holds. Proof. The proof follows as in Theorem 1.4.21, but we recall that in this case we do not need to impose hypothesis of small time step given in (1.100) in order to apply the Discrete Gronwall Lemma, since we use the form of the terms (un−1σn,∇¯u) and (un−1∇un,¯ σ) to bound them in a suitable way. Moreover, although in this linear scheme LC it is not clear if the relation σn=∇vn holds, it will be possible to obtain error estimates for vn. Theorem 1.5.7 (Error estimates for vn)Under hypothesis of Theorem 1.5.6. Let vnbe the solution of (1.49) (corresponding to the scheme LC), and assume the regularity: vtt ∈L2(0, T;L2(Ω)).(1.115) Then, the a priori error estimate (1.109) holds. Proof. Since in the scheme LC it is not clear the relation σn=∇vn, we will argue in a different way of Theorem 1.4.23. Indeed, we test (1.106) by Aen v, and using the H¨older and Young inequalities, we obtain 1 2δtken vk2 1+k 2kδten vk2 1+1 2kAen vk2 0≤C(ku(tn) + unk2 1ken uk2 1+kξn 3k2 0).(1.116) Observe that from (1.101) we have r P n=1ku(tn)−unk2 1≤C(T)k, which implies that kunk1≤C+ku(tn)k1.(1.117) Then, multiplying (1.116) by k, adding from n= 1 to n=rand using (1.117), we obtain (recall that e0 v= 0) ker vk2 1+k r X n=1kAen vk2 0≤C+Cku!‘k2 L∞H1k r X n=1ken uk2 1+Ck2Ztr 0kvtt(t)k2 0dt. (1.118) 36
Therefore, taking into account (1.99), (1.115) and (1.101), from (1.118) we conclude (1.109). 1.6 Numerical simulations The aim of this section is to compare the results of several numerical simulations that we have carried out using the schemes derived through the paper. We are considering a finite element discretization in space associated to the variational formulation of schemes US, LC and UV, where the P1-continuous approximation is taken for uh,σhand vh(where his the spatial parameter). Moreover, we have chosen the 2D domain Ω = [0,2]2using a structured mesh, and all the simulations are carried out using FreeFem++ software. The linear iterative method used to approach the nonlinear schemes US and UV is the Newton Method, and in all the cases, the iterative method stops when the relative error in L2-norm is less than ε= 10−6. 1.6.1 Positivity In this subsection, we compare the schemes US and LC in terms of positivity. For the fully discretization of both schemes is not clear the positivity of the variable uh. In fact, for the time-discrete scheme US the existence of nonnegative solution (un, vn) was proved (see Theorem 1.4.3 and Remark 1.4.1), but for the time-discrete scheme LC, although we can prove that vnis nonnegative, the nonnegativity of unis not clear. For this reason, in Figures 1.2-1.5, we compare the positivity of the variables uhand vhin both schemes taking meshes in space increasingly thinner (h=1 35,h=1 75 and h=1 150). In all the cases, we choose k= 10−5and the initial conditions are (see Figure 1.1) u0=−10xy(2 −x)(2 −y)exp(−10(y−1)2−10(x−1)2) + 10.0001 and v0=200xy(2 −x)(2 −y)exp(−30(y−1)2−30(x−1)2)+0.0001. In the case of the scheme US, we observe that although uhis negative for some x∈Ω in some times tn>0, when h→0 these values are closer to 0; while in the case of the scheme LC, when h→0 very negative cell densities uhare obtained for some x∈Ω in some times tn>0 (see Figures 1.2-1.4). On the other hand, the same behavior is observed for the minimum of vhin both schemes. In fact, independently of h, if v0is positive, then vhalso is positive (we show this behavior in Figure 1.5 for the case h=2 70, but the same holds for the cases h=2 150 and h=2 300). 37
(a) Initial cell density u0(b) Initial chemical concentration v0 Figure 1.1: Initial conditions. Remark 1.6.1 In Figures 1.3 and 1.4 there are also negative values of minimum of uhfor the scheme US, but those are of order 10−2and 10−4respectively. 1.6.2 Unconditional Stability In this subsection, we compare numerically the stability with respect to two energies E(u, σ) and E(u, v). Following line to line the proof of Lemma 1.4.6, we deduce the unconditional energy-stability for the fully discrete schemes corresponding to schemes US and LC (for the modified energy E(u, σ)). In fact, if (un,σn) is any solution of the fully discrete schemes corresponding to schemes US or LC, the following relation holds RE(un,σn) := δtE(un,σn) + k∇unk2 0+1 2kσnk2 1≤0,∀n. (1.119) However, considering the “exact” energy (1.21), in the case of fully discrete schemes, it is not clear how to prove unconditional energy-stability of schemes US,LC and UV with respect to this energy. Moreover, it is interesting to study the behaviour of the corresponding residual RE(un, vn) := δtE(un, vn) + k∇unk2 0+1 2k∆hvnk2 0+1 2k∇vnk2 0. Indeed, taking k= 10−6,h=2 50 and the initial conditions u0=−10xy(2 −x)(2 −y)exp(−10(y−1)2−10(x−1)2) + 10.0001 and v0=20xy(2 −x)(2 −y)exp(−30(y−1)2−30(x−1)2)+0.0001, we obtain that: 38
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 Time -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 Minimum of u h US LC Figure 1.2: Minimum values of uh, with h=1 35. (a) The schemes US and LC satisfy the energy decreasing in time property (1.59) for the modified energy E(u, σ), see Figure 1.6. (b) The schemes US and LC satisfy (1.119), see Figure 1.7. (c) The schemes US,LC and UV satisfy the energy decreasing in time property for the exact energy E(u, v), that is, E(un, vn)≤E(un−1, vn−1) for all n, see Figure 1.8. (d) The schemes US,LC and UV have RE(un, vn)>0 for some tn≥0. However, it is observed that the residues RE(un, vn) of the schemes US and LC in those times tn where each residue is positive, its values are less than the residues of the scheme UV, see Figure 1.9. 1.7 Appendix A In order to prove the solvability of (1.56), we will use the Leray-Schauder fixed point theorem. With this aim, we define the operator R:L4(Ω) ×L4(Ω) →L4(Ω) ×L4(Ω) by R(eu, e σ) = 39
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 Time ×10-3 -160 -140 -120 -100 -80 -60 -40 -20 0 20 Minimum of u h US LC Figure 1.3: Minimum values of uh, with h=1 75. (u, σ), such that (u, σ) solves the following linear decoupled problems 1 k(u, ¯u)+(∇u, ∇¯u) = 1 k(un−1,¯u)−(eu+e σ,∇¯u),∀¯u∈H1(Ω), 1 k(σ,¯ σ) + hBσ,¯ σi=1 k(σn−1,¯ σ)−(eu2,∇· ¯ σ),∀¯ σ∈H1 σ(Ω). (1.120) 1. Ris well defined. Let (eu, e σ)∈L4(Ω) ×L4(Ω) and consider the following bilinear forms ea:H1(Ω) ×H1(Ω) →R,eb:H1 σ(Ω) ×H1 σ(Ω) →R, and the linear operators l1:H1(Ω) →Rand l2:H1 σ(Ω) →Rgiven by ea(u, ¯u) = 1 k(u, ¯u)+(∇u, ∇¯u),eb(σ,¯ σ) = 1 k(σ,¯ σ) + hBσ,¯ σi, l1(¯u) = 1 kZΩ un−1¯u−ZΩeu+e σ·∇¯uand l2(¯ σ) = 1 kZΩ σn−1¯ σ−ZΩeu2∇· ¯ σ, for all u, ¯u∈H1(Ω) and σ,¯ σ∈H1 σ(Ω). Then, using the H¨older inequality and Sobolev embeddings, we can verify that eaand ebare continuous and coercive on H1(Ω) and H1 σ(Ω) respectively, and l1∈H1(Ω)0and l2∈H1 σ(Ω)0. Thus, from Lax-Milgram theorem, there exists a unique (u, σ)∈H1(Ω) ×H1 σ(Ω) ,→L4(Ω) ×L4(Ω) solution of (1.120). 40
0 0.5 1 1.5 2 2.5 3 Time ×10-3 -180 -160 -140 -120 -100 -80 -60 -40 -20 0 20 Minimum of u h US LC Figure 1.4: Minimum values of uh, with h=1 150. 2. Let us now prove that all possible fixed points of λR (with λ∈(0,1]) are bounded. In fact, observe that if (u, σ) is a fixed point of λR, then (u, σ) satisfies 1 λea(u, ¯u) = 1 k(un−1,¯u)−(u+σ,∇¯u),∀¯u∈H1(Ω), 1 λeb(σ,¯ σ) = 1 k(σn−1,¯ σ)−(u2,∇· ¯ σ),∀¯ σ∈H1 σ(Ω), (1.121) (because λR(u, σ) = (u, σ) implies R(u, σ) = (1 λu, 1 λσ)). Proceeding as in Part A of the proof of Theorem 1.4.3, we can prove that if (u, σ) is a solution of (1.121), then u≥0, which implies that u=u+. Then, multiplying (1.121)1and (1.121)2by λ, testing by ¯u=uand ¯ σ=1 2σand taking into account that λ∈(0,1], we obtain 1 4k(u, σ)k2 0+k 2k(∇u, σ)k2 L2×H1≤Cλ2k(un−1,σn−1)k2 0≡C(un−1,σn−1). Thus, we deduce that k(u, σ)kL4≤Ck(u, σ)k1≤C, where the constant Cdepends on data (Ω, un−1,σn−1). 3. We prove that Ris continuous. Let {(eul,e σl)}l∈N⊂L4(Ω) ×L4(Ω) be a sequence such that (eul,e σl)→(eu, e σ) in L4(Ω) ×L4(Ω),as l→+∞.(1.122) 41
Chapter 2 Energy-stable fully discrete approximation for a chemo-repulsion model with quadratic production 2.1 Introduction The aim of this paper is to study an unconditional energy-stable fully discrete scheme for the following parabolic-parabolic repulsive-productive chemotaxis model (with quadratic production term): ∂tu−∆u=∇·(u∇v) in Ω, t > 0, ∂tv−∆v+v=u2in Ω, t > 0, ∂u ∂n=∂v ∂n= 0 on ∂Ω, t > 0, u(x,0) = u0(x)≥0, v(x,0) = v0(x)≥0 in Ω, (2.1) where Ω is a n−dimensional open bounded domain, n= 2,3, with boundary ∂Ω. The unknowns for this model are u(x, t)≥0, the cell density, and v(x, t)≥0, the chemical concentration. Problem (2.1) is conservative in u, because the total mass RΩu(t) remains constant in time, as we can check integrating equation (2.1)1in Ω, d dt ZΩ u= 0,i.e. ZΩ u(t) = ZΩ u0:= m0,∀t > 0. In Chapter 1was proved that problem (2.1) is well-posed, because there exists global in time “weak-strong” solutions in the following sense: u≥0 and v≥0 a.e. (t, x)∈ 48
(0,+∞)×Ω, (u, v)∈L∞(0,+∞;L2(Ω) ×H1(Ω)) ∩L2(0, T;H1(Ω) ×H2(Ω)),∀T > 0, (∂tu, ∂tv)∈Lq0(0, T;H1(Ω)0×L2(Ω)),∀T > 0, (∇u, ∇v)∈L2(0,+∞;L2(Ω) ×H1(Ω)), (2.2) where q0= 2 in the 2-dimensional case (2D) and q0= 4/3 in the 3-dimensional case (3D) (q0 is the conjugate exponent of q= 2 in 2D and q= 4 in 3D), satisfying the u-equation (2.1)1in a variational sense, the v-equation (2.1)2a.e. (t, x)∈(0,+∞)×Ω, and the following energy inequality for a.e. t0, t1:t1≥t0≥0: E(u(t1), v(t1))−E(u(t0), v(t0))+Zt1 t0 (k∇u(s)k2 L2+1 2k∆v(s)k2 L2+1 2k∇v(s)k2 L2)ds ≤0,(2.3) where E(u, v) = 1 2kuk2 L2+1 4k∇vk2 L2. Moreover, assuming that the following regularity criteria is satisfied: (u, ∇v)∈L∞(0,+∞;H1(Ω) ×H1(Ω)),(2.4) (which, at least in 2D domains, is always true), it was proved that there exists a unique global in time strong solution of (2.1) satisfying (u, v)∈L∞(0,+∞;H2(Ω)2)∩L2(0, T;H3(Ω)2), (∂tu, ∂tv)∈L∞(0,+∞;L2(Ω) ×H1(Ω)) ∩L2(0,+∞;H1(Ω) ×H2(Ω)), (∂ttu, ∂ttv)∈L2(0,+∞;H1(Ω)0×L2(Ω)). (2.5) In particular, (2.5)1implies k(u, v)kL∞(0,+∞;L∞×L∞)≤C. (2.6) Therefore, it is desired to design numerical methods for the model (2.1) conserving at the discrete level the main properties of the continuous model, such as mass-conservation, energy-stability, positivity and regularity. There are only a few works about numerical analysis for chemotaxis models. For instance, for the Keller-Segel system (i.e. with chemo-attraction and linear production), Filbet studied in [4] the existence of discrete solutions and the convergence of a finite volume scheme. Saito, in [8,9], proved error estimates for a conservative Finite Element (FE) approximation. A mixed FE approximation is studied in [6]. In [3], some error estimates are proved for a fully discrete discontinuous FE method. However, as far as we know, there are not works studying 49
FE schemes satisfying the property of energy-stability related to the energy inequality (2.3). In this paper, we propose an unconditional energy-stable fully discrete scheme, which inherit some properties from the continuous model, such as mass-conservation, and weak and strong estimates analogues to (2.2) and (2.5)-(2.6). Moreover, with respect to the nonnegativity of the discrete cell and chemical variables, un hand vn h, we can deduce that vn h≥0 (see Remark 2.3.2), but the cell density nonnegativity un h≥0 can not be assured. In order to design the scheme, we follow the ideas presented in Chapter 1, and we reformulate (2.1) introducing a new variable σ=∇vinstead of v. Then, model (2.1) is rewritten as: ∂tu−∇·(∇u) = ∇·(uσ) in Ω, t > 0, ∂tσ−∇(∇·σ) + σ+ rot(rot σ) = ∇(u2) in Ω, t > 0, ∂u ∂n= 0 on ∂Ω, t > 0, σ·n= 0,[rot σ×n]tang = 0 on ∂Ω, t > 0, u(x,0) = u0(x)≥0,σ(x,0) = ∇v0(x) in Ω, (2.7) where (2.7)2has been obtained applying the gradient operator to equation (2.1)2and adding the term rot(rot σ) using the fact that rot σ= rot(∇v) = 0. Once system (2.7) is solved, we can recover vfrom u2by solving ∂tv−∆v+v=u2in Ω, t > 0, ∂v ∂n= 0 on ∂Ω, t > 0, v(x,0) = v0(x)>0 in Ω. (2.8) This chapter is organized as follows: In Section 2.2, we give the notation and some preliminary results that will be used along this paper. In Section 2.3, we study the FE Backward Euler scheme corresponding to formulation (2.7)-(2.8), including mass-conservation, unconditional energy-stability, solvability, weak and strong estimates, convergence towards weak solutions, and optimal error estimates of the scheme. In Section 2.4, we propose two different linear iterative methods in order to approach the nonlinear scheme proposed in Section 2.3, which are an energy-stable Picard’s method and the Newton’s method. We prove the solvability and the convergence of these methods to the nonlinear scheme. Finally, in Section 2.5, we present some numerical results in agreement with the theoretical analysis about the error estimates. 2.2 Notations and preliminary results We recall some functional spaces which will be used throughout this paper. We will consider the usual Sobolev spaces Hm(Ω) and Lebesgue spaces Lp(Ω),1≤p≤ ∞,with norms k·km 50
and k·kLp, respectively. In particular, the L2(Ω)-norm will be denoted by k·k0. We denote by H1 σ(Ω) := {u∈H1(Ω) : u·n= 0 on ∂Ω}and we will use the following equivalent norms in H1(Ω) and H1 σ(Ω), respectively (see [7] and [1, Corollary 3.5], respectively): kuk2 1=k∇uk2 0+ZΩ u2 ,∀u∈H1(Ω), kσk2 1=kσk2 0+krot σk2 0+k∇·σk2 0,∀σ∈H1 σ(Ω). If Zis a general Banach space, its topological dual will be denoted by Z0. Moreover, the letters C, K will denote different positive constants (independent of discrete parameters) which may change from line to line (or even within the same line). We define the linear elliptic operators Aw =g⇔ −∆w+w=gin Ω, ∂w ∂n= 0 on ∂Ω,(2.9) and Bσ=h⇔(−∇(∇·σ) + rot (rot σ) + σ=hin Ω, σ·n= 0,[rot σ×n]tang = 0 on ∂Ω,(2.10) which, in variational form, are given by A:H1(Ω) →H1(Ω)0and B:H1 σ(Ω) →H1 σ(Ω)0 such that hAw, ¯wi= (∇w, ∇¯w)+(w, ¯w),∀w, ¯w∈H1(Ω), hBσ,¯ σi= (σ,¯ σ)+(∇·σ,∇· ¯ σ) + (rot σ,rot ¯ σ),∀σ,¯ σ∈H1 σ(Ω). We assume the H2-regularity of problems (2.9) and (2.10). Consequently, we have the existence of some constants C > 0 such that kwk2≤CkAwk0,∀w∈H2(Ω),and kσk2≤CkBσk0,∀σ∈H2(Ω).(2.11) Along this paper, we will use repeatedly the classical 3D interpolation inequality kukL3≤Ckuk1/2 0kuk1/2 1∀u∈H1(Ω).(2.12) Finally, we will use the following results (see [5] and [10]): 51
Lemma 2.2.1 Assume that δ, β, k > 0and dn≥0satisfy dn+1 −dn k+δdn+1 ≤β, ∀n≥0. Then, for any n0≥0, dn≤(1 + δk)−(n−n0)dn0+δ−1β, ∀n≥n0. Lemma 2.2.2 (Uniform discrete Gronwall lemma) Let k > 0and dn, gn, hn≥0such that dn+1 −dn k≤gndn+hn,∀n≥0. If for any r∈N, there exist a1(tr),a2(tr)and a3(tr)depending on tr=kr, such that k n0+r−1 X n=n0 gn≤a1(tr), k n0+r−1 X n=n0 hn≤a2(tr), k n0+r−1 X n=n0 dn≤a3(tr), for any integer n0≥0, then dn≤a2(tr) + a3(tr) trexp {a1(tr)},∀n≥r. As consequence of Lemma 2.2.2 and Discrete Gronwall Lemma, we have the following result (see Corollary 1.2.4): Corollary 2.2.3 Under hypothesis of Lemma 2.2.2. Let k0>0be fixed, then the following estimate holds for all k≤k0 dn≤C(d0, k0)∀n≥0. 2.3 Fully Discrete Backward Euler Scheme in variables (u, σ) This section is devoted to design an unconditionally energy-stable scheme for model (2.1) (for a modified energy in variables (u, σ)), using a FE discretization in space and the backward Euler discretization in time (considered for simplicity on a uniform partition of [0, T] with time step k=T/N : (tn=nk)n=N n=0 ). Concerning the space discretization, we consider {Th}h>0be a family of shape-regular and quasi-uniform triangulations of Ω made up of simplexes K(triangles in two dimensions and tetrahedra in three dimensions), so that Ω = 52
∪K∈ThK, where h= maxK∈ThhK, with hKbeing the diameter of K. Further, let Nh={ai}i∈I denote the set of all the nodes of Th. We choose the following continuous FE spaces for u, σand v: (Uh,Σh, Vh)⊂H1×H1 σ×W1,6,generated by Pk,Pm,Prwith k, m, r ≥1. Now, let Ah:Uh→Uh,Bh:Σh→Σhand e Ah:Vh→Vhbe the linear operators defined, respectively, as follows: (Ahuh,¯uh) = (∇uh,∇¯uh)+(uh,¯uh),∀¯uh∈Uh, (Bhσh,¯ σh)=(∇·σh,∇· ¯ σh) + (rot σh,rot ¯ σh)+(σh,¯ σh),∀¯ σh∈Σh, (e Ahvh,¯vh)=(∇vh,∇¯vh)+(vh,¯vh),∀¯vh∈Vh. (2.13) Moreover, we choose the following interpolation operators: Ru h:H1(Ω) →Uh,Rσ h:H1 σ(Ω) →Σh,Rv h:H1(Ω) →Vh such that for all u∈H1(Ω), σ∈H1 σ(Ω) and v∈H1(Ω), Ru hu∈Uh,Rσ hσ∈Σhand Rv hv∈Vh satisfy (∇(Ru hu−u),∇¯uh)+(Ru hu−u, ¯uh) = 0,∀¯uh∈Uh,(2.14) (∇·(Rσ hσ−σ),∇· ¯ σh) + (rot(Rσ hσ−σ),rot ¯ σh) + (Rσ hσ−σ,¯ σh) = 0,∀¯ σh∈Σh,(2.15) (∇(Rv hv−v),∇¯vh)+(Rv hv−v, ¯vh) = 0,∀¯vh∈Vh,(2.16) respectively. Observe that, from Lax-Milgram Theorem, the interpolation operators Ru h,Rσ h and Rv hare well defined. Moreover, the following interpolation errors hold 1 hkRu hu−uk0+kRu hu−uk1≤Chkkukk+1 ∀u∈Hk+1(Ω),(2.17) 1 hkRσ hσ−σk0+kRσ hσ−σk1≤Chmkσkm+1 ∀σ∈Hm+1(Ω),(2.18) 1 hkRv hv−vk0+kRv hv−vk1≤Chrkvkr+1 ∀v∈Hr+1(Ω).(2.19) Also, the following stability properties will be used k(Ru hu, Rσ hσ,Rv hv)kW1,6≤Ck(u, σ, v)k2,(2.20) which can be obtained from (2.17)-(2.19), using the inverse inequality k(uh,σh, vh)kW1,6≤Ch−1k(uh,σh, vh)k1for all (uh,σh, vh)∈Uh×Σh×Vh, and comparing Ru,σ,v hwith an average interpolation of Clement or Scott-Zhang type (which is stable in W1,6-norm). 53
Lemma 2.3.1 Assume the H2-regularity for problems (2.9)-(2.10) given in (2.11). Then, the following estimates hold kuhkW1,6≤CkAhuhk0∀uh∈Uh,kvhkW1,6≤Cke Ahvhk0∀vh∈Vh,(2.21) kσhkW1,6≤CkBhσhk0∀σh∈Σh.(2.22) Proof. First, we consider regular functions associated to the discrete functions Ahuh, e Ahvhand Bhσh. We define u(h), v(h)∈H2(Ω) and σ(h)∈H2(Ω) as the solutions of problems −∆u(h) + u(h) = Ahuhin Ω, ∂u(h) ∂n= 0 on ∂Ω,(2.23) −∆v(h) + v(h) = e Ahvhin Ω, ∂v(h) ∂n= 0 on ∂Ω,(2.24) and −∇(∇·σ(h)) + rot (rot σ(h)) + σ(h) = Bhσhin Ω, σ(h)·n= 0,[rot σ(h)×n]tang = 0 on ∂Ω.(2.25) In particular, from (2.11), ku(h)k2≤CkAhunk0,kv(h)k2≤Cke Ahvhk0and kσ(h)k2≤CkBhσhk0.(2.26) Now, we decompose the W1,6-norm as: kuhkW1,6≤ kuh−Ru hu(h)kW1,6+kRu hu(h)−u(h)kW1,6+ku(h)kW1,6:= I1+I2+I3,(2.27) kvhkW1,6≤ kvh−Rv hv(h)kW1,6+kRv hv(h)−v(h)kW1,6+kv(h)kW1,6:= H1+H2+H3,(2.28) kσhkW1,6≤ kσh−Rσ hσ(h)kW1,6+kRσ hσ(h)−σ(h)kW1,6+kσ(h)kW1,6:= J1+J2+J3.(2.29) In order to bound Ji(i= 1,2), we test (2.25)1by ¯ σh∈Σhand using (2.13)2we have (∇·σh,∇· ¯ σh) +(rot σh,rot ¯ σh)+(σh,¯ σh) = (∇·σ(h),∇· ¯ σh) + (rot σ(h),rot ¯ σh)+(σ(h),¯ σh),∀¯ σh∈Σh.(2.30) By subtracting at both sides of equality (2.30) the terms (∇·Rσ hσ(h),∇·¯ σh), (rotRσ hσ(h),rot ¯ σh) and (Rσ hσ(h),¯ σh), testing by ¯ σh=σh−Rσ hσ(h)∈Σhand using the H¨older and Young inequalities, we deduce kσh−Rσ hσ(h)k1≤CkRσ hσ(h)−σ(h)k1≤Chkσ(h)k2,(2.31) 54
where in the last inequality interpolation error (2.18) was used. Then, using in (2.31) the inverse inequality kσhkW1,6≤Ch−1kσhk1for all σh∈Σh, we conclude that for i= 1,2 Ji≤Ch−1kRσ hσ(h)−σ(h)k1≤Ckσ(h)k2.(2.32) Finally, J3=kσ(h)kW1,6≤Ckσ(h)k2.(2.33) Therefore, using (2.32)-(2.33) in (2.29), and taking into account (2.26), we deduce (2.22). Proceeding analogously for Iiand Hi(i= 1,2,3), we deduce (2.21). 2.3.1 Definition of the scheme By taking into account the reformulation (2.7), we consider the following FE Backward Euler Scheme in variables (u, σ) (Scheme US, from now on) which is a first order in time, nonlinear and coupled scheme: •Initialization: We fix (u0 h,σ0 h) = (Ru hu0,Rσ hσ0)∈Uh×Σhand v0 h=Rv hv0∈Vh. Then, RΩu0 h=RΩu0=m0. Time step n: Given (un−1 h,σn−1 h)∈Uh×Σh, compute (un h,σn h)∈Uh×Σhsolving ((δtun h,¯uh)+(∇un h,∇¯uh)+(un hσn h,∇¯uh)=0,∀¯uh∈Uh, (δtσn h,¯ σh)+(Bhσn h,¯ σh)−2(un h∇un h,¯ σh)=0,∀¯ σh∈Σh,(2.34) where, in general, we denote δtan h=an h−an−1 h k. Once the scheme US is solved, given vn−1 h∈Vh, we can recover vn h=vn h((un h)2)∈Vh solving: (δtvn h,¯vh)+(e Ahvn h,¯vh) = ((un h)2,¯vh),∀¯vh∈Vh.(2.35) Given un h∈Uhand vn−1 h∈Vh, Lax-Milgram theorem implies that there exists a unique vn h∈Vhsolution of (2.35). Remark 2.3.2 By using the mass-lumping technique in all terms of (2.35) excepting the self-diffusion term (∇vn h,∇¯vh), approximating by P1-continuous FE and imposing a condition based on a geometrical property of the triangulation, related to the fact that the interior angles of the triangles or tetrahedra must be at most π/2, we can prove that if vn−1 h≥0then vn h≥0. However, in all numerical simulations that we have made without using mass-lumping, we have not found any example in which, beginning with v0 h≥0we obtain vn h(ai)<0, for some n > 0and ai. 55
2.3.2 Solvability, Energy-Stability and Convergence Assuming that the functions ¯uh= 1 ∈Uhand ¯vh= 1 ∈Vh, we deduce that the scheme US conserves in time the total mass RΩun h, that is, ZΩ un h=ZΩ un−1 h=···=ZΩ u0 h, and we have the following behavior for RΩvn h: δtZΩ vn h=ZΩ (un h)2−ZΩ vn h. Now, we establish some results concerning to the solvability and energy-stability of scheme US, but we will omit their proofs because those follow the same ideas given in Chapter 1(Theorem 1.4.3, Lemma 1.4.6 and Theorem 1.4.8, respectively). Theorem 2.3.3 (Unconditional existence and conditional uniqueness) There exists (un h,σn h)∈Uh×Σhsolution of the scheme US. Moreover, if kk(un h,σn h)k4 1is small enough, (2.36) then the solution is unique. Remark 2.3.4 In the case of 2D domains, from estimate (2.54) below, the uniqueness restriction (2.36) can be relaxed to kK2 0small enough, where K0is a constant depending on data (Ω, u0,σ0), but independent of (k, h)and n. Remark 2.3.5 In 3D domains, using the inverse inequality kuhk1≤C hkuhk0(see Lemma 4.5.3 in [2], p. 111) and estimate (2.41) below, we have that k(un h,σn h)k4 1≤C h4k(un h,σn h)k4 0≤C h4C2 0 and therefore, the uniqueness restriction (2.36) can be rewrite as kC1 h4small enough,(2.37) where C1is a positive constant depending on data (Ω, u0,σ0), but independent of n. 56
Definition 2.3.6 A numerical scheme with solution (un,σn)is called energy-stable with respect to the energy E(u, σ) = 1 2kuk2 0+1 4kσk2 0,(2.38) if this energy is time decreasing, that is E(un h,σn h)≤E(un−1 h,σn−1 h),∀n. (2.39) Lemma 2.3.7 (Unconditional stability) The scheme US is unconditionally energy-stable with respect to E(u, σ). In fact, if (un h,σn h)is a solution of the scheme US, then the following discrete energy law holds δtE(un h,σn h)+k 2kδtun hk2 0+k 4kδtσn hk2 0+k∇un hk2 0+1 2kσn hk2 1= 0.(2.40) Remark 2.3.8 Looking at (2.40), we can say that the scheme US introduces the following two first order “numerical dissipation terms”: k 2kδtun hk2 0and k 4kδtσn hk2 0. From the (local in time) discrete energy law (2.40), we deduce the following global in time estimates for (un h,σn h) solution of scheme US: Theorem 2.3.9 (Uniform Weak estimates of scheme US) Let (un h,σn h)be a solution of scheme US. Then, the following estimates hold k(un h,σn h)k2 0+k2 n X m=1k(δtum h, δtσm h)k2 0+k n X m=1k(∇um h,σm h)k2 L2×H1≤C0,∀n≥1,(2.41) k n+n0 X m=n0+1k(um h,σm h)k2 1≤C0+C1(nk),∀n≥1,(2.42) where the integer n0≥0is arbitrary, with positive constants C0, C1depending on the data (Ω, u0,σ0), but independent of (k, h)and (n, n0). 57
2.3.4 Error estimates at finite time In this subsection, we will obtain error estimates for any solution (un h,σn h) of the scheme US and vn hof (2.35), with respect to sufficiently regular solutions (u, σ) of (2.7) and vof (2.8) respectively. In our analysis, in order to obtain optimal error estimates we need to assume that both spaces Uh,Σhare generated by Pm-continuous FE and Vhis generated by Pm+1-continuous FE, with m≥1. This is a natural assumption taking into account that the energy norm for vin the continuous model has one order greater than the energy norms for u, σ. Error estimates for scheme US We start introducing the following notations for the errors at t=tn:en u=u(tn)−un hand en σ=σ(tn)−σn h, and for the discrete in time derivative of these errors: δten u=en u−en−1 u k and δten σ=en σ−en−1 σ k. Then, subtracting (2.7) at t=tnand the scheme US, we have that (en u, en σ) satisfies (δten u,¯uh)+(∇en u,∇¯uh)+(en uσ(tn) + un hen σ,∇¯uh)=(ξn 1,¯uh),∀¯uh∈Uh,(2.71) (δten σ,¯ σh) + hBen σ,¯ σhi= 2(en u∇u(tn) + un h∇en u,¯ σh)+(ξn 2,¯ σh),∀¯ σh∈Σh,(2.72) where ξn 1, ξn 2are the consistency errors associated to the scheme US, that is, ξn 1=δt(u(tn))− ut(tn) and ξn 2=δt(σ(tn)) −σt(tn). Now, considering the interpolation operators Ru hand Rσ h defined in (2.14)-(2.15), we decompose en uand en σas follows en u= (I−Ru h)u(tn) + Ru hu(tn)−un h=en u,i +en u,h,(2.73) en σ= (I−Rσ h)σ(tn) + Rσ hσ(tn)−σn h=en σ,i +en σ,h,(2.74) where en u,i is the interpolation error and en u,h is the discrete error of u. Then, taking into account (2.14)-(2.15), from (2.71)-(2.74) we have δten u,h,¯uh+(∇en u,h,∇¯uh)+(en u,hσ(tn) + un hen σ,h,∇¯uh) = (ξn 1,¯uh) −δten u,i,¯uh−(en u,iσ(tn) + un hen σ,i,∇¯uh)+(en u,i,¯uh),∀¯uh∈Uh,(2.75) δten σ,h,¯ σh+(Bhen σ,h,¯ σh) = (ξn 2,¯ σh) + 2(en u,h∇u(tn) + un h∇en u,h,¯ σh) +2(en u,i∇u(tn) + un h∇en u,i,¯ σh)−δten σ,i,¯ σh,∀¯ σh∈Σh.(2.76) Notice that RΩen u,h = 0 (since u0 h=Ru hu0and from (2.14)RΩRu hu(tn) = RΩu(tn) = m0), hence the following norms are equivalents: k∇en u,hk0≃ ken u,hk1. 64
Theorem 2.3.21 We assume that there exists (u, σ)an exact solution of (2.7) with the following regularity: ((u, σ)∈L∞(0, T;Hm+1(Ω)×Hm+1(Ω)),(ut,σt)∈L2(0, T;Hm+1(Ω)×Hm+1(Ω)), (utt,σtt)∈L2(0, T;H1(Ω)0×H1 σ(Ω)0). (2.77) Let (un h,σn h)be a solution of the scheme US. Then, if k(k(u, σ)k4 L∞(H1)+k(u, σ)k2 L∞(H2))is small enough,(2.78) the following a priori error estimate holds k(en u,h, en σ,h)kl∞L2∩l2H1≤C(T)(k+hm+1) (2.79) where C(T) = K1Texp(K2T), with K1, K2>0independent of (k, h). Recall that u, σare approximated by Pm-continuous FE. Proof. Taking ¯uh=en u,h n (2.75), ¯ σh=1 2en σ,h in (2.76) and adding, the terms (un h∇en u,h, en σ,h) cancel, and we obtain δt1 2ken u,hk2 0+1 4ken σ,hk2 0+1 2k(en u,h, en σ,h)k2 1= (ξn 1, en u,h) + 1 2(ξn 2, en σ,h)−δten u,i, en u,h −1 2δten σ,i, en σ,h−(en u,h,σ(tn)·∇en u,h −∇u(tn)·en σ,h)−(en u,i,σ(tn)·∇en u,h −∇u(tn)·en σ,h) −(un h, en σ,i ·∇en u,h −∇en u,i ·en σ,h)+(en u,i, en u,h) = 8 X m=1 Im.(2.80) Then, using the H¨older and Young inequalities, the 3D interpolation inequality (2.12), the interpolation errors (2.17)-(2.18), the stability property (2.20) and the hypothesis (2.77), we control the terms on the right hand side of (2.80) as follows I1+I2≤εk(en u,h, en σ,h)k2 1+Cεk(ξn 1, ξn 2)k2 (H1)0×(H1 σ)0 ≤εk(en u,h, en σ,h)k2 1+Ck Ztn tn−1k(utt(t),σtt(t))k2 (H1)0×(H1 σ)0dt, (2.81) I5≤ ken u,hkL3(k∇u(tn)k0ken σ,hkL6+k∇·σ(tn)k0ken u,hkL6) ≤εk(en u,h, en σ,h)k2 1+Cεk(∇u(tn),∇·σ(tn))k4 0ken u,hk2 0,(2.82) 65
I6≤ ken u,ik0(k∇en u,hk0kσ(tn)kL∞+k∇u(tn)kL3ken σ,hkL6) ≤εk(en u,h, en σ,h)k2 1+Cεk(∇u(tn),σ(tn))k2 L3×L∞ken u,ik2 0 ≤εk(en u,h, en σ,h)k2 1+Ch2(m+1)ku(tn)k2 m+1 ≤εk(en u,h, en σ,h)k2 1+Ch2(m+1),(2.83) I7≤ |(en u,h, en σ,i ·∇en u,h −∇en u,i ·en σ,h)|+|(Ru hu(tn), en σ,i ·∇en u,h −∇en u,i ·en σ,h)| ≤εk(en u,h, en σ,h)k2 1+Cεken u,hk2 0k(en u,i, en σ,i)k2 W1,3×L∞+CεkRu hu(tn)k2 W1,3∩L∞k(en u,i, en σ,i)k2 0 ≤εk(en u,h, en σ,h)k2 1+Ck(u(tn),σ(tn))k2 2ken u,hk2 0+Ch2(m+1)k(u(tn),σ(tn))k2 m+1ku(tn)k2 2 ≤εk(en u,h, en σ,h)k2 1+Ck(u(tn),σ(tn))k2 2ken u,hk2 0+Ch2(m+1),(2.84) I8≤ ken u,ik0ken u,hk0≤εken u,hk2 1+Ch2(m+1),(2.85) I3+I4≤ k(en u,h, en σ,h)k0k((I−Ru h)δtu(tn),(I−Rσ h)δtσ(tn))k0 ≤εk(en u,h, en σ,h)k2 1+Ch2(m+1)k(δtu(tn), δtσ(tn))k2 m+1 ≤εk(en u,h, en σ,h)k2 1+Ch2(m+1) kZtn tn−1k(ut,σt)k2 m+1dt, (2.86) where in the last inequality was used that (δtu(tn), δtσ(tn)) = 1 kZtn tn−1 (ut,σt). Therefore, taking εsmall enough, from (2.80)-(2.86) we obtain δt1 2ken u,hk2 0+1 4ken σ,hk2 0+k(en u,h, en σ,h)k2 1 ≤Ck Ztn tn−1k(utt(t),σtt(t))k2 (H1)0×(H1 σ)0dt +Ch2(m+1) kZtn tn−1k(ut,σt)k2 m+1dt +Ch2(m+1) +C(k(u(tn),σ(tn))k4 1+k(u(tn),σ(tn))k2 2)ken u,hk2 0.(2.87) Then, multiplying (2.87) by k, adding from n= 1 to n=r, recalling that e0 u,h =e0 σ,h = 0, taking into account (2.77), we obtain h1 4−k Cik(er u,h, er σ,h)k2 0+k r X n=1 k(en u,h, en σ,h)k2 1≤Ck2+Ch2(m+1) +C k r−1 X n=0 ken u,hk2 0. Therefore, assuming hypothesis (2.78) and using the Discrete Gronwall Lemma, we conclude (2.79). 66
Remark 2.3.22 Under hypothesis of Theorem 2.3.21, one has in particular k(un h,σn h)k2 1≤C+C(T)k+h2(m+1) k. Therefore, under the hypothesis h2(m+1) k≤C, (2.88) we have the estimate k(un h,σn h)k2 1≤C, (2.89) hence the hypothesis (2.36) providing uniqueness of the scheme is equivalent to ksmall enough. Finally, since for any choice of (k, h)either (2.37) (see Remark 2.3.5) or (2.88) holds, one has the uniqueness of (un h,σn h)solution of (2.34) only imposing ksmall enough. Error estimates for vn hsolution of (2.35) We introduce the following notation for the errors in t=tn:en v=v(tn)−vn h, and for the discrete in time derivative of this error: δten v=en v−en−1 v k. Then, subtracting (2.8) at t=tn and (2.35), we have that en vsatisfies (δten v,¯vh) + hAen v,¯vhi= ((u(tn) + un h)en u,¯vh)+(ξn 3,¯vh),∀¯vh∈Vh,(2.90) where ξn 3is the consistency error associated to (2.35), that is, ξn 3=δt(v(tn)) −vt(tn). Now, considering the interpolation operator Rv hdefined in (2.16), we decompose en vas follows en v= (I−Rv h)v(tn) + Rv hv(tn)−vn h=en v,i +en v,h,(2.91) where en v,i is the interpolation error and en v,h is the discrete error of v. Then, taking into account (2.16), from (2.90)-(2.91) we have δten v,h,¯vh+( e Ahen v,h,¯vh) = (ξn 3,¯vh) + ((u(tn) + un h)(en u,h+en u,i),¯vh)−(δten v,i,¯vh),∀¯vh∈Vh.(2.92) Theorem 2.3.23 Under hypothesis of Theorem 2.3.21. Let vn hbe the solution of (2.35), and assume the following regularity for vexact solution of (2.8): (vt, vtt)∈L2(0, T;Hm+2(Ω) ×H1(Ω)0).(2.93) Then, the a priori error estimate holds ken v,hkl∞L2∩l2H1≤C(T)(k+hm+1),(2.94) where C(T) = K1Texp(K2T), with K1, K2>0independent of (k, h). 67
Proof. Taking ¯vh=en v,h in (2.92) and using the H¨older and Young inequalities, we obtain δt1 2ken v,hk2 0+k 2kδten v,hk2 0+1 2ken v,hk2 1≤Ckξn 3k2 (H1)0 +Cku(tn) + un hk2 L3(ken u,hk2 0+ken u,ik2 0) + Ck(I−Rv h)δtv(tn)k2 0.(2.95) Using (2.17), (2.19) and proceeding as in (2.81) and (2.86), we bound the terms on the right hand side in (2.95) and we deduce δtken v,hk2 0+ken v,hk2 1≤Ck Ztn tn−1kvtt(t)k2 (H1)0dt +Cku(tn) + un hk2 L3ken u,hk2 0 +Cku(tn) + un hk2 L3h2(m+1)ku(tn)k2 m+1 +Ch2(m+2) kZtn tn−1kvtk2 m+2dt. (2.96) Then, multiplying (2.96) by k, adding from n= 1 to n=r, we obtain (recall e0 v,h = 0): ker v,hk2 0+k r X n=1ken v,hk2 1≤Ck2Ztr 0kvtt(t)k2 (H1)0dt +Cken u,hk2 l∞L2k r X n=1ku(tn) + un hk2 L3 +Ch2(m+1)k r X n=1ku(tn) + un hk2 L3+Ch2(m+2) Ztr 0kvtk2 m+2dt. Then, using (2.42), (2.77), (2.93) and (2.79), we conclude (2.94). Theorem 2.3.24 Under hypothesis of Theorem 2.3.23, but assuming the regularity: vtt ∈L2(0, T;L2(Ω)),(2.97) the a priori error estimate ken v,hkl∞H1∩l2W1,6≤C(T)(k+hm+1) (2.98) holds, where C(T) = K1Texp(K2T), with K1, K2>0independent of (k, h). Proof. Taking ¯vh=e Ahen v,h in (2.92) and using the H¨older and Young inequalities, we obtain δt1 2ken v,hk2 1+k 2kδten v,hk2 1+1 2ke Ahen v,hk2 0≤Ckξn 3k2 0+Cku(tn) + un hk2 L3ken u,hk2 L6 +Ck(u(tn) + un h)en u,ik2 0+Ck(I−Rv h)δtv(tn)k2 0.(2.99) 68
Using the H¨older inequality, the interpolation error (2.17), the stability property (2.20) and the hypothesis (2.77), we have k(u(tn) + un h)en u,ik2 0≤Cku(tn) + Ru hu(tn)k2 L∞ken u,ik2 0+Cken u,hk2 L6ken u,ik2 L3 ≤Ch2(m+1) +Cken u,hk2 L6.(2.100) Therefore, from (2.99), proceeding as in (2.81) and (2.86) and using (2.100), we deduce δtken v,hk2 1+ke Ahen v,hk2 0≤Ck Ztn tn−1kvtt(t)k2 0dt +(Cku(tn) + un hk2 L3+C)ken u,hk2 L6+Ch2(m+1) +Ch2(m+2) kZtn tn−1kvtk2 m+2dt. Now, in order to bound the term ku(tn) + un hk2 L3, we split the argument into two cases: 1. Estimates assuming h << f(k)(hsmall enough with respect to k): From (2.79) we have that k r P n=1ken u,hk2 1≤C(T)(k2+h2(m+1)), which implies that ken u,hk1≤C(T)(k1/2+hm+1 k1/2).(2.101) Moreover, using the interpolation inequality (2.12), (2.79), (2.20), (2.77) and (2.101), we obtain ku(tn) + un hk2 L3≤Cku(tn)k2 L3+CkRu hu(tn)k2 L3+Cken u,hk2 L3≤C+Cken u,hk0ken u,hk1 ≤C+C(T)(k+hm+1)(k1/2+hm+1 k1/2)≤C(2.102) under he hypothesis h2(m+1) k1/2≤C. (2.103) 2. Estimates assuming k << g(k)(ksmall enough with respect to h): Using the inverse inequality kuhkL3≤C h1/2kuhk0for all uh∈Uh, (2.20), (2.77) and (2.79), we have that ku(tn) + un hk2 L3≤Cku(tn)k2 L3+CkRu hu(tn)k2 L3+Cken u,hk2 L3 ≤C hken u,hk2 0+C≤C(T) h(k2+h2(m+1)) + C≤C under the hypothesis k2 h< C. (2.104) 69
Therefore, since for any choice of (k, h) either (2.103) or (2.104) holds, we arrive at δtken v,hk2 1+ke Ahen v,hk2 0≤Ck Ztn tn−1kvtt(t)k2 0dt +Cken u,hk2 L6+Ch2(m+1) +Ch2(m+2) kZtn tn−1kvtk2 m+2dt. (2.105) Multiplying (2.105) by k, adding from n= 1 to n=r, recalling that e0 v,h = 0 and using (2.77), (2.93), (2.97) and (2.79), we conclude (2.98). 2.4 Linear iterative methods to approach the Backward Euler scheme In this section, we propose two different linear iterative methods to approach the Backward Euler scheme US, which are an energy-stable Picard’s method and the Newton’s method. We prove the solvability and the convergence of these methods to the nonlinear scheme. 2.4.1 Picard Method In order to approximate the solution (un h,σn h) of the nonlinear scheme US, we consider the following Picard method: Let (un−1 h,σn−1 h)∈Uh×Σhbe fixed. Given ul−1 h∈Uh(assuming u0 h=un−1 hat the first iteration step), find (ul h,σl h)∈Uh×Σhsolving the linear coupled problem: 1 k(ul h,¯uh)+(∇ul h,∇¯uh)+(ul−1 hσl h,∇¯uh) = 1 k(un−1 h,¯uh),∀¯uh∈Uh, 1 k(σl h,¯ σh)+(Bhσl h,¯ σh)−2(ul−1 h∇ul h,¯ σh) = 1 k(σn−1 h,¯ σh),∀¯ σh∈Σh, (2.106) until the stopping criteria max kul h−ul−1 hk0 kul−1 hk0 ,kσl h−σl−1 hk0 kσl−1 hk0≤tol (with tol > 0 being a tolerance parameter). Theorem 2.4.1 (Unconditional Unique Solvability) There exists a unique (ul h,σl h) solution of (2.106). 70
Proof. Since (2.106) can be rewritten as a square linear algebraic system, it is sufficient to prove uniqueness. Suppose that there exist (ul h,1,σl h,1),(ul h,2,σl h,2)∈Uh×Σhtwo possible solutions of (2.106). Then defining ul h=ul h,1−ul h,2and σl h=σl h,1−σl h,2, we have that (ul h,σl h)∈Uh×Σhsatisfies 1 k(ul h,¯uh)+(∇ul h,∇¯uh)+(ul−1 hσl h,∇¯uh) = 0,∀¯uh∈Uh,(2.107) 1 k(σl h,¯ σh)+(Bhσl h,¯ σh)−2(ul−1 h∇ul h,¯ σh) = 0,∀¯ σh∈Σh.(2.108) Taking ¯uh=ul hand ¯ σh=1 2σl hin (2.107)-(2.108) and adding, the terms (ul−1 h∇ul h,σl h) cancel, and we obtain 1 2kk(ul h,σl h)k2 0+1 2k(∇ul h,σl h)k2 L2×H1≤0, and thus we conclude that k(ul h,σl h)k1= 0, which implies ul h,1=ul h,2and σl h,1=σl h,2. Theorem 2.4.2 (Local uniqueness of solution of scheme US and Convergence of Picard’s method) Given (un−1 h,σn−1 h), there exists r > 0(large enough) such that if kk(un−1 h,σn−1 h)k4 1and kr4are small enough,(2.109) then the scheme US has a unique solution (un h,σn h)in Br((un−1 h,σn−1 h)) := {(u, σ)∈Uh× Σh:k(u−un−1 h,σ−σn−1 h)k1≤r}. Moreover, the sequence of solutions {ul h,σl h}l≥0of the iterative algorithm (2.106) (assuming (u0 h,σ0 h) = (un−1 h,σn−1 h)at the first iteration step), converges to (un h,σn h)strongly in H1(Ω). Proof. We consider the operator R:Uh→Uh, given by R(eu) = u, where (u, σ) satisfies (2.106) with ul−1 h=euand (ul h,σl h)=(u, σ), that is, 1 k(u, ¯uh)+(∇u, ∇¯uh)+(euσ,∇¯uh) = 1 k(un−1 h,¯uh),∀¯uh∈Uh,(2.110) 1 k(σ,¯ σh)+(Bhσ,¯ σh)−2(eu∇u, ¯ σh) = 1 k(σn−1 h,¯ σh),∀¯ σh∈Σh.(2.111) Observe that from Theorem 2.4.1, we have that for any eu∈Uhthere exists a unique (u, σ)∈Uh×Σhsolution of (2.110)-(2.111). Thus, Ris well defined. Now, before to prove that Ris contractive, we will construct a ball Br(un−1 h) = {u∈Uh:ku−un−1 hk1≤r} ⊂ Uh 71
such that R(Br(un−1 h)) ⊆Br(un−1 h). In order to define r, we consider w=u−un−1 hand τ=σ−σn−1 h. Then, from (2.110)-(2.111) we have that (w, τ) verifies 1 k(w, ¯uh)+(∇w, ∇¯uh) = −(euτ,∇¯uh)−(∇un−1 h,∇¯uh)−(euσn−1 h,∇¯uh),∀¯uh∈Uh,(2.112) 1 k(τ,¯ σh)+(Bhτ,¯ σh) = 2(eu∇w, ¯ σh)−(Bhσn−1 h,¯ σh) + 2(eu∇un−1 h,¯ σh),∀¯ σh∈Σh.(2.113) Testing by ¯uh=wand ¯ σh=1 2τin (2.112)-(2.113) and adding, the terms (eu∇w, τ) cancel, and using the fact that ZΩ w= 0 as well as the 3D interpolation inequality (2.12), we obtain 1 2kk(w, τ)k2 0+1 2k(w, τ)k2 1≤1 8k(w, τ)k2 1+Ck(un−1 h,σn−1 h)k2 1 +1 8keu−un−1 hk2 1+1 8kun−1 hk2 1+1 8k(w, τ)k2 1+Ck(un−1 h,σn−1 h)k4 1k(w, τ)k2 0.(2.114) Therefore, from (2.114) we deduce 1 2k−Ck(un−1 h,σn−1 h)k4 1k(w, τ)k2 0+1 4k(w, τ)k2 1≤Ck(un−1 h,σn−1 h)k2 1+1 8keu−un−1 hk2 1.(2.115) Thus, if k < 1 2Ck(un−1 h,σn−1 h)k4 1 , from (2.115) we conclude k(w, τ)k2 1≤Ck(un−1 h,σn−1 h)k2 1+1 2keu−un−1 hk2 1.(2.116) Then, choosing r > 0 large enough such that Ck(un−1 h,σn−1 h)k2 1≤1 2r2,(2.117) from (2.116) we deduce that R(Br(un−1 h)) ⊆Br(un−1 h). Then, we take the restriction of R to Br(un−1 h), that is, Rr:Br(un−1 h)→Br(un−1 h). Let’s prove that Rris contractive. Let eu1,eu2∈Br(un−1 h), and (u1,σ1) and (u2,σ2) solutions of (2.110)-(2.111) corresponding to eu1 and eu2respectively (i.e., Rr(eu1) = u1and Rr(eu2) = u2). Then, from (2.110)-(2.111) we have that (u1−u2,σ1−σ2)∈Uh×Σhsatisfies 1 k(u1−u2,¯uh)+(∇(u1−u2),∇¯uh)+(eu1(σ1−σ2),∇¯uh)+((eu1−eu2)σ2,∇¯uh)=0,∀¯uh∈Uh, 72
1 k(σ1−σ2,¯ σh)+(Bh(σ1−σ2),¯ σh)−2(eu1∇(u1−u2),¯ σh)−2((eu1−eu2)∇u2,¯ σh) = 0,∀¯ σh∈Σh. Testing by ¯uh=u1−u2,¯ σh=1 2(σ1−σ2) and adding, the terms (eu1(σ1−σ2),∇(u1−u2)) cancel, and using the H¨older and Young inequalities, the 3D interpolation inequality (2.12) and taking into account that RΩu1−u2= 0, we obtain 1 2kk(u1−u2,σ1−σ2)k2 0+ku1−u2k2 1+1 2kσ1−σ2k2 1 ≤Ckeu1−eu2k1(kσ2k1ku1−u2kL3+ku2k1kσ1−σ2kL3) ≤1 4keu1−eu2k2 1+1 2ku1−u2k2 1+1 4kσ1−σ2k2 1+Ck(u1−u2,σ1−σ2)k2 0k(u2,σ2)k4 1 and thus, we deduce that 1 kk(u1−u2,σ1−σ2)k2 0+ku1−u2k2 1+1 2kσ1−σ2k2 1 ≤1 2keu1−eu2k2 1+Ck(u1−u2,σ1−σ2)k2 0k(u2,σ2)k4 1.(2.118) Therefore, since from (2.116) and (2.117) we have k(u2,σ2)k4 1≤C(r4+k(un−1 h,σn−1 h)k4 1), if 1 2k> Cr4and 1 2k> Ck(un−1 h,σn−1 h)k4 1), from (2.118) we have kRr(eu1)−Rr(eu2)k2 1≤1 2keu1−eu2k2 1, which implies that Rris contractive. Then, as a consequence of the Banach fixed point theorem, we conclude that there exists a unique fixed point of Rr,Rr(u) = u. Thus, (u, σ) is the unique solution of the scheme US in Br(un−1 h). Additionally, we conclude that the sequence of solutions {ul h,σl h}l≥0of the iterative algorithm (2.106), where (u0 h,σ0 h) = (un−1 h,σn−1 h), converges to the solution (un h,σn h). Remark 2.4.3 In the case of 2D Domains, from estimate (2.54), the restriction (2.109)1 can be relaxed to k≤K0, where K0is a constant depending on data (Ω, u0,σ0), but independent of (k, h)and n. Remark 2.4.4 We have that the restriction (2.109)1is equivalent to (2.36). Therefore, under hypothesis of Theorem 2.3.21 and arguing as in Remark 2.3.22, the conclusion of Theorem 2.4.2 remains true only assuming ksmall enough. 73
Bibliography [1] C. Amrouche and N.E.H. Seloula, Lp-theory for vector potentials and Sobolev’s inequalities for vector fields: application to the Stokes equations with pressure boundary conditions. Math. Models Methods Appl. Sci. 23 (2013), no. 1, 37–92. [2] S.C. Brenner and L.R. Scott, The mathematical theory of finite element methods, Third edition, Texts in Applied Mathematics, 15. Springer, New York (2008). [3] Y. Epshteyn and A. Izmirlioglu, Fully discrete analysis of a discontinuous finite element method for the Keller-Segel chemotaxis model. J. Sci. Comput. 40 (2009), no. 1-3, 211–256. [4] F. Filbet, A finite volume scheme for the Patlak-Keller-Segel chemotaxis model. Numer. Math. 104 (2006), no. 4, 457–488. [5] Y. He and K. Li, Asymptotic behavior and time discretization analysis for the nonstationary Navier-Stokes problem. Numer. Math. 98 (2004), no. 4, 647–673. [6] A. Marrocco, Numerical simulation of chemotactic bacteria aggregation via mixed finite elements. M2AN Math. Model. Numer. Anal. 37 (2003), no. 4, 617–630. [7] J. Necas, Les M´ethodes Directes en Th´eorie des Equations Elliptiques. Editeurs Academia, Prague (1967). [8] N. Saito, Conservative upwind finite-element method for a simplified Keller-Segel system modelling chemotaxis. IMA J. Numer. Anal. 27 (2007), no. 2, 332–365. [9] N. Saito, Error analysis of a conservative finite-element approximation for the Keller- Segel system of chemotaxis. Commun. Pure Appl. Anal. 11 (2012), no. 1, 339–364. [10] J. Shen, Long time stability and convergence for fully discrete nonlinear Galerkin methods. Appl. Anal. 38 (1990), 201–229. 80
Chapter 3 Asymptotic behaviour for a chemo-repulsion system with quadratic production: The continuous problem and fully discrete numerical schemes 3.1 Introduction The directed movement of cells in response to a chemical stimulus is known in biology as chemotaxis. More specifically, if the cells move towards regions of high chemical concentration, the motion is called chemoattraction, while if the cells move towards regions of lower chemical concentration, the motion is called chemorepulsion. Models for chemoattraction and chemorepulsion motion has been studied in literature (see [4,9,7,10] and references therein). One of the most important characteristics of chemoattractant models is that the finite blow up of solutions can happen in space dimension greater or equal to 2; while in chemorepulsion models this phenomenon is not expected. Many works have been devoted to study in what cases and how this phenomenon takes place. In those cases in which blow-up phenomenon does not happen, it is interesting to study the asymptotic behaviour of the solutions of the model. In fact, in [14], Osaki and Yagi studied the convergence of the solution of the Keller-Segel model to a stationary solution in the one-dimensional case. In [8], the convergence of the solution of the Keller-Segel model with an additional term of cross-diffusion to a steady state was shown. In [4] the authors proved the convergence to constant state for a chemorepulsion model with linear production. Therefore, taking into account the results above, the aim of this paper is to 81
study the asymptotic behaviour of the following parabolic-parabolic repulsive-productive chemotaxis model (with quadratic production term): ∂tu−∆u=∇·(u∇v) in Ω, t > 0, ∂tv−∆v+v=u2in Ω, t > 0, ∂u ∂n=∂v ∂n= 0 on ∂Ω, t > 0, u(x,0) = u0(x)≥0, v(x,0) = v0(x)≥0 in Ω, (3.1) where Ω is a n−dimensional open bounded domain, n= 2,3, with boundary ∂Ω; and the unknowns are u(x, t)≥0, the cell density, and v(x, t)≥0, the chemical concentration. This model has been studied in Chapter 1, and it was proved that model (3.1) is well-posed: there exists global in time weak-strong solution in the sense of Definition 3.2.1 (below), and, for 2Ddomains, there exists a unique global in time strong solution. On the other hand, another interesting topic is to study the asymptotic behaviour of fully discrete numerical schemes approximating (3.1). In fact, in [5] Guill´en-Gonz´alez and Samsidy studied and proved asymptotic convergence for a fully discrete finite element scheme for a Ginzburg-Landau model for nematic liquid crystal flow. In [12] Merlet and Pierre studied the asymptotic behaviour of the Backward Euler scheme applied to gradient flows. It is important to notice that, in chemotaxis models, there are few works studying large-time behaviour for fully discrete schemes. We refer to [2], where the authors shown convergence at infinite time of a finite volume scheme for a Keller-Segel model with an additional term of cross-diffusion. Meanwhile, the behavior at infinite time of a fully discrete scheme for model (3.1) seem to be still an open problem. For this reason, in this paper we also study the large-time behavior for two fully discrete numerical schemes associated to model (3.1). This chapter is organized as follows: In Section 3.2, we study the asymptotic behavior of the global weak-strong solutions for the model (3.1), and we prove the exponential convergence as time goes to infinity to a constant state. In Section 3.3, we analyze this same behavior for two fully discrete numerical schemes associated to system (3.1): the nonlinear backward Euler in the variables (u, v), and the nonlinear scheme defined in Chapter 2by introducing the auxiliary variable σ=∇v. Moreover, in order to analyze the asymptotic behaviour for the backward Euler scheme, we study its solvability and unconditional energystability. Finally, in Section 3.4, we compare the numerical schemes throughout several numerical simulations. 3.1.1 Notation We recall some functional spaces which will be used throughout this paper. We will consider the usual Sobolev spaces Hm(Ω) and Lebesgue spaces Lp(Ω),1≤p≤ ∞,with norms k·km 82
and k·kLp, respectively. In particular, the L2(Ω)-norm will be denoted by k·k0. We denote by H1 σ(Ω) := {u∈H1(Ω) : u·n= 0 on ∂Ω}and we will use the following equivalent norms in H1(Ω) and H1 σ(Ω), respectively (see [13] and [1, Corollary 3.5], respectively): kuk2 1=k∇uk2 0+ZΩ u2 ,∀u∈H1(Ω),(3.2) kσk2 1=kσk2 0+krot σk2 0+k∇·σk2 0,∀σ∈H1 σ(Ω).(3.3) In particular, (3.3) implies that k∇vk2 1=k∇vk2 0+k∆vk2 0,∀v:∇v∈H1 σ(Ω). If Zis a general Banach space, its topological dual will be denoted by Z0. Moreover, the letters C, K will denote different positive constants (independent of discrete parameters) which may change from line to line (or even within the same line). 3.2 Continuous problem First we give the following definition of weak-strong solutions for problem (3.1). Definition 3.2.1 (Weak-strong solutions) Given (u0, v0)∈L2(Ω)×H1(Ω) with u0, v0≥ 0a.e. x∈Ω, a pair (u, v)is called weak-strong solution of problem (3.1) in (0,+∞), if u≥0,v≥0a.e. (t, x)∈(0,+∞)×Ω, (u, v)∈L∞(0,+∞;L2(Ω) ×H1(Ω)) ∩L2(0, T;H1(Ω) ×H2(Ω)),∀T > 0, ∂tu∈Lq0(0, T; (H1(Ω))0)and ∂tv∈Lq0(0, T;L2(Ω)),∀T > 0,(3.4) where q0= 2 in 2D and q0= 4/3in 3D (q0is the conjugate exponent of q= 2 in 2D and q= 4 in 3D); the following variational formulation holds ZT 0h∂tu, ¯ui+ZT 0 (∇u, ∇¯u) + ZT 0 (u∇v, ∇¯u) = 0,∀¯u∈Lq(0, T;H1(Ω)),∀T > 0,(3.5) the following equation holds pointwisely ∂tv−∆v+v=u2,a.e. (t, x)∈(0,+∞)×Ω,(3.6) the initial conditions (3.1)4are satisfied and the following energy inequality (in integral version) holds for a.e. t0, t1:t1≥t0≥0: E(u(t1), v(t1)) −E(u(t0), v(t0)) + Zt1 t0k∇u(s)k2 0+1 2k∇v(s)k2 1ds ≤0,(3.7) where E(u(t), v(t)) = 1 2ku(t)k2 0+1 4k∇v(t)k2 0. 83
Remark 3.2.2 In particular, the energy inequality (3.7) is valid for t0= 0. Moreover, (3.7) shows the dissipative character of the model with respect to the total energy E(u(t), v(t)). Remark 3.2.3 (Positivity) u≥0in 2D domains and v≥0in any (2Dor 3D) dimension are a consequence of (3.4)-(3.6). Indeed, this follows from the fact that in these cases we can test (3.5) by u−:= min{u, 0} ∈ L2(0, T;H1(Ω)) and (3.6) by v−:= min{v, 0} ∈ L2(0, T;H2(Ω)) ,→L2(0, T;L2(Ω)). Observe that the problem (3.1) conserves in time the total mass RΩu, because taking ¯u= 1 in (3.5), d dt ZΩ u= 0,i.e. ZΩ u(t) = ZΩ u0,∀t > 0. Moreover, integrating (3.6) in Ω we deduce the following behavior of RΩv: d dt ZΩ v+ZΩ v=ZΩ u2. We recall that in Chapter 1it was proved the existence of weak-strong solutions of problem (3.1) (satisfying in particular the energy inequality (3.7)), through convergence of a time-discrete numerical scheme associated to model (3.1). 3.2.1 Convergence at infinite time In this subsection, we will prove the exponential convergence of any weak-strong solution (u, v) of problem (3.1) obtained by Galerkin approximations. First, we will prove exponential bounds for weak-strong norms a.e. t≥0. Theorem 3.2.4 Let (u, v)be any weak-strong solution of problem (3.1) obtained by Galerkin approximations. Then, the following estimates hold k(u(t)−m0,∇v(t))k2 0≤C0e−2t,a.e. t≥0.(3.8) kv(t)−(m0)2k2 0≤C0e−t,∀t≥0,(3.9) where m0:= 1 |Ω|ZΩ u0and C0is a positive constant depending on the data (u0, v0), but independent of t. 84
Proof. For each m≥1, we say that (um, vm) is a Galerkin solution if um: [0,+∞)→Um and vm: [0,+∞)→Vmare C1functions and satisfy (∂tum,¯u)+(∇um,∇¯u)+(um∇vm,∇¯u) = 0,∀¯u∈Um, t ≥0 (∂tvm,¯v)+(∇vm,∇¯v)+(vm,¯v)−((um)2,¯v) = 0,∀¯v∈Vm, t ≥0, um(0) = um 0:= Pm(u0), vm(0) = vm 0:= Qm(v0), (3.10) where Umand Vmare finite dimensional spaces generated by orthonormal eigenfunctions of the operator (−∆ + I), with ∆ and Ibeing the laplacian and identity operators; Pm: L2(Ω) →Umdenotes the projection from L2(Ω) onto Um, and Qm:H1(Ω) →Vmthe projection from H1(Ω) onto Vm. Then, (3.10) can be regarded as a Cauchy problem for a first order ordinary differential system in time, and the classical existence and uniqueness theory for ordinary differential systems implies that, for every m≥1, there exist Tm>0 and unique functions um: [0, Tm)→Umand Vm: [0, Tm)→Vmthat solve (3.10), with either Tm= +∞or lim sup t→Tmk(um(t), vm(t))k0= +∞. Now, we are going to deduce some estimates for (um, vm) showing that only Tm= +∞can be true. We define eum:= um−m0and taking ¯u=eumand ¯v=−1 2∆vmin (3.10), we arrive at 1 2 d dtkeum(t)k2 0+1 2k∇vm(t)k2 0+keum(t)k2 1+1 2k∇vm(t)k2 1= 0,(3.11) from which we deduce that (um,∇vm) is bounded in L∞(0,+∞;L2(Ω) ×L2(Ω)), (∇um,∇vm) is bounded in L2(0,+∞;L2(Ω) ×H1(Ω)).(3.12) Moreover, Moreover, we observe that the function ym(t) = RΩvm(x, t)dx2satisfies (ym)0(t)+ ym(t)≤wm(t), with wm(t) = kum(t)k4 0. In fact, it follows by taking ¯v= 1 in (3.10), multiplying the resulting equation by ZΩ vm(x, t)dxand using the Young inequality. Therefore, ym(t) = ym(0) e−t+Zt 0 e−(t−s)wm(s)ds, which implies that ZΩ vm(x, t)dx2≤ZΩ vm 0(x)dx2+kumk4 L∞(0,+∞;L2),∀t≥0.(3.13) Then, from (3.12) and (3.13), we deduce that vmis bounded in L∞(0,+∞;H1(Ω)) ∩L2(0, T;H2(Ω)),∀T > 0.(3.14) 85
Taking into account (3.12) and (3.14), we can deduce (∂tum, ∂t∇vm) is bounded in L4/3(0, T;H1(Ω)0×H1(Ω)0). Therefore, proceeding as in Theorem 1.4.11, we obtain that there exists a subsequence m0of m, and (u, v) weak-strong solution of (3.1), such that (um0, vm0) converges to (u, v) weakly-* in L∞(0,+∞;L2(Ω)×H1(Ω)), weakly in L2(0, T;H1(Ω)× H2(Ω)) and strongly in L2(0, T;L2(Ω) ×H1(Ω)) ∩C([0, T]; H1(Ω)0×Lp(Ω)), for any T > 0, 1≤p < 6. Therefore, in particular k(eum0(t),∇vm0(t))k2 0→ k(eu(t),∇v(t))k2 0,a.e. t≥0.(3.15) Moreover, from the equality (3.11), we deduce keum0(t)k2 0+1 2k∇vm0(t)k2 0≤ k(um0 0−m0,∇vm0 0)k2 0e−2t,∀t≥0.(3.16) Thus, from (3.15)-(3.16), we arrive at (3.8). Finally, testing (3.6) by ev:= v−(m0)2, one can obtain 1 2 d dtkev(t)k2 0+kev(t)k2 1=ZΩ ((u(x, t))2−(m0)2)ev(x, t)dx, which, using the H¨older and Young inequalities, implies that d dtkev(t)k2 0+kev(t)k2 1≤ ku(t) + m0k2 L3ku(t)−m0k2 0.(3.17) Therefore, from (3.17) and (3.8), we can deduce for all t≥0, kev(t)k2 0≤ kv0−(m0)2k2 0e−t+e−tZt 0ku(s)−m0k2 0ku(s) + m0k2 L3esds ≤C0e−t+C0e−tZt 0ku(s) + m0k2 L3e−2sesds ≤C0e−t+C0e−tZt 0k∇u(s)k2 0e−sds +C0e−tZt 0ku(s) + m0k2 0e−sds, from which, using (3.12), we conclude (3.9). In next theorem, we will show, for large times, exponential bounds for more regular norms. Theorem 3.2.5 Let ε > 0. Under hypothesis of Theorem 3.2.4, there exists a constant C1>0such that if ε2≤1 2C1it holds k(u(t)−m0,∇v(t))k2 1≤2εe−1 2(t−t2),a.e. t≥t2(ε),(3.18) with t2:= t2(ε)≥0a large enough time that will be obtained in the proof. 86
Proof. We define Fm(t) := keum(t)k2 1+1 2k∇vm(t)k2 1. Then, from (3.8) and (3.11), we have that Z+∞ t Fm(s)ds ≤ k(um(t)−m0,∇vm(t))k2 0≤C0e−2t, which, in particular, implies that for all δ > 0, there exists a large enough time t0=t0(δ)≥0 such that Z+∞ t0 Fm(s)ds ≤δ. (3.19) Then, taking into account that Fm(t) satisfies (3.19), proceeding as in [3, Lemma 2.1], we have that for all δ > 0, t≥t0(δ) and τ > 0, there exists a time ¯ t∈[t, t +τ] such that Fm(¯ t)≤2δ τ.(3.20) Indeed, the set of points ¯ t∈[t, t +τ] satisfying (3.20) has measure greater than τ/2. Now, in order to obtain strong estimates, we take ¯u=−∆umand ¯v=1 2∆2vmin (3.10), and proceeding as in (1.33), we arrive at d dtk∇eum(t)k2 0+1 2k∆vm(t)k2 0+k∆eum(t)k2 0+k∆vm(t)k2 1≤Ckeum(t)k6 1+Ck∇vm(t)k6 1.(3.21) Then, adding (3.11) and (3.21), we have d dtkeum(t)k2 1+1 2k∇vm(t)k2 1+keum(t)k2 2+k∇vm(t)k2 2≤C1keum(t)k2 1+1 2k∇vm(t)k2 13, or equivalently, Fm(t) satisfies F0 m(t) + Gm(t)≤C1Fm(t)3,(3.22) with Gm(t) = keum(t)k2 2+k∇vm(t)k2 2. Therefore, taking into account that Fm(t) satisfies (3.22), proceeding as in [3, Lemma 2.2], we can deduce that for any ε > 0 and t1≥0, Fm(t1)≤ε/2⇒Fm(t)≤ε, ∀t∈t1, t1+1 2C1ε2.(3.23) Thus, as consequence of (3.20) and (3.23), following the proof of [3, Theorem 2.3], we conclude that for any ε > 0, taking τ=1 4C1ε2,δ=1 16C1εand t0=t0(δ) such that Fm(t) satisfies (3.19), where C1is the constant in the estimate (3.22), it holds Fm(t)≤ε, ∀t≥t2(ε) := t0(δ) + 1 4C1ε2.(3.24) 87
Therefore, from (3.22) and (3.24), using the fact that Gm(t)≥Fm(t) and taking εsuch that ε2≤1 2C1, we deduce F0 m(t) + 1 2Gm(t)≤0,∀t≥t2:= t2(ε) and Fm(t)≤Fm(t2)e−1 2(t−t2)≤εe−1 2(t−t2),∀t≥t2.(3.25) Moreover, from (3.24), (3.22) and (3.19) we have that (eum,∇vm) is bounded in L∞(t2,+∞;H1(Ω) ×H1(Ω)) ∩L2(t2,+∞;H2(Ω) ×H2(Ω)). Then, using the fact that (∂teum, ∂t∇vm) is bounded in L2(t2,+∞;L2(Ω) ×L2(Ω)), a compactness result of Aubin- Lions type implies that (eum,∇vm) is relatively compact in L2(t2, t3;H1(Ω) ×H1(Ω)) for all t3≥t2. Therefore, in particular for some subsequence m0of m, we have k(eum0(t),∇vm0(t))k2 1→ k(eu(t),∇v(t))k2 1,a.e. t≥t2, and using (3.25) we arrive at (3.18). 3.3 Fully Discrete Schemes associated to system (3.1) In this section, we study the large-time behavior for two fully discrete schemes associated to model (3.1): the nonlinear backward Euler for model (3.1), and the nonlinear scheme defined in Chapter 2by introducing the auxiliary variable σ=∇v. Along this section we will use repeatedly the following result (see[6, Lemma 4.1]): Lemma 3.3.1 Assume that δ, k > 0and β, dn≥0satisfy (1 + δk)dn+1 −dn≤βk, ∀n≥0. Then, dn≤(1 + δk)−(n−n0)dn0+δ−1β, ∀n≥n0≥0. 3.3.1 Scheme UV The first scheme that will be studied in this paper is obtained by using FE backward Euler for the system (3.1) (considered for simplicity a uniform partition of [0, T] with time step k=T/N : (tn=nk)n=N n=0 ). Concerning the space discretization, we consider {Th}h>0be a family of shape-regular and quasi-uniform triangulations of Ω made up of simplexes (triangles 88
in two dimensions and tetrahedra in three dimensions), so that Ω = ∪K∈ThK, where h= maxK∈ThhK, with hKbeing the diameter of K. Further, let Nh={ai}i∈Idenote the set of all the nodes of Th. We choose finite element spaces for uand v, which we denote by (Uh, Vh)⊂H1×W1,6generated by (Pm,P2m)-continuous FE, with m≥1. Then, we consider the following first order in time, nonlinear and coupled scheme (Scheme UV, from now on): Initialization: Let (u0 h, v0 h)∈Uh×Vhbe a suitable approximation of (u0, v0)∈L2(Ω) × H1(Ω), as h→0, with 1 |Ω|ZΩ u0 h=1 |Ω|ZΩ u0=m0, and satisfying (3.27) below. Time step n: Given (un−1 h, vn−1 h)∈Uh×Vh, compute (un h, vn h)∈Uh×Vhsolving ((δtun h,¯uh)+(∇un h,∇¯uh)+(un h∇vn h,∇¯uh) = 0,∀¯uh∈Uh, (δtvn h,¯vh)+(∇vn h,∇¯vh)+(vn h,¯vh)−((un h)2,¯vh)=0,∀¯vh∈Vh,(3.26) where we denote in general δtan=an−an−1 k. For the initial approximation (u0 h, v0 h)∈Uh×Vh we assume that there exists a positive constant Cindependent of (k, h) such that k(u0 h, v0 h)kL2×H1≤Ck(u0, v0)kL2×H1.(3.27) Existence, energy-stability and convergence Assuming that the functions ¯uh= 1 ∈Uhand ¯vh= 1 ∈Vh, we deduce that the scheme UV conserves the total mass RΩun h, that is, ZΩ un h=ZΩ un−1 h=···=ZΩ u0 h,(3.28) and we have the following behavior for RΩvn h δtZΩ vn h=ZΩ (un h)2−ZΩ vn h.(3.29) Theorem 3.3.2 (Unconditional existence) There exists (un h, vn h)∈Uh×Vhsolution of the scheme UV. 89
Known results The scheme US also conserves the total mass RΩun h(satisfying (3.28)), and also has the behaviour for RΩvn hgiven in (3.29). The existence of (un h,σn h)∈Uh×Σhsolution of the scheme US,vn hsolution of (3.54), and the unconditional energy-stability of the scheme US was proved in Chapter 2. In fact, the following discrete energy law holds δtE(un h,σn h)+k 2kδtun hk2 0+k 4kδtσn hk2 0+k∇un hk2 0+1 2kσn hk2 1= 0,(3.55) where E(un h,σn h) = 1 2kun hk2 0+1 4kσn hk2 0. From the (local in time) discrete energy law (3.55), the following global in time weak estimates for (un h,σn h) are deduced (see Theorem 2.3.9): k(un h,σn h)k2 0+k2 n X m=1k(δtum h, δtσm h)k2 0+k n X m=1k(∇um h,σm h)k2 L2×H1≤C0,∀n≥1, k n0+n X m=n0+1kum hk2 1≤C0+C1(nk),∀n≥1, where n0≥0 is any integer and C0, C1are positive constants depending on the data (Ω, u0,σ0), but independent of (k, h) and (n, n0). Large-time behavior of scheme US Theorem 3.3.8 Let (un h,σn h)be a solution of the scheme US associated to an initial data (u0 h,σ0 h)∈Uh×Σhwhich is a suitable approximation of (u0,σ0)∈L2(Ω)×L2(Ω), as h→0, with 1 |Ω|ZΩ u0 h=1 |Ω|ZΩ u0=m0. Then, k(un h−m0,σn h)k2 0≤C0e−2 1+2kkn,∀n≥0,(3.56) kX m>nkeum hk2 1+1 2kσm hk2 1≤C0e−2 1+2kkn,∀n≥0,(3.57) where C0is a positive constant depending on the data (u0,σ0), but independent of (k, h)and n. Proof. Taking ¯uh=eun h:= un h−m0in (3.53)1,¯ σh=1 2σn hin (3.54) and using (3.28), we obtain δt1 2keun hk2 0+1 4kσn hk2 0+k 2kδteun hk2 0+k 4kδtσn hk2 0+keun hk2 1+1 2kσn hk2 1= 0.(3.58) 96
Then, from (3.58) we can obtain (1 + 2k)keun hk2 0+1 2kσn hk2 0−keun−1 hk2 0+1 2kσn−1 hk2 0≤0.(3.59) Then, applying Lemma 3.3.1 to (3.59), and proceeding as in (3.49), we arrive at (3.56). Finally, from (3.56) and (3.58), we have that for all n≥0, kX m>nkeum hk2 1+1 2kσm hk2 1≤ k(un h−m0,σn h)k2 0≤C0e−2 1+2kkn. Corollary 3.3.9 Let vn h=vn h((un h)2)be a solution of (3.54) associated to an initial data v0 h∈Vhwhich is a suitable approximation of v0∈H1(Ω), as h→0. Then (3.45) holds. Proof. The proof follows as in Theorem 3.3.7. Now, in order to obtain more regular estimates, we consider the linear operators e Ah: Uh→Uhand Bh:Σh→Σhdefined as follows ((e Ahuh,¯uh) = (∇uh,∇¯uh)+(uh,¯uh),∀¯uh∈Uh, (Bhσh,¯ σh) = (∇·σh,∇· ¯ σh) + (rot σh,rot ¯ σh)+(σh,¯ σh),∀¯ σh∈Σh. Then, we rewrite (3.53) as ((δtun h,¯uh)+(e Ahun h,¯uh)−(un h,¯uh)+(un hσn h,∇¯uh)=0,∀¯uh∈Uh, (δtσn h,¯ σh)+(Bhσn h,¯ σh)−2(un h∇un h,¯ σh)=0,∀¯ σh∈Σh.(3.60) Moreover, the following estimates hold (see for instance, Lemma 2.3.1): kuhkW1,6≤CkAhuhk0∀uh∈Uh,kσhkW1,6≤CkBhσhk0∀σh∈Σh. Theorem 3.3.10 Under hypothesis of Theorem 3.3.8, the following estimate holds kX m>nk(e Aheum h, Bhσm h)k2 0≤C1 ke−2 1+2kk(n−1) +1 k2e−4 1+2kk(n−1) + 1e−2 1+2kkn ,∀n≥1, (3.61) where Cis a positive constant independent of (k, h)and n. 97
Proof. We define Fn:= keun hk2 1+1 2kσn hk2 1. Then, from (3.57) we have that X m>n Fm≤1 kC0e−2 1+2kkn,∀n≥0, which, taking into account that Fm≥0 for all m∈N, implies that Fn≤1 kC0e−2 1+2kk(n−1),∀n≥1.(3.62) Now, taking ¯uh=e Aheun hand ¯ σh=Bhσn hin (3.60), we have 1 2δtk(eun h,σn h)k2 1+k 2k(δteun h, δtσn h)k2 1+k(e Aheun h, Bhσn h)k2 0≤ keun hk2 1+τk(e Aheun h, Bhσn h)k2 0 +Cτk∇eun hk2 L3kσn hk2 L6+Cτkeun hk2 L6k∇·σn hk2 L3+Cτk∇eun hk2 L3keun hk2 L6+Cτ(m0)2k(∇eun h,∇·σn h)k2 0 ≤ keun hk2 1+τk(e Aheun h, Bhσn h)k2 0+τk(∇eun h,∇·σn h)k2 L6+Cτkeun hk2 1kσn hk4 1+Cτkeun hk6 1 +Cτkσn hk2 1keun hk4 1+Cτ(m0)2k(eun h,σn h)k2 1.(3.63) Therefore, taking into account that k(uh,σh)k2 W1,6≤Ck(e Ahuh, Bhσh)k2 0for all (uh,σh)∈ Uh×Σh(see Lemma 2.3.1), from (3.63) (choosing τsmall enough) we deduce δtk(eun h,σn h)k2 1+k(e Aheun h, Bhσn h)k2 0≤C1k(eun h,σn h)k2 13+C2k(eun h,σn h)k2 1.(3.64) Then, from (3.64), taking into account (3.62) and (3.57), we deduce for all n≥1, kX m>nk(e Aheum h, Bhσm h)k2 0≤ k(eun h,σn h)k2 1+ ( 1 k2C2 0C1e−4 1+2kk(n−1) +C2)kX m>nk(eum h,σm h)k2 1 ≤1 kC0e−2 1+2kk(n−1) +C3 1 k2e−4 1+2kk(n−1) +C4e−2 1+2kkn, from which we conclude (3.61). Remark 3.3.11 In the case of the scheme UV it is not clear how to obtain one more regular estimate equivalent to the obtained in Theorem 3.3.10 for the scheme US. In fact, a key step in the proof of Theorem 3.3.10, is to integrate by parts in the term (un hσn h,∇(e Aheun h)) arriving at (∇un h·σn h,e Aheun h) + (un h∇·σn h,e Aheun h), which it is not possible for the scheme UV in the term (un h∇vn h,∇(e Aheun h)), because un h∇vn hdoes not have a derivative in L2(Ω). 98
3.4 Numerical Simulations In this section we will compare the results of several numerical simulations that we have carried out using the schemes studied in the paper. We are considering P1-continuous approximation for un h,σn hand P2-continuous approximation for vn h. Moreover, we have chosen the 2D domain Ω = [0,2]2using a structured mesh, and all the simulations are carried out using FreeFem++ software. The linear iterative method used to approach the nonlinear schemes US and UV is the Newton Method, and in all the cases, the iterative method stops when the relative error in L2-norm is less than ε= 10−6. 3.4.1 Positivity The aim of this subsection is to compare the fully discrete schemes UV and US in terms of positivity. Theoretically, for both schemes, is not clear the positivity of the variables un hand vn h. In fact, in some simulations, we obtain numerical results in which un his negative. For example, choosing k= 10−5, the initial conditions (see Figure 3.1): u0=−10xy(2 −x)(2 −y)exp(−10(y−1)2−10(x−1)2) + 10.0001 and v0=200xy(2 −x)(2 −y)exp(−30(y−1)2−30(x−1)2)+0.0001, and taking meshes in space increasingly thinner (h=1 10,h=1 20,h=1 35 and h=1 75), we (a) Initial cell density u0(b) Initial chemical concentration v0 Figure 3.1: Initial conditions. obtain that in both schemes, the discrete cell density un htakes negative values for some x∈Ω in some times tn>0 (see Figures 3.2-3.5). Moreover, as htends to 0, (a) the behaviour of both schemes is increasingly similar, and (b) the negative values taken for un hin both 99
schemes are closer to 0. This is in accordance with the results obtained in Chapter 1, where it was proved that the only time-discrete schemes corresponding to the schemes UV and US are equivalents and have nonnegative solution (un, vn). 0 0.005 0.01 0.015 Time -14 -12 -10 -8 -6 -4 -2 0 2 Min un h US UV Figure 3.2: Minimum values of un h, with h=1 10. Remark 3.4.1 In the computations, the execution time for the scheme UV is smaller than the execution time for the scheme US. In fact, the scheme UV is twice faster than the scheme US. On the other hand, with respect to the discrete chemical concentration vn h, we observe that the same behavior is obtained for the minimum of vn hin both schemes. In fact, independently of h, if v0is positive, then vhalso is positive (we show this behavior in Figure 3.6 for the case h=1 35, but the same holds for the cases h=1 10,h=1 20 and h=1 75). 3.4.2 Energy-Stability In Lemma 3.3.3, the unconditional energy-stability for the scheme UV with respect to the energy E(u, v) was proved. In fact, if (un h, vn h) is any solution of the scheme UV, the following relation holds δtE(un h, vn h) + k∇un hk2 0+1 2k(Ah−I)vn hk2 0+1 2k∇vn hk2 0≤0,∀n. 100
0 0.005 0.01 0.015 Time -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 Min un h US UV Figure 3.3: Minimum values of un h, with h=1 20. On the other hand, in Chapter 2it was proved the unconditional energy-stability for the scheme US with respect to the modified energy E(u, σ). Even more, if (un h,σn h) is any solution of the scheme US, it holds δtE(un h,σn h) + k∇un hk2 0+1 2kσn hk2 1≤0,∀n. Then, the aim of this subsection is to compare numerically the energy-stability of the schemes UV and US with respect to the energy E(u, v). Indeed, if we take k= 10−6,h=1 25 and the initial conditions u0=−10xy(2 −x)(2 −y)exp(−10(y−1)2−10(x−1)2) + 10.0001 and v0=20xy(2 −x)(2 −y)exp(−30(y−1)2−30(x−1)2)+0.0001, we obtain that: (a) The schemes UV and US satisfy the energy decreasing in time property for the energy E(u, v), that is, E(un h, vn h)≤E(un−1 h, vn−1 h) for all n, see Figure 3.7. (b) The schemes UV and US satisfy (see Figure 3.8) RE(un h, vn h) := δtE(un h, vn h) + k∇un hk2 0+1 2k(Ah−I)vn hk2 0+1 2k∇vn hk2 0≤0,∀n. 101
0 0.005 0.01 0.015 Time -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 Min un h US UV Figure 3.4: Minimum values of un h, with h=1 35. 0 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 Time -0.05 -0.04 -0.03 -0.02 -0.01 0 0.01 Min un h US UV Figure 3.5: Minimum values of un h, with h=1 75. 102
0 0.005 0.01 0.015 Time 0 0.5 1 1.5 Min vn h US UV Figure 3.6: Minimum values of vn h, with h=1 35. 0 0.005 0.01 0.015 0.02 0.025 Time 150 200 250 300 350 400 450 500 E(un h,vn h) US UV Figure 3.7: Energy E(un h, vn h) of schemes UV and US. 103
0 0.005 0.01 0.015 0.02 0.025 Time -8 -7 -6 -5 -4 -3 -2 -1 0 RE(u n h,vn h) ×104 US UV Figure 3.8: Residue RE(un h, vn h) of schemes UV and US. 104
Bibliography [1] C. Amrouche and N.E.H. Seloula, Lp-theory for vector potentials and Sobolev’s inequalities for vector fields: application to the Stokes equations with pressure boundary conditions. Math. Models Methods Appl. Sci. 23 (2013), no. 1, 37–92. [2] M. Bessemoulin-Chatard and A. J¨ungel, A finite volume scheme for a Keller-Segel model with additional cross-diffusion. IMA J. Numer. Anal. 34 (2014), no. 1, 96–122. [3] B. Climent-Ezquerra, F. Guill´en-Gonz´alez and M.A. Rodr´ıguez-Bellido, Stability for nematic liquid crystals with stretching terms. Internat. J. Bifur. Chaos Appl. Sci. Engrg. 20 (2010), no. 9, 2937–2942. [4] T. Cieslak, P. Lauren¸cot and C. Morales-Rodrigo, Global existence and convergence to steady states in a chemorepulsion system. Parabolic and Navier-Stokes equations. Part 1, 105–117, Banach Center Publ., 81, Part 1, Polish Acad. Sci. Inst. Math., Warsaw, 2008. [5] F. Guill´en-Gonz´alez and M. Samsidy, Stability and convergence at infinite time of several fully discrete schemes for a Ginzburg-Landau model for nematic liquid crystal flows. Discrete Contin. Dyn. Syst. 32 (2012), no. 12, 4229–4246. [6] Y. He and K. Li, Asymptotic behavior and time discretization analysis for the nonstationary Navier-Stokes problem. Numer. Math. 98 (2004), no. 4, 647–673. [7] T. Hillen and K.J. Painter, A user’s guide to PDE models for chemotaxis. J. Math. Biol. 58 (2009), no. 1-2, 183–217. [8] S. Hittmeir and A. J¨ungel, Cross diffusion preventing blow-up in the two-dimensional Keller-Segel model. SIAM J. Math. Anal. 43 (2011), no. 2, 997–1022. [9] D. Horstmann, From 1970 until present: the Keller-Segel model in chemotaxis and its consequences. I. Jahresber. Deutsch. Math.-Verein. 105 (2003), no. 3, 103–165. 105
Remark 4.3.1 The idea is to define a fully discrete scheme associated to (4.7), taking in general ε=ε(k, h), such that ε(k, h)→0as (k, h)→0, where kis the time step and hthe mesh size. Observe that multiplying (4.7)1by F0 ε(uε), (4.7)2by −∆vε, integrating over Ω and adding, again the chemotactic and production terms cancel, and we obtain the following energy law d dt ZΩFε(uε) + 1 2|∇vε|2dx+ZΩF00 ε(uε)|∇uε|2+|∆vε|2+|∇vε|2dx= 0. In particular, the modified energy Eε(u, v) = ZΩFε(u) + 1 2|∇v|2dx is decreasing in time. Then, we consider a fully discrete approximation using FE in space and backward Euler in time (considered for simplicity on a uniform partition of [0, T] with time step k=T/N : (tn=nk)n=N n=0 ). Let Ω be a polygonal domain. We consider a shaperegular and quasi-uniform family of triangulations of Ω, denoted by {Th}h>0, with simplices K,hK=diam(K) and h:= maxK∈ThhK, so that Ω = ∪K∈ThK. Moreover, in this case we will assume the following hypothesis: (H) The triangulation is structured in the sense that all simplices have a right angle. We choose the following continuous FE spaces for uand v: (Uh, Vh)⊂H1(Ω)2,generated by P1,Pmwith m≥1. Remark 4.3.2 The right angled requirement and the choice of P1-continuous FE for Uhare necessary in order to obtain the relation (4.10) below, which is essential in order to prove the energy-stability of the scheme UV (see Theorem 4.3.7 below). Let Jbe the set of vertices of Thand {aj}j∈Jthe coordinates of these vertices. We denote the Lagrange interpolation operator by Πh:C(Ω) →Uh, and we introduce the discrete semiinner product on C(Ω) (which is an inner product in Uh) and its induced discrete seminorm (norm in Uh): (u1, u2)h:= ZΩ Πh(u1u2),|u|h=p(u, u)h.(4.8) 112
Remark 4.3.3 In Uh, the norms |·|hand k·k0are equivalents uniformly with respect to h (see [4]). We consider also the L2-projection Qh:L2(Ω) →Uhgiven by (Qhu, ¯u)h= (u, ¯u),∀¯u∈Uh,(4.9) and the standard H1-projection Rh:H1(Ω) →Vh. Moreover, for each ε∈(0,1) we consider the construction of the operator Λε:Uh→L∞(Ω)d×dgiven in [3], satisfying that Λεuhis a symmetric and positive definite matrix for all uh∈Uhand a.e. xin Ω, and the following relation holds (Λεuh)∇Πh(F0 ε(uh)) = ∇uhin Ω.(4.10) Basically, Λεuhis a constant by elements matrix such that (4.10) holds by elements. We highlight that (4.10) is satisfied due to the right angled constraint (H) and the choice of P1-continuous FE for Uh. Moreover, the following stability estimate holds [3,10] kΛε(uh)kr Lr≤C(1 + kuhk2 1),∀uh∈Uh(for r= 2(d+ 1)/d),(4.11) where the constant C > 0 is independent of εand h. We recall the result below concerning to Λε(·) (see [3, Lemma 2.1]). Lemma 4.3.4 Let k·kdenote the spectral norm on Rd×d. Then for any given ε∈(0,1) the function Λε:Uh→[L∞(Ω)]d×dis continuous and satisfies εξTξ≤ξTΛε(uh)ξ≤ε−1ξTξ, ∀ξ∈Rd,∀uh∈Uh.(4.12) In particular, for all uh 1, uh 2∈Uhand K∈Thwith vertices {aK l}d l=0, it holds k(Λε(uh 1)−Λε(uh 2))|Kk ≤ ε−2max l=1,...,d{|uh 1(aK l)−uh 2(aK l))|+|uh 1(aK 0)−uh 2(aK 0))|},(4.13) where aK 0is the right-angled vertex. Let Ah:Vh→Vhbe the linear operator defined as follows (Ahvh,¯v)=(∇vh,∇¯v)+(vh,¯v),∀¯v∈Vh. Then, the following estimate holds (see for instance, Lemma 2.3.1): kvhkW1,6≤CkAhvhk0,∀vh∈Vh.(4.14) Taking into account the regularized problem (4.7), we consider the following first order in time, nonlinear and coupled scheme: 113
•Scheme UV : Initialization: Let (u0 h, v0 h) = (Qhu0, Rhv0)∈Uh×Vh. Time step n: Given (un−1 ε, vn−1 ε)∈Uh×Vh, compute (un ε, vn ε)∈Uh×Vhsolving ((δtun ε,¯u)h+ (∇un ε,∇¯u) + (Λε(un ε)∇vn ε,∇¯u) = 0,∀¯u∈Uh, (δtvn ε,¯v)+(Ahvn ε,¯v)−(un ε,¯v)=0,∀¯v∈Vh,(4.15) where, in general, we denote δtan:= an−an−1 k. 4.3.1 Mass conservation and Energy-stability Since ¯u= 1 ∈Uhand ¯v= 1 ∈Vh, we deduce that the scheme UV is conservative in un ε, that is, (un ε,1) = (un ε,1)h= (un−1 ε,1)h=···= (u0 h,1)h= (u0 h,1) = (Qhu0,1) = (u0,1) := m0,(4.16) and we have the following behavior for RΩvn ε: δtZΩ vn ε+ZΩ vn ε=ZΩ un ε=m0.(4.17) Lemma 4.3.5 (Estimate of RΩvn ε)The following estimate holds ZΩ vn ε≤(1 + k)−nZΩ v0+m0,∀n≥0.(4.18) Proof. From (4.17) we have (1 + k)RΩvn ε−RΩvn−1 ε≤k m0, and therefore, applying Lemma 4.2.1 (for δ= 1 and β=m0), we arrive at ZΩ vn ε≤(1 + k)−nZΩ v0 ε+m0= (1 + k)−nZΩ Rhv0+m0, which implies (4.18). Definition 4.3.6 A numerical scheme with solution (un ε, vn ε)is called energy-stable with respect to the energy Eh ε(u, v)=(Fε(u),1)h+1 2k∇vk2 0(4.19) if this energy is time decreasing, that is, Eh ε(un ε, vn ε)≤Eh ε(un−1 ε, vn−1 ε)for all n≥1. 114
Theorem 4.3.7 (Unconditional stability) The scheme UV is unconditional energy stable with respect to Eh ε(u, v). In fact, if (un ε, vn ε)is a solution of UV, then the following discrete energy law holds δtEh ε(un ε, vn ε) + εk 2kδtun εk2 0+k 2kδt∇vn εk2 0+εk∇un εk2 0+k(Ah−I)vn εk2 0+k∇vn εk2 0≤0.(4.20) Proof. Testing (4.15)1by ¯u= Πh(F0 ε(un ε)) and (4.15)2by ¯v= (Ah−I)vn ε, adding and taking into account that Λε(un ε) is symmetric as well as (4.10) (which implies that ∇Πh(F0 ε(un ε)) = Λ−1 ε(un ε)∇un ε), the terms −(Λε(un ε)∇vn ε,∇Πh(F0 ε(un ε))) = −(∇vn ε,Λε(un ε)∇Πh(F0 ε(un ε))) = −(∇vn ε,∇un ε) and (un ε,(Ah−I)vn ε)=(∇un ε,∇vn ε) cancel, and we obtain (δtun ε, F0 ε(un ε))h+ZΩ (∇un ε)T·Λ−1 ε(un ε)·∇un εdx +δt1 2k∇vn εk2 0+k 2kδt∇vn εk2 0+k(Ah−I)vn εk2 0+k∇vn εk2 0= 0.(4.21) Moreover, observe that from the Taylor formula we have Fε(un−1 ε) = Fε(un ε) + F0 ε(un ε)(un−1 ε−un ε) + 1 2F00 ε(θun ε+ (1 −θ)un−1 ε)(un−1 ε−un ε)2, and therefore, F0 ε(un ε)δtun ε=δtFε(un ε)+k 2F00 ε(θun ε+ (1 −θ)un−1 ε)(δtun ε)2.(4.22) Then, using (4.22) and taking into account that Πhis linear and F00 ε(s)≥εfor all s∈R, we have (δtun ε, F0 ε(un ε))h=ZΩ Πh(δtun ε·F0 ε(un ε)) =δtZΩ Πh(Fε(un ε))+k 2ZΩ Πh(F00 ε(θun ε+ (1 −θ)un−1 ε)(δtun ε)2) ≥δt(Fε(un ε),1)h+εk 2|δtun ε|2 h.(4.23) Thus, from (4.12), (4.21), (4.23) and Remark 4.3.3, we arrive at (4.20). Corollary 4.3.8 (Uniform estimates) Assume that (u0, v0)∈L2(Ω) ×H1(Ω). Let (un ε, vn ε)be a solution of scheme UV. Then, it holds (Fε(un ε),1)h+1 2kvn εk2 1+k n X m=1 εk∇um εk2 0+k(Ah−I)vm εk2 0+k∇vm εk2 0≤C0,∀n≥1,(4.24) 115
k n+n0 X m=n0+1kvm εk2 W1,6≤C1(1 + kn),∀n≥1,(4.25) where the integer n0≥0is arbitrary, with the constants C0, C1>0depending on the data (Ω, u0, v0), but independent of k, h, n and ε. Moreover, if ε∈(0, e−2), the following estimates hold ZΩ (Πh(un ε−))2≤C0ε, and ZΩ|un ε| ≤ m0+C√ε, ∀n≥1,(4.26) where un ε−:= min{un ε,0} ≤ 0and the constant C > 0depends on the data (Ω, u0, v0), but is independent of k, h, n and ε. Proof. First, using the inequality s(ln s −1) ≤s2for all s > 0 (which implies Fε(s)≤ C(s2+ 1) for all s≥0) and taking into account that (u0 h, v0 h)=(Qhu0, Rhv0), u0≥0 (and therefore, u0 h≥0), as well as the definition of Fε, we have Eh ε(u0 h, v0 h) = ZΩ Πh(Fε(u0 h)) + 1 2k∇v0 hk2 0≤CZΩ Πh((u0 h)2+ 1) + 1 2k∇v0 hk2 0 ≤C(ku0 hk2 0+k∇v0 hk2 0+ 1) ≤C(ku0k2 0+kv0k2 1+ 1) ≤C0,(4.27) with the constant C0>0 depending on the data (Ω, u0, v0), but independent of k, h, n and ε. Therefore, from the discrete energy law (4.20) and (4.27), we have Eh ε(un ε, vn ε) + k n X m=1 εk∇um εk2 0+k(Ah−I)vm εk2 0+k∇vm εk2 0≤Eh ε(u0 h, v0 h)≤C0.(4.28) Thus, from (4.18) and (4.28) we conclude (4.24). Moreover, adding (4.20) from m=n0+1 to m=n+n0, and using (4.14) and (4.24), we deduce (4.25). By other hand, if ε∈(0, e−2), from (4.6)2and taking into account that Fε(s)≥0 for all s∈R, we have 1 2ε(un ε−(x))2≤Fε(un ε(x)) for all un ε∈Uh; and therefore, using that (Πh(u))2≤Πh(u2) for all u∈C(Ω), we have 1 2εZΩ (Πh(un ε−))2≤1 2εZΩ Πh((un ε−)2)≤ZΩ Πh(Fε(un ε)) ≤C0, where in the last inequality (4.24) was used. Thus, we obtain (4.26)1. Finally, considering un ε+:= max{un ε,0} ≥ 0, taking into account that un ε=un ε++un ε−and |un ε|=un ε+−un ε−= un ε−2un ε−, using the H¨older and Young inequalities as well as (4.16) and (4.26)1, we have ZΩ|un ε| ≤ ZΩ Πh|un ε|=ZΩ un ε−2ZΩ Πh(un ε−)≤m0+CZΩ (Πh(un ε−))21/2≤m0+C√ε, which implies (4.26)2. 116
Remark 4.3.9 The l∞(L1)-norm is the only norm in which un εis bounded independently of (k, h)and ε(see (4.26)2). However, we can also obtain ε-dependent bounds for un ε. In fact, from (4.6) and taking into account that ε∈(0, e−2), we can deduce that ε 2s2≤Fε(s) + 2 for all s∈R, which together with (4.24), implies that (√ε un ε)is bounded in l∞(L2)∩l2(H1). Remark 4.3.10 (Approximated positivity) 1. From (4.26)1, the following estimate holds max n≥0kΠh(un ε−)k2 0≤C0ε. 2. Assuming Vhfurnished by P1-continuous FE and considering the following approximation for the v-equation: (δtvn ε,¯v)h+ ( e Ahvn ε,¯v)h−(un ε,¯v)h= 0,∀¯v∈Vh,(4.29) where e Ah:Vh→Vhis the operator defined by (e Ahvh,¯v)h= (∇vh,∇¯v)+(vh,¯v)hfor all ¯v∈Vh, then the unconditional energy-stability also holds and the following estimates are satisfied max n≥0kΠh(vn ε−)k2 0≤Cε and k n X m=1kΠh(vn ε−)k2 1≤Cε(kn),(4.30) where the constant Cis independent of k, h, n and ε. In fact, testing by ¯v= Πh(vn ε−)∈ Vhin (4.29), taking into account that (∇Πh(vn ε+),∇Πh(vn ε−)) ≥0(owing to the interior angles of the triangles or tetrahedra are less than or equal to π/2), and using again that (Πh(v))2≤Πh(v2)for all v∈C(Ω), we have 1 k+ 1kΠh(vn ε−)k2 0+k∇Πh(vn ε−)k2 0≤ZΩ Πhun ε+1 kvn−1 εvn ε− ≤ZΩ Πhun ε−+1 kvn−1 ε−vn ε− ≤1 21 k+ 1kΠh(vn ε−)k2 0+1 2kΠh(un ε−)k2 0+1 2kkΠh(vn−1 ε−)k2 0, from which, using (4.26)1, we arrive at 1 21 k+ 1kΠh(vn ε−)k2 0+k∇Πh(vn ε−)k2 0≤1 2C0ε+1 2kkΠh(vn−1 ε−)k2 0.(4.31) Therefore, if v0 h≥0(which holds for instance by considering v0 h=e Rhv0, where e Rh is an average interpolator of Clement or Scott-Zhang type, and using that v0≥0), using Lemma 4.2.1 in (4.31), we conclude (4.30)1. Finally, multiplying (4.31) by k and adding from m= 1 to m=n, and using again that v0 h≥0, we arrive at (4.30)2. 117
4.3.2 Well-posedness In this subsection, we will prove the well-posedness of the scheme UV. We recall that, taking into account that we remain in finite dimension, all norms are equivalents. Theorem 4.3.11 (Unconditional existence) There exists at least one solution (un ε, vn ε) of the scheme UV. Proof. We will use the Leray-Schauder fixed point theorem. With this aim, given (un−1 ε, vn−1 ε)∈Uh×Vh, we define the operator R:Uh×Vh→Uh×Vhby R(eu, ev) = (u, v), such that (u, v)∈Uh×Vhsolves the following linear decoupled problem u∈Uhs.t. 1 k(u, ¯u)h+ (∇u, ∇¯u) = 1 k(un−1 ε,¯u)h−(Λε(eu)∇ev, ∇¯u),∀¯u∈Uh,(4.32) v∈Vhs.t. 1 k(v, ¯v)+(Ahv, ¯v) = 1 k(vn−1 ε,¯v)+(eu, ¯v),∀¯v∈Vh.(4.33) 1. Ris well defined. Applying the Lax-Milgram theorem to (4.32) and (4.33), we can deduce that, for each (eu, ev)∈Uh×Vh, there exists a unique (u, v)∈Uh×Vhsolution of (4.32)-(4.33). 2. Let us now prove that all possible fixed points of λR (with λ∈(0,1]) are bounded. In fact, observe that if (u, v) is a fixed point of λR, then R(u, v) = (1 λu, 1 λv), and therefore (u, v) satisfies the coupled system 1 k(u, ¯u)h+ (∇u, ∇¯u) + λ(Λε(u)∇v, ∇¯u) = λ k(un−1 ε,¯u)h,∀¯u∈Uh, 1 k(v, ¯v)+(Ahv, ¯v)−λ(u, ¯v) = λ k(vn−1 ε,¯v),∀¯v∈Vh. (4.34) Then, testing (4.34)1and (4.34)2by ¯u= Πh(F0 ε(u)) and ¯v= (Ah−I)vrespectively, proceeding as in Theorem 4.3.7 and taking into account that λ∈(0,1], we obtain (Fε(u),1)h+1 2k∇vk2 0+kεk∇uk2 0+k(Ah−I)vk2 0+k∇vk2 0 ≤(Fε(λun−1 ε),1)h+λ2 2k∇vn−1 εk2 0≤C(un−1 ε, vn−1 ε),(4.35) where the last estimate is λ-independent (arguing as in (4.27)). Moreover, proceeding as in Lemma 4.3.5 and Corollary 4.3.8 (taking into account (4.35)), we deduce k(u, v)kL1×H1≤C, where the constant Cdepends on data (Ω, un−1 ε, vn−1 ε, ε), but it is independent of λand h. 118
3. We prove that Ris continuous. Let {(eul,evl)}l∈N⊂Uh×Vh,→W1,∞(Ω)2be a sequence such that (eul,evl)→(eu, ev) in Uh×Vhas l→+∞.(4.36) In particular, since we remain in finite dimension, {(eul,evl)}l∈Nis bounded in W1,∞(Ω)2. Then, if we denote (ul, vl) = R(eul,evl), we can deduce 1 2kk(ul, vl)k2 0+1 2k∇ulk2 0+1 2kvlk2 1 ≤1 2kk(un−1 ε, vn−1 ε)k2 0+ (1 + keulk2 1)2/rk∇evlk2 L∞+Ckeulk2 0≤C, where in the first inequality (4.11) was used and Cis a constant independent of l∈N. Therefore, {(ul, vl) = R(eul,evl)}l∈Nis bounded in Uh×Vh,→W1,∞(Ω)2. Then, there exists a subsequence of {R(eul,e σl)}l∈N, still denoted by{R(eul,e σl)}l∈N, such that R(eul,evl)→(u0, v0) in W1,∞(Ω)2,as l→+∞.(4.37) Then, from (4.36)-(4.37) and using Lemma 4.3.4, a standard procedure allows us to pass to the limit, as lgoes to +∞, in (4.32)-(4.33) (with (eul,evl) and (ul, vl) instead of (eu, ev) and (u, v) respectively), and we deduce that R(eu, ev)=(u0, v0). Therefore, we have proved that any convergent subsequence of {R(eul,evl)}l∈Nconverges to R(eu, ev) in Uh×Vh, and from uniqueness of R(eu, ev), we conclude that the whole sequence R(eul,evl)→R(eu, ev) in Uh×Vh. Thus, Ris continuous. Therefore, the hypotheses of the Leray-Schauder fixed point theorem (in finite dimension) are satisfied and we conclude that the map Rhas a fixed point (u, v), that is R(u, v) = (u, v), which is a solution of the scheme UV. Lemma 4.3.12 (Conditional uniqueness) If k g(h, ε)<1(where g(h, ε)↑+∞as h↓0 or ε↓0), then the solution (un ε, vn ε)of the scheme UV is unique. Proof. Suppose that there exist (un,1 ε, vn,1 ε),(un,2 ε, vn,2 ε)∈Uh×Vhtwo possible solutions of the scheme UV. Then, defining u=un,1 ε−un,2 εand v=vn,1 ε−vn,2 ε, we have that (u, v)∈Uh×Vhsatisfies, for all (¯u, ¯v)∈Uh×Vh, 1 k(u, ¯u)h+ (∇u, ∇¯u) + (Λε(un,1 ε)∇v, ∇¯u) + ((Λε(un,1 ε)−Λε(un,2 ε))∇vn,2 ε,∇¯u) = 0,(4.38) 1 k(v, ¯v)+(Ahv, ¯v)=(u, ¯v).(4.39) 119
Taking ¯u=u, ¯v=Ahvin (4.38)-(4.39), adding the resulting expressions and using the fact that ZΩ u= 0 and the equivalence of the norms k·k0and |·|hin Uhgiven in Remark 4.3.3, we obtain 1 kk(u, ∇v)k2 0+k(u, Ahv)k2 H1×L2≤ kuk0kAhvk0 +kΛε(un,1 ε)kL6k∇vkL3k∇uk0+kΛε(un,1 ε)−Λε(un,2 ε)kL∞k∇vn,2 εk0k∇uk0 ≤1 4kAhvk0+kuk2 0+1 4k∇uk2 0+1 4kAhvk2 0+CkΛε(un,1 ε)k4 L6k∇vk2 0 +1 4k∇uk2 0+kΛε(un,1 ε)−Λε(un,2 ε)k2 L∞k∇vn,2 εk2 0. Then, taking into account (4.11), (4.13), (4.24), (4.26)2and using the inverse inequalities: kuhk2 L6≤C1(h)kuhk2 Lr,kuhk2 1≤C2(h)kuhk2 L1and kuhk2 L∞≤C3(h)kuhk2 0for all uh∈Uh, we have k(u, ∇v)k2 0+k 2k(u, Ahv)k2 H1×L2≤k1 + CkΛε(un,1 ε)k4 L6k(u, ∇v)k2 0+kC0Cε−2kuk2 L∞ ≤k1 + C1(h)2(1 + C2(h))4/r +kC0C3(h)ε−2k(u, ∇v)k2 0:= k g(h, ε)k(u, ∇v)k2 0. Therefore, if k g(h, ε)<1, we conclude that u= 0, and therefore (from (4.39)) v= 0. 4.4 Scheme US In this section, we propose another energy-stable nonlinear fully discrete scheme associated to model (4.1), which is obtained by introducing the auxiliary variable σ=∇v. In fact, taking into account the functions λεand Fεand its derivatives (given in (4.4)-(4.5)), another regularized version of problem (4.1) reads: Find uε: Ω×[0, T]→Rand σε: Ω×[0, T]→Rd such that ∂tuε−∇·(λε(uε)∇(F0 ε(uε))) −∇·(uεσε) = 0 in Ω, t > 0, ∂tσε+ rot(rot σε)−∇(∇·σε) + σε=uε∇(F0 ε(uε)) in Ω, t > 0, ∂uε ∂n= 0 on ∂Ω, t > 0, σε·n= 0,[rot σε×n]tang = 0 on ∂Ω, t > 0, uε(x,0) = u0(x)≥0,σε(x,0) = ∇v0(x),in Ω. (4.40) This kind of formulation considering σ=∇vas auxiliary variable has been used in the construction of numerical schemes for other chemotaxis models (see for instance [17] and 120
Chapter 2of this PhD thesis). Once problem (4.40) is solved, we can recover vεfrom uε solving ∂tvε−∆vε+vε=uεin Ω, t > 0, ∂vε ∂n= 0 on ∂Ω, t > 0, vε(x,0) = v0(x)≥0 in Ω. (4.41) Observe that multiplying (4.40)1by F0 ε(uε), (4.40)2by σε, and integrating over Ω, we obtain the following energy law d dt ZΩFε(uε) + 1 2|σε|2dx+ZΩ λε(uε)|∇(F0 ε(uε))|2dx+kσεk2 1= 0. In particular, the modified energy Eε(u, σ) = ZΩFε(u) + 1 2|σ|2dx is decreasing in time. Then, we consider a fully discrete approximation of the regularized problem (4.40) using a FE discretization in space and the backward Euler discretization in time (again considered for simplicity on a uniform partition of [0, T] with time step k=T/N : (tn=nk)n=N n=0 ). Concerning the space discretization, we consider the triangulation as in the scheme UV, but in this case without imposing the constraint (H) related with the right-angles simplices. We choose the following continuous FE spaces for uε,σε, and vε: (Uh,Σh, Vh)⊂H1(Ω)3,generated by P1,Pm,Prwith m, r ≥1. Then, we consider the following first order in time, nonlinear and coupled scheme: •Scheme US: Initialization: Let (u0 h,σ0 h) = (Qhu0,e Qh(∇v0)) ∈Uh×Σh. Time step n: Given (un−1 ε,σn−1 ε)∈Uh×Σh, compute (un ε,σn ε)∈Uh×Σhsolving ((δtun ε,¯u)h+ (λε(un ε)∇Πh(F0 ε(un ε)),∇¯u) = −(λε(un ε)σn ε,∇¯u),∀¯u∈Uh, (δtσn ε,¯ σ)+(Bhσn ε,¯ σ) = (λε(un ε)∇Πh(F0 ε(un ε)),¯ σ),∀¯ σ∈Σh,(4.42) where Qhis the L2-projection on Uhdefined in (4.9), e Qhis the standard L2-projection on Σh, and the operator Bhis defined as (Bhσn ε,¯ σ) = (rot σn ε,rot ¯ σ)+(∇·σn ε,∇· ¯ σ)+(σn ε,¯ σ). 121
4.5.1 Mass conservation and Energy-stability Observe that the scheme UZSW is also conservative in u(satisfying (4.16)) and also has the behavior for RΩvngiven in (4.17). Definition 4.5.3 A numerical scheme with solution (un ε, zn ε,σn ε, wn ε)is called energy-stable with respect to the energy E(w, σ) = kwk2 0+1 2kσk2 0(4.57) if this energy is time decreasing, that is, E(wn ε,σn ε)≤E(wn−1 ε,σn−1 ε)for all n≥1. Theorem 4.5.4 (Unconditional stability) The scheme UZSW is unconditional energy stable with respect to E(w, σ). In fact, if (un ε, zn ε,σn ε, wn ε)is a solution of UZSW, then the following discrete energy law holds δtE(wn ε,σn ε) + kkδtwn εk2 0+k 2kδtσn εk2 0+ZΩ λε(un−1 ε)|∇zn ε|2+kσn εk2 1= 0.(4.58) Proof. The proof follows taking (¯z, ¯ σ,¯w, ¯u) = (zn ε,σn ε,2wn ε, δtun ε) in (4.56). From the (local in time) discrete energy law (4.58), we deduce the following global in time estimates for (un ε, zn ε,σn ε, wn ε) solution of the scheme UZSW: Corollary 4.5.5 (Uniform Weak estimates) Assume that (u0, v0)∈L2(Ω) ×H1(Ω). Let (un ε, zn ε,σn ε, wn ε)be a solution of scheme UZSW. Then, the following estimate holds kwn εk2 0+kσn εk2 0+k n X m=1 ZΩ λε(um−1 ε)|∇zm ε|2+kσm εk2 1≤C0,∀n≥1,(4.59) with the constant C0>0depending on the data (Ω, u0, v0), but independent of k, h, n and ε. Proof. Proceeding as in (4.27) and taking into account that u0≥0 and (u0 h,σ0 h, w0 h,ε) = (Qhu0,e Qh(∇v0),b Qh(pFε(u0(x)) + A)), we have that kw0 h,εk2 0+1 2kσ0 hk2 0=kb Qh(pFε(u0) + A)k2 0+1 2ke Qh(∇v0)k2 0≤ZΩ (Fε(u0) + A) + 1 2k∇v0k2 0 ≤CZΩ ((u0)2+ 1) + 1 2k∇v0k2 0≤C(ku0k2 0+kv0k2 1+ 1) ≤C0,(4.60) with the constant C0>0 depending on the data (Ω, u0, v0), but independent of k, h, n and ε. Therefore, multiplying the discrete energy law (4.58) by k, adding from m= 1 to m=n and using (4.60), we arrive at (4.59). 128
4.5.2 Well-posedness Theorem 4.5.6 (Unconditional unique solvability) There exists a unique (un ε, zn ε,σn ε, wn ε) solution of scheme UZSW. Proof. By linearity of the scheme UZSW, it suffices to prove uniqueness. Suppose that there exist (un ε,1, zn ε,1,σn ε,1, wn ε,1),(un ε,2, zn ε,2,σn ε,2, wn ε,2)∈Uh×Zh×Σh×Whtwo possible solutions of UZSW. Then defining u=un ε,1−un ε,2,z=zn ε,1−zn ε,2,σ=σn ε,1−σn ε,2and w=wn ε,1−wn ε,2, we have that (u, z, σ, w)∈Uh×Zh×Σh×Whsatisfies 1 k(u, ¯z)+(λε(un−1 ε)∇z, ∇¯z) = −(un−1 εσ,∇¯z),∀¯z∈Zh, 1 k(σ,¯ σ)+(Bhσ,¯ σ)=(un−1 ε∇z, ¯ σ),∀¯ σ∈Σh, 1 k(w, ¯w) = 1 2k1 √Fε(un−1 ε)+AF0 ε(un−1 ε)u, ¯w,∀¯w∈Wh, (z, ¯u) = 1 √Fε(un−1 ε)+AF0 ε(un−1 ε)w, ¯u,∀¯u∈Uh. (4.61) Taking (¯z, ¯ σ,¯w, ¯u) = (z, σ,2w, 1 ku) in (4.61) and adding, we obtain 2 kkwk2 0+1 kkσk2 0+ZΩ λε(un−1 ε)|∇z|2+kσk2 1= 0. Taking into account that λε(un−1 ε)≥ε, we deduce that (∇z, σ, w)=(0,0,0), hence z=C:= cte. Moreover, using the fact that w= 0 and z=C, from (4.61)4we conclude that z= 0. Finally, taking ¯z=uin (4.61)1(which is possible thanks to the choice Uh⊆Zh), since (∇z, σ) = (0,0) we conclude u= 0. 4.6 Numerical simulations The aim of this section is to compare the results of several numerical simulations using the schemes derived throughout the paper. We choose the spaces for (u, z, σ, w) generated by P1-continuous FE. Moreover, we have chosen the 2D domain [0,2]2with a structured mesh (then (H) holds and the scheme UV can be defined), and all the simulations are carried out using FreeFem++ software. In the comparison, we will also consider the classical Backward Euler scheme for model (4.1), which is given for the following first order in time, nonlinear and coupled scheme: •Scheme BEUV : Initialization: Let (u0, v0)=(Qhu0, Rhv0)∈Uh×Vh. 129
Time step n: Given (un−1, vn−1)∈Uh×Vh, compute (un, vn)∈Uh×Vhsolving ((δtun,¯u)+(∇un,∇¯u) = −(un∇vn,∇¯u),∀¯u∈Uh, (δtvn,¯v)+(Ahvn,¯v) = (un,¯v),∀¯v∈Vh. Remark 4.6.1 The scheme BEUV has not been analyzed in the previous sections because it is not clear how to prove its energy-stability. In fact, observe that the scheme UV (which is the “closest” approximation to the scheme BEUV considered in this paper) differs from the scheme BEUV in the use of the regularized functions Fε,F0 εand F00 ε(see (4.5) and Figure 4.1) and in the approximation of cross-diffusion term (u∇v, ∇¯u), which are crucial for the proof of the energy-stability of the scheme UV. The linear iterative methods used to approach the solutions of the nonlinear schemes UV,US and BEUV are the following Picard methods, in which, we denote (un ε, vn ε,σn ε) := (un, vn,σn). (i) Picard method to approach a solution (un, vn) of the scheme UV Initialization (l= 0): Set (u0, v0) = (un−1, vn−1)∈Uh×Vh. Algorithm: Given (ul, vl)∈Uh×Vh, compute (ul+1, vl+1)∈Uh×Vhsuch that (1 k(ul+1,¯u)h+ (∇ul+1,∇¯u) = 1 k(un−1,¯u)h−(Λε(ul)∇vl+1,∇¯u),∀¯u∈Uh, 1 k(vl+1,¯v)+(Ahvl+1,¯v) = 1 k(vn−1,¯v)+(ul,¯v),∀¯v∈Vh, until the stopping criteria max kul+1 −ulk0 kulk0 ,kvl+1 −vlk0 kvlk0≤tol. (ii) Picard method to approach a solution (un,σn) of the scheme US Initialization (l= 0): Set (u0,σ0)=(un−1,σn−1)∈Uh×Σh. Algorithm: Given (ul,σl)∈Uh×Σh, compute (ul+1,σl+1)∈Uh×Σhsuch that 1 k(ul+1,¯u)h+ (∇(ul+1,∇¯u)−(∇ul,∇¯u) =1 k(un−1,¯u)h−(λε(ul)∇Πh(F0 ε(ul)),∇¯u)−(λε(ul)σl+1,∇¯u),∀¯u∈Uh, 1 k(σl+1,¯ σ) + hBσl+1,¯ σi=1 k(σn−1,¯ σ)+(λε(ul)∇Πh(F0 ε(ul)),¯ σ),∀¯ σ∈Σh, until the stopping criteria max kul+1 −ulk0 kulk0 ,kσl+1 −σlk0 kσlk0≤tol. Note that a residual term (∇(ul+1 −ul),∇¯u) is considered. 130
(iii) Picard method to approach a solution (un, vn) of the scheme BEUV Initialization (l= 0): Set (u0, v0) = (un−1, vn−1)∈Uh×Vh. Algorithm: Given (ul, vl)∈Uh×Vh, compute (ul+1, vl+1)∈Uh×Vhsuch that (1 k(ul+1,¯u)+(∇ul+1,∇¯u) = 1 k(un−1,¯u)−(ul∇vl+1,∇¯u),∀¯u∈Uh, 1 k(vl+1,¯v)+(Ahvl+1,¯v) = 1 k(vn−1,¯v)+(ul,¯v),∀¯v∈Vh, until the stopping criteria max kul+1 −ulk0 kulk0 ,kvl+1 −vlk0 kvlk0≤tol. Remark 4.6.2 In all cases, first we compute vl+1 (resp. σl+1) solving the v-equation (resp. σ-system) and then, inserting vl+1 (resp. σl+1) in u-equation, we compute ul+1. 4.6.1 Positivity of un In this subsection, we compare the positivity of the variable un∈Uhin the four schemes. Here, we choose the space Vhgenerated by P2-continuous FE. We recall that for the three schemes studied in this paper, namely schemes UV,UZSW and US, it is not clear the positivity of the variable un. Moreover, for the schemes UV and US, it was proved that Πh(un ε−)→0 in L2(Ω) as ε→0 (see Remarks 4.3.10 and 4.4.4); while for the scheme UZSW this fact is not clear. For this reason, in Figures 4.3-4.5 we compare the positivity of the variable un εin the schemes, taking ε= 10−3,ε= 10−5and ε= 10−8. In the scheme UZSW we fix A= 1 (and thus, Fε(s) + A≥1 for all s∈R). We consider the time step k= 10−5, the tolerance parameter for the linear iterative methods tol = 10−4and the initial conditions (see Figure 4.2) u0=−10xy(2 −x)(2 −y)exp(−10(y−1)2−10(x−1)2) + 10.0001, v0=100xy(2 −x)(2 −y)exp(−30(y−1)2−30(x−1)2)+0.0001. Note that u0, v0>0 in Ω, min(u0) = u0(1,1) = 0.0001 and max(v0) = v0(1,1) = 100.0001. Moreover, for the schemes UV and UZSW we take the mesh size h=1 40, while for the scheme US it was necessary to take h=1 80, because for thicker meshes we had convergence problems of the iterative method. In the case of the schemes UV and US, we observe that although un εis negative for some x∈Ω in some times tn>0, when ε→0 these values are closer to 0; while in the case of 131
(a) Initial cell density u0(b) Initial chemical concentration v0 Figure 4.2: Initial conditions. the scheme UZSW, this same behavior is not observed (see Figures 4.3-4.5). Finally, in the case of the scheme BEUV (see Figure 4.6), we have also observed negative values for the minimum of unin some times tn>0, with more negative values than in the schemes UV and US. Remark 4.6.3 In Figures 4.3 and 4.5 there are also negative values of minimum of un εfor ε= 10−5and ε= 10−8, but those are of order 10−4and 10−7respectively in both figures. 4.6.2 Energy stability In this subsection, we compare numerically the stability of the schemes UV,UZSW,US and BEUV with respect to the “exact” energy Ee(u, v) = ZΩ F0(u(x))dx+1 2k∇vk2 0,(4.62) where F0(u) := F(u+) = 1,if u≤0, uln(u)−u+ 1,if u > 0. It was proved that the schemes UV,UZSW and US are unconditionally energy-stables with respect to modified energies obtained in terms of the variables of each scheme. Even more, some energy inequalities are satisfied (see Theorems 4.3.7,4.4.2 and 4.5.4). However, 132
0 0.005 0.01 0.015 Time -0.025 -0.02 -0.015 -0.01 -0.005 0 0.005 Min un ǫ ǫ=10-3 ǫ=10-5 ǫ=10-8 Figure 4.3: Minimum values of un εcomputed using the scheme UV. it is not clear how to prove the energy-stability of these schemes with respect to the “exact” energy Ee(u, v) given in (4.62), which comes from the continuous problem (4.1) (see (4.3)). Therefore, it is interesting to compare numerically the schemes with respect to this energy Ee(un, vn), and to study the behaviour of the corresponding discrete residual of the energy law (4.3): REe(un, vn) := δtEe(un, vn)+4ZΩ|∇pun +|2dx+k(Ah−I)vnk2 0+k∇vnk2 0.(4.63) 1. First test: We consider k= 10−3,h=1 40,tol = 10−4and the initial conditions (see Figure 4.7) u0= 7w+ 7.0001 and v0=−7w+ 7.0001, where w:= cos(2πx)cos(2πy). We choose Vhgenerated by P2-continuous FE. Then, we obtain that: (i) The scheme BEUV satisfies the energy decreasing in time property for the exact energy Ee(u, v), that is, Ee(un, vn)≤Ee(un−1, vn−1)∀n. (4.64) Its behaviour can be observed in Figure 4.8. The same behaviour is obtained for the schemes UV and US independently of the choice of ε. In the case of the scheme 133
0 0.005 0.01 0.015 Time -0.35 -0.3 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 Min un ǫ ǫ=10-3 ǫ=10-5 ǫ=10-8 Figure 4.4: Minimum values of un εcomputed using the scheme UZSW. UZSW, this property (4.64) is not satisfied for any value of ε. Indeed, increasing energies are obtained for different values of ε(see Figure 4.9). (ii) The scheme BEUV satisfies the discrete energy inequality REe(un, vn)≤0 for REe(un, vn) defined in (4.63) (see Figure 4.10). The same is observed for the schemes UV and US independently of the choice of ε. In the case of the scheme UZSW, it is observed that this discrete energy inequality is not satisfied for any value of ε. Indeed, the residual REe(un ε, vn ε) obtained for each εreaches very large positive values (see Figure 4.11). 2. Second test: We consider k= 10−5,h=1 20,tol = 10−4and the initial conditions u0= 14w+ 14.0001 and v0=−14w+ 14.0001, with the function was before. Now, we choose the space Vhgenerated by P1-continuous FE. Then, we obtain that: (i) The schemes BEUV,UV and US satisfy the energy decreasing in time property (4.64), independently of the choice of ε. (ii) The schemes UV and US satisfy the discrete energy inequality REe(un ε, vn ε)≤0, independently of the choice of ε; while the scheme BEUV have RE(un, vn)>0 for some n≥0 (see Figure 4.12). 134
0 0.005 0.01 0.015 Time -0.025 -0.02 -0.015 -0.01 -0.005 0 0.005 0.01 Min un ǫ ǫ=10-3 ǫ=10-5 ǫ=10-8 Figure 4.5: Minimum values of un εcomputed using the scheme US. 4.7 Conclusions In this paper we have developed three new mass-conservative and unconditionally energystable fully discrete FE schemes for the chemorepulsion production model (4.1), namely UV, US and UZSW. From the theoretical point of view we have obtained: (i) The well-posedness of the numerical schemes (with conditional uniqueness for the nonlinear schemes UV and US). (ii) The nonlinear scheme UV is unconditional energy-stable with respect to the energy Eh ε(u, v) given in (4.19), under the constraint (H) on the space triangulation related with the right-angles and assuming that Uhis approximated by P1-continuous FE. (iii) The nonlinear scheme US and the linear scheme UZSW are unconditional energystables with respect to the modified energies Eh ε(u, σ) (given in (4.44)) and E(w, σ) (given in (4.57)) respectively, without the constraint on the triangulation related with the right-angles simplices and assuming that Uhcan be approximated by P1-continuous and Pk-continuous FE respectively, for any k≥1. (iv) It is not clear how to prove the energy-stability of the nonlinear scheme BEUV with respect to the energy Ee(u, v) (given in (4.62)) or some modified energy (see Remark 4.6.1). 135
0 0.005 0.01 0.015 Time -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 Min un h Figure 4.6: Minimum values of uncomputed using the scheme BEUV. (v) In the schemes UV and US there is a control for Πh(un ε−) in L2-norm, which tends to 0 as ε→0. This allows to conclude the nonnegativity of the solution un εin the limit when ε→0. This property is not clear for the linear scheme UZSW. On the other hand, from the numerical simulations, we can conclude: (i) There are initial conditions for which the scheme UZSW is not energy stable with respect to the energy Ee(u, v), that is, the decreasing in time property (4.64) is not satisfied for any value of ε. Indeed, time increasing energies are obtained for different values of ε. (ii) For the three compared nonlinear schemes (UV,US and BEUV), only the scheme US has convergence problems for the linear iterative method. However, these problems are overcomed considering thinner meshes. (iii) The schemes UV and US have decreasing in time energy Ee(u, v), independently of the choice of ε. In fact, the discrete energy inequality REe(un ε, vn ε)≤0 is satisfied in all cases, for REe(un ε, vn ε) defined in (4.63). (iv) The scheme BEUV has decreasing in time energy Ee(u, v), but the discrete energy inequality REe(un, vn)≤0 is not satisfied for some n≥0. 136
(a) Initial cell density u0(b) Initial chemical concentration v0 Figure 4.7: Initial conditions. (v) Finally, it was observed numerically that, for the schemes UV and US,un ε−→0 as ε→0; while for the scheme UZSW this behavior was not observed. 137