Analysis of a fully discrete approximation for the classical Keller–Segel model: Lower and a priori bounds

Computers and Mathematics with Applications 85 (2021) 69–81 Contents lists available at ScienceDirect Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa Analysis of a fully discrete approximation for the classical Keller–Segel model: Lower and a priori bounds Juan Vicente Gutiérrez-Santacreu a,∗,1, José Rafael Rodríguez-Galvánb,1 aDpto. de Matemática Aplicada I, E. T. S. I. Informática, Universidad de Sevilla, Avda. Reina Mercedes, s/n., E-41012 Sevilla, Spain bDepartamento de Matemáticas. Facultad de Ciencias. Campus Universitario de Puerto Real, Universidad de Cádiz. E-11510 Puerto Real, Cádiz, Spain ARTICLE INFO Keywords: Keller–Segel equations Non-linear parabolic equations Finite-element approximation Lower bounds A priori bounds ABSTRACT This paper is devoted to constructing approximate solutions for the classical Keller–Segel model governing chemotaxis. It consists of a system of nonlinear parabolic equations, where the unknowns are the average density of cells (or organisms), which is a conserved variable, and the average density of chemoattractant. The numerical proposal is made up of a crude finite element method together with a mass lumping technique and a semi-implicit Euler time integration. The resulting scheme turns out to be linear and decouples the computation of variables. The approximate solutions keep lower bounds – positivity for the cell density and nonnegativity for the chemoattractant density–, are bounded in the 𝐿1(𝛺)-norm, satisfy a discrete energy law, and have a priori energy estimates. The latter is achieved by means of a discrete Moser–Trudinger inequality. As far as we know, our numerical method is the first one that can be encountered in the literature dealing with all of the previously mentioned properties at the same time. Furthermore, some numerical examples are carried out to support and complement the theoretical results. 1. Introduction 1.1. Aims In 1970/71 Keller and Segel [1,2] attempted to derive a set of equations for modeling chemotaxis – a biological process through which an organism (or a cell) migrates in response to a chemical stimulus being attractant or repellent. It is nowadays well-known that the work of Keller and Segel turned out to be somehow biologically inaccurate since their equations provide unrealistic solutions; a little more precisely, solutions that blow up in finite time. Such a phenomenon does not occur in nature. Even though the original Keller–Segel equations are less relevant from a biological point of view, they are mathematically of great interest. Much of work for the Keller–Segel equations has been carried out in developing purely analytical results, whereas there are very few numerical results in the literature. This is due to the fact that solving numerically the Keller–Segel equations is a challenging task because their solutions exhibit many interesting mathematical properties which are not easily adapted to a discrete framework. For instance, solutions to the Keller–Segel equations satisfy lower bounds (positivity and nonnegativity) and enjoy an energy law, which is obtained by testing ∗Corresponding author. E-mail addresses: [email protected] (J.V. Gutiérrez-Santacreu), [email protected] (J.R. Rodríguez-Galván). 1JVGS and JRRG were partially supported by the Spanish Grant No. PGC2018-098308-B-I00 from Ministerio de Ciencias e Innovación - Agencia Estatal de Investigación, Spain with the participation of FEDER. the equations against non linear functions. Cross-diffusion mechanisms governing the chemotactic phenomena are the responsible for the fact that the Keller–Segel equations are so difficult to analyze not only theoretically but also numerically. In spite of being a limited model, it is hoped that developing and analyzing numerical methods for the classical Keller–Segel equations may open new roads to deeper insights and better understandings for dealing with the numerical approximation of other chemotaxis models — biologically more realistic, but which are, however, inspired on the original Keller–Segel formulation. In a nutshell, these other chemotaxis models are modifications of the Keller–Segel equations in order to avoid the non-physical blow up of solutions and hence produce solutions being closer to chemotaxis phenomena. For these other chemotaxis models, it is recommended the excellent surveys of Hillen and Painter [3], Horstamann [4,5], and, more recently, Bellomo, Bellouquid, Tao, and Winkler [6]. In these surveys the authors reviewed to date as to when they were written the state of art of modeling and mathematical analysis for the Keller–Segel equations and their variants. It is our aim in this work to design a fully discrete algorithm for the classical Keller–Segel equations based on a finite element discretization whose discrete solutions satisfy lower and a priori bounds. https://doi.org/10.1016/j.camwa.2021.01.009 Received 22 September 2020; Received in revised form 5 January 2021; Accepted 8 January 2021 Available online 3 February 2021 0898-1221/©2021 Elsevier Ltd. All rights reserved. J.V. Gutiérrez-Santacreu and J.R. Rodríguez-Galván Computers and Mathematics with Applications 85 (2021) 69–81 1.2. The Keller–Segel equations Let 𝛺 ⊂ R2be a bounded domain, with 𝒏being its outwarddirected unit normal vector to 𝛺, and let [0, 𝑇 ]be a time interval. Take 𝑄= (0, 𝑇 ] × 𝛺and 𝛴= (0, 𝑇 ] × 𝜕𝛺. Then the boundary-value problem for the Keller–Segel equations reads a follows. Find 𝑢∶ 𝑄→(0,∞) and 𝑣∶ 𝑄→[0,∞) satisfying {𝜕𝑡𝑢−𝛥𝑢 = −∇ ⋅(𝑢∇𝑣)in 𝑄, 𝜕𝑡𝑣−𝛥𝑣 =𝑢−𝑣in 𝑄, (1) subject to the initial conditions 𝑢(0) = 𝑢0and 𝑣(0) = 𝑣0in 𝛺, (2) and the boundary conditions ∇𝑢⋅𝒏= 0 and ∇𝑣⋅𝒏= 0 on 𝛴. (3) Here 𝑢is the average density of organisms (or cells), which is a conserved variable, and 𝑣is the average density of chemical sign, which is a nonconserved variable. System (1) was motivated by Keller and Segel [1] describing the aggregation phenomena exhibited by the amoeba Dictyostelium discoideum due to an attractive chemical substance referred to as chemoattractant, which is generated by the own amoeba and is, nevertheless, degraded by living conditions. Moreover, diffusion is also presented in the motion of amebae and chemoattractant. The diffusion phenomena performed by cells and chemoattractant are modeled by the terms −𝛥𝑢 and −𝛥𝑣, respectively, whereas the aggregation mechanism is described by the term −∇ ⋅(𝑢∇𝑣). It is this nonlinear term that is the major difficulty in studying system (1). Further the production and degradation of chemoattractant are associated with the term 𝑢−𝑣. Concerning the mathematical analysis for system (1), Nagai, Senba, and Yoshida [7] proved existence, uniqueness and regularity of solutions under the condition ∫𝛺𝑢0(𝒙) d𝒙∈ (0,4𝜋). In proving this, a variant of Moser–Tridinguer’s inequality was used. In the particular case that 𝛺 be a ball, the above-mentioned condition becomes ∫𝛺𝑢0(𝒙) d𝒙∈ (0,8𝜋). Herrero and Velázquez [8] dealt with the first blow-up framework by constructing some radially symmetric two-dimensional solutions which blow up within finite time. The next progress in this sense with 𝛺 being non-radial and simply connected was the work of Horstmann and Wang [9] who found some unbounded solutions provided that ∫𝛺𝑢0(𝒙) d𝒙>4𝜋and ∫𝛺𝑢0(𝒙) d𝒙∉ {4𝑘𝜋 |𝑘∈N}. So far there is no supporting evidence as to whether such solutions may evolve to produce a blow-up phenomenon within finite time or whether, on the contrary, may increase to infinity with time. In three dimensions, Winkler [10] proved that there exist radially symmetric solutions blowing up in finite time for any value of ∫𝛺𝑢0(𝒙) d𝒙. The main tool [10] in proving blow-up solutions is the energy law which stems from system (1). Nevertheless, an inadequate approximation of lower bounds can trigger off oscillations of the variables, which can lead to spurious, blow-up solutions. Concerning the numerical analysis for system (1), very little is said about numerical algorithms which keep lower bounds, are 𝐿1(𝛺)- bounded and have a discrete energy law. Proper numerical treatment of these properties is made difficult by the fact that the non-linearity occurs in the highest order derivative. Numerical algorithms are mainly designed so as to keep lower bounds and to be mass-preserving. We refer the reader to [11–16]. We were pointed out by a referee the paper [17]. In it, the authors rewrote system (1) by using several ad hoc auxiliary variables so that the resulting discontinuous Galerkin method fulfilled a discrete energy law, but no lower bounds were proved. As far as we are concerned, there is no numerical method facing lower bounds as well as a discrete energy law. 1.3. Notation We collect here as a reference some standard notation used throughout the paper. For 𝑝∈ [1,∞], we denote by 𝐿𝑝(𝛺)the usual Lebesgue space, i.e., 𝐿𝑝(𝛺)={𝑣∶𝛺→R∶𝑣Lebesgue-measurable,∫𝛺|𝑣(𝒙)|𝑝d𝒙<∞}. or 𝐿∞(𝛺)={𝑣∶𝛺→R∶𝑣Lebesgue-measurable,ess sup 𝒙∈𝛺|𝑣(𝒙)|<∞}. This space is a Banach space endowed with the norm ‖𝑣‖𝐿𝑝(𝛺)= (∫𝛺|𝑣(𝒙)|𝑝d𝒙)1∕𝑝if 𝑝∈ [1,∞) or ‖𝑣‖𝐿∞(𝛺)= ess sup𝒙∈𝛺|𝑣(𝒙)|if 𝑝= ∞. In particular, 𝐿2(𝛺)is a Hilbert space. We shall use (𝑢, 𝑣)=∫𝛺𝑢(𝒙)𝑣(𝒙)d𝒙 for its inner product and ‖⋅‖for its norm. Let 𝛼= (𝛼1, 𝛼2) ∈ N2be a multi-index with |𝛼|=𝛼1+𝛼2, and let 𝜕𝛼 be the differential operator such that 𝜕𝛼=(𝜕 𝜕𝑥1)𝛼1(𝜕 𝜕𝑥2)𝛼2. For 𝑚≥0and 𝑝∈ [1,∞), we consider 𝑊𝑚,𝑝(𝛺)to be the Sobolev space of all functions whose 𝑚derivatives are in 𝐿𝑝(𝛺), i.e., 𝑊𝑚,𝑝(𝛺)={𝑣∈𝐿𝑝(𝛺) ∶ 𝜕𝑘𝑣∈𝐿2(𝛺) ∀ |𝑘|≤𝑚} associated to the norm ‖𝑓‖𝑊𝑚,𝑝(𝛺)=(∑ |𝛼|≤𝑚‖𝜕𝛼𝑓‖𝑝 𝐿𝑝(𝛺))1∕𝑝 for 1≤𝑝 < ∞, and ‖𝑓‖𝑊𝑚,𝑝(𝛺)= max |𝛼|≤𝑚‖𝜕𝛼𝑓‖𝐿∞(𝛺)for 𝑝= ∞. For 𝑝= 2, we denote 𝑊𝑚,2(𝛺) = 𝐻𝑚(𝛺). Moreover, we make of use the space 𝐻2 𝑁(𝛺) = {𝑣∈𝐻2(𝛺) ∶ ∫𝛺 𝑣(𝒙) d𝒙= 0 and 𝜕𝒏𝑣= 0 on 𝜕𝛺}, for which is known that ‖𝑣‖𝐻2 𝑁(𝛺)and ‖𝛥𝑣‖are equivalent norms. 1.4. Outline The remainder of this paper is organized in the following way. In the next section we state our finite element space and some tools. In particular, we prove a discrete version of a variant of Moser– Trudinger’s inequality. In Section 3, we apply our ideas to discretize system (1) in space and time for defining our numerical method and formulate our main result. Next is Section 4dedicated to demonstrating lower bounds, a discrete energy law, and a priori bounds all of which are local in time for approximate solutions. This is accomplished in a series of lemmas where the final argument is an induction procedure on the time step so as to obtain the above mentioned properties globally in time. Finally, in Section 5, we consider two numerical examples regarding blow-up and non blow-up scenarios. 2. Technical preliminaries This section is mainly devoted to setting out the hypotheses and some auxiliary results concerning the finite element space that will use throughout this work. 2.1. Hypotheses We construct the finite element approximation of (1) under the following assumptions on the domain, the mesh, and the finite element space. 70 J.V. Gutiérrez-Santacreu and J.R. Rodríguez-Galván Computers and Mathematics with Applications 85 (2021) 69–81 (H1) Let 𝛺be a convex, bounded domain of R2with a polygonal boundary, and let 𝜃𝛺be the minimum interior angle at the vertices of 𝜕𝛺. (H2) Let {ℎ}ℎ>0be a family of acute, shape-regular, quasi-uniform triangulations of 𝛺made up of triangles, so that 𝛺= ∪𝑇∈ℎ𝑇, where ℎ= max𝑇∈ℎℎ𝑇, with ℎ𝑇being the diameter of 𝑇. More precisely, we assume that (a) there exists 𝛼 > 0, independent of ℎ, such that min{diam 𝐵𝑇∶𝑇∈ℎ}≥𝛼ℎ, where 𝐵𝑇is the largest ball contained in 𝑇, and (b) there exists 𝛽 > 0such that every angle between two edges of a triangle 𝑇is bounded by 𝜋 2−𝛽. Further, let ℎ= {𝒂𝑖}𝑖∈𝐼be the coordinates of the nodes of ℎ. (H3) Associated with ℎis the finite element space 𝑋ℎ={𝑥ℎ∈𝐶0(𝛺) ∶ 𝑥ℎ|𝑇∈1(𝑇),∀𝑇∈ℎ}, where 1(𝑇)is the set of linear polynomials on 𝑇. Let {𝜑𝒂}𝒂∈ℎ be the standard basis functions for 𝑋ℎ. 2.2. Auxiliary results Our first result is concerned with the sign of the entries of the rigid matrix. Proposition 2.1. Let 𝛺be a polygonal. Consider 𝑋ℎto be constructed over ℎbeing acute. Then, for each 𝑇∈ℎwith vertices {𝒂1,𝒂2,𝒂3}, there exists a constant 𝐶neg >0, depending on 𝛽, but otherwise independent of ℎ and 𝑇, such that ∫𝑇 ∇𝜑𝒂𝑖⋅∇𝜑𝒂𝑗d𝒙≤−𝐶neg (4) for all 𝒂𝑖,𝒂𝑗∈𝑇with 𝑖≠𝑗, and ∫𝑇 ∇𝜑𝒂𝑖⋅∇𝜑𝒂𝑖d𝒙≥𝐶neg (5) for all 𝒂𝑖∈𝑇. A proof of (4) and (5) can be found in [18]. Both local and global finite element properties for 𝑋ℎwill be needed such as inverse estimates and bounds for the interpolation error. We first recall some local inverse estimates. See [19, Lem. 4.5.3] or [20, Lem. 1.138] for a proof. Proposition 2.2. Let 𝛺be polygonal. Consider 𝑋ℎto be constructed over ℎbeing quasi-uniform. Then, for each 𝑇∈ℎand 𝑝∈ [2,∞], there exists a constant 𝐶inv >0, independent of ℎand 𝑇, such that, for all 𝑥ℎ∈𝑋ℎ, ‖∇𝑥ℎ‖𝐿𝑝(𝑇)≤𝐶inv ℎ−1‖𝑥ℎ‖𝐿𝑝(𝑇)(6) and ‖∇𝑥ℎ‖𝐿∞(𝑇)≤𝐶inv ℎ−2 𝑝‖∇𝑥ℎ‖𝐿𝑝(𝑇).(7) Concerning global inverse inequalities, we need the following. Proposition 2.3. Let 𝛺be polygonal. Consider 𝑋ℎto be constructed over ℎbeing quasi-uniform. Then for each 𝑝∈ [2,∞], there exists a constant 𝐶inv >0, independent of ℎ, such that, for all 𝑥ℎ∈𝑋ℎ, ‖𝑥ℎ‖𝐿∞(𝛺)≤𝐶inv ℎ−1‖𝑥ℎ‖,(8) ‖∇𝑥ℎ‖𝐿𝑝(𝛺)≤𝐶inv ℎ−1‖𝑥ℎ‖𝐿𝑝(𝛺),(9) ‖∇𝑥ℎ‖𝐿𝑝(𝛺)≤𝐶inv ℎ−2( 1 2−1 𝑝)‖∇𝑥ℎ‖𝐿2(𝛺),(10) and ‖∇𝑥ℎ‖𝐿∞(𝛺)≤𝐶inv ℎ−2 𝑝‖∇𝑥ℎ‖𝐿𝑝(𝛺).(11) We introduce ℎ∶𝐶(𝛺)→𝑋ℎ, the standard nodal interpolation operator, such that ℎ𝜂(𝒂) = 𝜂(𝒂)for all 𝒂∈ℎ. Associated with ℎ, a discrete inner product on 𝑋ℎis defined by (𝑥ℎ, 𝑥ℎ)ℎ=∫𝛺 ℎ(𝑥ℎ(𝒙)𝑥ℎ(𝒙)) d𝒙. We also introduce ‖𝑥ℎ‖ℎ= (𝑥ℎ, 𝑥ℎ) 1 2 ℎ. Local and global error bounds for ℎare as follows (c.f. [19, Thm. 4.4.4] or [20, Thm. 1.103] for a proof). Proposition 2.4. Let 𝛺be polygonal. Consider 𝑋ℎto be constructed over ℎbeing quasi-uniform. Then, for each 𝑇∈ℎ, there exists 𝐶app >0, independent of ℎand 𝑇, such that ‖𝜑−ℎ𝜑‖𝐿1(𝑇)≤𝐶appℎ2‖∇2𝜑‖𝐿1(𝑇)∀𝜑∈𝑊2,1(𝑇).(12) Proposition 2.5. Let 𝛺be polygonal. Consider 𝑋ℎto be constructed over ℎbeing quasi-uniform. Then it follows that there exists 𝐶app >0, independent of ℎ, such that ‖∇(𝜑−ℎ𝜑)‖𝐿2(𝛺)≤𝐶appℎ‖𝛥𝜑‖𝐿2(𝛺)∀𝜑∈𝐻2(𝛺).(13) Corollary 2.6. Let 𝛺be polygonal. Consider 𝑋ℎto be constructed over ℎbeing quasi-uniform. Let 𝑛∈N. Then it follows that there exist three positive constants 𝐶app,𝐶com, and 𝐶sta, independent of ℎ, such that ‖𝑥𝑛 ℎ−ℎ(𝑥𝑛 ℎ)‖𝐿1(𝛺)≤𝐶app𝑛(𝑛− 1)ℎ2∫𝛺|𝑥ℎ(𝒙)|𝑛−2|∇𝑥ℎ(𝒙)|2d𝒙,(14) ‖𝑥ℎ𝑥ℎ−ℎ(𝑥ℎ𝑥ℎ)‖𝐿1(𝛺)≤𝐶comℎ‖𝑥ℎ‖𝐿2(𝛺)‖∇𝑥ℎ‖𝐿2(𝛺)(15) and ‖𝑥𝑛 ℎ‖𝐿1(𝛺)≤‖ℎ(𝑥𝑛 ℎ)‖𝐿1(𝛺)≤𝐶sta‖𝑥𝑛 ℎ‖𝐿1(𝛺),(16) where 𝐶sta depends on 𝑛. Proof. Let 𝑇∈ℎand compute ∇2(𝑥𝑛 ℎ) = 𝑛(𝑛− 1)𝑥𝑛−2 ℎ 𝑑 ∑ 𝑖,𝑗=1 𝜕𝑖𝑥ℎ𝜕𝑗𝑥ℎon 𝑇 . Then, from (12) and the above identity, we have ‖𝑥𝑛 ℎ−ℎ(𝑥𝑘 ℎ)‖𝐿1(𝑇)≤𝐶appℎ2 𝐾‖∇2(𝑥𝑛 ℎ)‖𝐿1(𝑇) ≤𝐶app𝑛(𝑛− 1)ℎ2∫𝑇|𝑥ℎ(𝒙)|𝑛−2|∇𝑥ℎ(𝒙)|2d𝒙. Summing over 𝑇∈ℎyields (14). The proof of (15) follows very closely the arguments of (14) for 𝑛= 2. The first part of assertion (16) is a simple application of Jensen’s inequality, whereas the second part follows from (14) on using Hölder’s inequality, (9) for 𝑝=𝑛and, later on, reverse Minkowski’s inequality. □ The proof of the following proposition can be found in [21]. It is a generalization of a Moser–Trudinger-type inequality. Proposition 2.7 (Moser–Trudinger).Let 𝛺be polygonal with 𝜃𝛺being the minimum interior angle at the vertices of 𝛺. Then there exists a constant 𝐶𝛺>0depending on 𝛺such that for all 𝑢∈𝐻1(𝛺)with ‖∇𝑢‖≤1and ∫𝛺𝑢(𝒙) d𝒙= 1, it follows that ∫𝛺 𝑒𝛼|𝑢(𝒙)|2d𝒙≤𝐶𝛺,(17) where 𝛼≤2𝜃𝛺. Corollary 2.8. Let 𝛺be polygonal with 𝜃𝛺being the minimum interior angle at the vertices of 𝜕𝛺. Consider 𝑋ℎto be constructed over ℎbeing 71 J.V. Gutiérrez-Santacreu and J.R. Rodríguez-Galván Computers and Mathematics with Applications 85 (2021) 69–81 quasi-uniform. Let 𝑢ℎ∈𝑋ℎwith 𝑢ℎ>0. Then it follows that there exists a constant 𝐶MT >0, independent of ℎ, such that ∫𝛺 ℎ(𝑒𝑢ℎ(𝒙)) d𝒙≤𝐶𝛺(1 + 𝐶MT‖∇𝑢ℎ‖2)𝑒 1 8𝜃𝛺‖∇𝑢ℎ‖2+1 |𝛺|‖𝑢ℎ‖𝐿1(𝛺). (18) Proof. From (14), we have ∫𝛺 ℎ(𝑒𝑢ℎ(𝒙)) d𝒙=∫𝛺 (1 + 𝑢ℎ(𝒙)) d𝒙+ ∞ ∑ 𝑛=2 1 𝑛!∫𝛺 ℎ(𝑢𝑛 ℎ(𝒙))d𝒙 ≤∞ ∑ 𝑛=0 1 𝑛!∫𝛺 𝑢𝑛 ℎ(𝒙) d𝒙 + ∞ ∑ 𝑛=2 𝐶app𝑛(𝑛− 1)ℎ2 𝑛!∫𝛺|∇𝑢ℎ(𝒙)|2𝑢𝑛−2 ℎ(𝒙) d𝒙 =∫𝛺 (1 + 𝐶appℎ2|∇𝑢ℎ(𝒙)|2)𝑒𝑢ℎ(𝒙)d𝒙. (19) Let 𝑣ℎ=𝑢ℎ−𝑚 ‖∇𝑢ℎ‖with 𝑚=1 |𝛺|∫𝛺𝑢ℎ(𝒙) d𝒙. Young’s inequality gives 𝑢ℎ=‖∇𝑢ℎ‖𝑣ℎ+𝑚≤1 8𝜃𝛺‖∇𝑢ℎ‖2+ 2𝜃𝛺|𝑣ℎ|2+𝑚. (20) Thus, combining (19) and (20) yields, on noting (11) for 𝑝= 2 and (17), that ∫𝛺 ℎ(𝑒𝑢ℎ(𝒙)) d𝒙≤𝑒 1 8𝜃𝛺‖∇𝑢ℎ‖2+𝑚∫𝛺 (1 + 𝐶appℎ2‖∇𝑢ℎ(𝒙)‖2|∇𝑣ℎ(𝒙)|2) 𝑒2𝜃𝛺|𝑣ℎ(𝒙)|2d𝒙 ≤𝑒 1 8𝜃𝛺‖∇𝑢ℎ‖2+𝑚(1 + 𝐶app𝐶inv‖∇𝑢ℎ(𝒙)‖2) ×∫𝛺 𝑒2𝜃𝛺|𝑣ℎ(𝒙)|2d𝒙 ≤𝐶𝛺(1 + 𝐶app𝐶inv‖∇𝑢ℎ(𝒙)‖2)𝑒 1 8𝜃𝛺‖∇𝑢ℎ‖2+𝑚.□ An (average) interpolation operator into 𝑋ℎwill be required in order to properly initialize our numerical method. We refer to [22,23]. Proposition 2.9. Let 𝛺be polygonal. Consider 𝑋ℎto be constructed over ℎbeing quasi-uniform. Then there exists an (average) interpolation operator ℎfrom 𝐿1(𝛺)to 𝑋ℎsuch that ‖ℎ𝜓‖𝑊𝑠,𝑝(𝛺)≤𝐶sta‖𝜓‖𝑊𝑠,𝑝(𝛺)for 𝑠= 0,1and 1≤𝑝≤∞,(21) and ‖ℎ(𝜓) − 𝜓‖𝑊𝑠,𝑝(𝛺)≤𝐶appℎ1+𝑚−𝑠‖𝜓‖𝑊𝑚+1,𝑝(𝛺)for 0≤𝑠≤𝑚≤1.(22) Moreover, let − 𝛥ℎbe defined from 𝑋ℎto 𝑋ℎas − (  𝛥ℎ𝑥ℎ, 𝑥ℎ)ℎ= (∇𝑥ℎ,∇𝑥ℎ)for all 𝑥ℎ∈𝑋ℎ,(23) and let 𝑥(ℎ) ∈ 𝐻2 𝑁(𝛺)be such that {−𝛥𝑥(ℎ)=− 𝛥ℎ𝑥ℎin 𝛺, 𝜕𝒏𝑥(ℎ)=0 on 𝜕𝛺. (24) From elliptic regularity theory, the well-posedness of (24) is ensured by the convexity assumption stated in (H1) and ‖𝑥(ℎ)‖𝐻2 𝑁(𝛺)≤𝐶‖− 𝛥ℎ𝑥ℎ‖.(25) See [24] for a proof. Proposition 2.10. Let 𝛺be a convex polygon. Consider 𝑋ℎto be constructed over ℎbeing quasi-uniform. Then there exists a constant 𝐶Lap >0, independent of ℎ, such that ‖∇(𝑥(ℎ) − 𝑥ℎ)‖𝐿2(𝛺)≤𝐶Lapℎ‖ 𝛥ℎ𝑥ℎ‖𝐿2(𝛺).(26) Proof. We refer the reader to [18] for a proof which uses (13) and (15).□ Corollary 2.11. Let 𝛺be a convex polygon. Consider 𝑋ℎto be constructed over ℎbeing quasi-uniform. Then, for each 𝑝∈ [2,∞], there exists a constant 𝐶sta >0, independent of ℎ, such that ‖∇𝑥ℎ‖𝐿𝑝(𝛺)≤𝐶sta‖− 𝛥ℎ𝑥ℎ‖.(27) Proof. The triangle inequality gives ‖∇𝑥ℎ‖𝐿𝑝(𝛺)≤‖∇𝑥ℎ− ∇ℎ𝑥(ℎ)‖𝐿𝑝(𝛺)+‖∇ℎ𝑥(ℎ)‖𝐿𝑝(𝛺) and hence applying (10),(26),(25),(22),(21), and Sobolev’s inequality yields (27).□ 3. Presentation of main result We now define our numerical approximation of system (1). Assume that (𝑢0, 𝑣0) ∈ 𝐻1(𝛺) × 𝐻2(𝛺)with 𝑢0>0and 𝑣0≥0a. e. in 𝛺. We begin by approximating the initial data (𝑢0, 𝑣0)by (𝑢0 ℎ, 𝑣0 ℎ) ∈ 𝑋2 ℎ as follows. Define 𝑢0 ℎ=ℎ𝑢0,(28) which satisfies 𝑢0 ℎ>0a. e. in 𝛺, ‖𝑢0 ℎ‖𝐿1(𝛺)≤𝐶sta‖𝑢0‖𝐿1(𝛺)and ‖𝑢0 ℎ‖≤𝐶sta‖𝑢0‖ (29) and 𝑣0 ℎ=ℎ𝑣0,(30) which satisfies 𝑣0 ℎ≥0a. e. in 𝛺, ‖𝑣0 ℎ‖𝐻1(𝛺)≤𝐶‖𝑣0‖𝐻1(𝛺)and ‖ 𝛥ℎ𝑣0 ℎ‖≤𝐶‖𝛥𝑣0‖.(31) Given 𝑁∈N, we let 0 = 𝑡0< 𝑡1<⋯< 𝑡𝑁−1 < 𝑡𝑁=𝑇be a uniform partitioning of [0,T] with time step 𝑘=𝑇 𝑁. To simplify the notation we define the time-increment operator 𝛿𝑡𝜙𝑛+1 ℎ=𝜙𝑛+1 ℎ−𝜙𝑛 ℎ 𝑘. Known (𝑢𝑛 ℎ, 𝑣𝑛 ℎ) ∈ 𝑋ℎ×𝑋ℎ, find (𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ) ∈ 𝑋ℎ×𝑋ℎsuch that (𝛿𝑡𝑢𝑛+1 ℎ, 𝑥ℎ)ℎ+ (∇𝑢𝑛+1 ℎ,∇𝑥ℎ) = (∇𝑣𝑛 ℎ, 𝑢𝑛+1 ℎ∇𝑥ℎ)(32) and (𝛿𝑡𝑣𝑛+1 ℎ, 𝑥ℎ)ℎ+ (∇𝑣𝑛+1 ℎ,∇𝑥ℎ)+(𝑣𝑛+1 ℎ, 𝑥ℎ)ℎ= (𝑢𝑛+1 ℎ, 𝑥ℎ)ℎ(33) for all 𝑥ℎ∈𝑋ℎ. It should be noted that scheme (32)–(33) combines a finite element method together a mass-lumping technique to treat some terms and a semi-implicit time integrator. The resulting scheme is linear and decouples the computation of 𝑢𝑛+1 ℎand 𝑣𝑛+1 ℎ. In order to carry out our numerical analysis we must rewrite the chemotaxis term by using a barycentric quadrature rule as follows. Let 𝑇∈ℎand consider 𝒃𝑇∈𝑇to be the barycenter of 𝑇. Then let 𝑢𝑛+1 ℎ be the interpolation of 𝑢𝑛+1 ℎinto 𝑋0 ℎ, with 𝑋0 ℎbeing the space of all piecewise constant functions over ℎ, defined by 𝑢𝑛+1 ℎ|𝑇=𝑢𝑛+1 ℎ(𝒃𝑇).(34) As a result, one has (∇𝑣𝑛 ℎ, 𝑢𝑛+1 ℎ∇𝑥ℎ) = ∑ 𝑇∈𝐾|𝑇|∇𝑣𝑛 ℎ⋅∇𝑥ℎ𝑢𝑛+1 ℎ(𝒃𝑇) = (∇𝑣𝑛 ℎ,𝑢𝑛+1 ℎ∇𝑥ℎ).(35) Let us define 0(𝑢ℎ, 𝑣ℎ) = 1 2‖𝑣ℎ‖2 ℎ+1 2‖∇𝑣ℎ‖2− (𝑢ℎ, 𝑣ℎ)ℎ+ (log 𝑢ℎ, 𝑢ℎ)ℎ,(36) 1(𝑢ℎ, 𝑣ℎ) = ‖𝑢ℎ‖2 ℎ+‖ 𝛥ℎ𝑣ℎ‖2.(37) 72 J.V. Gutiérrez-Santacreu and J.R. Rodríguez-Galván Computers and Mathematics with Applications 85 (2021) 69–81 and, for each 𝜀, 𝛿 ∈ (0,1), 𝜀,𝛿 0(𝑢0 ℎ, 𝑣0 ℎ) ∶= 1 𝛿𝑒 +‖𝑢0 ℎ‖𝐿1(𝛺) 𝛿 ×(𝐶𝛺 𝜀+𝜀+(1 + 𝛿) |𝛺|(‖𝑣0 ℎ‖𝐿1(𝛺)+‖𝑢0 ℎ‖𝐿1(𝛺))).(38) Associated with the above definitions, consider 0(𝑢ℎ, 𝑣ℎ) = 1 𝛿0(𝑢ℎ, 𝑣ℎ) + 𝛿,𝜀 0(𝑢ℎ, 𝑣ℎ), 1(𝑢ℎ, 𝑣ℎ) = (1 + 1 𝛿)0(𝑢ℎ, 𝑣ℎ) + 𝛿,𝜀 0(𝑢ℎ, 𝑣ℎ)+2|𝛺| 𝑒, and 2(𝑢ℎ, 𝑣ℎ) = 0(𝑢ℎ, 𝑣ℎ) + 0(𝑢ℎ, 𝑣ℎ) + 1(𝑢ℎ, 𝑣ℎ). Finally, define (𝑢ℎ, 𝑣ℎ) = 𝑒2(𝑢ℎ,𝑣ℎ)+𝑇 1 21 2 2(𝑢ℎ,𝑣ℎ)(0(𝑢ℎ, 𝑣ℎ)+𝐶𝑇 3 1(𝑢ℎ, 𝑣ℎ)+𝐶𝑇 ‖𝑢ℎ‖𝐿1(𝛺)), The definition of the above quantities will be apparent later. We are now prepared to state the main result of this paper. Theorem 3.1. Assume that hypotheses (H1)–(H3) are satisfied. Let (𝑢0, 𝑣0) ∈ 𝐻1(𝛺) × 𝐻2 𝑁(𝛺)with 𝑢0>0such that ‖𝑢0‖𝐿1(𝛺)∈ (0,4𝜃𝛺)and 𝑣0≥0, and take 𝑢0 ℎ>0and 𝑣0 ℎ≥0defined by (28) and (30), respectively. Assume that (ℎ, 𝑘)fulfill 𝑘 ℎ2(𝑢0 ℎ, 𝑣0 ℎ)<1 2𝐶(39) and ℎ1− 2 𝑝1 2(𝑢0 ℎ, 𝑣0 ℎ)<𝐶neg 𝐶.(40) Then the sequence {(𝑢𝑚 ℎ, 𝑣𝑚 ℎ)}𝑁 𝑚=0 computed via (32) and (33) satisfies the following properties, for all 𝑚∈ {0,…, 𝑁}: •Lower bounds: 𝑢𝑚 ℎ(𝒙)>0(41) and 𝑣𝑚 ℎ(𝒙)≥0(42) for all 𝒙∈𝛺, •𝐿1(𝛺)-bounds: ‖𝑢𝑚 ℎ‖𝐿1(𝛺)=‖𝑢0 ℎ‖𝐿1(𝛺)(43) and ‖𝑣𝑚 ℎ‖𝐿1(𝛺)≤‖𝑣0 ℎ‖𝐿1(𝛺)+‖𝑢0 ℎ‖𝐿1(𝛺).(44) •A discrete energy law: 0(𝑢𝑚 ℎ, 𝑣𝑚 ℎ) + 𝑘 𝑚 ∑ 𝑟=1 (‖𝛿𝑡𝑣𝑟 ℎ‖2 ℎ+𝑘‖−1 2 ℎ(𝑢𝑟 ℎ)∇𝑢𝑟 ℎ−1 2 ℎ(𝑢𝑟 ℎ)∇𝑣𝑟−1 ℎ‖2) ≤0(𝑢0 ℎ, 𝑣0 ℎ), (45) where ℎis defined in (61). Moreover, if we are given ℎsuch that 𝐶ℎ1− 2 𝑝1(𝑢0 ℎ, 𝑣0 ℎ)≤5 12 ,(46) it follows that 1(𝑢𝑚 ℎ, 𝑣𝑚 ℎ) + 𝑘 2 𝑚 ∑ 𝑟=1 (‖∇𝑢𝑟 ℎ‖2+‖∇ 𝛥ℎ𝑣𝑟 ℎ‖2)≤(𝑢0 ℎ, 𝑣0 ℎ).(47) Remark 3.2. The constant 𝐶is not easy to compute in practice. Hence (39) and (40) should only be seen as theoretical conditions, meaning that 𝑘∕ℎ2and ℎ2− 1 𝑝have to be sufficiently small to reach (41) and (42) on (0, 𝑇 ]. As system (32)–(33) is linear, existence follows from uniqueness. The latter is an immediate outcome of a priori bounds for (𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ). 4. Proof of main result In this section we address the proof of Theorem 3.1. Rather than prove en masse the estimates in Theorem 3.1, because all of them are connected, we have divided the proof into various subsections for the sake of clarity. The final argument will be an induction procedure on 𝑛relied on the semi-explicit time discretization employed in (32). 4.1. Lower bounds and a discrete energy law We first demonstrate lower bounds for (𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ)and, as a consequence of this, a discrete local-in-time energy law is established. Lemma 4.1 (Lower Bounds).Assume that (H1)–(H3) are satisfied. Let 𝑢𝑛 ℎ>0and 𝑣𝑛 ℎ≥0and let ‖ 𝛥ℎ𝑣𝑛 ℎ‖2≤(𝑢0 ℎ, 𝑣0 ℎ).(48) Then if one chooses (ℎ, 𝑘)satisfying (39) and (40), it follows that the solution pair (𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ) ∈ 𝑋2 ℎcomputed via (32) and (33) is lower bounded, i.e, for all 𝒙∈𝛺, 𝑢𝑛+1 ℎ(𝒙)>0(49) and 𝑣𝑛+1 ℎ(𝒙)≥0.(50) Proof. Since 𝑢𝑛+1 ℎand 𝑣𝑛+1 ℎare piecewise linear polynomial functions, it will suffice to prove that (49) and (50) hold at the nodes. To do this, let 𝑇∈ℎbe a fixed triangle with vertices {𝒂1,𝒂2,𝒂3}, and choose two of them, i.e. 𝒂𝑖,𝒂𝑗∈𝑇with 𝑖≠𝑗. Then, from (6) for 𝑝= ∞,(7),(27), and (48), we have on noting (35) that ∫𝑇 𝜑𝒂𝑖∇𝑣𝑛 ℎ⋅∇𝜑𝒂𝑗d𝒙=∫𝑇 𝜑𝒂𝑖(𝒃𝑇)∇𝑣𝑛 ℎ⋅∇𝜑𝒂𝑗d𝒙 ≤|𝑇|‖𝜑𝒂𝑖‖𝐿∞(𝑇)‖∇𝑣𝑛 ℎ‖𝐿∞(𝑇)‖∇𝜑𝒂𝑗‖𝐿∞(𝑇) ≤𝐶ℎ‖∇𝑣𝑛 ℎ‖𝐿∞(𝑇)≤𝐶ℎ1− 2 𝑝‖∇𝑣𝑛 ℎ‖𝐿𝑝(𝑇) ≤𝐶ℎ1− 2 𝑝‖ 𝛥ℎ𝑣𝑛 ℎ‖𝐿2(𝛺)≤𝐶ℎ1− 2 𝑝1 2(𝑢0 ℎ, 𝑣0 ℎ). (51) If we now compare (4) with (51), we find on recalling (40) that ∫𝑇 ∇𝜑𝒂𝑖⋅∇𝜑𝒂𝑗d𝒙−∫𝑇 𝜑𝒂𝑖∇𝑣𝑛 ℎ⋅∇𝜑𝒂𝑗d𝒙 ≤−𝐶neg +𝐶ℎ1− 2 𝑝1 2 0(𝑢0 ℎ, 𝑣0 ℎ)<0 and on summing over 𝑇∈ supp 𝜑𝒂𝑖∩ supp 𝜑𝒂𝑗that (∇𝜑𝒂𝑖,∇𝜑𝒂𝑗)−(𝜑𝒂𝑖∇𝑣𝑛 ℎ,∇𝜑𝒂𝑗)<0.(52) Analogously, we have, from (5), that (∇𝜑𝒂𝑖,∇𝜑𝒂𝑖)−(𝜑𝒂𝑖∇𝑣𝑛 ℎ,∇𝜑𝒂𝑖)>0(53) holds under assumption (40). Let 𝑢min ℎ∈𝑋ℎbe defined as 𝑢min ℎ=∑ 𝒂∈ℎ 𝑢− ℎ(𝒂)𝜑𝒂, where 𝑢− ℎ(𝒂) = min{0, 𝑢𝑛+1 ℎ(𝒂)}. Analogously, one defines 𝑢max ℎ∈𝑋ℎas 𝑢max ℎ=∑ 𝒂∈ℎ 𝑢+ ℎ(𝒂)𝜑𝒂, where 𝑢+ ℎ(𝒂) = max{0, 𝑢𝑛+1 ℎ(𝒂)}. It is easy to check that 𝑢𝑛+1 ℎ=𝑢min ℎ+𝑢max ℎ. 73 J.V. Gutiérrez-Santacreu and J.R. Rodríguez-Galván Computers and Mathematics with Applications 85 (2021) 69–81 Set 𝑥ℎ=𝑢min ℎin (32) using (35) to get (𝛿𝑡𝑢𝑛+1 ℎ, 𝑢min ℎ)ℎ+ (∇𝑢𝑛+1 ℎ,∇𝑢min ℎ)−(𝑢𝑛+1 ℎ∇𝑣𝑛 ℎ,∇𝑢min ℎ)=0.(54) Our goal is to show that 𝑢min ℎ≡0. Indeed, note that (𝑢𝑛+1 ℎ, 𝑢min ℎ)ℎ= (𝑢min ℎ+𝑢max ℎ, 𝑢min ℎ)ℎ=‖𝑢min ℎ‖2 ℎ. and hence (𝛿𝑡𝑢𝑛+1 ℎ, 𝑢min ℎ)ℎ=1 𝑘‖𝑢min ℎ‖2 ℎ−1 𝑘(𝑢𝑛 ℎ, 𝑢min ℎ)≥‖𝑢min ℎ‖2 ℎ,(55) where we have used the fact that 𝑢𝑛 ℎ>0. One can further write (∇𝑢𝑛+1 ℎ,∇𝑢min ℎ)−(𝑢𝑛+1 ℎ∇𝑣𝑛 ℎ,∇𝑢min ℎ) =(∇𝑢max ℎ,∇𝑢min ℎ)−(𝑢max ℎ∇𝑣𝑛 ℎ,∇𝑢min ℎ) + (∇𝑢min ℎ,∇𝑢min ℎ)−(𝑢min ℎ∇𝑣𝑛 ℎ,∇𝑢min ℎ), whereupon we deduce from (52) and (53) that (∇𝑢max ℎ,∇𝑢min ℎ)−(𝑢max ℎ∇𝑣𝑛 ℎ,∇𝑢min ℎ) =∑ 𝒂≠ 𝒂∈ℎ 𝑢max ℎ(𝒂)𝑢min ℎ( 𝒂)[(∇𝜑𝒂,∇𝜑 𝒂)−(𝜑𝒂∇𝑣𝑛 ℎ,∇𝜑 𝒂)] +∑ 𝒂∈ℎ 𝑢max ℎ(𝒂)𝑢min ℎ(𝒂)[(∇𝜑𝒂,∇𝜑𝒂)−(𝜑𝒂∇𝑣𝑛 ℎ,∇𝜑𝒂)]≥0, since 𝑢max ℎ(𝒂)𝑢min ℎ( 𝒂)≤0and 𝑢max ℎ(𝒂)𝑢min ℎ(𝒂)=0. Therefore, (∇𝑢min ℎ,∇𝑢min ℎ)−(𝑢min ℎ∇𝑣𝑛 ℎ,∇𝑢min ℎ)≤(∇𝑢𝑛+1 ℎ,∇𝑢min ℎ)−(𝑢𝑛+1 ℎ∇𝑣𝑛 ℎ,∇𝑢min ℎ)(56) As a result, we infer on applying (55) and (56) into (54) that ‖𝑢min ℎ‖2 ℎ+𝑘‖∇𝑢min ℎ‖2≤𝑘(𝑢min ℎ∇𝑣𝑛 ℎ,∇𝑢min ℎ). But inequalities (8),(16) for 𝑛= 2,(27) for 𝑝= 2,(48), and (39) allow us to estimate ‖𝑢min ℎ‖2 ℎ+𝑘‖∇𝑢min ℎ‖2≤𝑘‖𝑢min ℎ‖𝐿∞(𝛺)‖∇𝑣𝑛 ℎ‖𝐿2(𝛺)‖∇𝑢min ℎ‖2 ≤𝐶𝑘 ℎ2(𝑢0 ℎ, 𝑣0 ℎ)‖𝑢min ℎ‖2 ℎ+𝑘 2‖∇𝑢min ℎ‖2 ≤1 2‖𝑢min ℎ‖2 ℎ+𝑘 2‖∇𝑢min ℎ‖2; thereby ‖𝑢min ℎ‖2 ℎ≤0, which, in turn, implies that 𝑢min ℎ≡0and hence 𝑢𝑛+1 ℎ≥0. It remains to prove that indeed 𝑢𝑛+1 ℎ>0. We proceed by contradiction. Let  𝒂∈ℎ be such 𝑢𝑛+1 ℎ( 𝒂)=0. Substitute 𝑥ℎ=𝜑 𝒂into (32) to arrive at 0< 𝑢𝑛 ℎ( 𝒂)∫𝛺 𝜑 𝒂= (∇𝑢𝑛+1 ℎ,∇𝜑 𝒂) − (∇𝑣𝑛 ℎ, 𝑢𝑛+1 ℎ∇𝜑 𝒂) =∑  𝒂≠𝒂∈ℎ 𝑢𝑛+1 ℎ(𝒂){(∇𝜑𝒂,∇𝜑 𝒂)−(𝜑𝒂∇𝑣𝑛 ℎ,∇𝜑 𝒂)}≤0. In the last line we have utilized (52) and the fact that 𝑢𝑛+1 ℎ≥0. This gives a contradiction. It is now a simple matter to show that (50) holds. It completes the proof. □ We are now concerned with obtaining 𝐿1(𝛺)bounds for (𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ). In particular, we will see that Eq. (32) is mass-preserving. Lemma 4.2 (𝐿1(𝛺)-Bounds).Under the conditions of Lemma 4.1, the discrete solution pair (𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ) ∈ 𝑋2 ℎcomputed via (32) and (33) fulfills ‖𝑢𝑛+1 ℎ‖𝐿1(𝛺)=‖𝑢0 ℎ‖𝐿1(𝛺)(57) and ‖𝑣𝑛+1 ℎ‖𝐿1(𝛺)≤‖𝑣0 ℎ‖𝐿1(𝛺)+‖𝑢0 ℎ‖𝐿1(𝛺).(58) Proof. On choosing 𝑥ℎ= 1 in (32), it follows immediately after a telescoping cancellation that ∫𝛺 𝑢𝑛+1 ℎ(𝒙) d𝒙=∫𝛺 𝑢0 ℎ(𝒙) d𝒙.(59) Consequently, we get that (57) holds from (49) and (29). Now let 𝑥ℎ= 1 in (33) to get ∫𝛺 𝑣𝑛+1 ℎ(𝒙) d𝒙+𝑘∫𝛺 𝑣𝑛+1 ℎ(𝒙) d𝒙=∫𝛺 𝑣𝑛 ℎ(𝒙) d𝒙+𝑘∫𝛺 𝑢𝑛+1 ℎ(𝒙) d𝒙. A simple calculation shows that ∫𝛺 𝑣𝑛+1 ℎ(𝒙) d𝒙=1 (1 + 𝑘)𝑛+1 ∫𝛺 𝑣0 ℎ(𝒙) d𝒙+(∫𝛺 𝑢0 ℎ(𝒙) d𝒙)𝑛+1 ∑ 𝑗=1 𝑘 (1 + 𝑘)𝑗, where we have used (59). Inequality (58) is proved by applying (50).□ Once the positivity of 𝑢𝑛+1 ℎhas been proved, we are in a position to reformulate equation (32) so as to be able to obtain a discrete energy law, which exactly mimics its counterpart at the continuous level. Lemma 4.3. Under the conditions of Lemma 4.1. Eq. (32) can be equivalently written as (𝛿𝑡𝑢𝑛+1 ℎ, 𝑥ℎ)ℎ+ (∇𝑢𝑛+1 ℎ,∇𝑥ℎ) = (∇𝑣𝑛 ℎ,ℎ(𝑢𝑛+1 ℎ)∇𝑥ℎ),(60) where ℎ(𝑢𝑛+1 ℎ) ∈ R2×2 is a piecewise constant, diagonal matrix defined as follows. Let 𝑇∈ℎ. Then there exist two pairs (𝒂𝑇 𝑢𝑖,𝒂𝑇 𝑢𝑖) ∈ 𝑇2,𝑖= 1,2, such that [ℎ(𝑢𝑛+1 ℎ)|𝑇]𝑖𝑖 =⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖) − 𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖) 𝑓′(𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖)) − 𝑓′(𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖)) if 𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖) − 𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖)≠0, 𝑢𝑛+1 ℎ(𝒃𝑇)if 𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖) − 𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖)=0, (61) where 𝑓(𝑠) = 𝑠log 𝑠−𝑠. Proof. We must identify ℎ(𝑢𝑛+1 ℎ)≡𝐼2𝑢𝑛+1 ℎ, where 𝐼2is the 2 ×2 identity matrix. To do so, first observe that 𝑓′(𝑠) = log 𝑠. Then it is easy to see that, for each 𝜀 > 0and 𝑐∈ (0,∞), there exist two points 0< 𝑢 < 𝑢 such that |𝑢−𝑢|< 𝜀, with 𝑢≤𝑐≤𝑢, so that 𝑓′(𝑢) − 𝑓′(𝑢) 𝑢−𝑢=1 𝑐. Let 𝑇∈ℎand choose 𝑐=𝑢𝑛+1 ℎ(𝒃𝑇). We are allowed to choose 𝜀small enough such that 𝑢≥min 𝑡∈(−1,1) 𝑢𝑛+1 ℎ(𝒃𝑇+𝑟𝒃𝑇𝒆𝑖𝑡)≠max 𝑡∈(−1,1) 𝑢𝑛+1 ℎ(𝒃𝑇+𝑟𝒃𝑇𝒆𝑖𝑡)≥𝑢, where 𝑟𝒃𝑇=𝑑𝑖𝑠𝑡(𝒃𝑇, 𝜕𝑇 )and 𝒆𝑖is the i𝑡ℎ vector of the canonical basis of R2. Therefore, there exists a pair (𝒂𝑇 𝑢𝑖,𝒂𝑇 𝑢𝑖)such that 𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖) = 𝑢and 𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖) = 𝑢and hence one defines [ℎ(𝑢𝑛+1 ℎ)|𝑇]𝑖𝑖 = 𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖) − 𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖) 𝑓′(𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖)) − 𝑓′(𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖)) . In the case that min𝑡∈(−1,1) 𝑢𝑛+1 ℎ(𝒃𝑇+𝑟𝒃𝑇𝒆𝑖𝑡) = max𝑡∈(−1,1) 𝑢𝑛+1 ℎ(𝒃𝑇+𝑟𝒊𝑇𝒆𝑖𝑡), one defines [ℎ(𝑢𝑛+1 ℎ)|𝑇]𝑖𝑖 =𝑢𝑛+1 ℎ(𝒃𝑇). This completes the proof. □ It is now shown that system (32)–(33) enjoys a discrete energy law locally in time. Lemma 4.4 (A Discrete Energy Law).Under the conditions of Lemma 4.1, the discrete solution (𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ) ∈ 𝑋2 ℎcomputed via (32) and (33) satisfies 0(𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ) − 0(𝑢𝑛 ℎ, 𝑣𝑛 ℎ) + 𝑘‖𝛿𝑡𝑣𝑛+1 ℎ‖2 ℎ +𝑘‖−1 2 ℎ(𝑢𝑛+1 ℎ)∇𝑢𝑛+1 ℎ−1 2 ℎ(𝑢𝑛+1 ℎ)∇𝑣𝑛 ℎ‖2≤0 (62) where 0(⋅,⋅)is defined in (36). 74 J.V. Gutiérrez-Santacreu and J.R. Rodríguez-Galván Computers and Mathematics with Applications 85 (2021) 69–81 Proof. First of all, recall that 𝑓(𝑠) = 𝑠log 𝑠−𝑠; therefore, we are allowed to compute 𝑓′(𝑢𝑛+1 ℎ)due to (49). Select 𝑥ℎ=ℎ𝑓′(𝑢𝑛+1 ℎ)−𝑣𝑛 ℎin (60) and 𝑥ℎ=𝛿𝑡𝑣𝑛+1 ℎin (33) to get (𝛿𝑡𝑢𝑛+1 ℎ, 𝑓′(𝑢𝑛+1 ℎ) − 𝑣𝑛 ℎ)ℎ+ (∇𝑢𝑛+1 ℎ,∇(ℎ𝑓′(𝑢𝑛+1 ℎ) − 𝑣𝑛 ℎ)) −(∇𝑣𝑛 ℎ,ℎ(𝑢𝑛+1 ℎ)∇(ℎ𝑓′(𝑢𝑛+1 ℎ) − 𝑣𝑛 ℎ)) = 0 (63) and 1 2𝑘(‖𝑣𝑛+1 ℎ‖2 ℎ+‖∇𝑣𝑛+1 ℎ‖2) − 1 2𝑘(‖𝑣𝑛 ℎ‖2 ℎ+‖∇𝑣𝑛 ℎ‖2) +1 2𝑘‖𝑣𝑛+1 ℎ−𝑣𝑛 ℎ‖2 ℎ+1 2𝑘‖∇(𝑣𝑛+1 ℎ−𝑣𝑛 ℎ)‖2+‖𝛿𝑡𝑣𝑛+1 ℎ‖2 ℎ− (𝑢𝑛+1 ℎ, 𝛿𝑡𝑣𝑛+1 ℎ)ℎ= 0. (64) We next pair some terms from (63) and (64) in order to handle them together. It is not hard to see that − (𝛿𝑡𝑢𝑛+1 ℎ, 𝑣𝑛 ℎ)ℎ− (𝑢𝑛+1 ℎ, 𝛿𝑡𝑣𝑛+1 ℎ)ℎ= − 1 𝑘(𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ)ℎ−1 𝑘(𝑢𝑛 ℎ, 𝑣𝑛 ℎ)ℎ.(65) In view of (61), there holds 𝜕𝑥𝑖ℎ(𝑓′(𝑢𝑛+1 ℎ)) = 𝑓′(𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖)) − 𝑓′(𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖)) [𝒂𝑇 𝑢𝑖−𝒂𝑇 𝑢𝑖]𝑖 , since ℎis piecewise linear. As a result, one deduces from (61) that ∇ℎ𝑓′(𝑢𝑛+1 ℎ) = −1 ℎ(𝑢𝑛+1 ℎ)∇𝑢𝑛+1 ℎ. Therefore, (∇𝑢𝑛+1 ℎ,∇(ℎ𝑓′(𝑢𝑛+1 ℎ) − 𝑣𝑛 ℎ)) = (∇𝑢𝑛+1 ℎ,−1 ℎ(𝑢𝑛+1 ℎ)∇𝑢𝑛+1) − (∇𝑢𝑛+1 ℎ,∇𝑣𝑛 ℎ) and (∇𝑣𝑛 ℎ,ℎ(𝑢𝑛+1 ℎ)∇(ℎ𝑓′(𝑢𝑛+1 ℎ) − 𝑣𝑛 ℎ)) = (∇𝑣𝑛 ℎ,∇𝑢𝑛+1 ℎ) − (∇𝑣𝑛 ℎ,ℎ(𝑢𝑛+1 ℎ)∇𝑣𝑛 ℎ), which imply that (∇𝑢𝑛+1 ℎ,∇(ℎ𝑓′(𝑢𝑛+1 ℎ) − 𝑣𝑛 ℎ)) − (∇𝑣𝑛 ℎ,ℎ(𝑢𝑛+1 ℎ)∇(ℎ𝑓′(𝑢𝑛+1 ℎ) − 𝑣𝑛 ℎ)) =‖−1 2 ℎ(𝑢𝑛+1 ℎ)∇𝑢𝑛+1 ℎ−−1 2 ℎ(𝑢𝑛+1 ℎ)∇𝑣𝑛 ℎ‖2.(66) A Taylor polynomial of 𝑓round 𝑢𝑛+1 ℎevaluated at 𝑢𝑛 ℎyields 𝑓(𝑢𝑛 ℎ) = 𝑓(𝑢𝑛+1 ℎ) − 𝑓′(𝑢𝑛+1 ℎ)(𝑢𝑛+1 ℎ−𝑢𝑛 ℎ) + 𝑓′′(𝑢𝑛+𝜃 ℎ) 2(𝑢𝑛+1 ℎ−𝑢𝑛 ℎ)2, where 𝜃∈ (0,1) such that 𝑢𝑛+𝜃 ℎ=𝜃𝑢𝑛+1 ℎ+ (1 − 𝜃)𝑢𝑛 ℎ. Hence, (𝜕𝑡𝑢𝑛+1 ℎ, 𝑓′(𝑢𝑛+1 ℎ))ℎ=1 𝑘(𝑓(𝑢𝑛+1 ℎ),1)ℎ−1 𝑘(𝑓(𝑢𝑛+1 ℎ),1)ℎ +𝑘 2(𝑓′′(𝑢𝑛+𝜃 ℎ),(𝛿𝑡𝑢𝑛+1 ℎ)2)ℎ. In fact, one can write the above expression as (𝜕𝑡𝑢𝑛+1 ℎ, 𝑓′(𝑢𝑛+1 ℎ))ℎ= (𝑢𝑛+1 ℎ,log 𝑢𝑛+1 ℎ)ℎ−(𝑢𝑛 ℎ,log 𝑢𝑛 ℎ)ℎ+𝑘(𝑓′′(𝑢𝑛+𝜃 ℎ),(𝛿𝑡𝑢𝑛+1 ℎ)2)ℎ (67) owing to ∫𝛺 𝑢𝑛+1 ℎ(𝒙) d𝒙=∫𝛺 𝑢𝑛 ℎ(𝒙) d𝒙. On adding (63) and (64), we verify (62) from (65),(66), and (67).□ 4.2. A priori bounds Now that we have accomplished the discrete energy law (62) for system (32)–(33), our goal is to derive a priori energy bounds. It will be no means obvious since 0(𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ)does not provide directly any control over 𝑢𝑛+1 ℎand 𝑣𝑛+1 ℎ. The key ingredient will be the discrete Moser–Trudinger inequality (18). Lemma 4.5 (Control of (𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ)ℎ).Assume that the conditions of Lemma 4.1 are satisfied. Let 𝑢0∈𝐻1(𝛺)such that ‖𝑢0‖𝐿1(𝛺)∈ (0,4𝜃𝛺). Then there exist 𝛿, 𝜀 ∈ (0,1) such that (𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ)ℎ≤1 𝛿0(𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ) + 𝜀,𝛿 0(𝑢0 ℎ, 𝑣0 ℎ),(68) where 0(⋅,⋅)and 𝜀,𝛿 0(⋅,⋅)are given in (36) and (38), respectively. Proof. Let 𝛿, 𝜀 ∈ (0,1) such that (1 + 𝛿)2[8𝜃𝛺𝐶MT𝜀+ 1]‖𝑢0‖𝐿1(𝛺) 8𝜃𝛺 ≤1 2.(69) Using Jensen’s inequality and invoking (57), one finds − log ∫𝛺 ℎ(𝑒(1+𝛿)𝑣𝑛+1 ℎ(𝒙)) ‖𝑢0 ℎ‖𝐿1(𝛺) d𝒙= − log ∑ 𝒂∈ℎ 𝑒(1+𝛿)𝑣𝑛+1 ℎ(𝒂) 𝑢𝑛+1 ℎ(𝒂) 𝑢𝑛+1 ℎ(𝒂) ‖𝑢0 ℎ‖𝐿1(𝛺) ×∫𝛺 𝜑𝒂(𝒙) d𝒙 ≤∑ 𝒂∈ℎ − log (𝑒(1+𝛿)𝑣𝑛+1 ℎ(𝒂) 𝑢𝑛+1 ℎ(𝒂))𝑢𝑛+1 ℎ(𝒂) ‖𝑢0 ℎ‖𝐿1(𝛺) ×∫𝛺 𝜑𝑎(𝒙) d𝒙 = − 1 + 𝛿 ‖𝑢0 ℎ‖𝐿1(𝛺) (𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ)ℎ +1 ‖𝑢0 ℎ‖𝐿1(𝛺) (log 𝑢𝑛+1 ℎ, 𝑢𝑛+1 ℎ)ℎ and hence −(𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ)ℎ+ (log 𝑢𝑛+1 ℎ, 𝑢𝑛+1 ℎ)ℎ≥𝛿(𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ)ℎ +‖𝑢0 ℎ‖𝐿1(𝛺)log ‖𝑢0 ℎ‖𝐿1(𝛺) −‖𝑢0 ℎ‖𝐿1(𝛺)log ×∫𝛺 ℎ(𝑒(1+𝛿)𝑣𝑛+1 ℎ(𝒙)) d𝒙. In virtue of (18), we can bound log ∫𝛺 ℎ(𝑒(1+𝛿)𝑣𝑛+1 ℎ(𝒙)) d𝒙≤log(𝐶𝛺(1 + 𝐶MT(1 + 𝛿)2‖∇𝑣𝑛+1 ℎ‖2)) +(1 + 𝛿)2 8𝜃𝛺‖∇𝑣𝑛+1 ℎ‖2+1 + 𝛿 |𝛺|‖𝑣𝑛+1 ℎ‖𝐿1(𝛺) = log( 𝐶𝛺 𝜀) + log(𝜀(1 + 𝐶MT(1 + 𝛿)2‖∇𝑣𝑛+1 ℎ‖2)) +(1 + 𝛿)2 8𝜃𝛺‖∇𝑣𝑛+1 ℎ‖2+1 + 𝛿 |𝛺|‖𝑣𝑛+1 ℎ‖𝐿1(𝛺) ≤𝐶𝛺 𝜀+𝜀+𝐶MT𝜀(1 + 𝛿)2‖∇𝑢𝑛+1 ℎ‖2 +(1 + 𝛿)2 8𝜃𝛺‖∇𝑣𝑛+1 ℎ‖2+1 + 𝛿 |𝛺|‖𝑣𝑛+1 ℎ‖𝐿1(𝛺). Therefore, −(𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ)ℎ+ (log 𝑢𝑛+1 ℎ, 𝑢𝑛+1 ℎ)ℎ≥𝛿(𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ)ℎ+‖𝑢0 ℎ‖𝐿1(𝛺) log ‖𝑢0 ℎ‖𝐿1(𝛺) −‖𝑢0 ℎ‖𝐿1(𝛺)(𝐶𝛺 𝜀+𝜀) −𝐶MT𝜀(1 + 𝛿)2‖𝑢0 ℎ‖𝐿1(𝛺)‖∇𝑣𝑛+1 ℎ‖2 −(1 + 𝛿)2‖𝑢0 ℎ‖𝐿1(𝛺) 8𝜃𝛺‖∇𝑣𝑛+1 ℎ‖2 −(1 + 𝛿)‖𝑢0 ℎ‖𝐿1(𝛺) |𝛺|‖𝑣𝑛+1 ℎ‖𝐿1(𝛺). 75 J.V. Gutiérrez-Santacreu and J.R. Rodríguez-Galván Computers and Mathematics with Applications 85 (2021) 69–81 On recalling (36) and on noting (69), it follows from 𝑥log 𝑥 > −1 𝑒for 𝑥 > 0that 0(𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ)≥1 2‖𝑣𝑛+1 ℎ‖2 ℎ+1 2‖∇𝑣𝑛+1 ℎ‖2+𝛿(𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ)ℎ +‖𝑢0 ℎ‖𝐿1(𝛺)(log ‖𝑢0 ℎ‖𝐿1(𝛺)−𝐶𝛺 𝜀−𝜀) −(1 + 𝛿)2[8𝜃𝛺𝐶MT𝜀+ 1]‖𝑢0‖𝐿1(𝛺) 8𝜃𝛺‖∇𝑣𝑛+1 ℎ‖2 −(1 + 𝛿)‖𝑢0 ℎ‖𝐿1(𝛺) |𝛺|‖𝑣𝑛+1 ℎ‖𝐿1(𝛺) ≥1 2‖𝑣𝑛+1 ℎ‖2 ℎ+𝛿(𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ)ℎ−1 𝑒 −‖𝑢0 ℎ‖𝐿1(𝛺)(𝐶𝛺 𝜀+𝜀+(1 + 𝛿) |𝛺|‖𝑣𝑛+1 ℎ‖𝐿1(𝛺)). Finally, we have, from (58), that 𝛿(𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ)ℎ≤ 0(𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ) + 1 𝑒 +‖𝑢0 ℎ‖𝐿1(𝛺)(𝐶𝛺 𝜀+𝜀 +(1 + 𝛿) |𝛺|(‖𝑣0 ℎ‖𝐿1(𝛺)+‖𝑢0 ℎ‖𝐿1(𝛺))); thus, proving the result. □ The following is an immediate consequence of Lemma 4.5. Corollary 4.6 (Control of ‖𝑢𝑛+1 ℎlog 𝑢𝑛+1 ℎ‖𝐿1(𝛺)).Under the conditions of Lemma 4.1, there holds ‖𝑢𝑛+1 ℎlog 𝑢𝑛+1 ℎ‖𝐿1(𝛺)≤(1 + 1 𝛿)0(𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ) + 𝜀,𝛿 0(𝑢0 ℎ, 𝑣0 ℎ)+2|𝛺| 𝑒.(70) Proof. Write ‖𝑢𝑛+1 ℎlog 𝑢𝑛+1‖𝐿1(𝛺)= − ∫𝛺 𝑢𝑛+1 ℎ(𝒙) log−𝑢𝑛+1(𝒙) d𝒙 +∫𝛺 𝑢𝑛+1 ℎ(𝒙) log+𝑢𝑛+1(𝒙) d𝒙. (71) Clearly, from −1 𝑒≤𝑥log 𝑥≤0for 𝑥∈ [0,1], one gets −∫𝛺 𝑢𝑛+1 ℎ(𝒙) log−𝑢𝑛+1(𝒙) d𝒙≤|𝛺| 𝑒.(72) By the definition of 0(𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ), one can easily deduce from (68) that ∫𝛺 𝑢𝑛+1 ℎ(𝒙) log+𝑢𝑛+1(𝒙) d𝒙≤ 0(𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ) −∫𝛺 𝑢𝑛+1 ℎ(𝒙) log−𝑢𝑛+1(𝒙) d𝒙 + (𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ)ℎ ≤(1 + 1 𝛿)0(𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ) +𝜀,𝛿 0(𝑢0 ℎ, 𝑣0 ℎ) + |𝛺| 𝑒. (73) Thus, inserting (72) and (73) into (71) yields (70).□ At this point a local-in-time, a priori bound for 𝑢𝑛+1 ℎand 𝑣𝑛+1 ℎon which an induction procedure will be applied is derived. Lemma 4.7 (A Priori Bounds).Suppose that the conditions of Lemma 4.1 are fulfilled. Let (𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ) ∈ 𝑋2 ℎbe the discrete solution computed via (32) and (33). Then there holds 1(𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ) − 1(𝑢𝑛 ℎ, 𝑣𝑛 ℎ)+2𝑘‖ 𝛥ℎ𝑣𝑛+1 ℎ‖2 +𝑘(𝛾,ℎ 1(𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ)‖∇𝑢𝑛+1 ℎ‖2+1 2‖∇ 𝛥ℎ𝑣𝑛+1 ℎ‖2) ≤𝐶 𝑘(‖𝛿𝑡𝑣𝑛+1 ℎ‖ℎ+‖𝛿𝑡𝑣𝑛+1 ℎ‖2 ℎ)‖𝑢𝑛+1 ℎ‖2 ℎ +𝐶 𝑘‖𝑢𝑛+1 ℎlog 𝑢𝑛+1 ℎ‖𝐿1(𝛺)+𝐶𝑘‖𝑢0 ℎ‖𝐿1(𝛺), (74) where 1(⋅,⋅)is defined in (37) and 𝛾,ℎ 1(𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ) ∶= 4 3−𝛾−𝛾3‖𝑢𝑛+1 ℎlog 𝑢𝑛+1 ℎ‖𝐿1(𝛺)−𝐶ℎ1− 2 𝑝(𝑢0 ℎ, 𝑣0 ℎ). Proof. Set 𝑥ℎ=𝑢𝑛+1 ℎin (32) to obtain ‖𝑢𝑛+1 ℎ‖2 ℎ−‖𝑢𝑛 ℎ‖2 ℎ+‖𝑢𝑛+1 ℎ−𝑢𝑛 ℎ‖2 ℎ+ 2𝑘‖∇𝑢𝑛+1 ℎ‖2= 2 𝑘(∇𝑣𝑛 ℎ, 𝑢𝑛+1 ℎ∇𝑢𝑛+1 ℎ).(75) Consider 𝑇∈ℎand let 𝒃𝑇be its barycenter to write 𝑢𝑛+1 ℎ(𝒃𝑇) = 𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖)+∇𝑢𝑛+1 ℎ|𝑇⋅(𝒃𝑇−𝒂𝑇 𝑢𝑖) and 𝑢𝑛+1 ℎ(𝒃𝑇) = 𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖)+∇𝑢𝑛+1 ℎ|𝑇⋅(𝒃𝑇−𝒂𝑇 𝑢𝑖); thereby, 𝑢𝑛+1 ℎ(𝒃𝑇) = 1 2(𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖)+𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖))+ 1 2∇𝑢𝑛+1 ℎ|𝑇⋅((𝒃𝑇−𝒂𝑇 𝑢𝑖)+(𝒃𝑇−𝒂𝑇 𝑢𝑖)).(76) To deal with the right-hand side of (75), we proceed as follows. Let us write (∇𝑣𝑛 ℎ, 𝑢𝑛+1 ℎ∇𝑢𝑛+1 ℎ) = 2 ∑ 𝑖=1 ∑ 𝑇∈ℎ∫𝑇 𝜕𝒙𝑖𝑣𝑛 ℎ𝑢𝑛+1 ℎ(𝒃𝑇)𝜕𝒙𝑖𝑢𝑛+1 ℎd𝒙.(77) Thus, on substituting (76) into (77), we arrive at ∫𝑇 𝜕𝒙𝑖𝑣𝑛 ℎ𝑢𝑛+1 ℎ(𝒃𝑇)𝜕𝒙𝑖𝑢𝑛+1 ℎd𝒙=∫𝑇 𝜕𝒙𝑖𝑣𝑛 ℎ𝑢𝑛+1 ℎ(𝒃𝑇) 𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖) − 𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖) (𝒂𝑇 𝑢𝑖−𝒂𝑇 𝑢𝑖)𝑖 d𝒙 =1 2∫𝑇 𝜕𝒙𝑖𝑣𝑛 ℎ(𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖) + 𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖)) 𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖) − 𝑢𝑛+1 ℎ(𝒂𝑇 𝑢𝑖) (𝒂𝑇 𝑢𝑖−𝒂𝑇 𝑢𝑖)𝑖 d𝒙 +∫𝑇 𝜕𝒙𝑖𝑣𝑛 ℎ∇𝑢𝑛+1 ℎ⋅((𝒃𝑇−𝒂𝑇 𝑢𝑖) + (𝒃𝑇−𝒂𝑇 𝑢𝑖))𝜕𝒙𝑖𝑢𝑛+1 ℎd𝒙 =1 2∫𝑇 𝜕𝒙𝑖𝑣𝑛 ℎ𝜕𝒙𝑖ℎ(𝑢𝑛+1 ℎ)2d𝒙 +∫𝑇 𝜕𝒙𝑖𝑣𝑛 ℎ∇𝑢𝑛+1 ℎ⋅((𝒃𝑇−𝒂𝑇 𝑢𝑖) + (𝒃𝑇−𝒂𝑇 𝑢𝑖))𝜕𝒙𝑖𝑢𝑛+1 ℎd𝒙, which combined with (77) shows that (∇𝑣𝑛 ℎ,𝑢𝑛+1 ℎ∇𝑢𝑛+1 ℎ) = 1 2(∇𝑣𝑛+1 ℎ,∇ℎ(𝑢𝑛+1 ℎ)2) +1 2(∇(𝑣𝑛 ℎ−𝑣𝑛+1 ℎ),∇ℎ(𝑢𝑛+1 ℎ)2) + 2 ∑ 𝑖=1 ∑ 𝑇∈ℎ∫𝑇 𝜕𝒙𝑖𝑣𝑛 ℎ∇𝑢𝑛+1 ℎ⋅(𝒃𝑇−𝒂𝑇 𝑢𝑖)𝜕𝒙𝑖𝑢𝑛+1 ℎd𝒙 + 2 ∑ 𝑖=1 ∑ 𝑇∈ℎ∫𝑇 𝜕𝒙𝑖𝑣𝑛 ℎ∇𝑢𝑛+1 ℎ⋅(𝒃𝑇−𝒂𝑇 𝑢𝑖)𝜕𝒙𝑖𝑢𝑛+1 ℎd𝒙 ∶= 𝐼1+𝐼2+𝐼3+𝐼4. (78) It remains to bound each term of (78). We first proceed with 𝐼1. Choose 𝑥ℎ=1 2ℎ(𝑢𝑛+1 ℎ)2in (33) to write 𝐼1= − 1 2(𝛿𝑡𝑣𝑛+1 ℎ,ℎ(𝑢𝑛+1 ℎ)2)ℎ−1 2(𝑣𝑛+1 ℎ,ℎ(𝑢𝑛+1 ℎ)2)ℎ +1 2(𝑢𝑛+1 ℎ,ℎ(𝑢𝑛+1 ℎ)2)ℎ ∶= 𝐽1+𝐽2+𝐽3. An estimate for 𝐽1is easily computed from (16) for 𝑛= 2,4and the Gagliardo–Nirenberg interpolation ‖𝑥ℎ‖𝐿4(𝛺)≤𝐶‖𝑥ℎ‖1 2‖𝑥ℎ‖1 2 𝐻1(𝛺). It is given by 𝐽1≤𝐶‖𝛿𝑡𝑣𝑛+1 ℎ‖ℎ‖𝑢𝑛+1 ℎ‖2 𝐿4(𝛺) ≤𝐶‖𝛿𝑡𝑣𝑛+1 ℎ‖ℎ(‖∇𝑢𝑛+1 ℎ‖‖𝑢𝑛+1 ℎ‖ℎ+‖𝑢𝑛+1 ℎ‖2 ℎ) ≤𝐶(‖𝛿𝑡𝑣𝑛+1 ℎ‖ℎ+‖𝛿𝑡𝑣𝑛+1 ℎ‖2 ℎ)‖𝑢𝑛+1 ℎ‖2 ℎ+𝛾 4‖∇𝑢𝑛+1 ℎ‖,2 76 J.V. Gutiérrez-Santacreu and J.R. Rodríguez-Galván Computers and Mathematics with Applications 85 (2021) 69–81 where 𝛾 > 0is a constant to be adjusted later on. From (49) and (50), we know that 𝐽2≤0. For 𝐽3, we use the interpolation ‖𝑥ℎ‖𝐿3(𝛺)≤ 𝛾‖∇𝑥ℎ‖2 3‖𝑥ℎlog |𝑥ℎ|‖1 3 𝐿1(𝛺)+𝐶𝛾 ‖𝑥ℎlog |𝑥ℎ|‖𝐿1(𝛺)+𝐶𝛾 ‖𝑥ℎ‖1 3 𝐿1(𝛺)for 𝛾 > 0 (see [7, Lemma 3.5]) and (16) for 𝑛= 3 to obtain 𝐽3≤𝐶‖𝑢𝑛+1 ℎ‖3 𝐿3(𝛺) ≤𝛾3 2‖𝑢𝑛+1 ℎlog 𝑢𝑛+1 ℎ‖𝐿1(𝛺)‖∇𝑢𝑛+1 ℎ‖2 +𝐶‖𝑢𝑛+1 ℎlog 𝑢𝑛+1 ℎ‖3 𝐿1(𝛺)+𝐶‖𝑢𝑛+1 ℎ‖𝐿1(𝛺). Inequality (9) for 𝑝= 2 shows that 𝐼2= − 𝑘 2(∇𝛿𝑡𝑣𝑛+1 ℎ,∇ℎ(𝑢𝑛+1 ℎ)2) ≤𝐶𝑘 ℎ2‖𝛿𝑡𝑣𝑛+1 ℎ‖ℎ‖𝑢𝑛+1 ℎ‖2 𝐿4(𝛺) ≤𝐶𝑘 ℎ2(‖𝛿𝑡𝑣𝑛+1 ℎ‖ℎ+‖𝛿𝑡𝑣𝑛+1 ℎ‖2 ℎ)‖𝑢𝑛+1 ℎ‖2 ℎ+𝛾 4‖∇𝑢𝑛+1 ℎ‖2. We treat 𝐼3and 𝐼4together. Thus, 𝐼3+𝐼4≤𝐶ℎ‖∇𝑣𝑛 ℎ‖𝐿∞(𝛺)‖∇𝑢𝑛+1 ℎ‖2≤𝐶ℎ1− 2 𝑝‖∇𝑣𝑛 ℎ‖𝐿𝑝(𝛺)‖∇𝑢𝑛+1 ℎ‖2 ≤𝐶ℎ1− 2 𝑝‖ 𝛥ℎ𝑣𝑛 ℎ‖‖∇𝑢𝑛+1 ℎ‖2≤𝐶ℎ1− 2 𝑝(𝑢0 ℎ, 𝑣0 ℎ)‖∇𝑢𝑛+1 ℎ‖2. In the above we used (11),(27), and (48). The estimates for the 𝐼𝑖’s applied to (75) lead to ‖𝑢𝑛+1 ℎ‖2 ℎ−‖𝑢𝑛 ℎ‖2 ℎ+‖𝑢𝑛+1 ℎ−𝑢𝑛 ℎ‖2 ℎ+ 2𝑘‖∇𝑢𝑛+1 ℎ‖2 ≤𝐶(1 + 𝑘 ℎ2)𝑘(‖𝛿𝑡𝑣𝑛+1 ℎ‖ℎ+‖𝛿𝑡𝑣𝑛+1 ℎ‖2 ℎ)‖𝑢𝑛+1 ℎ‖2 ℎ +𝑘(𝛾+𝛾3‖𝑢𝑛+1 ℎlog 𝑢𝑛+1 ℎ‖𝐿1(𝛺)+𝐶ℎ1− 2 𝑝(𝑢0 ℎ, 𝑣0 ℎ)) ‖∇𝑢𝑛+1 ℎ‖2 +𝐶𝑘‖𝑢𝑛+1 ℎlog 𝑢𝑛+1 ℎ‖3 𝐿1(𝛺)+𝐶𝑘‖𝑢𝑛+1 ℎ‖𝐿1(𝛺). (79) Choose 𝑥ℎ= −  𝛥2𝑣𝑛+1 ℎin (33) to get ‖ 𝛥ℎ𝑣𝑛+1 ℎ‖2 ℎ−‖ 𝛥ℎ𝑣𝑛 ℎ‖2 ℎ+‖ 𝛥ℎ(𝑣𝑛+1 ℎ−𝑣𝑛 ℎ)‖2 ℎ +2𝑘‖∇ 𝛥ℎ𝑣𝑛+1 ℎ‖2 ℎ+ 2𝑘‖ 𝛥ℎ𝑣𝑛+1 ℎ𝑣𝑛+1 ℎ‖2 ℎ= 2𝑘(∇𝑢𝑛+1 ℎ,∇ 𝛥ℎ𝑣𝑛+1 ℎ) ≤𝑘(4 3‖∇𝑢𝑛+1 ℎ‖2 +3 2‖∇ 𝛥ℎ𝑢𝑛+1 ℎ‖2). (80) The proof follows by use of (79) and (80).□ 4.3. Induction argument The essential step to finishing up the proof of Theorem 3.1 is an induction argument on 𝑛. We need to verify that the overall sequence {𝑢𝑚 ℎ}𝑁 𝑚=0 provided by system (32)–(33) accomplishes the estimates from Theorem 3.1. Observe first that (𝑢0 ℎ, 𝑣0 ℎ)is uniformly bounded with regard to ℎ, because of (29) and (31), and hence we are allowed to choose (ℎ, 𝑘) satisfying (39) and (40). ∙Case (𝑚= 1). We want to prove Theorem 3.1 for 𝑚= 1. Inequality (48) holds trivially, since (𝑢0 ℎ, 𝑣0 ℎ)is bounded independently of (ℎ, 𝑘); thereby, from (49) and (50), we obtain, for 𝑛= 0, that, for all 𝒙∈𝛺, 𝑢1 ℎ(𝒙)>0(81) and 𝑢1 ℎ(𝒙)≥0.(82) Likewise, we have, by (57) and (58) for 𝑛= 0, that ‖𝑢1 ℎ‖𝐿1(𝛺)=‖𝑢0 ℎ‖𝐿1(𝛺) and ‖𝑣1 ℎ‖𝐿1(𝛺)≤‖𝑣0 ℎ‖𝐿1(𝛺)+‖𝑢0 ℎ‖𝐿1(𝛺). In view of (81) and (82), inequality (62) for 𝑛= 0 shows that 0(𝑢1 ℎ, 𝑣1 ℎ) − 0(𝑢0 ℎ, 𝑣0 ℎ) + 𝑘‖𝛿𝑡𝑣1 ℎ‖2 ℎ +𝑘‖−1 2 ℎ(𝑢1 ℎ)∇𝑢1 ℎ−1 2 ℎ(𝑢1 ℎ)∇𝑣0 ℎ‖2≤0 which, in turn, gives 0(𝑢1 ℎ, 𝑣1 ℎ)≤0(𝑢0 ℎ, 𝑣0 ℎ).(83) Applying (83) to (68) and (70) for 𝑛= 0 yields that (𝑢1 ℎ, 𝑣1 ℎ)ℎ≤1 𝛿0(𝑢0 ℎ, 𝑣0 ℎ) + 𝛿,𝜀 0(𝑢0 ℎ, 𝑣0 ℎ) ∶= 0(𝑢0 ℎ, 𝑣0 ℎ), (84) and ‖𝑢1 ℎlog 𝑢1 ℎ‖𝐿1(𝛺)≤(1 + 1 𝛿)0(𝑢0 ℎ, 𝑣0 ℎ) + 𝜀,𝛿 0(𝑢0 ℎ, 𝑣0 ℎ)+2|𝛺| 𝑒 ∶= 1(𝑢0 ℎ, 𝑣0 ℎ). (85) As a result of applying (84) and (85) to (62) for 𝑛= 0, we find 𝑘‖𝛿𝑡𝑣1 ℎ‖2 ℎ≤ 0(𝑢0 ℎ, 𝑣0 ℎ) − 0(𝑢1 ℎ, 𝑣1 ℎ)ℎ ≤ 0(𝑢0 ℎ, 𝑣0 ℎ) + 0(𝑢0 ℎ, 𝑣0 ℎ) + 1(𝑢0 ℎ, 𝑣0 ℎ) ∶= 2(𝑢0 ℎ, 𝑣0 ℎ). Selecting 𝛾to be sufficiently small such that 𝛾+𝛾32(𝑢0 ℎ, 𝑣0 ℎ)≤5 12 (86) and recalling (46), this implies, from (83), that 𝛾,ℎ 1(𝑢1 ℎ, 𝑣1 ℎ)≥4 3−𝛾−𝛾32(𝑢0 ℎ, 𝑣0 ℎ) − 𝐶ℎ1− 1 𝑝1(𝑢0 ℎ, 𝑣0 ℎ)≥1 2; thus, one can find upon using (74) for 𝑛= 0 that 1(𝑢1 ℎ, 𝑣1 ℎ) − 1(𝑢0 ℎ, 𝑣0 ℎ)+2𝑘‖ 𝛥ℎ𝑣1 ℎ‖2+𝑘 2(‖∇𝑢1 ℎ‖2+‖∇ 𝛥ℎ𝑣1 ℎ‖2) ≤𝐶𝑘(‖𝛿𝑡𝑣1 ℎ‖ℎ+‖𝛿𝑡𝑣1 ℎ‖2 ℎ)‖𝑢1 ℎ‖2+𝐶𝑘3 1(𝑢0 ℎ, 𝑣0 ℎ) + 𝐶𝑘‖𝑢0 ℎ‖𝐿1(𝛺). (87) Grönwall’s inequality now provides the bound 1(𝑢1 ℎ, 𝑣1 ℎ) + 𝑘 2(‖∇𝑢1 ℎ‖2+‖∇ 𝛥ℎ𝑣1 ℎ‖2)≤(𝑢0 ℎ, 𝑣0 ℎ). Theorem 3.1 is therefore verified for 𝑚= 1. ∙Case 𝑚=𝑛+ 1. Assume that the bounds in Theorem 3.1 are valid for all 𝑚∈ {1,…, 𝑛}. Consequently, it follows that 1(𝑢𝑛 ℎ, 𝑣𝑛 ℎ)+ 𝑘 2 𝑛 ∑ 𝑟=1 (‖∇𝑢𝑟 ℎ‖2+‖∇ 𝛥ℎ𝑣𝑟 ℎ‖2) ≤1(𝑢0 ℎ, 𝑣0 ℎ) + 𝐶 𝑘 𝑛 ∑ 𝑟=1(‖𝛿𝑡𝑣𝑟 ℎ‖ℎ+‖𝛿𝑡𝑣𝑟 ℎ‖2 ℎ)‖𝑢𝑟 ℎ‖2(88) +𝐶 𝑘 𝑛 ∑ 𝑟=1(3 1(𝑢0 ℎ, 𝑣0 ℎ) + ‖𝑢0 ℎ‖𝐿1(𝛺)) holds on the basis of 𝛾,ℎ 1(𝑢𝑚 ℎ, 𝑣𝑚 ℎ)≥1 2for all 𝑚∈ {0,…, 𝑛}. Then we want to prove Theorem 3.1 for 𝑚=𝑛+ 1. Indeed, by the induction hypothesis (47) for 𝑚=𝑛, it is clear that (48) holds; therefore, one has (49) and (50). That inequalities (43) and (44) are satisfied for 𝑚=𝑛+ 1 is simply by noting (57) and (58). Combining (62) and the induction hypothesis (45) for 𝑚=𝑛, we deduce (45) for 𝑚=𝑛+ 1, which implies max 𝑚∈{0,…,𝑛+1} 0(𝑢𝑚 ℎ, 𝑣𝑚 ℎ)≤0(𝑢0 ℎ, 𝑢0 ℎ).(89) As a result of this, we have, by (68) and (70), that (𝑢𝑛+1 ℎ, 𝑣𝑛+1 ℎ)ℎ≤0(𝑢0 ℎ, 𝑣0 ℎ)(90) and ‖𝑢𝑛+1 ℎlog 𝑢𝑛+1 ℎ‖𝐿1(𝛺)≤1(𝑢0 ℎ, 𝑣0 ℎ).(91) 77