1. Introduction
1.1. Aims
In 1970/71 Keller and Segel [15, 16] 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 wellknown 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 is 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 the equations against non linear functions. Crossdiffusion mechanisms governing the chemotactic phenomena are the responsible for the Keller–Segel equations to be so difficult 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 nonphysical blow up of solutions and hence produce solution being closer to chemotaxis phenomena. For these other chemotaxis models, it is recommended the excellent surveys of Hillen and Painter [11], Horstamann [12, 13], and, more recently, Bellomo, Bellouquid, Tao, and Winkler [1]. 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.
1.2. The Keller–Segel equations
Let be a bounded domain, with
being its outwarddirected unit normal vector to
, and let be a time interval. Take and . Then the boundaryvalue problem for the Keller–Segel equations reads a follows. Find and satisfying(1) 
subject to the initial conditions
(2) 
and the boundary conditions
(3) 
Here denotes is the average density of organisms (or cells), which is a conserved variable, and is the average density of chemioattranct, which is a nonconserved variable.
System (1) was motivated by Keller and Segel [15] 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 ameobae and chemoattractant.
The diffusion phenomena performed by cells and chemoattractant are modelled by the terms and , respectively, while 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 [17] proved existence, uniqueness and regularity of solutions under the condition by means of one particular instance of Moser–Tridinguer’s inequality; in the particular case that be a ball, such a condition becomes , and Herrero and Velázquez [10] dealt with the first blow up framework by constructing some radially symmetric solutions which blow up within finite time. The next progress in this sense with being nonradial and simply connected was the work of Horstmann and Wang [14] who found some unbounded solutions provided that and . So far there is no supporting evidence as to whether such solutions may evolve to produce a blowup phenomenon within finite time or whether, on the contrary, may increase to infinity with time. The main tool in proving blowup solutions is the energy law which stems from system (1). Observe that an inadequate approximation of lower bounds can trigger oscillations of the variables, which can lead to spurious, blowup solutions.
Concerning the numerical analysis for system (1), very little is said about numerical algorithms which keep lower bounds, are mass conservative and have a discrete energy law. Proper numerical treatment of these properties is made difficult by the fact that the nonlinearity occurs in the highest order derivative. Numerical algorithms are mainly designed so as to keep lower bounds. We refer the reader to [18, 5, 22, 3, 21] which mainly deal with lower bounds and mass conservation. As far as we are concerned, there is no numerical method coping with a discrete energy law.
1.3. Notation
We collect here as a reference some standard notation used throughout the paper. For , we denote by the usual Lebesgue space, i.e.,
or
This space is a Banach space endowed with the norm if or if . In particular, is a Hilbert space. We shall use for its inner product and for its norm.
Let be a multiindex with , and let be the differential operator such that
For and , we consider to be the Sobolev space of all functions whose derivatives are in , i.e.,
associated to the norm
For , we denote . Moreover, we make of use the space
for which is known that .
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 4, 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 5 dedicated 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.
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 consider the finite element approximation of (1) under the following assumptions on the domain, the mesh and the finite element space. Moreover, these assumptions will be mentioned on stating our main result.

Let be a convex, bounded domain of with a polygonal boundary, and let be the minimum interior angle at the vertices of .

Let be a family of acute, shaperegular, quasiuniform triangulations of made up of triangles, so that , where , with being the diameter of . More precisely, we assume that

there exists , independent of , such that
where is the largest ball contained in , and

there exists such that every angle between two edges of a triangle is bounded by .
Further, let be the coordinates of the nodes of .


