Numerical solution for an aggregation equation with degenerate diffusion
Roberto Carlos Cabrales Juan Vicente Gutiérrez-Santacreu José Rafael Rodríguez-Galván


Numerical solution for a degenerate aggregation equation]Numerical solution for an aggregation equation with degenerate diffusion

R. C. Cabrales]Roberto Carlos Cabrales^
 Instituto de Investigación Multidisciplinaria 
en Ciencia y Tecnología, Universidad de la Serena, 
Benavente 980, La Serena, Chile. E-mail: rcabrales@userena.cl
RCC was partially supported by project PR18511 from Dirección de 
Investigación of Universidad de La Serena
J. V. Gutiérrez-Santacreu]Juan Vicente Gutiérrez-Santacreu^†
† Dpto. de Matemtica Aplicada I, E. T. S. I. Informtica, Universidad de Sevilla. Avda. Reina Mercedes, s/n. E-41012 Sevilla, Spain.  E-mail: juanvi@us.es. 
JVGS and JRRG were partially supported by the Spanish Grant No. PGC2018-098308-B-I00 from Ministerio de Ciencias, Innovacin y Universidades with the participation of FEDER

J. R. Rodrguez-Galvn]Jos Rafael Rodrguez-Galvn^
 Departamento de Matemticas. Facultad de Ciencias. Campus Universitario de Puerto Real, Universidad de Cdiz. 11510 Puerto Real, Cdiz E-mail: rafael.rodriguez@uca.es






 A numerical method for approximating weak solutions of an aggregation equation with degenerate diffusion is introduced. The numerical method consists of a stabilized finite element method together with a mass lumping technique and an extra stabilizing term plus a semi–implicit Euler time integration. Then we carry out a rigorous passage to the limit as the spatial and temporal discretization parameters tend to zero, and show that the sequence of finite element approximations converges toward the unique weak solution of the model at hands. In doing so, nonnegativity is attained due to the stabilizing term and the acuteness on partitions of the computational domain, and hence a priori energy estimates of finite element approximations are established. As we deal with a nonlinear problem, some form of strong convergence is required. The key compactness result is obtained via an adaptation of a Riesz–Frchet–Kolmogorov criterion by perturbation. A numerical example is also presented.