Associated with is the finite element space
where is the set of linear polynomials on . Let be the standard basis functions for .
2.2. Auxiliary results
In the subsequent sections, we will make use of some technical results concerning the abovementioned hypotheses.
A key tool for proving lower bounds for the finite element approximation of (1) is the acuteness of .
Proposition 2.1.
Let be a polygonal. Consider to be constructed over being acute. Then, for each with vertices , there exists a constant , depending on , but otherwise independent of and , such that
(4) 
for all with , and
(5) 
for all .
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
[2, Lem. 4.5.3] or [6, Lem. 1.138] for a proof.Proposition 2.2.
Let be polygonal. Consider to be constructed over being quasiuniform. Then, for each and , there exists a constant , independent of and , such that, for all ,
(6) 
and
(7) 
Concerning global inverse inequalities, we need the following.
Proposition 2.3.
Let be polygonal. Consider to be constructed over being quasiuniform. Then for each , there exists a constant , independent of , such that, for all ,
(8) 
(9) 
(10) 
and
(11) 
We introduce , the standard nodal interpolation operator, such that for all . Associated with , a discrete inner product on is defined by
We also introduce
Local and global error bounds for are as follows (c.f. [2, Thm. 4.4.4] or [6, Thm. 1.103] for a proof).
Proposition 2.4.
Let be polygonal. Consider to be constructed over being quasiuniform. Then, for each , there exists , independent of and , such that
(12) 
Proposition 2.5.
Let be polygonal. Consider to be constructed over being quasiuniform. Then it follows that there exists , independent of , such that
(13) 
Corollary 2.6.
Let be polygonal. Consider to be constructed over being quasiuniform. Let . Then it follows that there exist two positive constants and , independent of , such that
(14) 
(15) 
and
(16) 
where depends on .
Proof.
Let and compute
Then, from (12) and the above identity, we have
Summing over yields (14). The proof of (15) follows very closely the arguments of (14) for . The first part of assertion (16) is a simple application of Jensen’s inequality, while the second part follows from (14) on using Hölder’s inequality, (9) for and, later on, Minkowski’s inequality. ∎
The proof of the following proposition can be found in [4]. It is a generalization of a Moser–Trudingertype inequality.
Proposition 2.7 (MoserTrudinger).
Let be polygonal with being the minimum interior angle at the vertices of . Then there exists a constant depending on such that for all with and , it follows that
(17) 
where .
Corollary 2.8.
Let be polygonal with being the minimum interior angle at the vertices of . Consider to be constructed over being quasiuniform. Let with . Then it follows that there exists a constant , independent of , such that
(18) 
Proof.
From (14), we have
(19) 
Let with . Young’s inequality gives
(20) 
Thus, combining (19) and (20) yields, on noting (7) for and (17), that
∎
An (average) interpolation operator into will be required in order to properly initialize our numerical method. We refer to [19, 7].
Proposition 2.9.
Let be polygonal. Consider to be constructed over being quasiuniform. Then there exists an (average) interpolation operator from to such that
(21) 
and
(22) 
Moreover, let be defined from to as
(23) 
and let be such that
(24) 
From elliptic regularity theory, the wellposedness of (24) is ensured by the convexity assumption stated in and
See [8] for a proof.
Proposition 2.10.
Let be a convex polygon. Consider to be constructed over being quasiuniform. Then there exists a constant , independent of , such that
(25) 
Corollary 2.11.
Let be a convex polygon. Consider to be constructed over being quasiuniform. Then, for each , there exists a constant , independent of , such that
(26) 
3. Presentation of main result
We now define our numerical approximation of system (1). Assume that with and a. e. in .
We begin by approximating the initial data by as follows. Define
(27) 
which satisfies
(28) 
and
(29) 
which satisfies
(30) 
Given , we let be a uniform partitioning of [0,T] with time step . Given , find such that
(31) 
and
(32) 
It should be noted that scheme (31)(32) combines a finite element method together a masslumping technique to treat some terms and a semiimplicit time integrator. The resulting scheme is linear and decouples the computation of .
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 barycentre of . Then let be the interpolation of into , with being the space of all piecewise constant functions over , defined by
(33) 
As a result, one has
(34) 
From now on we will use to denote a generic constant independent of the discretization parameters .
Let us define
(35) 
(36) 
and, for each ,
(37) 
Associated with the above definitions, consider
and
Finally, define
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 – are satisfied. Let with and , and take and defined by (27) and (29), respectively. Assume that fulfill
(38) 
and
(39) 
Then the sequence computed via (31) and (32) satisfies the following properties, for all :

Lower bounds:
(40) and
(41) for all ,

bounds:
(42) and
(43)
Moreover, if we are given such that
(45) 
it follows that
(46) 
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 semiexplicit time discretization employed in (31).
4.1. Lower bounds and a discrete energy law
We first demonstrate lower bounds for and, as a consequence of this, a discrete localintime energy law is established.
Comments
There are no comments yet.