[
    [
    December 30, 2023
=====================



2010 Mathematics Subject Classification. 65M60, 35K55, 45K05, 35K20.

Keywords. Finite-element approximation; Aggregation equation; Nonlinear diffusion.





§ INTRODUCTION



 §.§ The model

 Let Ω⊂^d, d=2 or 3, be a bounded domain and T>0 be a fixed time. We consider an aggregation equation with degenerate diffusion term which reads as follows. Find ρ: Ω̅×[0,T]→ [0,∞) such that

    ∂_t ρ-Δ A(ρ)+∇·(ρ∇ K*ρ)=0     Q:=Ω× (0,T],

subject to the boundary condition

    (-∇ A(ρ) + ρ∇ K * ρ) · = 0    Σ:=∂Ω×(0,T]

and the initial condition

    ρ(0)=ρ^0    Ω,

where * stands for the convolution operator and  is the outward-pointing unit vector to ∂Ω.

Equation (<ref>) arises in many models in biology, where ρ represents the population density,  K * ρ stands for the density of the chemo-attractant, and A(ρ) models the local repulsion. Patlak–Keller–Segel models <cit.> governing the movement of species by chemotaxis are a particular instance, which correspond to considering A(ρ)=ρ^m and  K()=-1/2πlog || for d=2 or K()=1/3 ω_d d(2-d)1/||^d-2 for d≥ 3, with ω_d being the volume of the unit ball in ^d.  Pure aggregation equations modeling biological swarming <cit.> result from ruling  out the diffusion term -Δ A(ρ) and from selecting K(·) to be the Newtonian potential, repulsive-attractive Morse potential, or power law potential.



While there is a rich body of literature on the mathematical analysis of equation (<ref>) supported by numerical simulations, very few results on numerical analysis are available for the situation considered here. Carrillo, Chertock, and Huang <cit.> introduced a positivity-preserving entropy-decreasing finite volume scheme for (<ref>) which takes into account a confinement potential term as well.



The existence and uniqueness of a weak solution to equation (<ref>) was established by Bertozzi and Slepčev <cit.> for A(ρ) being degenerate and K(·) satisfying  some regularity assumptions.  It is this degeneracy of A that is the major source of difficulties in studying equation (<ref>). The existence proof consists of three steps: (a) introducing a regularized problem via the diffusion term A(ρ), (b) establishing a maximum principle and a priori energy bounds independent of the regularizing parameter, and (c) proving compactness for the regularized problem. In particular, the compactness of the regularized solutions is obtained by using some results borrowed from <cit.> based on the Riesz–Frchet–Kolmogorov criterion on Lebesgue spaces.

Our aim in this work is to construct a sequence of fully discrete approximations and analyze its convergence toward the unique solution to (<ref>)-(<ref>). Our algorithm uses a stabilized finite element method combined with a mass lumping technique plus a semi–implicit Euler time integration. This resulting scheme is conditionally solvable and mass conserving, and preserves nonnegativity under acute partitions of the computational domain. A priori energy bounds are obtained in a different way from those in <cit.> since a discrete maximum principle does not hold. The lack of such a discrete maximum principle is overcame with the use of a nodal truncating operator <cit.>. A version of the Riesz–Frchet–Kolmogorov compactness criterion on Lebesgue spaces by perturbation <cit.> allows the passage to the limit in the nonlinear terms as the spatial and temporal discretization parameters tend to zero in order to reach the unique weak solution of (<ref>)-(<ref>).






















 §.§ Notation

For p∈[1,∞], we denote by L^p(Ω) the usual Lebesgue space, i.e.,

    L^p(Ω) = {v : Ω→  :  v , ∫_Ω |v()|^p d<∞}.

or

    L^∞(Ω) = {v : Ω→  :  v ,  esssup_∈Ω |v()|<∞}.



This space is a Banach space endowed with the norm
v_L^p(Ω)=(∫_Ω|v()|^p  d)^1/p if p∈[1, ∞) or v_L^∞(Ω)= esssup_∈Ω|v()| if p=∞. In particular,  L^2(Ω) is a Hilbert space.  We shall use
(u,v)=∫_Ωu()v() d for its inner product and · for its norm.

Let α = (α_1, α_2, . . . , α_d)∈N^d be a
multi-index with |α|=α_1+α_2+...+α_d, and let
∂^α be the differential operator such that

    ∂^α=
    (∂/∂x_1)^α_1...(∂/∂x_d)^α_d.



For m≥ 0 and p∈[1, ∞), we define W^m,p(Ω) to be the Sobolev space of all functions whose m derivatives are in L^p(Ω), i.e.,

    W^m,p(Ω) = {v ∈ L^p(Ω) : ∂^k v ∈ L^2(Ω) ∀  |k|≤ m}
 associated to the norm

    f_W^m,p(Ω)   =(∑_|α|≤ m∂^α f^p_L^p(Ω))^1/p        for 1 ≤ p < ∞, 
    f_W^m,p(Ω)   =max_|α|≤ m∂^α f_L^∞(Ω),         for  p = ∞.

For p=2, we denote W^m,2(Ω)=H^m(Ω) and its dual as (H^m(Ω))'. The dual pairing between H^1(Ω) and (H^1(Ω))' is denoted by <·, ·>.

Let X be a Banach space. Thus, L^p(0,T;X) denotes the space of Bochner-measurable, X-valued functions on (0, T) such that ∫_0^Tf(s)^p_X d s<∞ for p∈ [1,∞) or  esssup_s∈(0,T)f(s)_X<∞ for
p=∞.



 §.§ Outline of the paper 
 The layout of the paper is as follows. In section 2 we introduce the hypotheses for constructing the finite element approximation of (<ref>) as  well as some auxiliary results. In section 3 we present our finite element method which includes a stabilizing term and combines a semi-implicit time integration. Afterwards we state our main theorem which is proved in the subsequent sections. The well-posedness of our algorithm is carried out in section 4. Non-negativity under the acuteness of the mesh and a priori energy estimates are obtained in section 5.  Section 6 deals with the compactness of the finite element approximations. The passage to the limit toward the unique weak solution of (<ref>) is reported in section 7. To finish off, we present a numerical example in section 8.


§ THE DISCRETE SETTING

This section is mainly devoted to the numerical tools for approximating the solution to problem (<ref>)-(<ref>).


 §.§ Hypotheses
 Herein we set out the hypotheses that will be required for the domain, the mesh, and the finite element space.


  (H1) Let Ω be a convex, bounded domain of ^d with polygonal (d=2) or polyhedral (d=3) Lipschitz-continuous boundary.

  (H2) Let {ℰ_h}_h>0  be a family of simplicial partitions of  Ω that is acute, shape-regular, and quasi-uniform, so that Ω=∪_E∈ℰ_hE, where h=max_E∈ℰ_h h_E, with h_E being the diameter of E. More precisely, we assume that


  (a)  there exists α>0, independent of h, such that

    min{ diam   B_E  :   E∈ℰ_h }≥α h,

where B_E is the largest ball contained in E, and

  (b) there exists β > 0 such that  every angle between two edges (or faces)  of a triangle (or a tetrahedron) is bounded by π/2 - β.

Further, let 𝒩_h = {_i}_i∈ I denote the set of all the nodes of ℰ_h.

  (H3) A conforming finite element space associated with
ℰ_h is assumed for approximating H^1(Ω). Let  𝒫_1(E) be the set of linear polynomials on  E; the space of continuous, piecewise polynomial
functions on ℰ_h  is then denoted as

    D_h = {ρ̅_h ∈C^0(Ω)  : ρ̅_h|_E ∈𝒫_1(E),  ∀ E ∈ℰ_h },

whose shape functions are {φ_}_∈𝒩_h.







 §.§ Technical preliminaries
 Under hypotheses (H1)–(H3) we collect some properties that will be used in the subsequent analysis.

To start with, we state a consequence of the acuteness of the mesh needed for proving non-negativity of the finite element approximation.
 Let E∈ℰ_h with vertices {a_0,⋯a_d}. Then there exists a constant C_ neg>0, depending on β, but otherwise independent of h and E, such that

    ∫_E ∇φ__i·∇φ__j≤ - C_ neg h^d-2
 for all _i,_j∈ E with i≠j, and

    ∫_E∇φ__i·∇φ__i≥ C_ neg h^d-2

for all _i∈ E.

 For every d-simplex E∈ℰ_h and for every vertex a_i∈ E, we denote by F_a_i the opposite face  to a_i and by n_a_i the exterior (to the d-simplex E) unit normal vector to the face F_a_i. Write

    ∇φ__i|_E=-1/h_F__i__i,

where h_F__i is the distance of _i to the hyperplane which contains F__i. Then we have

    ∇φ__i|_E·∇φ__j|_E=1/h_F__i1/h_F__j__i·__j.

Note that __i·__j =cos(__i __j)=cos(F__i F__j-π)=-cos(F__i F__j).
Integrating over E gives

    [                                             ∫_E∇φ__i·∇φ__j                                                          =                        -  |E|1/h_F__i1/h_F__jcos(F__iF__j);                                                                                                                     ≤                         -  |B_E|1/h_F__i1/h_F__jcos(π/2-β);                                                                                                                     = - π^d/2/2^dΓ(d/2+1)( diam B_E)^d1/h_F__i1/h_F__jcos(π/2-β);                                                                                                                     ≤       - α^d π^d/2/2^dΓ(d/2+1)h^d1/h_F__i1/h_F__jcos(π/2-β);                                                                                                                     ≤                   - α^d π^d/2/2^dΓ(d/2+1)cos(π/2-β) h^d-2, ]

where we have used that the fact that |B_E|=π^d/2/2^dΓ(d/2+1)( diam B_E)^d with Γ(·) being  Euler's gamma function.

The same argument as in the proof of (<ref>) yields (<ref>).

Some inverse inequalities are provided in the following proposition.
 Let E∈ℰ_h. There exists a constant C_ inv>0, independent of h and E, such that, for all ρ̅_h∈𝒫_1(E),

    ρ̅_h_H^1(E)≤C_ inv/hρ̅_h_L^2(E)

and

    ρ̅_h_L^2(E)≤C_ inv/hρ̅_h_(H^1(E))'.


 The proof of (<ref>) can be found in <cit.> or <cit.>.

To obtain (<ref>), we use a duality argument. Let π_h be the L^2(E) orthogonal interpolation operator from L^2(E) into 𝒫_1(E). Then, from (<ref>), we find

    [                                                                                  ρ̅_h_L^2(E)                                                                                            =                    sup_0≠ρ∈ L^2(E)(ρ̅_h, ρ)/ρ_L^2(E)≤sup_0≠ρ∈ L^2(E)(ρ̅_h, π_hρ)/π_hρ_L^2(E);                                                                                                                                                                                         ≤ C_ inv/hsup_0≠ρ∈ L^2(Ω)(ρ̅_h, π_hρ)/π_hρ_H^1(E)≤C_ inv/hsup_0≠ρ_h∈ D_h(ρ̅_h, ρ_h)/ρ_h_H^1(E);                                                                                                                                                                                         ≤                            C_ inv/hsup_0≠ρ∈ H^1(Ω)(ρ̅_h, ρ)/ρ_H^1(E)=C_ inv/hρ̅_h_(H^1(E))'. ]


 There holds

    ρ̅_h_H^1(Ω)≤C_ inv/hρ̅_h_L^2(Ω)

and

    ρ̅_h_L^2(Ω)≤C_ inv/hρ̅_h_(H^1(Ω))'.



Let ℐ_h be the nodal interpolation operator from C^0(Ω̅) to D_h and  consider

    (ρ_h,ρ̅_h)_h=∫_Ωℐ_h( ρ_hρ̅_h)=∑_∈𝒩_hρ_h()·ρ̅_h() ∫_Ωφ_

for all ρ_h,ρ̅_h∈  D_h, with the induced norm ρ_h_h = √(( ρ_h, ρ_h)_h). It is well-known that there exists a constant C_ eq>1, independent of h, such that

    ρ_h_h≤ρ_h_L^2(Ω)≤ C_ eqρ_h_h.

From the definition of ℐ_h, one can straightforwardly check the following.
 Let E∈ℰ_h. It follows that

    ℐ_hφ_L^∞(E)≤φ_L^∞(E)    φ∈ L^∞(E)

and

    ∇ℐ_hφ_L^∞(E)≤∇φ_L^∞(E)    φ∈ W^1,∞(E).


 There holds

    ℐ_hφ_L^∞(Ω)≤φ_L^∞(Ω)    φ∈ L^∞(Ω)

and

    ∇ℐ_hφ_L^∞(Ω)≤∇φ_L^∞(Ω)    φ∈ W^1,∞(Ω).





 Let E∈ℰ_h. There exists a constant C_ app>0, independent of h and E, such that

    φ-ℐ_hφ_L^∞(E)+ h ∇(φ-ℐ_hφ)_L^∞(E)≤ C_  app h^2 ∇^2 φ_L^∞(E)    φ∈ W^2,∞(E).

and

    ∇(φ-ℐ_hφ)_L^2(E)≤ C_  app h ∇^2 φ_L^2(E)    φ∈ H^2(E).


  The proof of (<ref>) and (<ref>) can be found in <cit.> or <cit.>.

 There holds

    φ-ℐ_hφ_L^∞(Ω)+ h ∇(φ-ℐ_hφ)_L^∞(Ω)≤ C_  app h^2 ∇^2 φ_L^∞(Ω)    φ∈ W^2,∞(Ω)

and

    ∇(φ-ℐ_hφ)_L^2(Ω)≤ C_  app h ∇^2 φ_L^2(Ω)    φ∈ H^2(Ω).


 There exists a constant C_ com>0, independent of h and E, such that

    ρ_h ρ̅_h-ℐ_h(ρ_h ρ_h)_L^1(Ω)≤ C_ com h^1/2 ρ_h_(H^1(Ω))' ∇ρ̅_h_L^∞(Ω)

and

    ρ_h ρ̅_h-ℐ_h(ρ_h ρ_h)_L^1(Ω)≤ C_ com h ρ_h_L^2(Ω) ∇ρ̅_h_L^2(Ω).



 On each element E∈ℰ_h, combine (<ref>), (<ref>), and (<ref>) to obtain

    [                    ρ_h ρ̅_h-ℐ_h(ρ_h ρ_h)_L^1(E)                                               ≤             C_ app h^2 |E|  ∇^2(ρ_hρ̅_h)_L^∞(E);                                                                                               ≤          C_ app h^2 |E| ∇ρ_h_L^∞(E)∇ρ̅_h_L^∞(E);                                                                                               ≤              C_ app h^2 ∇ρ_h_L^1(E)∇ρ̅_h_L^∞(E);                                                                                               ≤             C_ app h^5/2∇ρ_h_L^2(E)∇ρ̅_h_L^∞(E);                                                                                               ≤ C_ app C_ inv^2 h^1/2ρ_h_(H^1(E))'∇ρ̅_h_L^∞(E). ]

Estimate (<ref>) follows by summing up this last estimate over all the elements E ∈ℰ_h.

One can prove estimate (<ref>) in a similar fashion.



 Let f∈ C^0,1(R)  be monotonically increasing with Lipschitz constant C_ Lip. Then  it follows that, for all ρ_h∈ D_h,

    |∇ℐ_h f(ρ_h)|^2≤ C_ Lip∇ρ_h·∇ℐ_h f(ρ_h).


 On each element E∈ℰ_h, consider Ẽ to be an oriented, right element with vertices {_0^Ẽ, ⋯, _d^Ẽ}, where  _0^Ẽ is the vertex supporting the right angle, such that Ẽ⊂ E. Observe that

    ∂_x_iρ_h=ρ_h(_i^Ẽ)-ρ_h(_0^Ẽ)/(_i^Ẽ-_0^Ẽ)_i,

where ()_i is the ith component of . Since

    (f(x)-f(y))^2≤ C_ Lip (f(x)-f(y))(x-y)     x,y∈,

we have

    [                                                            C_ Lip∇ρ_h|_E·∇ℐ_hf(ρ_h)|_E                                                                                      = C_ Lip∑_i=1^d ρ_h(_i^Ẽ)-ρ_h(_0^Ẽ)/(_i^Ẽ-_0^Ẽ)_if(ρ_h(_i^Ẽ))-f(ρ_h(_0^Ẽ))/(_i^Ẽ-_0^Ẽ)_i;                                                                                                                                                                             ≥             ∑_i=1^d (f(ρ_h(_i^Ẽ))-f(ρ_h(_0^Ẽ)))^2/((_i^Ẽ-_0^Ẽ)_i)^2=|∇ℐ_hf(ρ_h)|_E|^2, ]

where e_i is the ith vector of the  Cartesian basis. We deduce (<ref>) upon summing over all the element E∈ℰ_h.

For each element E∈ℰ_h with vertices {a_0,⋯a_d},  we associate once and for all a vertex _E∈ E. Thus we define 𝒫_h(ρ_h)(x)= ρ_h(a_E) for all x∈ E.
 There exists a constant C_ int>0, independent of h, such that

    ρ_h-𝒫_hρ_h_L^2(Ω)≤ C_ int h ∇ρ_h_L^2(Ω)    ρ_h∈ D_h.


 Let ∈ E and write

    ρ_h(x)-𝒫_h(ρ_h)()=ρ_h()-ρ_h(_E)= ∇ρ_h|_E·(-_E).

Squaring and integrating over E gives

    ρ_h-𝒫_h(ρ_h)_L^2(E)≤ C  h ∇ρ_h_L^2(E)

 and  hence summing over E∈ℰ_h yields the desired result.

Moreover, let Δ̃_h be defined from D_h to D_h as

    -(Δ̃_h ϕ_h, ρ̅_h)_h=(∇ϕ_h, ∇ρ̅_h)  ρ̅_h∈ D_h,

and let ϕ(h)∈ H^2(Ω) be such that

    {[    -Δϕ(h)         = -Δ̃_h ϕ_h                  Ω,;    ∂_ϕ(h)         =         0                 ∂Ω. ].

The H^2(Ω)-regularity of ϕ(h) is ensured by the convexity assumption stated in (H1). See <cit.> for a proof.
 There exists a constant C_ Lap>0, independent of h, such that, for all ϕ_h∈ D_h,

    ∇(ϕ(h)-ϕ_h)_L^2(Ω)≤ C_ Lap h Δ̃_hϕ_h_L^2(Ω).


 Testing (<ref>) with ρ̅_h∈ D_h yields

    (∇ϕ(h), ∇ρ̅_h)=-(Δ̃_hϕ_h, ρ̅_h).

Combining the above equation and (<ref>), we write

    (∇ (ϕ(h)-ϕ_h), ∇ρ̅_h)=(Δ̃_hϕ_h, ρ̅_h)_h-(Δ̃_hϕ_h, ρ̅_h)

and hence

    (∇ (ℐ_hϕ(h)-ϕ_h), ∇ρ̅_h)=(Δ̃_hϕ_h, ρ̅_h)_h-(Δ̃_hϕ_h, ρ̅_h)+(∇ (ℐ_hϕ(h)-ϕ(h)), ∇ρ̅_h).

We now choose ρ̅_h=ℐ_hϕ(h)-ϕ_h to get

    [                         ∇(ℐ_hϕ(h)-ϕ_h)^2_L^2(Ω)                                               = (Δ̃_hϕ_h, ℐ_hϕ(h)-ϕ_h)_h-(Δ̃_hϕ_h, ℐ_hϕ(h)-ϕ_h);                                                           +(∇ (ℐ_hϕ(h)-ϕ(h)), ∇ (ℐ_hϕ(h)-ϕ_h)). ]

By (<ref>) and (<ref>), we have

    |(Δ̃_hϕ_h, ℐ_hϕ(h)-ϕ_h)_h-(Δ̃_hϕ_h, ℐ_hϕ(h)-ϕ_h)|≤ C_ com h Δ̃_hϕ_h_L^2(Ω)∇(ℐ_hϕ(h)-ϕ_h)_L^2(Ω)

and

    (∇ (ℐ_hϕ(h)-ϕ(h)), ∇ (ℐ_hϕ(h)-ϕ_h))≤ C_ app h Δ̃_hϕ_h_L^2(Ω)∇ (ℐ_hϕ(h)-ϕ(h))_L^2(Ω).

Consequently, estimate (<ref>) is satisfied by inserting (<ref>) and (<ref>) into (<ref>).


 There holds

    -Δ̃_hϕ_h_ L^2(Ω)≤C_ inv/h∇ϕ_h_L^2(Ω)

and

    ∇ϕ(h)_L^2(Ω)≤ C_ sta∇ϕ_h_ L^2(Ω).


 Select ρ̅_h=-Δ̃_hϕ_h in (<ref>) and use  (<ref>) to have

    Δ̃ϕ_h_L^2(Ω)^2≤∇ϕ_h_L^2(Ω)∇Δ̃_hϕ_h_L^2(Ω)≤C_ inv/h∇ϕ_h_L^2(Ω)Δ̃_hϕ_h_L^2(Ω),

which implies (<ref>). Inequality (<ref>) is obtained by using  (<ref>) and (<ref>), so we find that

    ∇ϕ(h)_L^2(Ω)≤∇ϕ_h_L^2(Ω)+C_ LaphΔ̃_hϕ_h_L^2(Ω)≤ (1+ C_ Lap C_ inv)∇ϕ_h_L^2(Ω).


In order to construct a proper sequence of initial approximations we need an interpolation operator that preserves non-negativity and has L^p-stability.  Let 𝒮𝒵_h be the variant of the Scott-Zhang interpolation operator defined in <cit.>, which satisfies the following.
 For p∈[1,∞], s=0,1, and m=0,1, there exist two constants C_ sta, C_ app>0, independent of h, such that

    𝒮𝒵_hφ_W^s,p(Ω)≤ C_ staφ_W^s,p(Ω)    φ∈ W^s,p(Ω).

and

    φ-𝒮𝒵_hφ_W^s,p(Ω)≤ C_  app h^m+1-sφ_W^m+1,p(Ω)    φ∈ W^m+1,p(Ω).

Moreover,

    φ≥0    Ω,    𝒮𝒵_hφ≥ 0    Ω.




Henceforth C denotes a generic constant whose value may change at each occurrence. This constant may depend on the data problem and the constants C_ neg, C_ inv, C_ eq,  C_ app, C_ Lip,  C_ com, and C_ Lap.


§ STATEMENT OF THE MAIN RESULT

Let ρ^0∈ L^∞(Ω) be nonnegative and consider ρ^0_h=𝒮𝒵_hρ^0. From (<ref>) and (<ref>),  we see that

    ρ^0_h≥0    Ω

and

    ρ_h^0_L^p(Ω)≤ C_ staρ^0_L^p(Ω)     p∈[1,∞].

Moreover, a regularization argument together with (<ref>) provides

    ρ^0_h→ρ^0     L^p(Ω)   h→0.


Let k=T/N with N∈N and consider {t_n}_n=0^N with t_n=k n.
Given ρ^n_h∈ D_h, compute ρ^n+1_h∈ D_h satisfying

    (δ_tρ^n+1_h,ρ̅_h)_h+ h^γ (∇ρ^n+1_h,∇ρ̅_h)+(∇ℐ_h  A([ρ^n+1_h]_T),∇ρ̅_h)-(ρ^n+1_h∇ℐ_h( K*[ρ^n_h]_T),∇ρ̅_h)=0,

where 0<γ<1  and [·]_T:D_h→ D_h is a nodal truncating operator defined as

    [ρ̅_h()]_T={[                   0                         ρ̅_h()∈(-∞, 0),;              ρ̅_h()                      ρ̅_h()∈ [0,B_L^∞],;               B_L^∞                     ρ̅_h()∈(B_L^∞, +∞), ].

with ∈𝒩_h and B_L^∞= e^ TΔ K_L^∞(R^d)ρ_0_L^1(Ω)ρ_0_L^∞(Ω). By an abuse of notation, the convolution K*[ρ_h,k]_T must be understood for [ρ_h,k]_T χ_Ω, where χ_Ω is the characteristic function. Moreover, the definition of B_L^∞ will be explained later on.

Due to the embedding W^2,∞(^d) into C^0(^d), the convolution K*[ρ_h,k]_T belongs to C^0(^d); therefore, we are allowed to consider ℐ_h ((K*[ρ_h,k]_T)|_Ω) that we write ℐ_h (K*[ρ_h,k]_T)  to simplify notation. Further, we have introduced δ_tρ^n+1_h=ρ^n+1_h-ρ^n_h/k.


For future references, note that ℐ_hA([ρ_h^n+1]_T)=ℐ_h A_T(ρ^n+1_h), where

    A_T(s)={[             0                   s∈(-∞,0),;          A(s)               s∈[0, B_L^∞],;      A(B_L^∞)               s∈(B_L^∞,+∞). ].


A weak solution for (<ref>) will be understood in the following sense <cit.>.
 A function ρ : Q→ [0, ∞) is a weak solution to (<ref>) with (<ref>) and (<ref>) if

    ρ∈ L^∞(Q ),      A(ρ) ∈ L^2(0, T, H^1(Ω)),   ∂_tρ∈ L^2(0,T; (H^1(Ω))'),

and

    {[ ∂_tρ -Δ A(ρ)+∇·(ρ(∇ K*ρ))                         =                         0                                L^2(0,T; (H^1(Ω))'),;                      ρ(0)                         =                       ρ_0                                          (H^1(Ω))'. ].






To establish convergence of the discrete solutions constructed via scheme (<ref>) toward the weak solution to (<ref>), we need to assume that


  (A1) A is  C^1([0,∞); [0,∞))  with A'>0  on (0, ∞) and A(0) = 0,

and


  (K1)  K ∈ W^2,∞(^d) is  such that K()=r(||) with r being nonincreasing.







Let us define ρ_h,k, ρ_h,k^-, ρ_h,k^+ :[0,T]→ D_h such that

    ρ_h,k=t-t_n+1/kρ^n+1_h+t_n-t/kρ^n_h,    t∈[t_n, t_n+1],


    ρ_h,k^-=ρ_h^n,   ρ_h,k^+=ρ^n+1_h,    t∈(t_n,t_n+1].


Our main result is summarized in the following theorem.
 Suppose that (A1), (K1), and (H1)-(H3) are satisfied. Then


  * there is a unique solution, ρ_h^n+1, to scheme (<ref>) provided that

    Ck (1 + 1/h) K_W^2,∞(^d)ρ_h^n_L^1(Ω)≤1/2.


  * ρ^n+1_h≥ 0
provided that

    C h^1-γK_W^2,∞(^d)ρ_h^n_L^1(Ω)< C_ neg.


  * and

    ∫_Ωρ^n+1_h()=∫_Ωρ^n_h().



Thus, the sequences of approximate solutions {ρ_h,k}_h,k>0 and {ρ_h,k^±}_h,k>0 constructed via scheme (<ref>)


  * are well-defined if

    Ck (1 + 1/h) K_W^2,∞(^d)ρ_h^0_L^1(Ω)≤1/2


  * satisfy

    ρ_h,k, ρ_h,k^±≥0       Q,

if

    C h^1-γ∇ K_W^2,∞(^d)ρ_h^0_L^1(Ω)< C_ neg.


  * and

    ∫_Ωρ_h,k(t,)=∫_Ωρ_h,k^±(t,)=∫_Ωρ_0h()     t∈[0,T].


Furthermore, the sequences of approximate solutions {ρ_h,k}_h,k>0 and {ρ_h,k^±}_h,k>0 converge toward the unique weak solution ρ of (<ref>), as (h,k)→ (0,0), in the sense that

    [ρ_h,k]_T, [ρ_h,k^±]_T→ρ L^2(0,T; L^p(Ω))

and

    A([ρ_h,k^+]_T)→ A(ρ)  L^2(0,T; L^p(Ω)) L^2(0,T; H^1(Ω))

with 1<p<∞.




From now on we assume that assumptions (A1), (K1), and (H1)-(H3) hold  without further comment on the statement of the results.




§ EXISTENCE AND UNIQUENESS OF DISCRETE SOLUTIONS 

In this section we prove the unique solvability of  scheme (<ref>). To simplify notation we suppress the superscript in ρ^n+1_h since there will be no ambiguity in setting ρ_h=ρ_h^n+1.

Before proceeding, we need an auxiliary result concerning the sign of ∇(K*[ρ^n_h]_T)·n on ∂Ω.
 It follows that

    ∇ (K*[ρ^n_h]_T)·≤ 0    ∂Ω.


 Let ∈∂Ω be such that () is well-defined at ∈∂Ω as being the outward unit normal vector and let s>0.  Write

    [      K*[ρ_h^n]_T(+  s )-K*[ρ_h^n]_T()                                     =     ∫_Ω (K(+  s -)-K (-))[ρ^n_h]_T()d;                                                                           = ∫_Ω (r(|+  s -|)-r(|-|))[ρ^n_h]_T()d. ]

In virtue of the decreasing property from (K1) and the convexity from (H1), we find that r(|+  s -|)-r(|-|)≤ 0 since |+  s -|≥ |-| and [ρ^n_h]_T≥ 0. Then

    ∂_ (K*[ρ^n_h]_T)()=lim_s→ 0^+(K*[ρ^n_h]_T)(+  s )-(K*[ρ^n_h]_T)()/s≤ 0.





The next lemma shows that scheme (<ref>) has at least one solution.  In doing so, we make use of Brouwer's theorem.


 Let ρ^n_h∈ L^1(Ω) be such that ρ_h^n≥0 in Ω.  Then scheme (<ref>) has at least one solution provided that

    C k K_W^2,∞(^d)ρ_h^n_L^1(Ω)≤1/2.

for a certain constant C>0 independent of h.

 
Let Φ: D_h→ D_h be defined  by  Φ(ρ̃_h)=ρ_h such that

    [ (ρ_h-ρ_h^n,ρ̅_h)_h+k h^γ(∇ρ_h,∇ρ̅_h)+k(∇ℐ_h  A_T( ρ̃_h),∇ρ̅_h)-k(ρ_h∇ℐ_h(K*[ρ^n_h]_T),∇ρ̅_h)=0. ]

Pick ρ̅_h=ρ_h to get

    [ ρ_h^2_h+k  h^γ∇ρ_h^2_L^2(Ω)= (ρ_h^n,ρ_h)_h -k (∇ℐ_h A_T(ρ̃_h), ∇ρ_h)+ k (ρ_h ∇ℐ_h (K*[ρ_h^n]_T), ∇ρ_h). ]

Cauchy-Schwarz' and Young's inequalities give

    (ρ_h^n,ρ_h)_h≤ 2 ρ_h^n^2_h+1/4ρ_h^2_h.

The second term on the right-hand of (<ref>) can be estimated on noting (<ref>) as

    [                                k (∇ℐ_h  A_T(ρ̃_h), ∇ρ_h)                                                        ≤                        k∇ℐ_h A_T(ρ̃_h)_L^2(Ω)∇ρ_h_L^2(Ω);                                                                                                                 ≤ 1/2k/h^1+γℐ_h  A_T(ρ̃_h)^2_L^2(Ω)+ 1/2k h^γ∇ρ_h^2_L^2(Ω);                                                                                                                 ≤      1/2 |Ω| k/h^1+γ  A(B_L^∞)+ 1/2 k  h^γ∇ρ_h^2_L^2(Ω). ]

The third term on the right-hand side of (<ref>) can be rewritten as

    k (ρ_h ∇ℐ_h (K*[ρ_h^n]_T), ∇ρ_h)= k (ρ_h ∇ K*[ρ_h^n]_T, ∇ρ_h) +k (ρ_h ∇(ℐ_h-ℐ)(K*[ρ_h^n]_T), ∇ρ_h),

where ℐ is the identity operator. Integrating by parts and using (<ref>) leads to

    [                                k (ρ_h ∇ K*[ρ_h^n]_T, ∇ρ_h)                                                          = - k/2(Δ K*[ρ_h^n]_T,  ρ_h^2)+ k((∇ K*[ρ^n_h]_T)·,ρ_h^2)_∂Ω;                                                                                                                     ≤                   C k Δ K_L^∞(^d)ρ_h^n_L^1(Ω)ρ_h^2_L^2(Ω), ]

where we have used the fact that [ρ_h^n]_T_L^1(Ω)≤ρ_h^n_L^1(Ω) since ρ_h^n≥0. Now note that (<ref>) and (<ref>) imply that

    [                       k (ρ_h ∇(ℐ_h-ℐ)(K*[ρ_h^n]_T), ∇ρ_h)                                                         ≤ C k h 
    ∇^2 K_L^∞(^d)[ρ_h^n]_T_L^1(Ω)ρ_h_L^2(Ω)∇ρ_h_L^2(Ω);                                                                                                                   ≤                 C k∇^2 K_L^∞(^d)ρ_h^n_L^1(Ω)ρ_h^2_L^2(Ω). ]

Combining (<ref>) and (<ref>), we find

    k (ρ_h ∇ℐ_h(K*[ρ_h^n]_T), ∇ρ_h)≤ C k K_W^2,∞(^d)ρ_h^n_L^1(Ω)ρ_h^2_L^2(Ω).

Putting (<ref>), (<ref>), (<ref>), and (<ref>) together, we arrive at the estimate

    3/2ρ_h^2_h+k  h^γ∇ρ_h^2_L^2(Ω)≤4ρ^n_h^2_h+ |Ω| k/h^1+γ  A(B_L^∞)+ C k K_W^2,∞(^d)ρ_h^n_L^1(Ω)ρ_h^2_L^2(Ω)

and, in view of (<ref>) and (<ref>),

    ρ_h_*=:ρ_h^2_h+k h^γ∇ρ_h^2_L^2(Ω)≤ 4ρ^n_h^2_h+ k/h^1+γ|Ω| A(B_L^∞):=R.

As a result, if we choose r>R, we find that ρ̃_h_*<r implies that Φ(ρ̃_h)_*<r.


Let us see that Φ is a continuous mapping from D_h into D_h with respect to the ·_*-norm. Suppose that ρ̃_h,m→ρ̃_h in the ·_*-norm as m→∞. Then we want to prove that Φ(ρ̃_h,m)→Φ(ρ̃_h) in the ·_*-norm as m→ 0. To do this, we compare (<ref>) and (<ref>) for ρ_h=ρ̃_h,m, and test against ρ_h=ρ̃_h,m-ρ_h to get

    [                                              Φ(ρ̃_h,m)-Φ(ρ̃_h)^2_*                                                                  =           -k(∇ℐ_h A_T(ρ̃_h,m)-∇ℐ_h A_T(ρ̃_h),∇(Φ(ρ̃_h,m)-Φ(ρ̃_h)));                                                                    +k ((Φ(ρ̃_h,m)-Φ(ρ̃_h)) ∇ℐ_h (K*[ρ_h^n]_T), ∇(Φ(ρ̃_h,m)-Φ(ρ̃_h))). ]


It is straightforward to see that

    -k(∇ℐ_h A_T(ρ̃_h,m)   -∇ℐ_h A_T(ρ̃_h),∇(Φ(ρ̃_h,m)-Φ(ρ̃_h)))
    
       ≤ k∇ℐ_h A_T(ρ̃_h,m)-∇ℐ_h A_T(ρ̃_h)_L^2(Ω)∇(Φ(ρ̃_h,m)-Φ(ρ̃_h))_L^2(Ω)
       ≤1/2k/h^γ∇ℐ_h A_T(ρ̃_h,m)-∇ℐ_h A_T(ρ̃_h)^2_L^2(Ω)+1/2 k  h^γ∇(Φ(ρ̃_h,m)-Φ(ρ̃_h))^2_L^2(Ω).

and, from (<ref>),

    k ((Φ(ρ̃_h,m)   -Φ(ρ̃_h)) ∇ℐ_h (K*[ρ_h^n]_T), ∇(Φ(ρ̃_h,m)-Φ(ρ̃_h)))
    
    ≤   C k  K_W^2,∞(^d)ρ_h^n_L^1(Ω)Φ(ρ̃_h,m)-Φ(ρ̃_h)^2_L^2(Ω).

Therefore, under (<ref>), we finally get

    Φ(ρ̃_h,m)-Φ(ρ̃_h)^2_*≤k/h^γ∇ℐ_h A_T(ρ̃_h,m)-∇ℐ_h A_T(ρ̃_h)^2_L^2(Ω).

As we are dealing with a finite-dimensional space, all norms are equivalent in D_h; and therefore we infer that ρ̃_h,m→ρ̃_h in  C(Ω̅) as m→∞. Since A_T is a continuous operator, we obtain that A_T(ρ̃_h,m)→  A_T(ρ̃_h) in C(Ω̅) as m→∞. This gives that  ℐ_h A_T(ρ̃_h,m)→ℐ_h A_T(ρ̃_h) in H^1(Ω) as m→∞. Now the continuity of Φ is obvious.

Next apply the Brouwer fixed-point theorem to conclude the proof.


Once we have proved existence, we turn to the question of uniqueness.

 Let ρ^n_h∈ L^1(Ω) such that ρ_h^n≥0 in Ω. Then  scheme (<ref>) possesses at most one solution provided that
(<ref>) holds.

 Suppose that there are two solutions ρ^1_h and ρ_h^2, respectively. Define ρ_h=ρ_h^1-ρ_h^2 which satisfies

    1/k(ρ_h,ρ̅_h)_h+ h^γ (∇ρ_h,∇ρ̅_h)+(∇ℐ_h(A([ρ_h^1]_T)-A([ρ_h^2]_T)), ∇ρ̅_h)-(ρ_h ∇ℐ_h(K*[ρ^n_h]_T),∇ρ̅_h)=0.

Let us define ϕ_h∈ D_h such that

    (∇ϕ_h,∇ρ̅_h)=(ρ_h,ρ̅_h)_h    ρ̅_h∈ D_h.

Select ρ̅_h=ϕ_h in (<ref>)  to get

    ∇ϕ_h^2_L^2(Ω)+ k h^γρ_h_h^2=-k(∇ℐ_h(A([ρ_h^1]_T)-A([ρ_h^2]_T)), ∇ϕ_h)+k(ρ_h ∇ℐ_h (K*[ρ_h^n]_T), ∇ϕ_h).

The first term on the right-hand side has negative sign. Indeed, by (<ref>), the mean-value theorem and (A1),

    -k(∇ℐ_h(A([ρ_h^1]_T)-A([ρ_h^2]_T)), ∇ϕ_h)   =-k(ℐ_h (A([ρ_h^1]_T)-A([ρ_h^2]_T)), ρ_h )_h
    
       =
    
    -k∑_∈𝒩_h A'(ξ_)([ρ_h^1()]_T-[ρ_h^2()]_T) ρ_h()  ∫_Ωφ_a() 
       ≤
    -k∑_∈𝒩_h A'(ξ_) ([ρ_h^1()]_T-[ρ_h^2()]_T)^2  ∫_Ωφ_a() ≤ 0,

where ξ_∈ ([ρ_h^1()]_T,[ρ_h^2()]_T) or ([ρ_h^2()]_T,[ρ_h^1()]_T).
For the second term, we proceed as follows. Combing  (<ref>)  and (<ref>), we have -Δ̃_hϕ_h=ρ_h. Thus, by
(<ref>), we write

    [              k(ρ_h∇ℐ_h( K*[ρ^n_h]_T), ∇ϕ_h)                                           =               -k(Δϕ(h)∇ K*[ρ^n_h]_T, ∇ϕ(h));                                                  -k(Δ̃ϕ_h∇ (ℐ_h-ℐ)(K*[ρ^n_h]_T), ∇ϕ(h));                                             -k(Δ̃_hϕ_h∇ℐ_h (K*[ρ^n_h]_T), ∇(ϕ_h-ϕ(h))). ]

Integration by parts shows that

    [         -k(Δϕ(h)∇ K*[ρ^n_h]_T, ∇ϕ(h))                                     =      k((∇ϕ(h)·∇)∇ K*[ρ_h^n]_T, ∇ϕ(h));                                            +k((∇ K*[ρ_h^n]_T·∇)∇ϕ(h),∇ϕ(h));                                                                           =      k((∇ϕ(h)·∇)∇ K*[ρ_h^n]_T, ∇ϕ(h));                                             - k/2(Δ K*[ρ_h^n]_T, |∇ϕ(h)|^2);                                       +k/2(|∇ϕ(h)|^2, (∇ K*[ρ_h^n]_T)·)_∂Ω, ]

which, from (<ref>) and (<ref>), gives

    -k(Δϕ(h)∇ℐ_h(K*[ρ^n_h]_T), ∇ϕ(h))≤ C k Δ K_L^∞(^d)ρ_h^n_L^1(Ω)∇ϕ_h^2_L^2(Ω).

In view of (<ref>), (<ref>) (<ref>), we have

    -k(Δϕ(h)∇ (ℐ_h-ℐ)(K*[ρ^n_h]_T), ∇ϕ(h))≤ C k ∇^2 K_L^∞(^d)ρ_h^n_L^1(Ω)∇ϕ_h^2_L^2(Ω).

Using (<ref>), (<ref>), and (<ref>),  leads to the estimate

    [   -k(Δ̃_hϕ_h∇ℐ_h( K*[ρ^n_h]_T), ∇(ϕ_h-ϕ(h)))                                            ≤ k  h Δ̃_hϕ_h^2_L^2(Ω)∇ K_L^∞(^d)ρ^n_h_L^1(Ω);                                                                                         ≤    Ck/h∇ K_L^∞(^d)ρ_h^n_L^1(Ω)∇ϕ_h^2_L^2(Ω). ]

Therefore,

    k(ρ_h∇ℐ_h (K*[ρ^n_h]_T), ∇ϕ_h)≤ Ck (1 + 1/h)  K_W^2,∞(^d)ρ_h^n_L^1(Ω)∇ϕ_h^2_L^2(Ω).


From (<ref>) and (<ref>), we estimate (<ref>) as

    ∇ϕ_h^2_L^2(Ω)+k h^γρ_h_h^2≤ Ck (1 + 1/h) K_W^2,∞(^d)ρ_h^n_L^1(Ω)∇ϕ_h_L^2(Ω)^2.

The result follows by taking k/h small enough so that (<ref>) holds.


It should be noted that condition (<ref>) is indeed more demanding than condition (<ref>) concerning the space and time parameters.



§ NON-NEGATIVITY AND A PRIORI ESTIMATES

In this section we show that the discrete solution ρ^n+1_h computed by (<ref>)  is nonnegative. Moreover, we derive some a priori energy estimates.

 Let ρ^n_h∈ L^1(Ω) be such that ρ_h^n≥0 in Ω. Assume that (<ref>) is satisfied. Then the solution ρ^n+1_h to scheme (<ref>) is nonnegative provided that (<ref>) holds.

 First of all, note that, for all E∈ℰ_h and for all _i,_j∈ E with i≠j,

    [                     ∫_Eφ__i∇ℐ_h (K*[ρ_h^n]_T)·∇φ__j                                                   ≤ |E| φ__j_L^∞(E)∇ℐ_h(K*[ρ_h^n]_T)_L^∞(E)∇φ__j_L^∞(E);                                                                                                       ≤                     C h^d-1∇ K_L^∞(^d)ρ_h^n_L^1(E), ]

where we used (<ref>) and (<ref>).  Comparing (<ref>) with (<ref>), we find that

    h^γ∫_E∇φ__i·∇φ__j   - ∫_Eφ__i∇ℐ_h(K*[ρ_h^n]_T)·∇φ__j
       ≤ h^d-2 h^γ (-C_ neg + C h^1-γ∇ K_L^∞(^d)ρ_h^n_L^1(E))<0

holds if we let C h^1-γ∇ K_L^∞(^d)ρ_h^n_L^1(E)< C_ neg, which is a consequence of (<ref>). As a result, summing over E∈ supp  φ__i∩ supp  φ__j yields

    h^γ (∇φ__i, ∇φ__j) - (φ__i∇ℐ_h (K*[ρ_h^n]_T),∇φ__j)<0.

Analogously, we have, from (<ref>), that

    h^γ(∇φ__i, ∇φ__i) - (φ__i∇ℐ_h (K*[ρ_h^n]_T),∇φ__i)>0

holds if we let C h^1-γ∇ K_L^∞(^d)ρ_h^n_L^1(E)< C_ neg, which is globally imposed in (<ref>).

Now let ρ_h^ min∈ D_h be defined as

    ρ_h^ min=∑_∈𝒩_hρ_h^-()φ_,

where ρ_h^-()=min{0, ρ^n+1_h()}. Analogously, one defines ρ^ max_h∈ D_h as

    ρ_h^ max=∑_∈𝒩_hρ_h^+()φ_,

where ρ_h^+()=max{0, ρ^n+1_h()}. Notice that ρ^n+1_h=ρ_h^ min+ ρ_h^ max. Set ρ̅_h=ρ^ min_h in (<ref>) to get

    [                                                     (δ_tρ^n+1_h,ρ^ min_h)_h                                                  + h^γ (∇ρ^n+1_h,∇ρ^ min_h);                                                                             +(∇ℐ_h  A([ρ^n+1_h]_T),∇ρ^ min_h)-(ρ^n+1_h∇ℐ_h( K*[ρ^n_h]_T),∇ρ^ min_h )=0. ]

We will handle each term of (<ref>) in order to show that ρ_h^ min≡0. Indeed, in virtue of the equality

    (ρ^n+1_h, ρ^ min_h)_h=(ρ^ min_h+ρ^ max_h, ρ^ min_h)_h=ρ^ min_h_h^2,

it follows that

    (δ_tρ^n+1_h,ρ^ min_h)_h= 1/k(ρ_h^ min_h^2-(ρ^n_h, ρ^ min_h))≥1/kρ_h^ min_h^2.

Observe that we have

    [                 (∇ℐ_h  A([ρ^n+1_h]_T),∇ρ^ min_h)                                                =   ∑_≠∈𝒩_h A([ρ_h^n+1()]_T) ρ^ min_h() (∇φ_, ∇φ_);                                                    +∑_∈𝒩_h A([ρ_h^n+1()]_T) ρ^ min_h() (∇φ_, ∇φ_);                                                                                                 = ∑_≠∈𝒩_h A([ρ_h^n+1()]_T) ρ^ min_h() (∇φ_, ∇φ_)>0 ]

from  (<ref>) and A([ρ_h^n+1()]_T) ρ^ min_h()≤0. By the decomposition

    h^γ (∇ρ^n+1_h,∇ρ^ min_h)   -(ρ^n+1_h∇ℐ_h( K*[ρ^n_h]_T),∇ρ^ min_h )
    
    
    =   h^γ (∇ρ^ max_h,∇ρ^ min_h)-(ρ^ max_h∇ℐ_h(K*[ρ^n_h]_T),∇ρ^ min_h )
    
       +h^γ (∇ρ^ min_h,∇ρ^ min_h)-(ρ^ min_h∇ℐ_h(K*[ρ^n_h]_T),∇ρ^ min_h ),

we deduce from (<ref>) and (<ref>)  that

    h^γ (∇ρ^ max_h,   ∇ρ^ min_h)-(ρ^ max_h∇ℐ_h (K*[ρ^n_h]_T),∇ρ^ min_h )
    
    
    =   ∑_≠∈𝒩_hρ^ max_h() ρ^ min_h()[ (h^γ (∇φ_,∇φ_ )-(φ_∇ℐ_h(K*[ρ^n_h]_T),∇φ_)]
    
       
    +∑_∈𝒩_hρ^ max_h() ρ^ min_h()[ (h^γ (∇φ_,∇φ_)-(φ_∇ℐ_h(K*[ρ^n_h]_T),∇φ_)]≥0

since ρ^ max_h() ρ^ min_h()≤0 and ρ^ max_h() ρ^ min_h()=0. Therefore,

    [                                      h^γ (∇ρ^ min_h,∇ρ^ min_h)                                                              -                         (ρ^ min_hℐ_h(∇ K*[ρ^n_h]_T),∇ρ^ min_h);                                                                                                                             ≤ h^γ (∇ρ^n+1_h,∇ρ^ min_h)-(ρ^n+1_h∇ℐ_h(K*[ρ^n_h]_T),∇ρ^ min_h). ]

As a result, we infer on applying (<ref>)-(<ref>)  into (<ref>)  that

    ρ^ min_h^2_h+k h^γ∇ρ^ min_h^2≤ k (ρ^ min_h∇ℐ_h(K*[ρ^n_h]_T),∇ρ^ min_h ).

We know from (<ref>) and (<ref>) that

    ρ^ min_h^2_h+k h^γ∇ρ^ min_h^2≤ C k  K_W^2,∞(^d)ρ^n_h_L^1(Ω)ρ^ min_h^2_h.

Thus, from (<ref>),

    ρ_h^ min_h^2≤0,

which implies that ρ^ min_h≡0 and hence  ρ^n+1_h≥ 0. It completes the proof.




Since we do not have a pointwise upper bound for ρ^n+1_h, we must slightly modify the argument leading to a priori energy estimates from <cit.>, which uses the maximum principle. 
 Assume that (<ref>) and (<ref>) are satisfied. Then  the sequence {ρ^n_h}_n=1^N computed via scheme (<ref>) satisfies

    ρ^n+1_h_L^1(Ω)=ρ_h^0_L^1(Ω):=B_L^1

and

    [ ρ^n+1_h^2_h+ ∑_m=0^n(k^2 ∂_tρ^n+1_h^2_h+k  h^γ∇ρ^m+1_h^2_L^2(Ω)+k∇ℐ_hA_T(ρ^m+1_h)^2_L^2(Ω));                                                   ≤ e^T B_L^1K_W^2,∞(^d)ρ^0_h^2_h:=B^2_L^2. ]


 We proceed by induction on n to  prove (<ref>). From (<ref>), we know that ρ^1_h≥ 0 is true by Lemma <ref> for (<ref>) and (<ref>). On selecting ρ̅_h=1 in (<ref>), we obtain (<ref>) for n=0. The same argument leads us to proving that (<ref>) holds from ρ^n_h≥ 0 and ρ^n_h_L^1(Ω)=ρ_0h_L^1(Ω) by induction hypothesis. At this point, it should be noted that (<ref>) and (<ref>)  combined with (<ref>) imply (<ref>) and (<ref>).

Now let ρ̅_h=ρ^n+1_h in (<ref>) to get

    ρ^n+1_h_h^2   +ρ_h^n+1-ρ_h^n_h^2+ 2 k h^γ∇ρ^n+1_h^2_L^2(Ω)+ 2 C_ Lip^-1 k∇ℐ_hA_T(ρ^n+1_h)^2_L^2(Ω)
       ≤ρ^n_h_h^2+2 k (ρ^n+1_h∇ℐ_h(K*[ρ^n_h]_T),∇ρ_h^n+1),

where we have used (<ref>) for f=A_T being non-decreasing and Lipschitzian. Repeating the argument that led to estimating (<ref>) and noting (<ref>) and (<ref>) yields

    ρ^n+1_h_h^2   +ρ_h^n+1-ρ_h^n_h^2+ h^γ∇ρ^n+1_h_L^2(Ω)^2+C_ Lip^-1∇ℐ_hA_T(ρ^n+1_h)^2_L^2(Ω)
       ≤ρ^n_h_h^2+ C k K_W^2,∞(^d)ρ^0_h_L^1(Ω)ρ_h^n+1^2_h.

By a discrete Grnwall lemma, we conclude that (<ref>) holds under condition (<ref>).


The constants B_L^1 and B_L^2 can be estimated uniformly with respect to h in term of ρ^0 from (<ref>).
 It follows that

    k∑_n=0^N-1δ_tρ^n+1_h^2_(H^1(Ω))'≤ C,

where C>0 is a constant independent of h.

 Apply a standard duality technique to obtain  (<ref>) from (<ref>) and (<ref>).

We end this section by summarizing the results of Lemmas <ref> and <ref>. Bounds (<ref>) and (<ref>) yield that

    ρ_h,k, ρ_h,k^±≥0      Q,


    {ρ_h,k}_h,k, {ρ_h,k^±}_h,k  L^∞(0,T; L^2(Ω)),


    { h^γ/2ρ_h,k^+}_h,k, {ℐ_hA([ρ_h,k^+]_T)}_h,k L^2(0,T; H^1(Ω)),

and

    {ρ_h,k}_h,k H^1(0,T; (H^1(Ω))')

and, by passing to the limit in a subsequence, denoted by (k,h) for convenience,  that there exists ρ∈ L^∞(0,T; L^2(Ω)) such that

    ρ_h,k, ρ_h,k^±→ρ     L^∞(0,T; L^2(Ω))*,

and

    ρ_h,k→ρ     H^1(0,T; (H^1(Ω))')

as (h,k)→ (0,0). Moreover, there exists χ∈ L^2(0,T, H^1(Ω)) such that

    ℐ_h A([ρ_h,k^+]_T)→χ     L^2(0,T; H^1(Ω))

as (h,k)→ (0,0).



§ COMPACTNESS

As we are dealing with a nonlinear equation, the key ingredient in passing to the limit is obtaining  compactness of the discrete solutions computed using (<ref>). Since we do not have control of the gradient of the discrete solutions due to the degenerate diffusion term, compactness turns out to be more complicated to achieve than in the non-degenerate case. We have split the proof into a series of four lemmas.
 There exists a nonincreasing function F_1: [0,∞)→ [0,∞) with F_1(z)→ 0 as z→ 0 such that for any sequence of  discrete solutions {ρ_h,k^+}_h,k computed via scheme (<ref>) satisfies

    [ρ^+_h,k]_T(t+δ)-[ρ^+_h,k]_T(t)_L^2(Ω)≤ F_1(ℐ_hA([ρ^+_h,k]_T(t+δ))-ℐ_hA([ρ^+_h,k]_T(t))_L^2(Ω)),

for all δ∈(0,T) and t∈[0,T-δ], and

    ∫_0^T∫_ω   |𝒫_h[ρ_h,k^+]_T(t,+δe_i)-𝒫_h[ρ_h,k^+]_T(t,)|^2 
       ≤ F_1[∫_0^T∫_ω|𝒫_hℐ_hA([ρ_h,k^+]_T(t,+δe_i))-𝒫_hℐ_hA([ρ_h,k^+]_T(t,))|^2 ],

for all ω⊂⊂Ω and 0<δ< dist(ω,∂Ω) with {_i}_i=1^d being the Cartesian basis of ^d.

 For x≥0 and y≥0,  define the following continuous function

    σ(x,y)={[ A(x)-A(y)/x-y                        x≠y,;         A'(x)                        x=y. ].
 
Let η>0 and consider  f(η)=min{σ (x,y): (x,y) ∈ [η, B_L^∞]×[0, B_L^∞]}. 
Then we have that f(η)>0 since A(x)'>0 for x>0 owing to (A1). Let

    𝒩^1_h={∈𝒩_h: |[ρ_h,k^+]_T(t+δ, )|<η |[ρ_h,k^+]_T(t, )|<η}

and
𝒩^2_h=𝒩_h\𝒩^1_h. Then, from (<ref>), we get

    [                           [ρ_h,k^+]_T(t+δ)-[ρ_h,k^+]_T(t)^2_L^2(Ω)                                                                  ≤                           C_ eq[ρ_h,k^+]_T(t+δ)-[ρ_h,k^+]_T(t)_h^2;                                                                                                                                     =          C_ eq∑_∈𝒩^1_h ([ρ_h,k^+]_T(t+δ,)-[ρ_h,k^+]_T(t, ))^2∫_Ωφ_;                                                                            +C_ eq∑_∈𝒩^2_h ([ρ_h,k^+]_T(t+δ,)-[ρ_h,k^+]_T(t, ))^2∫_Ωφ_;                                                                                                                                     ≤                                                       C_ eqη^2 |Ω|;                                                                    +C_ eqf^-2(η)ℐ_hA([ρ_h,k^+]_T(t+δ))-ℐ_hA([ρ_h,k^+]_T(t))^2_L^2(Ω). ]

Consider  F_μ(z)=inf_η>0{C_ eq|Ω|η^2+C_ eqμf^-2(η)z^2} with μ=1 to complete the proof of (<ref>).

For (<ref>), we reason along the same line as before. Define

    Q^1_ω={ (t,)∈ [0, T]×ω : |𝒫_h [ρ_h,k^+]_T(t,+δ_i) |< η | 𝒫_h[ρ_h,k^+]_T(t,) |< η}

and Q^2_ω=Q\ Q^1_ω. Estimating as before, we find

    ∫_0^T ∫_ω    |𝒫_h [ρ_h,k^+]_T(t,+δ_i)-𝒫_h[ρ_h,k^+]_T(t, )|^2
    
    
    =   ∫_Q^1_ω |𝒫_h [ρ_h,k^+]_T(t,+δ_i)-𝒫_h[ρ_h,k^+]_T(t, )|^2 
       +∫_Q^2_ω |𝒫_h [ρ_h,k^+]_T(t,+δ_i)-𝒫_h[ρ_h,k^+]_T)(t, )|^2  
    
    
    ≤   η |Q|+1/f^2(η)∫_0^T ∫_ω |𝒫_hℐ_hA([ρ_h,k^+]_T(t,+δ_i))-𝒫_hℐ_hA([ρ_h,k^+]_T(t, ))|^2 ,















which implies (<ref>). 

 Let δ∈ (0,T) and t∈[0,T-δ].  Assume that there exists B>0 such that the sequence of discrete solutions {ρ_h,k^+}_h,k computed via (<ref>) satisfies

    ℐ_hA ([ρ_h,k^+]_T(t+δ))_H^1(Ω)≤ B    ∇ℐ_hA ([ρ_h,k^+]_T(t))_H^1(Ω)≤ B.

and

    h^γ∇ρ_h,k^+(t+δ)_L^2(Ω)≤ B     h^γ∇ρ_h,k^+(t))_L^2(Ω)≤ B.

Then there exists a function G_B :[0,∞)→ [0,∞) being nondecreasing and satisfying G_B(ε)→ 0 as ε→ 0 such that

    [ρ_h,k^+(t+δ)]_T-[ρ_h,k^+(t)]_T^2_L^2(Ω)≤ G_B(ε)

providing that 
    ([ρ_h,k^+]_T(t+δ)-[ρ_h,k^+]_T(t), ℐ_h A([ρ_h,k^+]_T(t+δ))-ℐ_hA([ρ_h,k^+]_T(t)))_h≤ε

holds


We establish the lemma by contradiction. Assume that there exist κ>0 and two sequences {ρ_h_n,k_n^+(t+δ)}_n=1^∞ and {ρ_h_n,k_n^+(t)}_n=1^∞ such that

    ([ρ_h_n,k_n^+]_T(t+δ)-[ρ_h_n,k_n^+]_T(t), ℐ_h_n A([ρ_h_n,k_n^+]_T(t+δ))-ℐ_h_nA([ρ_h_n,k_n^+]_T(t)))_h_n≤1/n

and

    [ρ_h_n,k_n^+]_T(t+δ)-[ρ_h_n,k_n^+]_T(t)^2_L^2(Ω)>κ.

From (<ref>), we know that there exist w_1,w_2∈ L^2(Ω) and a subsequence of {[ρ_h_n,k_n^+]_T(t+δ)}_n=0^∞ and  { [ρ_h_n,k_n^+]_T(t)}_n=0^∞, still denoted by itself, such that

    ℐ_h_n A([ρ_h_n,k_n^+]_T(t+δ)) → A(ρ_1)     L^2(Ω)     n→∞

and

    ℐ_h_n A([ρ_h_n,k_n^+]_T(t)) → A(ρ_2)      L^2(Ω)     n→∞,

where ρ_1=A^-1(w_1) and ρ_2=A^-1(w_2). It is not hard to see from (<ref>) and (<ref>) that

    𝒫_h_nℐ_h_n A([ρ_h_n,k_n^+]_T(t+δ)) → A(ρ_1)     L^2(Ω)     n→∞

and

    𝒫_h_nℐ_h_n A([ρ_h_n,k_n^+]_T(t)) → A(ρ_2)      L^2(Ω)     n→∞.

Lebesgue's Dominated Convergence Theorem  implies that

    A^-1𝒫_h_nℐ_h_n A([ρ_h_n,k_n^+]_T(t+δ)) →ρ_1     L^2(Ω)     n→∞

and

    A^-1𝒫_h_nℐ_h_nA([ρ_h_n,k_n^+]_T(t))→ρ_2     L^2(Ω)     n→∞.


It is clear that

    𝒫_h_n[ρ_h_n,k_n^+]_T(t+δ)= A^-1𝒫_h_nℐ_h A([ρ_h_n,k_n^+]_T(t+δ))

and

    𝒫_h_n [ρ_h_n,k_n^+]_T(t)= A^-1𝒫_h_nℐ_h_n A([ρ_h_n,k_n^+]_T(t)).

In view of (<ref>) and (<ref>), we get

    [ρ_h_n,k_n^+]_T(t+δ)-𝒫_hρ_h_n,k^+]_T(t+δ)_L^2(Ω)≤ C h_n ∇ρ_h_n,k_n^+(t+δ)_L^2(Ω)≤ C h_n^1-γ B

and

    ρ_h_n,k_n^+(t)-𝒫_hρ_h_n,k_n^+(t)_L^2(Ω)≤ C h_n ∇ρ_h_n,k_n^+(t) _L^2(Ω)≤ C h_n^1-γ B.

Here, we used the fact that |∇ [ρ^+_h_n,k_n]_T(·)|≤ |∇ [ρ^+_h_n,k_n]_T(·)|. Therefore,

    [ρ_h_n,k_n^+]_T(t+δ)→ρ_1     [ρ_h_n,k_n^+]_T(t)→ρ_2     L^2(Ω)     n→∞.

On noting  (<ref>), we have

    [                        ([ρ_h_n,k_n^+]_T(t+δ)-[ρ_h_n,k_n^+]_T(t), ℐ_h_n A([ρ_h_n,k_n^+]_T(t+δ))-ℐ_h_nA([ρ_h_n,k_n^+]_T(t))); ≤1/n+ C h_n [ρ_h_n,k_n^+]_T(t+δ)-[ρ_h_n,k_n^+]_T(t)_L^2(Ω)∇( ℐ_h_n (A([ρ_h_n,k_n^+]_T(t+δ))-A([ρ_h_n,k_n^+]_T(t))))_L^2(Ω);                                                                                                ≤1/n+C |Ω|^1/2 h_n B_L^∞ B. ]

Passing to the limit in this last estimate yields

    (ρ_1-ρ_2, A(ρ_1)-A(ρ_1))=0,

which implies that ρ_1=ρ_2. As a result, we have ℐ_h_n A([ρ_h_n,k_n]_T(t+δ))-ℐ_hA([ρ_h_n,k_n]_T(t))_L^2(Ω)→ 0 as n→ +∞. But then [ρ_h_n,k_n]_T(t+δ)-[ρ_h_n,k_n]_T(t)_L^2(Ω)→ 0 as n→ +∞ from Lemma <ref>, which is a contradiction from (<ref>).


In order to prove the following lemma, we draw on <cit.>.
 Let δ∈(0,T) and t∈[0,T-δ]. Then it follows that

    ∫_0^T-δ([ρ_h,k^+]_T(t+δ)-[ρ_h,k]_T^+(t),  A([ρ_h,k^+]_T(t+δ))- A([ρ_h,k^+]_T(t)))_h ≤ C δ.


 Since ρ_h,k^+ is a time–stepping function, we only need to consider δ= r k, with r=1,⋯, N, and prove

    k∑_m=0^N-r ([ρ_h^m+r]_T-[ρ_h^m]_T,  A([ρ_h^m+r]_T)- A([ρ_h^m]_T))_h ≤ C (r k)^1/2.

Let us test (<ref>) against ρ̅_h=ℐ_hA([ρ_h^m+r]_T)-ℐ_hA([ρ_h^m]_T) to obtain

    (ρ^n+1_h-ρ_h^n   , A([ρ_h^m+r]_T)-A([ρ_h^m]_T))_h
    
       =-h^γ k (∇ρ^n+1_h,∇ℐ_hA([ρ_h^m+r]_T)-∇ℐ_hA([ρ_h^m]_T))
    
       -k(∇ℐ_h A([ρ_h^n+1]_T), ∇ℐ_hA([ρ_h^m+r]_T)-∇ℐ_hA([ρ_h^m]_T))
    
       + k (ρ^n+1_h∇ℐ_h (K*[ρ_h^n]_T), ∇ℐ_hA([ρ_h^m+r]_T)-∇ℐ_hA([ρ_h^m]_T).

Summing for n = m, ..., m - 1 + r, multiplying by k and summing for m = 0, ⋯, N - r  yields

    k∑_n=0^N-r(ρ^m+r_h-ρ_h^m,    A([ρ_h^m+r]_T)-A([ρ_h^m]_T))_h
     
        =-h^γ k∑_n=0^N-rk∑_n=m^m-1+r (∇ρ^n+1_h,∇ℐ_hA([ρ_h^m+r]_T)-∇ℐ_hA([ρ_h^m]_T))
    
       - k∑_n=0^N-rk∑_n=m^m-1+r(∇ℐ_h A([ρ_h^n+1]_T), ∇ℐ_hA([ρ_h^m+r]_T)-∇ℐ_hA([ρ_h^m]_T))
    
       + k∑_n=0^N-rk∑_n=m^m-1+r (ρ^n+1_h∇ℐ_h(A*[ρ_h^n]_T), ∇ℐ_hA([ρ_h^m+r]_T)-∇ℐ_hA([ρ_h^m]_T).

We now proceed to bound each term on the right-hand side. In doing so, we first apply a Fubini discrete rule to write

    [    h^γ k∑_n=0^N-rk∑_n=m^m-1+r (∇ρ^n+1_h,∇ℐ_hA([ρ_h^m+r]_T)-∇ℐ_hA([ρ_h^m]_T)); =h^γ k∑_n=0^N-1k∑_m=n-1+r^n̅ (∇ρ^n+1_h,∇ℐ_hA([ρ_h^m+r]_T)-∇ℐ_hA([ρ_h^m]_T)), ]

where

    n̅={[           0                    n<0,;           n             0≤ n ≤ N-r,;         N-r                  n>N-r. ].

Therefore, using |n̅-n-r+1|≤ r, we have, by (<ref>), that

    [                                                           h^γ k∑_n=0^N-rk∑_n=m^m-1+r (∇ρ^n+1_h,∇ℐ_hA([ρ_h^m+r]_T)-∇ℐ_hA([ρ_h^m]_T)); ≤ h^γ^1/2 k∑_n=0^N-r h^γ^1/2 k ∇ρ^n+1_h_L^2(Ω)(∑_m=n-1+r^n̅ k∇(ℐ_hA([ρ_h^m+r]_T)-ℐ_hA([ρ_h^m]_T))^2_L^2(Ω))^1/2(∑_m=n-1+r^n̅ k)^1/2;                                                                                                        ≤ C B^2_L^2  T^1/2 (rk)^1/2. ]

Analogously, we bound

    k∑_n=0^N-rk∑_n=m^m-1+r(∇ℐ_hA([ρ_h^n+1]_T),ℐ_hA([ρ_h^m+r]_T)-ℐ_hA([ρ_h^m]_T))≤ C T^1/2 B^2_L^2 (rk)^1/2

and

    k∑_n=0^N-rk∑_n=m^m-1+r    (ρ^n+1_h∇ℐ_h(K*[ρ_h^n]_T), ∇ℐ_hA([ρ_h^m+r]_T)-ℐ_hA([ρ_h^m]_T)
    
       ≤ C T K_W^2,∞(^d)ρ_h^0_L^1(Ω) B^2_L^2 (rk)^1/2.

Combining these above estimates gives

    k∑_n=0^N-r(ρ^m+r_h-ρ_h^m, A([ρ_h^m+r]_T)-A([ρ_h^m]_T))_h≤ C(rk)^1/2.

The proof is now completed on noting that |[ρ^m+r_h]_T()-[ρ_h^m()]_T|≤ |ρ^m+r_h()-ρ_h^m()| for all ∈𝒩_h.

In order to set out that  the sequence of {[ρ_h,k^+]_T}_h,k is precompact in L^2(Q), we will use the Riesz-Frchet-Kolmogorov compactness criterion.
 It follows that

    [ρ_h,k^+]_T→ρ_T L^2(Q) (h,k)→ (0,0),

where ρ_T is the truncating of the limiting function ρ obtained from the weak convergences.

 The proof should be understood for the subsequence obtained in (<ref>)  and  (<ref>).  We divide the proof into two parts:

Part I: We claim that for each ε>0 there exists 0<δ_0≤ T such that for all (h,k)>0 and all 0 < δ < δ_0

    ∫_0^T-ε[ρ_h,k^+]_T(·+δ)-[ρ_h,k]_T^+(·)^2_L^2(Ω)< ε.


By Lemma <ref>, we know that

    A([ρ_h,k^+]_T)_L^2(0,T; H^1(Ω))≤ ( B_L^∞^2 |Ω|+ B_L^2^2)^1/2:=B.

Consider 0<δ<ε and θ>1 and define

    E_θ(δ)={t∈ [0,T-ε]:   A([ρ_h,k^+]_T(t))_H^1(Ω)≤ B θ^1/2,   A([ρ_h,k^+]_T(t+δ))^2_H^1(Ω)≤ B θ^1/2,
    
       h^γ∇ρ_h,k^+(t)_L^2(Ω)≤ Bθ^1/2,   h^γ∇ρ_h,k^+(t+δ)^2_L^2(Ω)≤ Bθ^1/2,
    
       ([ρ_h,k^+]_T(t+δ)-[ρ_h,k^+]_T(t), A([ρ_h,k^+]_T(t+δ))-A([ρ_h,k^+]_T(t)))_h≤ C  θδ}.

By Chebyshev's inequality, we deduce that |E^c_θ(δ)|≤5/θ, where E^c_θ(δ) is the complementary set of E_θ(δ). Therefore, by Lemma <ref> combined with (<ref>),

    ∫_0^T-ε[ρ_h,k^+]_T(·+δ)-[ρ_h,k^+]_T(·)^2_L^2(Ω)≤ T G_B θ^1/2(Cθδ)+2B^2_L^∞5/θ.

On choosing θ=max{20B^2_L^∞/ε, 1} and δ_0>0 such that T G_ B θ^1/2(Cθδ)<ε/2, this leads to (<ref>).

Part II: We claim that for each ε>0 and each ω⊂⊂Ω there exists 0<δ_0<  dist(ω,∂Ω) such that

    ∫_0^T∫_ω|[ρ_h,k^+]_T(t,+δe_i))-[ρ_h,k^+]_T(t,))|^2 < ε,

for all (h,k)>0, and all 0 < δ≤δ_0 and i=1,⋯, d.


Using Minkowski's inequality, we have

    ∫_0^T ∫_ω |[ρ_h,k^+]_T(t,+δ_i)   -[ρ_h,k^+]_T(t, )|^2 
       ≤ C  ∫_0^T  [ρ_h,k^+]_T-𝒫_h[ρ_h,k^+]_T^2_L^2(Ω)
       + C ∫_0^T ∫_ω |𝒫_h [ρ_h,k^+]_T(t,+δ_i)-𝒫_h[ρ_h,k^+]_T(t, )|^2 dt.

We estimate each term on the right-hand side separately.  We have, by (<ref>) and (<ref>),  that

    ∫_0^T  [ρ_h,k^+]_T-𝒫_h[ρ_h,k^+]_T^2_L^2(Ω)≤ C B_L^2^2 h^1-γ,

where we have used that that fact |∇[ρ_h]_T|≤ |∇ρ_h| for all ρ_h∈ D_h.

Now we want to use (<ref>) to control the second term. Observe that

    ∫_0^T ∫_ω   |𝒫_hℐ_h A([ρ_h,k^+]_T(t,+δ_i))-𝒫_hℐ_h A([ρ_h,k^+]_T(t, ))|^2_L^2(Ω) 
    ≤    C  ∫_0^T  ℐ_h A([ρ_h,k^+]_T)-𝒫_hℐ_h A([ρ_h,k^+]_T)^2_L^2(Ω)
       + C ∫_0^T ∫_ω |ℐ_h A ([ρ_h,k^+]_T(t,+δ_i))-ℐ_h A([ρ_h,k^+]_T(t, ))|^2 dt.

It is easily to check, from (<ref>), that

    ∫_0^T  ℐ_h A([ρ_h,k^+]_T)-𝒫_hℐ_h A([ρ_h,k^+]_T)^2_L^2(Ω)≤ C B_L^2^2 h^2

and, by the Mean-Value Theorem,  that

    ∫_0^T∫_ω|ℐ_hA([ρ_h,k^+]_T(t,+δe_i)   -ℐ_hA([ρ_h,k^+]_T(t,))|^2 
       ≤δ^2 ∫_0^T ∇ℐ_hA([ρ_h,k^+(t)]_T)^2_L^2(Ω)≤δ^2 B^2_L^2.

Thus, by Lemma <ref>,

    ∫_0^T ∫_ω |𝒫_h [ρ_h,k^+]_T(t,+δ_i)-𝒫_h[ρ_h,k^+]_T(t, )|^2 ≤ F_1( B^2_L^2 (C h^2+δ^2)).


Therefore,

    ∫_0^T∫_ω|ρ_h,k^+(t,+δe_i)-ρ_h,k^+(t,)|^2 ≤ C B_L^2 h^1-γ+ F_2(δ^2 B^2_L^2)+F_2(C B^2_L^2  h^2).

Following the proof of <cit.> we infer that (<ref>) holds.


Finally, inequalities (<ref>) and (<ref>) are sufficient to prove that the sequence of {[ρ_h,k]_T}_h,k is precompact via the Riesz-Frchet-Kolmogorov compactness criterion. It is not hard to see that the limiting function ρ_T is the truncating of ρ.

We further infer that 

    [ρ_h,k^-]_T→ρ_T L^2(Q)   (h,k)→ (0,0).


As a consequence of Lemma <ref>, we have the following.
 There holds

    ℐ_hA([ρ_h,k^+]_T)→ A(ρ_T)     L^2(0,T; H^1(Ω))   (h,k)→ (0,0).


 Using Minkowski's inequality a few times, we see that

    [                 ∫_0^Tℐ_hA([ρ_h,k^+]_T)-A(ρ_T)^2_L^2(Ω)                                                      ≤ C ∫_0^Tℐ_hA([ρ_h,k^+]_T)-𝒫_hℐ_hA([ρ_h,k^+]_T)^2_L^2(Ω);                                                           +C∫_0^T𝒫_hℐ_hA([ρ_h,k^+]_T)-A([ρ_h,k^+]_T)^2_L^2(Ω);                                                                       +C ∫_0^TA([ρ_h,k^+]_T)-A(ρ_T)^2_L^2(Ω). ]

In view of (<ref>) and (<ref>), we obtain

    ∫_0^Tℐ_hA([ρ_h,k^+]_T)-𝒫_hℐ_hA([ρ_h,k^+]_T)^2_L^2(Ω)≤ C h^2 ∫_0^T∇ℐ_hA([ρ_h,k^+]_T)^2_L^2(Ω)→ 0

and

    [        ∫_0^T𝒫_hℐ_hA([ρ_h,k^+]_T)-A([ρ_h,k^+]_T)^2_L^2(Ω)                                                        =           ∫_0^T𝒫_hA([ρ_h,k^+]_T)-A([ρ_h,k^+]_T)^2_L^2(Ω);                                                                                                                 ≤                    C h^2 ∫_0^T ∇ A([ρ_h,k^+]_T)^2_L^2(Ω);                                                                                                                 ≤ C h^2 A'([ρ_h,k^+]_T)_L^∞(Q)∫_0^T ∇ [ρ_h,k^+]_T^2_L^2(Ω);                                                                                                                 ≤          C h^2-γ A'(B_L^∞) ∫_0^T h^γ∇ρ_h,k^+^2_L^2(Ω)→ 0 ]

as (h,k)→ (0,0). Lebesgue's Dominated Convergence
Theorem combined with (<ref>) provides

    ∫_0^TA([ρ_h,k^+]_T)-A(ρ_T)^2_L^2(Ω)→ 0

as (h,k)→ (0,0). Thus, the above convergence gives

    ℐ_hA([ρ_h,k^+]_T)→ A(ρ_T)     L^2(Q)   (h,k)→ (0,0).

Furthermore, it follows from (<ref>) that (<ref>) is satisfied; thus completing the proof.



§ PASSAGE TO THE LIMIT 

We briefly outline the main steps of the passage to the limit since the arguments are quite standard.

Let ρ̅∈ L^2(0,T;W^1,∞(Ω)). We know that 𝒮𝒵_hρ̅→ρ̅ in L^2(0,T;W^1,∞(Ω)) as h→ 0 from (<ref>). Then  selecting ρ̅_h=𝒮𝒵_hρ̅ in (<ref>), multiplying by k, and summing over n yields

    ∫_0^T (∂_tρ_h,k,𝒮𝒵_hρ̅)_h + ∫_0^T h^γ (∇ρ_h,k^+,∇𝒮𝒵_hρ̅) 
       +∫_0^T (∇ℐ_h  A([ρ_k,h^+]_T),∇𝒮𝒵_hρ̅)-∫_0^T(ρ_h,k^+∇ℐ_h((K*[ρ_h,k^-]_T),∇𝒮𝒵_hρ̅) =0.



  * For the time derivative, we have:

    ∫_0^T (∂_tρ_h,k,𝒮𝒵_hρ)_h=∫_0^T [(∂_tρ_k,h,𝒮𝒵_hρ)_h-(∂_tρ_h,k,𝒮𝒵_hρ)]+∫_0^T(∂_tρ_h,k,𝒮𝒵_hρ).

It is clear from (<ref>)  that

    ∫_0^T(∂_tρ_h,k, 𝒮𝒵_hρ)→∫_0^T<∂_tρ,ρ̅>.

To control the residual term, we use (<ref>) combined with (<ref>) to see

    ∫_0^T    [(δ_tρ^n+1_h,𝒮𝒵_hρ)_h-(δ_tρ^n+1_h,𝒮𝒵_hρ)]
       ≤ C h^1/2(∫_0^T ∂_tρ_h,k^2_(H^1(Ω))')^1/2(∫_0^T∇𝒮𝒵_hρ_L^∞(Ω)^2)^1/2→ 0.


  * For the dissipation terms, we have by (<ref>) and (<ref>) that

    h^γ∫_0^T (∇ρ_h,k^+,∇𝒮𝒵_hρ̅)→0

and

    ∫_0^T (∇ℐ_hA([ρ_h,k^+]_T),∇𝒮𝒵_hρ̅)→∫_0^T (∇ℐ_hA(ρ_T),∇ρ̅)


  * For the convolution term, we proceed as follows.

    [      ∫_0^T(ρ_h,k^+∇ℐ_h( K*[ρ_h,k^-]_T),∇𝒮𝒵_hρ̅)                                               = ∫_0^T(ρ_h,k^+∇  (ℐ_h-ℐ)(K*[ρ_h,k^-]_T,)∇𝒮𝒵_hρ̅);                                                         +∫_0^T(ρ_h,k^+∇ K*[ρ_h,k^-]_T,∇𝒮𝒵_hρ̅). ]

For the first term, we have

    ∫_0^T   (ρ_h,k^+∇ (ℐ_h-ℐ)(K*[ρ_h,k^-]_T),∇𝒮𝒵_hρ̅) 
       ≤ C h  ∇^2 K_L^∞(^d)ρ^0_h_L^1(Ω)(∫_0^T ρ_h,k^+^2_L^2(Ω))^1/2(∫_0^T∇𝒮𝒵_hρ_L^2(Ω)^2)^1/2→ 0.

For the second term, we apply (<ref>) to show

    ∫_0^T∇ K*[ρ_h,k^-]_T-K*ρ_T^2_L^∞(Ω)≤∇ K^2_L^∞(^d)∫_0^T [ρ_h,k^-]_T-ρ_T^2_L^1(Ω)→ 0,

which implies on recalling (<ref>)  that

    ∫_0^T(ρ_h,k^+∇ K*[ρ_h,k^-]_T,∇𝒮𝒵_hρ̅)→∫_0^T(ρ∇ K*ρ_T,∇𝒮𝒵_hρ̅) .

Therefore,

    ∫_0^T(ρ_h,k^+∇ℐ_h( K*[ρ_h,k^-]_T),∇𝒮𝒵_hρ̅)→∫_0^T(ρ∇ K*ρ_T,∇𝒮𝒵_hρ̅) .

as (h,k)→ (0,0).


The continuous assimilation of the initial datum is ensured by the compact embedding H^1(0,T; (H^1(Ω))') into C([0,T]; (H^1(Ω))') and (<ref>). Moreover, one can show that ρ(t)→ρ^0 in L^p(Ω) as t→ 0. For more details, see <cit.>.

Since ρ̅∈ L^2(0,T, W^1,∞(Ω)) is dense in L^2(0,T; H^1(Ω)), we have found ρ: Q→ [0,∞) such that

    ρ∈ L^∞(0,T; L^2(Ω))∩ H^1(0,T; (H^1(Ω))'),


    ∫_Ωρ(t)=∫_Ωρ_0     t∈[0,T],

and

    {[ ∂_tρ-Δ  A(ρ_T) +∇· (ρ∇ K*ρ_T)                             =                             0                                        L^2(0,T; (H^1(Ω))'),;                          ρ(0)                             =                           ρ_0                                                  (H^1(Ω))'. ].

To complete with the proof of Theorem <ref> we show the equivalence of problems (<ref>) and (<ref>).
 Problems (<ref>) and (<ref>) are equivalent.

 At this point the only thing we need to show is that ρ defined by (<ref>) satisfies ρ≤ B_L^∞ in Q. Indeed, define  ρ_ aux= e^ tΔ K_L^∞(R^d)ρ_0_L^1(Ω)ρ_0_L^∞(Ω) for t∈[0,T], and observe that, by (<ref>),

    (∂_t ρ_ aux, ρ̅)-(ρ_ aux∇ K*ρ_T,∇ρ̅)≤0

holds for all ρ̅∈ H^1(Ω) with ρ̅≤0. Substracting (<ref>) from (<ref>), and testing the resulting equation against ρ=(ρ_ aux-ρ)_-∈ H^1(Ω) yields

    1/2d/dtρ^2_L^2(Ω)-(∇ A_T(ρ), ∇ρ)-(ρ∇ K*ρ_T, ∇ρ)≤ 0,

or equivalently,

    1/2d/dtρ^2_L^2(Ω)+(A^'_T)^1/2(ρ)∇ρ^2_L^2(Ω)-(ρ∇ K*ρ_T, ∇ρ)≤ 0.

It follows again from (<ref>) and integration by parts that

    d/dtρ^2_L^2(Ω)≤Δ K_L^∞(^d)ρ_0_L^1(Ω)ρ^2_L^2(Ω)

and so ρ≤ρ_ aux≤ B_L^∞ by Grnwall's lemma. Therefore, ρ_T=ρ.




§ SIMULATION OF AGGREGATION PHENOMENA


In this section we illustrate how scheme (<ref>) can be used to approximate the unique 
weak solution to (<ref>) with (<ref>)-(<ref>). Moreover, we compare our numerical solution 
to that computed in <cit.>.



 §.§ Computational performance
 At this point we shall make two comments regarding scheme  (<ref>). 
Firstly, we need not use the truncating operator [·]_T to compute A([ρ^n+1_h]_T) and K*[ρ_h^n]_T because  the discrete approximations are non-negative and the unique weak solution being approximated is not expected to blow up. Furthermore the convolution term K*ρ^n_h cannot be exactly computed at the nodes in order to construct its nodal interpolation, so a quadrature formula must be utilized on simplexes. 

Then, our numerical method remains as:  Given ρ^n_h∈ D_h, compute ρ^n+1_h∈ D_h satisfying

    (δ_tρ^n+1_h,ρ̅_h)_h+ h^γ (∇ρ^n+1_h,∇ρ̅_h)+(∇ℐ_h  A(ρ^n+1_h),∇ρ̅_h)-(ρ^n+1_h∇ℐ_h( 𝒬_1(K*ρ^n_h))),∇ρ̅_h)=0,

where 𝒬_1 is the midpoint quadrature formula. The term (δ_tρ^n+1_h,ρ̅_h)_h can be computed by using a closed-nodal quadrature formula, and the term ℐ_h(𝒬_1(K*ρ^n_h)) can be rewritten as follows. Let ∈𝒩_h, then

    Q_1(K*ρ^n_h)()=∫_Ω𝒫_h(K(-) ρ^n_h()) d=∑_E∈ℰ_h K(-_̱E) |E|,

where 𝒫_h is a piecewise constant interpolation taking its value on each E∈ℰ_h at the barycenter _̱E.

We see no obstacle to analyzing algorithm (<ref>) using (<ref>) and the truncating [·] in A(·) as well, but we did not consider such a formulation in our analysis because it is tedious.


Scheme (<ref>), and its modification (<ref>), require the solution of  nonlinear algebraic systems at each time step, which can be approximately solved using fixed-point iterations. In doing so, we first observe that A(ρ^n+1_h)=𝒟(ρ^n+1_h)∇ρ^n+1_h, where 𝒟(ρ^n+1_h) is a piecewise constant, d× d diagonal matrix function over the mesh 𝒯_h constructed as follows. Let E∈ℰ_h and consider Ẽ⊂ E to be a right simplex (see Figure <ref>) with vertices
{a_j}_j=0,⋯,d with a_0 supporting
the right angle. Then

    [𝒟(ρ̅_h)|_E]_jj=
    {[ A(ρ̅_h(a_j))-A(ρ̅_h(a_0))ρ̅_h(a_j)-ρ̅_h(a_0)                    if ρ̅_h(a_j)-ρ̅_h(a_0)≠0,;                                            0                    if ρ̅_h(a_j)-ρ̅_h(a_0)=0. ].

Particularly, we choose a_0  to be the incenter of E. Thus, we take a_i=a_0+ r_E/2_i, where r_E is the inradius of the inscribed ball.



We linearize as follows. For i=0, select
ρ^n+1_h,i=ρ^n_h, then compute ρ^n+1_h,i+1 using
in (<ref>) the expression

    (ρ^n+1_h,i+1,ρ̅_h)_h+ h^γ (∇ρ^n+1_h,i+1,∇ρ̅_h)+   (𝒟(ρ^n+1_h,i)∇ρ^n+1_h,i+1,∇ρ̅_h)
    
       -(ρ^n+1_h,i+1∇ℐ_h( 𝒬_1(K*ρ^n_h))) =(ρ^n_h,ρ̅_h)_h.

As a stopping criterion for the iterations, we choose  ρ^n+1_h,i+1-ρ^n+1_h,i_L^2(Ω)<tol, with tol being the prescribed tolerance.

Finally, the computation of (<ref>) for each ∈𝒩_h constitutes the bottleneck in running scheme (<ref>). To make it possible in an acceptable amount of time, a parallel procedure on a high-performance cluster can be invoked since all the nodes ∈𝒩_h are independent of each other.



 §.§ A numerical experiment
 As the domain we take the square Ω= [-4,4]^2⊂^2.  The evolution starts from the initial datum ρ^0=1/4 χ_[-3,3]^2 being a rescaled characteristic function supported in the square [-3,3]^2, which is shown in Figure (<ref>).  The local repulsion term is chosen as A(ρ)=ν/mρ^m with ν=0.1 and m=3, and the kernel is set as K(x)=exp(-|x|^2)/π.

From an N_ square× N_ square uniform grid, obtained by dividing Ω into macroelements consisting of squares, we construct the mesh 𝒯_h by splitting each macroelement  into 14 acute triangles as indicated in Figure <ref>.

This way, for N_ square=120, we define a mesh consisting of 201600 acute triangles and 101281 vertices with h=8/(2 N_𝐱)≃ 0.033. For the time discretization, we carry out 1500 iterations with time step k=10^-1. Moreover we select γ=0.99 to be as large as possible in orden to reduce the impact of the stabilizing term, which has a smoothing effect on the dynamics of aggregation phenomena. It should be noticed that (h,k) do not fulfill (<ref>) what makes us believe that such a restriction is superfluous. We also performed some numerical tests with k=10^-2, obtaining quite similar results, which are omitted for brevity.

Specifically, we run our test on a machine with 16 Intel Xeon E5 2670 processors (2,6 GHz, 8-core), in a distributed memory architecture; thus using a total of 256 parallel threads. The MPI library on the FreeFem++ PDE solver <cit.> was selected as a software framework.

Using the above-described parallel computing environment for the computation of (<ref>) for each ∈𝒩_h, the converged solution is obtained by about 3 iterations with tolerance tol=10^-3 in the L^2(Ω)-norm. To be more precise, the average number of iterations is 2.81, with minimum and maximum equal to 2 and 11, respectively. So, each time step takes an average time of 88.85 seconds, of which 35.95 seconds (on average) are due to the parallel computation of (<ref>). The remaining time is occupied in solving the associated linear system, which spends  3.06 seconds (on average) for each iteration, and data I/O. 

Figure <ref> shows how the initial state changes into four peaks that are aggregated into a single component until reaching a final steady state. This result is in good agreement with that in <cit.>. The dynamics regarding the ·_L^∞(Ω)- and ·_L^1(Ω)-norms is reported in Figure <ref>. The ·_L^∞(Ω)- norm approaches the value 16 as of t=14, while the ·_L^1(Ω)-norm takes values around 8.92, which is comparable to ρ^0_L^1(Ω)=9.
















































































99
Alt_Luckhaus_1983Alt, H. W. Luckhaus, S., Quasilinear elliptic-parabolic differential equations. Math. Z. 183 (1983), no. 3, 311–341.

Arzela_GG Azrad, P.; Guilln-Gonzlez, F., Mathematical justification of the hydrostatic approximation in the primitive equations of geophysical fluid dynamics. SIAM J. Math. Anal. 33 (2001), no. 4, 847–859.

Bertozzi_Slepcev_2010Bertozzi, A. L.; Slepčev, D.,
Existence and uniqueness of solutions to an aggregation equation with degenerate diffusion. Commun. Pure Appl. Anal. 9 (2010), no. 6, 1617–1637.

Boi_Cappaso_Morale_2000Boi, S.; Capasso, V.; Morale, D., Modeling the aggregative behavior of ants of the species Polyergus rufescens. Spatial heterogeneity in ecological models (Alcal de Henares, 1998). Nonlinear Anal. Real World Appl. 1 (2000), no. 1, 163–176.



Brenner_Scott_2008Brenner, S. C.; Scott, L. R., The mathematical theory of finite element methods. Third edition. Texts in Applied Mathematics, 15. Springer, New York, 2008.

Carrillo_Chertock_Huang_2015Carrillo, J. A.; Chertock, A.; Huang, Y., A finite-volume method for nonlinear nonlocal equations with a gradient flow structure. Commun. Comput. Phys. 17 (2015), no. 1, 233–258.







Ern_Guermond_2004Ern, A.; Guermond, J.-L., Theory and practice of finite elements. Applied Mathematical Sciences, 159. Springer-Verlag, New York, 2004.

Girault_Lions_2001Girault, V.; Lions, J.-L., Two-grid finite-element schemes for the transient Navier-Stokes problem. M2AN Math. Model. Numer. Anal. 35 (2001), no. 5, 945–980.

GG_GS_2009Guilln-Gonzlez, F.; Gutirrez-Santacreu, J. V., Stability and convergence of two discrete schemes for a degenerate solutal non-isothermal phase-field model. M2AN Math. Model. Numer. Anal. 43 (2009), no. 3, 563?589.

GG_GS_2008Guilln-Gonzlez, F.; Gutirrez-Santacreu, J. V.,Unconditional stability and convergence of fully discrete schemes for 2D viscous fluids models with mass diffusion. Math. Comp. 77 (2008), no. 263, 1495–1524.

Gurtin_McCamy_1977 Gurtin, M. E.; MacCamy, R. C., On the diffusion of biological populations. Math. Biosci. 33 (1977), no. 1-2, 35–49.

Grisvard_1985Grisvard, P., Elliptic problems in nonsmooth domains. Monographs and Studies in Mathematics, 24. Pitman (Advanced Publishing Program), Boston, MA, 1985.

Hecht_2012 Hecht, F., New development in freefem++. J. Numer. Math. 20 (2012), no. 3-4, 251–265.

Hillen_Painter_2009Hillen, T.; Painter, K. J., A user's guide to PDE models for chemotaxis. J. Math. Biol. 58 (2009), no. 1-2, 183–217.

Horstmann_1970_IHorstmann, D., From 1970 until present: the Keller-Segel model in chemotaxis and its consequences. I. Jahresber. Deutsch. Math.-Verein. 105 (2003), no. 3, 103–165.

Horstmann_1970_IIHorstmann, D., From 1970 until present: the Keller-Segel model in chemotaxis and its consequences. II. Jahresber. Deutsch. Math.-Verein. 106 (2004), no. 2, 51–69.



Keller_Segel_1971Keller, Evelyn F.; Segel, Lee A.,Model for chemotaxis. J. Theor. Biol., 30 (, 1971), 225–234.

Milewski_Yang_2008 Milewski, P. A.; Yang, X., A simple model for biological aggregation with asymmetric sensing. Commun. Math. Sci. 6 (2008), no. 2, 397–416.

Mogilner_Edelstein-Keshet_1999Mogilner, A.; Edelstein-Keshet, L., A non-local model for a swarm, J. Math. Biol. 38 (1999), no. 6, 534–570.

Mogilner_Edelstein-Keshet_Bent_Spiros_2003Mogilner,  A. I.; Edelstein-Keshet,  L..; Bent, L.; Spiros, A., Mutual interactions, potentials, and individual distance in a social aggregation. J. Math. Biol. 47 (2003), no. 4, 353–389.

Patlak_1953Patlak, C. S., Random walk with persistence and external bias. Bull. Math. Biophys. 15, (1953) 311–338.




Topaz_Bertozzi_Lewis_2004 Topaz, C. M.; Bertozzi, A. L.; Lewis, M. A., Swarming patterns in a two-dimensional kinematic model for biological groups. SIAM J. Appl. Math. 65 (2004), no. 1, 152–174.

Topaz_Bertozzi_Lewis_2006 Topaz, C. M.; Bertozzi, A. L.; Lewis, M. A., A nonlocal continuum model for biological aggregation. Bull. Math. Biol. 68 (2006), no. 7, 1601–1623.




