New Interior Penalty Discontinuous Galerkin Methods for ... - CMU Math

c

2007 001

.... ..... ..... .... .... .... .... ....

New Interior Penalty Discontinuous Galerkin Methods for the Keller-Segel Chemotaxis Model Yekaterina Epshteyn∗ and Alexander Kurganov† Abstract We develop a family of new interior penalty discontinuous Galerkin methods for the Keller-Segel chemotaxis model. This model is described by a system of two nonlinear PDEs: a convection-diffusion equation for the cell density coupled with a reaction-diffusion equation for the chemoattractant concentration. It has been recently shown that the convective part of this system is of a mixed hyperbolic-elliptic type, which may cause severe instabilities when the studied system is solved by straightforward numerical methods. Therefore, the first step in the derivation of our new methods is made by introducing the new variable for the gradient of the chemoattractant concentration and by reformulating the original Keller-Segel model in the form of a convection-diffusion-reaction system with a hyperbolic convective part. We then design interior penalty discontinuous Galerkin methods for the rewritten Keller-Segel system. Our methods employ the central-upwind numerical fluxes, originally developed in the context of finite-volume methods for hyperbolic systems of conservation laws. In this paper, we consider Cartesian grids and prove error estimates for the proposed high-order discontinuous Galerkin methods. Our proof is valid for pre-blow-up times since we assume boundedness of the exact solution. We also show that the blow-up time of the exact solution is bounded from above by the blow-up time of our numerical solution. In the numerical tests presented below, we demonstrate that the obtained numerical solutions have no negative values and are oscillation-free, even though no slope limiting technique has been implemented.

AMS subject classification: 65M60, 65M12, 65M15, 92C17, 35K57 Key words: Keller-Segel chemotaxis model, convection-diffusion-reaction systems, discontinuous Galerkin methods, NIPG, IIPG, and SIPG methods, Cartesian meshes.

1

Introduction

The goal of this work is to design new Discontinuous Galerkin (DG) methods for the twodimensional (2-D) Keller-Segel chemotaxis model, [13, 28, 29, 30, 35, 37]. The DG methods have recently become increasingly popular thanks to their attractive features such as: • local, element-wise mass conservation; ∗

Department of Mathematical Sciences, Carnegie Mellon University, Pittsburgh, PA, 15213, [email protected] † Mathematics Department, Tulane University, New Orleans, LA 70118; [email protected]

2

Y. Epshteyn and A. Kurganov

• flexibility to use high-order polynomial and non-polynomial basis functions; • ability to easily increase the order of approximation on each mesh element independently; • ability to achieve almost exponential convergence rate when smooth solutions are captured on appropriate meshes; • block diagonal mass matrices, which are of great computational advantage if an explicit time integration is used; • suitability for parallel computations due to (relatively) local data communications; • applicability to problems with discontinuous coefficients and/or solutions; • The DG methods have been successfully applied to a wide variety of problems ranging from the solid mechanics to the fluid mechanics (see, e.g., [3, 7, 14, 15, 17, 20, 22, 40] and references therein). In this paper, we consider the most common formulation of the Keller-Segel system [13], which can be written in the dimensionless form as  ρt + ∇ · (χρ∇c) = ∆ρ, (x, y) ∈ Ω, t > 0, (1.1) ct = ∆c − c + ρ, subject to the Neumann boundary conditions: ∇ρ · n = ∇c · n = 0,

(x, y) ∈ ∂Ω.

Here, ρ(x, y, t) is the cell density, c(x, y, t) is the chemoattractant concentration, χ is a chemotactic sensitivity constant, Ω is a bounded domain in R2 , ∂Ω is its boundary, and n is a unit normal vector. It is well-known that solutions of this system may blow up in finite time, see, e.g., [26, 27] and references therein. This blow-up represents a mathematical description of a cell concentration phenomenon that occurs in real biological systems, see, e.g., [1, 8, 10, 11, 16, 38]. Capturing blowing up solutions numerically is a challenging problem. A finite-volume, [21], and a finite-element, [34], methods have been proposed for a simpler version of the Keller-Segel model,  ρt + ∇ · (χρ∇c) = ∆ρ, ∆c − c + ρ = 0, in which the equation for concentration c has been replaced by an elliptic equation using an assumption that the chemoattractant concentration c changes over much smaller time scales than the density ρ. A fractional step numerical method for a fully time-dependent chemotaxis system from [41] has been proposed in [42]. However, the operator splitting approach may not be applicable when a convective part of the chemotaxis system is not hyperbolic, which is a generic situation for the original Keller-Segel model as it was shown in [12], where the finitevolume Godunov-type central-upwind scheme was derived for (1.1) and extended to some other chemotaxis and haptotaxis models. The starting point in the derivation of the central-upwind scheme in [12] was rewriting the original system (1.1) in an equivalent form, in which the concentration equation is replaced with the corresponding equation for the gradient of c:  ρt + ∇·(χρw) = ∆ρ, w ≡ (u, v) := ∇c. wt − ∇ρ = ∆w − w,

DG Methods for the Keller-Segel System

3

This form can be considered as a convection-diffusion-reaction system Ut + f(U)x + g(U)y = ∆U + r(U),

(1.2)

where U := (ρ, u, v)T , f(U) := (χρu, −ρ, 0)T , g(U) := (χρv, 0, −ρ)T , and r(U) := (0, −u, −v)T . The system (1.2) is an appropriate form of the chemotaxis system if one wants to solve it numerically by a finite-volume method. Even though the convective part of the system (1.2) is not hyperbolic, some stability of the resulting central-upwind scheme was ensured by proving its positivity preserving property, see [12]. A major disadvantage of the system (1.2) is a mixed type of its convective part. When a high-order numerical method is applied to (1.2), a switch from a hyperbolic region to an elliptic one may cause severe instabilities in the numerical solution since the propagation speeds in the elliptic region are infinite. Therefore, in order to develop high-order DG methods for (1.1), we rewrite it in a different form, which is suitable for DG settings: ρt + (χρu)x + (χρv)y = ∆ρ, ct = ∆c − c + ρ, u = cx , v = cy ,

(1.3) (1.4) (1.5) (1.6)

where the new unknowns ρ, c, u, v satisfy the following boundary conditions: ∇ρ · n = ∇c · n = (u, v)T · n = 0,

(x, y) ∈ ∂Ω.

(1.7)

The new system (1.3)–(1.6) may also be considered as a system of convection-diffusion-reaction equations kQt + F(Q)x + G(Q)y = k∆Q + R(Q), (1.8) where Q := (ρ, c, u, v)T , the fluxes are F(Q) := (χρu, 0, −c, 0)T and G(Q) := (χρv, 0, 0, −c)T , the reaction term is R(Q) := (0, ρ − c, −u, −v), the constant k = 1 in the first two equations in (1.8), and k = 0 in the third and the fourth equations there. As we show in §3, the convective part of the system (1.8) is hyperbolic. In this paper, we develop a family of high-order DG methods for the system (1.8). The proposed methods are based on three primal DG methods: the Nonsymmetric Interior Penalty Galerkin (NIPG), the Symmetric Interior Penalty Galerkin (SIPG), and the Incomplete Interior Penalty Galerkin (IIPG) methods, [4, 18, 19, 39]. The numerical fluxes in the proposed DG methods are the fluxes developed for the semidiscrete finite-volume central-upwind schemes in [32] (see also [31, 33] and references therein). These schemes belong to the family of non-oscillatory central schemes, which are highly accurate and efficient methods applicable to general multidimensional systems of conservation laws and related problems. Like other central fluxes, the central-upwind ones are obtained without using (approximate) Riemann problem solver, which is unavailable for the system under consideration. At the same time, a certain upwinding information—one-sided speeds of propagation—is incorporated into the central-upwind fluxes. We consider Cartesian grids and prove the error estimates for the proposed high-order DG methods under the assumption of boundedness of the exact solution. We also show that the blow-up time of the exact solution is bounded from above by the blow-up time of the solution of our DG methods. In numerical tests presented in §6, we demonstrate that the obtained numerical

4

Y. Epshteyn and A. Kurganov

solutions have no negative values and are oscillation-free, even though no slope limiting technique has been implemented. We also demonstrate a high order of numerical convergence, achieved even when the final computational time gets close to the blowup time and the spiky structure of the solution is well developed. The paper is organized as follows. In §2, we introduce our notations and assumptions, and state some standard results. The new DG methods are presented in §3. The consistency and error analysis of the proposed methods are established in Sections 4 and 5 (some proof details are postponed to Appendix A). Finally, in §6, we perform several numerical experiments.

2

Assumptions, Notations, and Standard Results

We denote by Eh a nondegenerate quasi-uniform rectangular subdivision of the domain Ω (the quasi-uniformity requirement will only be used in §5 for establishing the rate of convergence with respect to the polynomials degree). The maximum diameter over all mesh elements is denoted by h and the set of the interior edges is denoted by Γh . To each edge e in Γh , we associate a unit normal vector ne = (nx , ny ). We assume that ne is directed from the element E 1 to E 2 , where E 1 denotes a certain element and E 2 denotes an element that has a common edge with the element E 1 and a larger index (this simplified element notation will be used throughout the paper). For a boundary edge, ne is chosen so that it coincides with the outward normal. The discrete space of discontinuous piecewise polynomials of degree r is denoted by  Wr,h (Eh ) = w ∈ L2 (Ω) : ∀E ∈ Eh , w|E ∈ Pr (E) , where Pr (E) is a space of polynomials of degree r over the element E. For any function w ∈ Wr,h , we denote the jump and average operators over a given edge e by [w] and {w}, respectively: for an interior edge e = ∂E 1 ∩ ∂E 2 , for a boundary edge e = ∂E 1 ∩ ∂Ω, 1

1

2

[w] := weE − weE , 1

[w] := weE ,

1

2

{w} := 0.5weE + 0.5weE , 1

{w} := weE ,

2

where weE and weE are the corresponding polynomial approximations from the elements E 1 and E 2 . We also recall that the following identity between the jump and the average operators is satisfied: [w1 w2 ] = {w1 }[w2 ] + {w2 }[w1 ]. (2.1) For the finite-element subdivision Eh , we define the broken Sobolev space  H s (Eh ) = w ∈ L2 (Ω) : w|E j ∈ H s (E j ), j = 1, . . . , Nh

with the norms

|||w|||0,Ω =

X

E∈Eh

kwk20,E

! 21

and |||w|||s,Ω =

X

E∈Eh

kwk2s,E

! 21

, s > 0,

where k · ks,E denotes the Sobolev s-norm over the element E. We now recall some well-known facts that will be used in the error analysis in §5. First, let us state some approximations properties and inequalities for the finite-element space.

5

DG Methods for the Keller-Segel System

Lemma 2.1 (hp Approximation, [5, 6]) Let E ∈ Eh and ψ ∈ H s (E), s ≥ 0. Then there exist a positive constant C independent of ψ, r, and h, and a sequence ψerh ∈ Pr (E), r = 1, 2, . . . , such that for any q ∈ [0, s]



ψ − ψerh

q,E

≤C

hµ−q kψks,E , r s−q

µ := min(r + 1, s).

(2.2)

Lemma 2.2 (Trace Inequalities, [2]) Let E ∈ Eh . Then for the trace operators γ0 and γ1 , there exists a constant Ct independent of h such that   1 (2.3) ∀w ∈ H s (E), s ≥ 1, kγ0 wk0,e ≤ Ct h− 2 kwk0,E + hk∇wk0,E ,   1 ∀w ∈ H s (E), s ≥ 2, kγ1 wk0,e ≤ Ct h− 2 k∇wk0,E + hk∇2 wk0,E , (2.4)

where e is an edge of the element E.

Lemma 2.3 ([39]) Let E be a mesh element with an edge e. Then there is a constant Ct independent of h and r such that ∀w ∈ Pr (E),

1

kγ0wk0,e ≤ Ct h− 2 rkwk0,E .

(2.5)

Lemma 2.4 ([4, 9]) There exists a constant C independent of h and r such that ∀w ∈ Wr,h (Eh ),

kwk20,Ω ≤ C

X

E∈Eh

k∇wk20,E

X 1 + k[w]k20,e |e| e∈Γ h

! 21

,

where |e| denotes the measure of e. Lemma 2.5 (Inverse Inequalities) Let E ∈ Eh and w ∈ Pr (E). Then there exists a constant C independent of h and r such that kwkL∞ (E) ≤ Ch−1 r 2 kwk0,E , kwk1,E ≤ Ch−1 r 2 kwk0,E .

(2.6) (2.7)

We also recall the following form of Gronwall’s lemma: Lemma 2.6 (Gronwall) Let ϕ, ψ, and φ be continuous nonnegative functions R tdefined on the interval a ≤ t ≤ b, and the function φ is nondecreasing. If ϕ(t) + ψ(t) ≤ φ(t) + a ϕ(s) ds for all t ∈ [a, b], then ϕ(t) + ψ(t) ≤ et−a φ(t). In the analysis below we also make the following assumptions: • Ω is a rectangular domain with the boundary ∂Ω = ∂Ωver ∪ ∂Ωhor , where ∂Ωver and ∂Ωhor denote the vertical and horizontal pieces of the boundary ∂Ω, respectively. We also split the set hor if interior edges, Γh , into two sets of vertical, Γver h , and horizontal, Γh , edges, respectively; • The degree of basis polynomials is r ≥ 2 and the maximum diameter of the elements is h < 1 (the latter assumption is only needed for simplification of the error analysis).

6

3

Y. Epshteyn and A. Kurganov

Description of the Numerical Scheme

We consider the Keller-Segel  χu   0 ∂F  = ∂Q  0  0

system (1.8). First, notice that the Jacobians of F and G are    0 χρ 0 χv 0 0 χρ    0 0 0  0 0 0  ∂G   0   =  and ,  ∂Q  0 −1 0 0  0 0 0    0 0 0 0 −1 0 0

and their eigenvalues are

F F F G λF 1 = χu, λ2 = λ3 = λ4 = 0 and λ1 = χv,

G G λG 2 = λ3 = λ4 = 0,

(3.1)

respectively. Hence, the convective part of (1.8) is hyperbolic. We now design semidiscrete interior penalty Galerkin methods for this system. We assume that at any time level t ∈ [0, T ] the solution, (ρ, c, u, v)T is approximated by (discontinuous) piecewise polynomials of the corresponding degrees rρ , rc , ru , and rv , which satisfy the following relation: rmax ≤ a, rmax := max{rρ , rc , ru , rv }, rmin := min{rρ , rc , ru , rv }, (3.2) rmin where a is a constant independent of rρ , rc , rp , and rq . Our new DG methods are formulated as follows. Find a continuous in time solution (ρDG (·, t), cDG(·, t), uDG(·, t), v DG (·, t)) ∈ Wrρρ ,h × Wrcc ,h × Wruu ,h × Wrvv ,h , which satisfies the following weak formulation of the chemotaxis system (1.3)–(1.6): Z XZ XZ XZ DG ρ DG ρ DG ρ {∇w ρ · ne }[ρDG ] {∇ρ · ne }[w ] + ε ∇ρ ∇w − ρt w + Ω

+σρ

E∈Eh E X rρ2 Z

e∈Γh

− Z

XZ

DG

[ρ

c cDG t w

+σc

ρ

][w ] −

e

+

E∈Eh E X r2 Z c

uDG w u +

|e|

∇c

c

∇w −

[cDG ][w c ] +

e XZ

−

e∈∂Ωver e

c

DG

Z Ω

cDG (w u )x +

u

u

nx w + σu

ρ

(w )x +

e∈Γh e

{∇c

cDG w c − X Z

(χρDG uDG )∗ nx [w ρ]

(χρDG v DG )∗ ny [w ρ] = 0,

e

XZ

X

X Z

e∈Γver h e

DG

Z

c

· ne }[w ] + ε

XZ

e∈Γh e

ρDG w c = 0,

(3.3) {∇w c · ne }[cDG ] (3.4)

Ω

(−cDG )∗u nx [w u ]

e∈Γver h e

E∈Eh E

X Z

χρ

X Z

e∈Γhor h DG

DG DG

E∈Eh E

χρDG v DG (w ρ )y +

XZ

e∈Γh

Ω

XZ

E∈Eh E

Ω

Z

|e|

e∈Γh e

e∈Γh e

e∈Γh ∪∂Ωver

ru2 |e|

Z e

[uDG ][w u ] = 0,

(3.5)

7

DG Methods for the Keller-Segel System Z Ω

−

v

DG

v

w +

X Z

XZ

c

DG

v

DG

v

(w )y +

E∈Eh E

c

X Z

e∈Γhor e h 2 X rv

ny w + σv

e∈Γh ∪∂Ωhor

e∈∂Ωhor e

|e|

and the initial conditions: Z Z DG ρ ρ (·, 0)w = ρ(·, 0)w ρ, Ω Z

u

DG

u

(·, 0)w =

Ω

Ω Z

u

u(·, 0)w ,

Ω

(−cDG )∗v ny [w v ] Z

[v DG ][w v ] = 0,

(3.6)

e

Z

ΩZ

c

DG

v

c

(·, 0)w =

DG

Z

ΩZ

v

(·, 0)w =

Ω

c(·, 0)w c, v(·, 0)w v .

(3.7)

Ω

Here, (w ρ , w c , w u , w v ) ∈ Wrρρ ,h × Wrcc ,h × Wruu ,h × Wrvv ,h are the test functions, σρ , σc , σu and σv are real positive penalty parameters. The parameter ε is equal to either −1, 0, or 1: these values of ε correspond to the SIPG, IIPG, or NIPG method, respectively. To approximate the convective terms in (3.3) and (3.5)–(3.6), we use the central-upwind fluxes from [32]: 1

2

in DG DG E aout ain aout (χρDG uDG )E u )e e − a (χρ − [ρDG ], (χρ u ) = out in out in a −a a −a out DG DG E 1 in DG DG E 2 b (χρ v )e − b (χρ v )e bout bin (χρDG v DG )∗ = − [ρDG ], bout − bin bout − bin 1 in DG E 2 aout (cDG )E )e aout ain e − a (c DG ∗ (−c )u = − − [uDG ], aout − ain aout − ain 1 in DG E 2 bout bin bout (cDG )E )e e − b (c DG ∗ − out [v DG ]. (−c )v = − out in in b −b b −b DG DG ∗

(3.8)

Here, aout , ain , bout , and bin are the one-sided local speeds in the x- and y- directions. Since the convective part of the system (1.3)–(1.6) is hyperbolic, these speeds can be estimated using the ∂F largest and the smallest eigenvalues of the Jacobian ∂Q and ∂G (see (3.1)): ∂Q     1 DG E 2 in DG E 1 DG E 2 aout = max (χuDG )E , (χu ) , 0 , a = min (χu ) , (χu ) , 0 , e e e e     (3.9) 1 1 DG E 2 DG E 2 bout = max (χv DG )E )e , 0 , bin = min (χv DG )E )e , 0 . e , (χv e , (χv Remark. If aout − ain = 0 at a certain element edge e, we set 1

1

2

2

DG DG E DG DG E (χρDG v DG )E v )e (χρDG uDG )E u )e e + (χρ e + (χρ , (χρDG v DG )∗ = , (χρ u ) = 2 2 1 1 DG E 2 DG E 2 (cDG )E )e (cDG )E )e e + (c e + (c DG ∗ DG ∗ , (−c )v = − , (−c )u = − 2 2 DG DG ∗

there. Notice that in any case, the following inequalities, aout ≤ 1, aout − ain

−ain ≤ 1, aout − ain

bout ≤ 1, bout − bin

and

−bin ≤ 1, bout − bin

(3.10)

8

Y. Epshteyn and A. Kurganov

are satisfied. From now on we will assume that aout −ain > 0 and bout −bin > 0 throughout the computational domain.

4

Consistency of the Numerical Scheme

In this section, we show that the proposed DG methods (3.3)–(3.6) are strongly consistent with the Keller-Segel system (1.3)–(1.6). Lemma 4.1 If the solution of (1.3)–(1.6) is sufficiently regular, namely, if (ρ, c) ∈ H 1 ([0, T ]) ∩ H 2 (Eh ) and (u, v) ∈ L2 ([0, T ]) ∩ H 2 (Eh ), then it satisfies the formulation (3.3)–(3.6). Proof: We first multiply equation (1.3) by w ρ ∈ Wrρρ ,h and integrate by parts on one element E to obtain Z Z Z Z Z Z Z ρ ρ ρ ρ ρ ρ ρt w + ∇ρ∇w − ∇ρ · ne w − χρu(w )x + χρunx w − χρv(w )y + χρvny w ρ = 0. E

E

∂E

E

E

∂E

Notice that continuity of ρ and u implies that at the edge e, Therefore, [ρ] = 0 and

1 ρE e

=

∂E

2 ρE e

and

1 (χρu)E e

(4.1) 2 = (χρu)E e .

1 aout − ain 1 1 1 E2 E1 (χρu)E (χρu) (χρu)E + = (χρu) = e e e e out in 2 2 a −a 1 in in E2 out a aout (χρu)E a e − a (χρu)e E1 E2 (χρu)e − out (χρu)e = = (χρu)∗ . = out a − ain a − ain aout − ain

{χρu} =

Summing now equation (4.1) over all elements E ∈ Eh , using the jump-average identity (2.1), R P P r2 R adding the penalty terms ε e∈Γh e {∇w ρ · ne }[ρ] and σρ e∈Γh |e|ρ e [ρ][w ρ ], and using the Neumann boundary conditions (1.7), we obtain that the solution of the system (1.3)–(1.6) satisfies equation (3.3). A similar procedure can be applied to show that the solution of (1.3)–(1.6) satisfies equations (3.4)–(3.6) as well. This concludes the consistency proof. 

5

Error Analysis

In this section, we prove the existence and show the convergence of the numerical solution using the Schauder’s fixed point theorem, [24]. In the analysis below we will assume that the exact solution of the system (1.3)–(1.6) is sufficiently regular for t ≤ T , where T is a pre-blow-up time. In particular we will assume that (ρ, c, u, v) ∈ H s1 ([0, T ]) ∩ H s2 (Ω),

s1 > 3/2, s2 ≥ 3,

(5.1)

which is needed for the h-analysis (convergence rate with respect to the mesh size), or (ρ, c, u, v) ∈ H s1 ([0, T ]) ∩ H s2 (Ω),

s1 > 3/2, s2 ≥ 6,

(5.2)

9

DG Methods for the Keller-Segel System

which is needed for the r-analysis (convergence rate with respect to the polynomial degree). Notice that these assumptions are reasonable since classical solutions of the Keller-Segel system (1.1) are regular (before the blow-up time) provided the initial data are sufficiently smooth, see [26] and references therein. We denote by ρe, e c, u e, and ve the piecewise polynomial interpolants of the exact solution components ρ, c, u, and v of the Keller-Segel system (1.3)–(1.6) and assume that these interpolants satisfy the approximation property (2.2). We then use the idea similar to [36] and define the following subset of the broken Sobolev space:  S = (φρ , φc , φu , φv ) ∈ L2 ([0, T ]) ∩ L∞ ([0, T ]) ∩ Wrρρ ,h × Wrcc ,h × Wruu ,h × Wrvv ,h : sup kφρ −

t∈[0,T ]

≤ Cρ



h

2sρ −4

rρ

sup kφ −

≤ Cc



ZT   X rρ2 2 ρ + |||∇(φ − ρe)|||0,Ω + k[(φρ − ρe)]k20,e |e| e∈Γ

e ck20,Ω

+

h

rc2sc −4

+

h

ru2su −4

h2 min(rv +1,sv )−2 + rv2sv −4

h

0 2 min(rρ +1,sρ )−2

h

2 min(ru +1,su )−2

2 min(rc +1,sc )−2

ZT   X r2 2 2 c c c + |||∇(φ − e c)|||0,Ω + k[(φ − e c)]k0,e |e| e∈Γ

2sρ −4

rρ

sup kφu − u ek0,Ω ≤ Ch [0,T ]

h

0 2 min(rρ +1,sρ )−2

c

t∈[0,T ]

ρek20,Ω

2 min(rc +1,sc )−2

h

2 min(ru +1,su )−2

h

+ rc2sc −4 ru2su −4  1 1 1 1 + + + , rρ2 rc2 ru2 rv2

+ 

h2 min(rv +1,sv )−2 + rv2sv −4



,



,

ZT   X r2 u kφu − u ek20,Ω + k[(φu − u e)]k20,e |e| e∈Γh 0  2 min(rρ +1,sρ )−2  h h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 ≤ Cu + + , + 2s −4 rc2sc −4 ru2su −4 rv2sv −4 rρ ρ   1 1 1 1 v sup kφ − e v k0,Ω ≤ Ch 2 + 2 + 2 + 2 , rρ rc ru rv [0,T ]

ZT   X r2 2 2 v v v k[(φ − ve)]k0,e kφ − vek0,Ω + |e| e∈Γ h 0  2 min(rρ +1,sρ )−2  h h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 ≤ Cv + + , + 2s −4 rc2sc −4 ru2su −4 rv2sv −4 rρ ρ where C, Cρ, Cc , Cu , and Cv are positive constants (which will be defined later) independent of h and the polynomial degrees (rρ , rc , ru , rv ), and the parameters sρ , sc , su , and sv denote the regularity of the corresponding components of the exact solution. Clearly the subset S is a closed convex subset of the broken Sobolev space and it is not empty since it contains the element (e ρ, e c, u e, e v). We first show that the functions in S are bounded.

10

Y. Epshteyn and A. Kurganov

Lemma 5.1 For any (φρ , φc , φu , φv ) ∈ S, there exist positive constants Mρ , Mc , Mu , and Mv independent of h, rρ , rc , ru , and rv , such that sup kφρ k∞,Ω ≤ Mρ ,

sup kφc k∞,Ω ≤ Mc ,

t∈[0,T ]

sup kφu k∞,Ω ≤ Mu ,

t∈[0,T ]

sup kφv k∞,Ω ≤ Mv . (5.3)

t∈[0,T ]

t∈[0,T ]

Proof: From the definition of the subset S, we have  2 min(rρ +1,sρ )−2  h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 h 2 ρ + sup kφ − ρek0,Ω ≤ Cρ + + . 2s −4 rc2sc −4 ru2su −4 rv2sv −4 rρ ρ t∈[0,T ]

Hence,

sup kφρ − ρek0,Ω ≤ M

t∈[0,T ]

h

2 rmin

.

Using the inverse inequality (2.6), we obtain

sup kφρ − ρek∞,Ω ≤ M1 rρ2 h−1 sup kφρ − ρek0,Ω ≤

t∈[0,T ]

t∈[0,T ]

rmax ∗ M ≤ M. rmin

This estimate implies that

sup kφρ k∞,Ω ≤ M + sup ke ρk∞,Ω ,

t∈[0,T ]

[0,T ]

which, together with the hp approximation property (see Lemma 2.1), yields the first bound in (5.3). The remaining three estimates in (5.3) are obtained in a similar manner.  We now define the solution operator A on S as follows: ∀(φρ , φc , φu , φv ) ∈ S,

A(φρ , φc , φu , φv ) = (φρL, φcL , φuL , φvL),

c,0 u,0 v,0 where the initial conditions are (φρ,0 ρ0 , e c0 , u e0 , e v 0 ), and the functions L , φL , φL , φL ) = (e

φρL ∈ Wrρρ ,h,t := H s ([0, T ]) ∩ Wrρρ ,h ,

φcL ∈ Wrcc ,h,t := H s ([0, T ]) ∩ Wrcc ,h ,

s > 3/2

φuL ∈ Wruu ,h,t := L2 ([0, T ]) ∩ L∞ ([0, T ]) ∩ Wruu ,h , φvL ∈ Wrvv ,h,t := L2 ([0, T ]) ∩ L∞ ([0, T ]) ∩ Wrvv ,h

are such that Z XZ XZ XZ ρ ρ ρ ρ ρ ρ {∇w ρ · ne }[φρL ] {∇φL · ne }[w ] + ε (φL )t w + ∇(φL )∇w − E∈Eh E

Ω

e∈Γh e

e∈Γh e

X Z X rρ2 ρ u ρ ρ ρ (χφρL φu )∗ nx [w ρ ] χφL φ (w )x + + σρ [φL ][w ] − |e| e∈Γver e∈Γh E∈Eh E h e e Z Z X X (χφρL φv )∗ ny [w ρ ] = 0, ∀w ρ ∈ Wrρρ ,h , (5.4) χφρL φv (w ρ)y + − XZ

Z

E∈Eh E

Z Ω

(φcL )t w c

+

e∈Γhor e h

XZ

E∈Eh E

+ σc

∇φcL ∇w c

X r2 c |e| e∈Γ h

Z e

−

XZ

e∈Γh e

[φcL ][w c ] +

Z Ω

{∇φcL

c

· ne }[w ] + ε

φcL w c −

Z Ω

XZ

e∈Γh e

φρL w c = 0,

{∇w c · ne }[φcL ] ∀w c ∈ Wrcc ,h ,

(5.5)

11

DG Methods for the Keller-Segel System Z

φuLw u

+

XZ

φcL (w u )x

E∈Eh E

Ω

X

+ σu

e∈Γh ∪∂Ωver

Z

φvL w v

+

XZ

ru2 |e|

φcL (w v )y

E∈Eh E

Ω

X

+ σv

e∈Γh ∪∂Ωhor

rv2 |e|

X Z

+ Z

(−φcL )∗u nx [w u ]

e∈Γver h e

[φuL ][w u ] = 0,

−

X Z

φcL nx w u

e∈∂Ωver e

∀w u ∈ Wruu ,h ,

(5.6)

e

X Z

+ Z

(−φcL )∗v ny [w v ]

e∈Γhor e h

[φvL ][w v ] = 0,

−

X Z

φcL ny w v

e∈∂Ωhor e

∀w v ∈ Wrvv ,h .

(5.7)

e

As before, the central-upwind numerical fluxes are utilized in (5.4)–(5.7): ρ u E ρ u E in in aout aout L aL L (χφL φ )e − aL (χφL φ )e = − [φρ ], out in out in L aL − aL aL − aL 1

(χφρL φu )∗

2

ρ v E ρ v E in in bout bout L bL L (χφL φ )e − bL (χφL φ )e − [φρ ], in out in L bout − b b − b L L L L 1

(χφρL φv )∗ =

1

2

2

in c E in aout (φcL)E aout e − aL (φL )e L aL (−φcL )∗u = − L − [φu ], in out in L aout − a a − a L L L L 1

(−φcL )∗v

(5.8)

2

in c E in bout (φcL )E bout e − bL (φL )e L bL =− L − [φv ], in out in L bout − b b − b L L L L

where the one-sided local speeds are:     u E1 u E2 in u E1 u E2 aout := max (χφ ) , (χφ ) , 0 , a := min (χφ ) , (χφ ) , 0 , e e e e L L     v E1 v E2 in v E1 v E2 bout L := max (χφ )e , (χφ )e , 0 , bL := min (χφ )e , (χφ )e , 0 .

(5.9)

Notice that the inequalities similar to (3.10), aout L ≤ 1, in aout L − aL

−ain L ≤ 1, aout − ain L L

bout L ≤ 1, in bout L − bL

and

−bin L ≤ 1, bout − bin L L

(5.10)

which are needed in our convergence proof, are satisfied for the local speeds defined in (5.9) as well (for simplicity, we assume that aout −ain 6= 0 and bout −bin 6= 0 throughout the computational domain). We now show that the operator A is well-defined by proving existence and uniqueness of (φρL , φcL, φuL , φvL ). Lemma 5.2 There exists a unique solution (φρL , φcL , φuL, φvL ) ∈ Wrρρ ,h,t × Wrcc ,h,t × Wruu ,h,t × Wrvv ,h,t of (5.4)–(5.7). Proof: First, notice that equations (5.4)–(5.5) can be rewritten as the explicit linear differential equations for φρL and φcL . Hence, there exists a unique solution (φρL , φcL ) ∈ Wrρρ ,h,t × Wrcc ,h,t .

12

Y. Epshteyn and A. Kurganov

Equations (5.6)–(5.7) can be rewritten as Z Z X Z XZ X ru2 c u u u u u (−φcL )∗u nx [w u ] φL (w )x − [φL ][w ] = − φL w + σu |e| e∈Γver E∈Eh E e∈Γh ∪∂Ωver h e e Ω Z X ∀w u ∈ Wruu ,h (5.11) φcL nx w u , + e∈∂Ωver e

Z Ω

φvL w v

+ σv

X

e∈Γh ∪∂Ωhor

rv2 |e|

Z

[φvL ][w v ]

=−

e

+

XZ

φcL (w v )y

E∈Eh

E X Z

φcL ny w v ,

e∈∂Ωhor e

−

X Z

(−φcL )∗v ny [w v ]

e∈Γhor e h

∀w v ∈ Wrvv ,h .

(5.12)

The bilinear form on the left-hand side (LHS) of equation (5.11) is coercive since for all ϕ ∈ Wruu ,h , Z Z X ru2 [ϕ][ϕ] ≥ kϕk20,Ω . ϕϕ + σu |e| e∈Γh ∪∂Ωver

Ω

e

It is also continuous on Wruu ,h ×Wruu ,h , while the linear form on the right-hand side (RHS) of (5.11) is continuous on Wruu ,h . Hence, there exists a unique solution of (5.11). The same argument is true for equation (5.12). This concludes the proof of the lemma.  Our next goal is to show that the operator A maps S into itself and that A is compact. By the second Shauder fixed-point theorem, [24], this will imply that the nonlinear mapping (φρ , φc , φu , φv ) ∈ S → A(φρ , φc , φu , φv ) has a fixed point denoted by (ρDG , cDG , uDG , v DG ). Theorem 5.3 Let the solution of (1.3)–(1.6) satisfy the assumption (5.1). Then for any (φρ , φc , φu , φv ) ∈ S, A(φρ , φc , φu , φv ) ∈ S.

Proof: Let (φρ , φc , φu , φv ) ∈ S and (φρL , φcL, φuL , φvL) = A(φρ , φc , φu , φv ). We introduce the following notation: τ ρ := φρL − ρe, ξ ρ := ρ − ρe, τ c := φcL − e c, ξ c := c − e c, (5.13) u u u v v τ := φL − u e, ξ = u − u e, τ := φL − ve, ξ v := v − ve. It follows from the consistency Lemma 4.1 that the exact solution of (1.3)–(1.6) satisfies not only equation (3.3) but also the similar equation Z X rρ2 Z XZ XZ XZ ρ ρ ρ ρ {∇w · ne }[ρ] + σρ {∇ρ · ne }[w ] + ε ∇ρ∇w − ρt w + [ρ][w ρ ] |e| e∈Γh e∈Γh e e∈Γh e E∈Eh E e Ω X Z XZ X Z XZ (χρv)∗∗ ny [w ρ ] = 0,(5.14) χρv(w ρ )y + (χρu)∗∗ nx [w ρ] − χρu(w ρ)x + − E∈Eh E

e∈Γver h e

E∈Eh E

e∈Γhor e h

where 1

(χρu)

∗∗

1

(χρv)

∗∗

2

in in E aout aout (χρu)E L aL e − aL (χρu)e − [ρ], := L in in aout aout L − aL L − aL 2

in in E bout bout (χρv)E e − bL (χρv)e L bL − [ρ], := L in in bout bout L − bL L − bL

13

DG Methods for the Keller-Segel System

in out in and the local speeds aout L , aL , bL , and bL are given by (5.9). Using (5.13), equation (5.14) can be rewritten as: Z X rρ2 Z XZ XZ XZ ρ ρ ρ ρ {∇w · ne }[e ρ] + σρ {∇e ρ · ne }[w ] + ε [e ρ][w ρ ] ∇e ρ∇w − ρet w + |e| e∈Γh e∈Γh e e∈Γh e E∈Eh E e Ω Z Z Z Z X X X X (χρv)∗∗ ny [w ρ ] χe ρφv (w ρ )y + (χρu)∗∗ nx [w ρ] − χe ρφu (w ρ)x + − E∈Eh E

=−

Z

ξtρ w ρ −

Ω

− σρ +

X

e∈Γh

XZ

rρ2 |e|

XZ

E∈Eh E

Z

e∈Γver h e

∇ξ ρ ∇w ρ +

[ξ ρ ][w ρ ] +

E∈Eh E

e

χξ ρ v(w ρ )y −

E∈Eh E

XZ

XZ

XZ

e∈Γh e

E∈Eh E

{∇ξ ρ · ne }[w ρ ] − ε

χξ ρ u(w ρ)x −

XZ

E∈Eh E

XZ

e∈Γh e

e∈Γhor e h

{∇w ρ · ne }[ξ ρ]

χe ρ(φu − u)(w ρ)x

χe ρ(φv − v)(w ρ )y .

E∈Eh E

(5.15)

Subtracting equation (5.15) from (5.4) and choosing w ρ = τ ρ , we obtain  X rρ2 1d  ρ kτ k0,Ω + |||∇τ ρ |||20,Ω + σρ k[τ ρ ]k20,e 2 dt |e| e∈Γh X Z XZ XZ ρ u ρ ρ ρ ((χφρLφu )∗ − (χρu)∗∗ ) nx [τ ρ ] χτ φ (τ )x − {∇τ · ne }[τ ] + = (1 − ε) XZ

+

E∈Eh E

+

Z

χτ ρ φv (τ ρ )y −

ξtρτ ρ

+

X

rρ2

e∈Γh

XZ

+

XZ

E∈Eh E

Ω

+ σρ

E∈Eh E

e∈Γver h e

E∈Eh E

e∈Γh e

|e|

Z e

X Z

e∈Γhor e h

ρ

ρ

∇ξ ∇τ −

[ξ ρ ][τ ρ ] −

XZ

E∈Eh E

χe ρ(φu − u)(τ ρ )x +

((χφρL φv )∗ − (χρv)∗∗ ) ny [τ ρ ] XZ

e∈Γh e

ρ

{∇ξ · ne }[τ ] + ε

χξ ρu(τ ρ )x −

XZ

E∈Eh E

ρ

XZ

XZ

e∈Γh e

{∇τ ρ · ne }[ξ ρ]

χξ ρ v(τ ρ )y

E∈Eh E

ρ . χe ρ(φv − v)(τ ρ )y =: T1ρ + T2ρ + ... + T14

(5.16)

Next, we bound each term on the RHS of (5.16) using standard DG techniques. The quantities εi in the estimates below are positive real numbers, which will be defined later. We begin with the first term on the RHS of (5.16). The Cauchy-Schwarz inequality yields: X k{∇τ ρ }k0,e k[τ ρ ]k0,e . |T1ρ| ≤ (1 − ε) e∈Γh

As before, we denote by E 1 and E 2 the two elements sharing the edge e. Then, using the inequality (2.5), we obtain

 X 1  X

ρ ρ ρ E1 ρ E2 k{∇τ }k0,e k[τ ]k0,e ≤ k[τ ρ ]k0,e

(∇τ )e + (∇τ )e 2 0,e 0,e e∈Γ e∈Γ h

h

14

Y. Epshteyn and A. Kurganov  C t rρ X  ≤ √ k∇τ ρ k0,E 1 + k∇τ ρ k0,E 2 k[τ ρ ]k0,e , 2 h e∈Γh

and hence, using the fact that |e| ≤ |T1ρ|

≤

ερ1

X

E∈Eh

k∇τ ρ k20,E

+

C1ρ

√

h, we end up with the following bound on T1ρ :

X rρ2 X rρ2 ρ ρ ρ 2 ρ 2 k[τ ]k0,e = ε1 |||∇τ |||0,Ω + C1 k[τ ρ ]k20,e . |e| |e| e∈Γ e∈Γ h

(5.17)

h

Consider now the second term on the RHS of (5.16). From Lemma 5.1 we know that φu is a bounded function, hence T2ρ can be bounded as follows: X (5.18) k(τ ρ )x k20,E + C2ρ kτ ρ k20,Ω ≤ ερ2 |||(τ ρ )x |||20,Ω + C2ρ kτ ρ k20,Ω . |T2ρ | ≤ ερ2 E∈Eh

Next, we bound the third term on the RHS of (5.16) as   out X  Z a 1 1 ρ ρ L u E E ρ |T3 | ≤ (χφ φ ) − (χρu) n [τ ] x e e L aout − ain ver L L e∈Γh e Z   in −a 2 2 ρ L u E E ρ + (χφ φ ) − (χρu) n [τ ] x e e L in aout − a L L e  Z in out −a a ρ L L ρ [φ − ρ]n [τ ] + x in L =: I + II + III. aout L − aL

(5.19)

e

Using (5.10) and (5.13), the first term on the RHS of (5.19) can be estimated by  X Z  ρ u E1 E1 ρ (φ φ ) − (ρu) n [τ ] I≤χ x e e L e∈Γver h

e Z X  Z ρ u E1 ρ ρ u E1 ρ (τ φ ) n [τ ] + (ξ φ ) n [τ ] ≤χ x x e e e∈Γver h

Ze e Z  u E1 ρ u E1 ρ e)ρ)e nx [τ ] + (ξ ρ)e nx [τ ] =: eI. + ((φ − u e

e

We now use the Cauchy-Schwarz inequality, the trace inequality (2.3), the inequality (2.5), the assumption (3.2), the approximation inequality (2.2), and the bound on φu from Lemma 5.1 to obtain the bound on eI:   2 min(rρ +1,sρ ) 2 min(ru +1,su ) X rρ2 1 h h 2 2 ρ ρ ∗ eI ≤ kτ k + K + C ∗∗ kφu − u ek20,Ω . k[τ ]k0,e + C + 0,Ω 2sρ 2su 2 |e| r rρ u e∈Γh

A similar bound can be derived for the second term II on the RHS of (5.19). To estimate the last term on the RHS of (5.19), we first use (5.13) and the definition of the one-sided local speeds (5.9) to obtain Z  X  ρ 2 ρ ρ f k[τ ]k0,e + [ξ ][τ ] := III. III ≤ C e∈Γver h

e

DG Methods for the Keller-Segel System

15

Then, using the Cauchy-Schwarz inequality, the trace inequality (2.3), and the approximation f as follows: inequality (2.2), we bound III X 2  rρ h2 min(rρ +1,sρ ) K h 1 ρ 2 f≤ + K k[τ ]k + C . III 2 0,e 2s rρ2 |e| rρ ρ e∈Γ h

Combining the above bounds on I, II, and III, we arrive at   2 min(rρ +1,sρ ) X rρ2 h h2 min(ru +1,su ) ρ ρ ∗ ρ 2 ρ 2 + C ∗∗ kφu − u ek20,Ω . k[τ ]k0,e + C |T3 | ≤ kτ k0,Ω + C3 + 2sρ 2su |e| r rρ u e∈Γh (5.20) The terms T4ρ and T5ρ are bounded in the same way as the terms T2ρ and T3ρ , respectively, and the bounds are: |T4ρ | ≤ ερ2 |||(τ ρ )y |||20,Ω + C4ρ kτ ρ k20,Ω . (5.21) and |T5ρ |

≤

kτ ρ k20,Ω

+

C5ρ

  2 min(rρ +1,sρ ) X rρ2 h h2 min(rv +1,sv ) ∗ ρ 2 + C ∗∗ kφv − e v k20,Ω . k[τ ]k0,e + C + 2sρ 2sv |e| r rρ v e∈Γ h

(5.22) The term T6ρ is bounded using the Cauchy-Schwarz inequality and the approximation inequality (2.2): h2 min(rρ +1,sρ ) |T6ρ | ≤ kτ ρ k20,Ω + C ∗ . (5.23) 2s rρ ρ Using the Cauchy-Schwarz inequality, Young’s inequality, and the approximation inequality (2.2) for ρ, we obtain the following bound for the term T7ρ : |T7ρ |

≤

ερ7 |||∇τ ρ |||20,Ω

+

2 min(rρ +1,sρ )−2 ∗h C . 2s −2 rρ ρ

(5.24)

The term T8ρ is bounded using the Cauchy-Schwarz inequality, the trace inequality (2.4), and the approximation inequality (2.2): |T8ρ | ≤ C8ρ

X rρ2 h2 min(rρ +1,sρ )−2 k[τ ρ ]k20,e + C ∗ . 2sρ −2 |e| r ρ e∈Γ

(5.25)

h

To bound the term T9ρ we use the trace inequality (2.5), inequality (2.3), the Cauchy-Schwarz inequality and Young’s inequality: |T9ρ | ≤ ερ9 |||∇τ ρ |||20,Ω + C ∗

h2 min(rρ +1,sρ )−2 2sρ −4

rρ

.

(5.26)

ρ Similarly, we bound the term T10 by: ρ ρ |T10 | ≤ C10

X rρ2 h2 min(rρ +1,sρ )−2 . k[τ ρ ]k20,e + C ∗ 2sρ −4 |e| r ρ e∈Γ h

(5.27)

16

Y. Epshteyn and A. Kurganov

ρ ρ For the terms T11 and T12 , we use our assumption on the smoothness of the exact solution together with the Cauchy-Schwarz inequality and the approximation inequality (2.2) to obtain the following bounds:

ρ | |T11

≤

ερ11 |||(τ ρ )x |||20,Ω

+

2 min(rρ +1,sρ ) ∗h , C 2s rρ ρ

ρ | |T12

≤

ερ11 |||(τ ρ )y |||20,Ω

+

2 min(rρ +1,sρ ) ∗h C . 2s rρ ρ

(5.28)

ρ . We first use (5.13) to obtain Consider now the term T13 ρ |T13 |

 Z X  Z u ρ u ρ (φ − u e)(τ )x + ξ (τ )x . ≤C E∈Eh

E

E

Then we apply the Cauchy-Schwarz inequality and the approximation inequality (2.2), which result in h2 min(ru +1,su ) ρ |T13 | ≤ ερ13 |||(τ ρ )x |||20,Ω + C ∗ + C ∗∗ kφu − u ek20,Ω . (5.29) ru2su ρ ρ : is obtained in the same way as the bound on T13 The bound on the term T14

ρ |T14 |

≤

ερ13 |||(τ ρ )y |||20,Ω

+

2 min(rv +1,sv ) ∗h C rv2sv

+ C ∗∗ kφv − vek20,Ω .

(5.30)

Finally, we plug the estimates (5.17)–(5.18) and (5.20)–(5.30) into (5.16) and use the assumption that h < 1 to obtain 1d ρ 2 kτ k0,Ω + (1 − ερ1 − ερ2 − ερ7 − ερ9 − ερ11 − ερ13 )|||∇τ ρ |||20,Ω 2 dt X rρ2 ρ k[τ ρ ]k20,Ω ) + (σρ − C1ρ − C3ρ − C5ρ − C8ρ − C10 |e| e∈Γh  2 min(rρ +1,sρ )−2  h h2 min(ru +1,su ) h2 min(rv +1,sv ) ρ ρ ρ 2 ∗ + ≤ (3 + C2 + C4 )kτ k0,Ω + Cρ + 2s −4 ru2su rv2sv rρ ρ +C ∗∗ (kφu − u ek20,Ω + kφv − vek20,Ω ).

(5.31)

We now choose ερi and the penalty parameter σρ so that the coefficients of the |||∇τ ρ |||20,Ω and P rρ2 ρ 2 e∈Γh |e| k[τ ]k0,Ω on the LHS of (5.31) are equal to 1/2. We then multiply equation (5.31) by 2 and integrate it in time from 0 to t. Taking into account that (φu , φv ) ∈ S and using the fact that τ 0 = 0, we obtain:  Zt Zt  2 X r ρ e ρ kτ ρ k2 k[τ ρ k20,e ≤ C kτ ρ k20,Ω + |||∇τ ρ |||20,Ω + 0,Ω |e| e∈Γh 0 0  2 min(rρ +1,sρ )−2  h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 uv h +C + + + . (5.32) 2s −4 rc2sc −4 ru2su −4 rv2sv −4 rρ ρ

17

DG Methods for the Keller-Segel System

Next, we apply Gronwall’s Lemma 2.6 and take the supremum with respect to t of the both sides of (5.32): sup kτ ρ k20,Ω [0,T ] ≤C

+ I

ZT 



0

|||∇τ

ρ

|||20,Ω

2 min(rρ +1,sρ )−2

h

2sρ −4

rρ

 X rρ2 ρ 2 k[τ ]k0,e + |e| e∈Γ h

2 min(rc +1,sc )−2

+

h

rc2sc −4

h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 + + ru2su −4 rv2sv −4



, (5.33)

where C I is a constant that depends on kρk(L∞ ([0,T ]);H 2 (Ω)) , kρt k(L∞ ([0,T ]);L2 (Ω)) , kuk(L∞ ([0,T ]);L2 (Ω)) , kvk(L∞ ([0,T ]);L2 (Ω)) , and T only. According to the definition on page 9, the estimate (5.33) implies that φρL ∈ S. Using similar techniques, it can be shown that (φcL , φuL, φvL ) ∈ S as well (see Appendix A for the detailed proof). Therefore, we have shown that A(S) ⊂ S, and the proof of Theorem 5.3 is now complete.  Let us recall that our goal is to show that the operator A has a fixed point. Equipped with Theorem 5.3, it remained to prove that A is compact. To this end, we need to show that A is continuous and equicontinuous. Lemma 5.4 The operator A is continuous and equicontinuous. Proof: We consider the sequence {(φρn , φcn , φun , φvn )} and assume that sup (k(φρn , φcn , φun , φvn ) − (φρ , φc , φu , φv )kS ) → 0 as n → ∞.

t∈[0,T ]

Let (φρL,n , φcL,n , φuL,n , φvL,n ) = A(φρn , φcn , φun , φvn )

(5.34)

and (φρL , φcL , φuL, φvL ) = A(φρ , φc , φu , φv ) (5.35) be two solutions of (5.4)–(5.7). We denote by (φbρL, φbcL , φbuL , φbvL) the difference between these two bc,0 bu,0 bv,0 solutions (note that (φbρ,0 L , φL , φL , φL ) = (0, 0, 0, 0)), subtract (5.35) from (5.34), and choose the test function in the resulting equation for ρ to be w ρ = φbρ . This yields: L

X rρ2 1 d bρ 2 kφL k0,Ω + |||∇φbρL|||20,Ω + σρ k[φbρL ]k20,e 2 dt |e| e∈Γh XZ XZ XZ ρ ρ ρ u bρ b b b {∇φL · ne }[φL] + χφL φ (φL )x + χφρL,n (φun − φu )(φbρL )x = (1 − ε) E∈Eh E

e∈Γh e

−

+

X Z

e∈Γver h e

XZ

E∈Eh E

+

X Z

e∈Γhor h

e

(χφbρL φu )∗ nx [φbρL ] +

χφbρL φv (φbρL)y +

X Z

e∈Γver h e

XZ

E∈Eh E

E∈Eh E


(χφρL,n φu )∗ − (χφρL,n φun )∗ nx [φbρL ]

χφρL,n (φvn

− φ )(φbρL )y − v

X Z

e∈Γhor e h

(χφbρL φv )∗ ny [φbρL ]


(χφρL,n φv )∗ − (χφρL,n φvn )∗ ) ny [φbρL ] =: R1 + R2 + ... + R9 .

(5.36)

18

Y. Epshteyn and A. Kurganov

We now bound each term on the RHS of (5.36). The term R1 can be bounded using the Cauchy-Schwarz inequality, Young’s inequality, and the inequality (2.5): X rρ2 1 k[φbρL ]k20,e . (5.37) |R1 | ≤ |||∇φbρL|||20,Ω + C1 6 |e| e∈Γ h

Next, applying the Cauchy-Schwarz and Young’s inequalities and using the boundedness of ||φu ||∞,Ω , established in Lemma 5.1, we obtain the following bound on R2 : 1 |R2 | ≤ |||(φbρL)x |||20,Ω + C2 kφbρLk20,Ω . 6

(5.38)

Using the Cauchy-Schwarz and Young’s inequalities and the fact that φρL,n ∈ S, we bound the term R3 by 1 (5.39) |R3 | ≤ |||(φbρL)x |||20,Ω + C3 kφun − φu k20,Ω 6 We then use the Cauchy-Schwarz inequality, the inequality (2.5), and the first numerical flux formula in (5.8) to estimate R4 : |R4 | ≤ kφbρL k20,Ω + C4

X rρ2 k[φbρL ]k20,e . |e| e∈Γ

(5.40)

h

We now consider the term R5 . It follows from formulae (5.8)–(5.9) that the numerical fluxes 1 (χφρL,n φu )∗ is the composition of the continuous functions with respect to the variables (φu )E e 2 and (φu )E e . Hence, we can apply the Cauchy-Schwarz inequality and the inequality (2.5) to R5 so that it is bounded by |R5 | ≤

k(χφρL,n φu )∗

−

2 (χφρL,n φun )∗ k0,Ω

X rρ2 + C5 k[φbρL ]k20,e . |e| e∈Γ

(5.41)

h

The terms R6 , R7 , R8 , and R9 are similar to the terms R2 , R3 , R4 , and R5 estimated in (5.38), (5.39), (5.40), and (5.41), respectively. Therefore, we obtain 1 |R6 | ≤ |||(φbρL)y |||20,Ω + C6 kφbρL k20,Ω , 6 1 bρ |R7 | ≤ |||(φL)y |||20,Ω + C7 kφvn − φv k20,Ω , 6 X rρ2 ρ 2 b k[φbρL]k20,e , |R8 | ≤ kφL k0,Ω + C8 |e| e∈Γ

(5.42) (5.43) (5.44)

h

|R9 | ≤ k(χφρL,n φv )∗ − (χφρL,n φvn )∗ k20,Ω + C9

Substituting the estimates (5.37)–(5.45) into (5.36) yields:

X rρ2 k[φbρL ]k20,e . |e| e∈Γ

(5.45)

h

 X rρ2 1 1 d bρ 2 ρ 2 b kφL k0,Ω + |||∇φL|||0,Ω + (σρ − C) ||[φbρL]||20,e ≤ C ∗ kφbρL k20,Ω + C ∗∗ kφun − φu k20,Ω 2 dt 2 |e| e∈Γh  ρ ρ ρ ρ v v 2 u ∗ u ∗ 2 v ∗ v ∗ 2 + kφn − φ k0,Ω + k(χφL,n φ ) − (χφL,n φn ) k0,Ω + k(χφL,n φ ) − (χφL,n φn ) k0,Ω ,

19

DG Methods for the Keller-Segel System

where the penalty parameter σρ is chosen sufficiently large so that the coefficient (σρ − C) is nonnegative. We now integrate the latter inequality with respect to time from 0 to t and apply Gronwall’s Lemma 2.6 to obtain kφbρLk20,Ω

  Zt  Zt  X rρ2 ρ,0 ρ ρ ||[φbL ]||20,e dt ≤ M kφL k20,Ω + kφun − φu k20,Ω + |||∇φbL|||20,Ω + (σρ − C) |e| e∈Γh 0 0  ρ ρ ρ ρ v v 2 u ∗ u ∗ 2 v ∗ v ∗ 2 + kφn − φ k0,Ω + k(χφL,n φ ) − (χφL,n φn ) k0,Ω + k(χφL,n φ ) − (χφL,n φn ) k0,Ω .

Finally, taking the supremum over t and since φbρ,0 L = 0, we arrive at sup t∈[0,T ]

kφbρLk20,Ω

 ZT  ZT  X rρ2 ρ ρ 2 ∗ 2 kφun − φu k20,Ω + kφvn − φv k20,Ω ||[φbL ]||0,e dt ≤ M + |||∇φbL|||0,Ω + |e| e∈Γh 0 0  ρ ρ ρ ρ u ∗ u ∗ 2 v ∗ v ∗ 2 + k(χφL,n φ ) − (χφL,n φn ) k0,Ω + k(χφL,n φ ) − (χφL,n φn ) k0,Ω .

This inequality together with the similar inequalities for φbc , φbu , and φbv , which can be obtained in an analogous way, imply continuity of the operator A. ρ c u v Applying similar techniques to the difference (φL , φL , φL , φL ) := (φρL, φcL , φuL, φvL )(t1 , x1 , y1 ) − (φρL , φcL, φuL , φvL )(t2 , x2 , y2) and using the fact that (φu , φv ) ∈ S, one can show that the operator A is equicontinuous.  Equipped with Lemma 5.4, we conclude that the operator A is compact. Hence, by the second Schauder fixed-point theorem, [24], it has at least one fixed point (ρDG , cDG , uDG, v DG ), which is the DG solution of (3.3)–(3.6). For this solution, we establish the convergence rate results, stated in the following theorem. Theorem 5.5 (L2 (H 1 )- and L∞ (L2 )-Error Estimates) Let the solution of the Keller-Segel system (1.3)–(1.6) satisfies the smoothness assumption (5.2). If the penalty parameters σρ , σc , σu , and σv in the DG method (3.3)–(3.9) are sufficiently large and rmin ≥ 2, then there exist constants Cρ and Cc , independent of h, rρ , rc , ru , and rv such that the following two error estimates hold: DG

kρ kc

DG

DG

− ρkL∞ ([0,T ];L2 (Ω)) + |||∇(ρ − ckL∞ ([0,T ];L2(Ω)) + |||∇(c

where E :=



hmin(rρ +1,sρ )−1 s −2

rρρ

DG

 21  ZT X 2 rρ 2 DG k[ρ − ρ]k0,e − ρ)|||L2 ([0,T ];L2(Ω)) + ≤ Cρ E, |e| e∈Γ − c)|||L2 ([0,T ];L2 (Ω)) +

0  ZT 0

h

2 X r2 2 c DG k[c − c]k0,e ≤ Cc E, |e| e∈Γ 1

h

hmin(rc +1,sc )−1 hmin(ru +1,su )−1 hmin(rv +1,sv )−1 + + + rcsc −2 rusu −2 rvsv −2



.

20

Y. Epshteyn and A. Kurganov

Proof: The result follows from the definition of space S, the fact that the DG solution is a fixed point of the compact operator A (defined above), the hp Approximation Lemma 2.1, and the triangle inequality.  Remark. The obtained error estimates are h-optimal, but only suboptimal for r. Finally, equipped with the results established in Theorem 5.5, we obtain the following bound for the blow-up time of the exact solution of the Keller-Segel system. Theorem 5.6 Let us denote by tb the blow-up time of the exact solution of the Keller-Segel system (1.1) and by tDG the blow-up time of the DG solution of (3.3)–(3.9). Then tb ≤ tDG b b . Proof: The solution ρ of the Keller-Segel model blows up if kρkL∞ (Ω) becomes unbounded in either finite or infinite time (see, e.g., [26, 27]). Therefore, in order to prove the theorem we need to establish an L∞ -error bound. Consider interpolant ρe From Theorem 5.5 we have the following L2 -error bound:   min(rρ +1,sρ )−1 hmin(rc +1,sc )−1 hmin(ru +1,su )−1 hmin(rv +1,sv )−1 h DG + + , + kρ − ρekL2 (Ω) ≤ Cρ s −2 rcsc −2 rusu −2 rvsv −2 rρρ which together with the inverse inequality (2.6) leads to the desired L∞ -error bound,  min(rρ +1,sρ )−2  h hmin(rc +1,sc )−2 hmin(ru +1,su )−2 hmin(rv +1,sv )−2 DG kρ − ρekL∞ (Ω) ≤ Cρ + + , + s −3 rcsc −3 rusu −3 rvsv −3 rρρ

which, in turn, implies that DG

kρ

kL∞ (Ω) ≤ ke ρkL∞ (Ω) +Cρ



hmin(rρ +1,sρ )−2 s −3

rρρ

hmin(rc +1,sc )−2 hmin(ru +1,su )−2 hmin(rv +1,sv )−2 + + + rcsc −3 rusu −3 rvsv −3

From the last estimate the statement of the theorem follows.

6



.



Numerical Example

In this section, we demonstrate the performance of the proposed DG method. In all our numerical experiments, we have used the third-order strong stability preserving Runge-Kutta method for the time discretization, [23]. No slope limiting technique has been implemented. The values of the penalty parameters used are σρ = σc = 1 and σu = σv = 0.01. We note that no instabilities have been observed when the latter two parameters were taken zero, however, since our convergence proof requires σu and σv to be positive, we only show the results obtained with positive σu and σv , which are almost identical to the ones obtained with σu = σv = 0. We consider the initial-boundary value problem for the Keller-Segel system in the square domain [− 21 , 21 ] × [− 21 , 12 ]. We take the chemotactic sensitivity χ = 1 and the bell-shaped initial data 2 2 2 2 ρ(x, y, 0) = 1200e−120(x +y ) , c(x, y, 0) = 600e−60(x +y ) . According to the results in [25], both components ρ and c of the solution are expected to blow up at the origin in finite time. This situation is especially challenging since capturing blowing up solution with shrinking support is extremely hard.

21

DG Methods for the Keller-Segel System

In Figures 6.1–6.4, we plot the logarithmically scaled density, ln(1+ρDG ), computed at different times on two different uniform grids with h = 1/51 (Figures 6.1 and 6.3) and h = 1/101 (Figures 6.2 and 6.4). The results shown in Figures 6.1–6.2 have been obtained with quadratic polynomials (i.e., rρ = rc = ru = rv = r = 2), while the solution shown in Figures 6.3–6.4 have been computed with the help of cubic polynomials (i.e., rρ = rc = ru = rv = r = 3). Numerical convergence of the scheme is verified by refining the mesh and by increasing the polynomial degree. As one can see, the computed solutions in a very good agreement at the smaller times (t = 1.46 · 10−5 , 2.99 · 10−5 , and 6.03 · 10−5 ). However, at time close to the blowup time (t = 1.21 · 10−4 ) the maximum value of ρDG grows while its support shrinks, and no mesh-refinement convergence is observed: the numerical solution keeps increasing when the mesh is refined. Using Theorem 5.6, we can conclude that in this example, the blow-up time of the exact solution is less or equal to the blow-up time of the DG solution, which is approximately tDG ≈ 1.21 · 10−4 . b We note that even though no slope limiting or any other positivity preserving techniques have been implemented, the computed solutions have never developed negative values and are oscillation-free.

12

12

10

10

8

8

6

6

4

4

2

2

0 0.5

0 0.5 0.5

0.5

0

0 0 −0.5

0 −0.5

−0.5

12

12

10

10

8

8

6

6

4

4

2

2

0 0.5

0 0.5

−0.5

0.5 0

0.5 0

0 −0.5

−0.5

0 −0.5

−0.5

Figure 6.1: h = 1/51, r = 2. Logarithmically scaled density computed at t = 1.46 · 10−5 (top left), t = 2.99 · 10−5 (top right), t = 6.03 · 10−5 (bottom left), and t = 1.21 · 10−4 ≈ tDG (bottom right). b Finally, we check the numerical order of the convergence of the proposed DG method. We first consider the smooth solution at a very small time t = 1.0 · 10−7 and test the convergence with

22

Y. Epshteyn and A. Kurganov

12

12

10

10

8

8

6

6

4

4

2

2

0 0.5

0 0.5 0.5

0.5

0

0 0 −0.5

0 −0.5

−0.5

12

12

10

10

8

8

6

6

4

4

2

2

0 0.5

−0.5

0 0.5 0.5 0

0.5 0

0 −0.5

−0.5

0 −0.5

−0.5

Figure 6.2: The same as in Figure 6.1 but with h = 1/101, r = 2.

respect to the mesh size h for the fixed r = 2 (piecewise quadratic polynomials). Since the exact solution for the Keller-Segel system is unavailable, we compute the reference solution by the proposed DG method on a fine mesh with h = 1/128 and using the fifth-order (r = 5) piecewise polynomials. We then use the obtained reference solution to compute the relative L2 - and relative H 1 -errors. These errors are presented in Table 6.1. From this table, one can see that the solution numerically converges to the reference solution with the (optimal) second order in the H 1 -norm which confirms the theoretical results predicted by our convergence analysis. Moreover, the achieved third order of convergence in the L2 -norm is optimal for quadratic piecewise polynomials. We then test the convergence of the proposed DG method with respect to the degree r of piecewise polynomials for the fixed h = 1/32. The obtained results, reported in Table 6.2, show that the error decreases almost exponentially when the polynomial degree increases (this is a typical situation when DG methods capture smooth solutions). We also compute the L2 -errors with respect to the reference solution, for the solutions plotted on Figures 6.1 and 6.2 at times t = 2.99 · 10−5 and t = 6.03 · 10−5. These times are close to the blowup time and the solutions develop a pick at the origin. The obtained errors are reported in Table 6.3. As one can see, even for the spiky solutions, the convergence rate is very high though it, as expected, deteriorates as t approaches tDG b .

23

DG Methods for the Keller-Segel System

12

12

10

10

8

8

6

6

4

4

2

2

0 0.5

0 0.5 0.5

0.5

0

0 0 −0.5

0 −0.5

−0.5

12

12

10

10

8

8

6

6

4

4

2

2

0 0.5

0 0.5

−0.5

0.5 0

0.5 0

0 −0.5

0 −0.5

−0.5

−0.5

Figure 6.3: The same as in Figures 6.1–6.2 but with h = 1/51, r = 3. h 1/4 1/8 1/16 1/32 1/64

L2 -error 3.0578 1.0290 0.0796 0.0075 0.0006

Rate – 1.6 3.7 3.4 3.6

H 1 -error

Rate

1.5591 1.2348 0.5206 0.0937 0.0157

– 0.35 1.3 2.5 2.6

Table 6.1: Relative errors as functions of the mesh size h; r = 2 is fixed. r

L2 -error

Rate

H 1 -error

2 3 4 5

7.5e-03 9.0e-04 8.0e-05 6.9e-06

– 5.2 8.4 11.0

9.4e-02 2.2e-02 2.6e-03 2.9e-04

Rate – 3.6 7.4 9.8

Table 6.2: Relative errors as functions of the piecewise polynomial degree r; h = 1/32 is fixed.

24

Y. Epshteyn and A. Kurganov

12

12

10

10

8

8

6

6

4

4

2

2

0 0.5

0 0.5 0.5

0.5

0

0 0 −0.5

0 −0.5

−0.5

12

12

10

10

8

8

6

6

4

4

2

2

0 0.5

−0.5

0 0.5 0.5

0.5

0

0 0 −0.5

0 −0.5

−0.5

−0.5

Figure 6.4: The same as in Figures 6.1–6.3 but with h = 1/101, r = 3.

Appendix A: Proof of Theorem 5.3 — Continuation In this appendix, we complete the proof of Theorem 5.3 by proving that (φcL , φuL , φvL) ∈ S. We begin with φcL and show that φcL ∈ S in a way similar to the proof of the fact that φρL ∈ S given in §5. First, from the consistency Lemma 4.1 we obtain that the exact solution of (1.3)–(1.6) satisfies equation (3.4), which may be rewritten as Z Ω

XZ

XZ

XZ

X r2 Z c {∇w · ne }[e c] + σc {∇e c · ne }[w ] + ε ∇e c∇w − e ct w + [e c][w c ] |e| e∈Γh e∈Γh e e∈Γh e E∈Eh E e Z Z Z XZ XZ + e cw c − ρew c = − ξtc w c − {∇ξ c · ne }[w c ] ∇ξ c ∇w c + c

Ω

−ε

XZ

e∈Γh e

c

Ω

Ω

c

c

{∇w · ne }[ξ ] − σc

c

E∈Eh E X r2 Z c

e∈Γh

|e|

c

c

[ξ ][w ] −

e

Z Ω

c

e∈Γh e c

c

ξw +

Z

ξ ρwc.

Ω

We then subtract equation (A.1) from equation (5.5) and set w c = τ c to obtain X r2 1d c 2 c kτ k0,Ω + kτ c k20,Ω + |||∇τ c |||20,Ω + σc k[τ c ]k20,e 2 dt |e| e∈Γ h

(A.1)

25

DG Methods for the Keller-Segel System

t = 2.99 · 10−5

h 1/51 1/101

L2 -error 5.5e-02 5.2e-03

t = 6.03 · 10−5

L2 -error 5.0e-02 1.1e-02

Rate – 3.4

Rate – 2.2

Table 6.3: Relative L2 -errors at two different times; r = 2 is fixed. =

Z Ω

+ε

ρ c

τ τ + (1 − ε) XZ

e∈Γh e

XZ

e∈Γh e

c

c

{∇τ · ne }[τ ] +

{∇τ c · ne }[ξ c ] + σc

X r2 c |e| e∈Γ h

Z

Z

ξtc τ c

Ω

[ξ c ][τ c ] +

e

+

XZ

E∈Eh E

Z

ξ cτ c −

Ω

c

c

∇ξ ∇τ − Z

XZ

e∈Γh e

{∇ξ c · ne }[τ c ]

ξ ρ τ c =: T1c + ... + T9c .

(A.2)

Ω

Next, we bound each term on the RHS of (A.2). We begin with the term T1c . We first bound it using the Cauchy-Schwarz and Young’s inequalities, and then apply the estimate (5.33). This results in 1 |T1c | ≤ kτ c k20,Ω + kτ ρ k20,Ω 4 2 min(rρ +1,sρ )−2  I C h h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 c 2 + + . (A.3) ≤ kτ k0,Ω + + 2s −4 4 rc2sc −4 ru2su −4 rv2sv −4 rρ ρ ρ The terms T2c , T3c , T4c , T5c , T6c , and T7c are similar to the terms T1ρ , T6ρ , T7ρ , T8ρ , T9ρ , and T10 estimated in (5.17), (5.23), (5.24), (5.25), (5.26), and (5.27), respectively. Hence, they can be bounded as follows: X r2 c (A.4) k[τ c ]k20,e , |T2c | ≤ εc2 |||∇τ c |||20,Ω + C2c |e| e∈Γ h

|T3c | ≤ kτ c k20,Ω + C ∗

2 min(rc +1,sc )

h

, rc2sc h2 min(rc +1,sc )−2 , |T4c | ≤ εc4 |||∇τ c |||20,Ω + C ∗ rc2sc −2 2 min(rc +1,sc )−2 X r2 c c c ∗h c 2 , |T5 | ≤ C5 k[τ ]k0,e + C |e| rc2sc −2 e∈Γ

(A.5) (A.6) (A.7)

h

h2 min(rc +1,sc )−2 , rc2sc −4 X τ2 h2 min(rc +1,sc )−2 c |T7c | ≤ C6c . k[τ c ]k20,e + C ∗ |e| rc2sc −4 e∈Γ

|T6c | ≤ εc6 |||∇τ c |||20,Ω + C ∗

(A.8) (A.9)

h

Finally, the last two terms on the RHS of (A.2), T8c and T9c , are bounded using the CauchySchwarz inequality, Young’s inequality, and the approximation inequality (2.2): |T8c | ≤ kτ c k20,Ω + C ∗

h2 min(rc +1,sc ) , rc2sc

|T9c | ≤ kτ c k20,Ω + C ∗

h2 min(rρ +1,sρ ) 2sρ

rρ

.

(A.10)

26

Y. Epshteyn and A. Kurganov

Now substituting (A.3)–(A.10) into (A.2) and using the assumption that h < 1, we obtain the following estimate for τ c : X r2 1d c 2 c kτ k0,Ω + (1 − εc2 − εc4 − εc6 )|||∇τ c |||20,Ω + (σc − C2c − C5c − C6c ) k[τ c ]k20,e 2 dt |e| e∈Γh  2 min(rρ +1,sρ )−2  2 min(ru +1,su )−2 2 min(rc +1,sc )−2 h h h2 min(rv +1,sv )−2 h c c 2 ≤ 3kτ k0,Ω + C∗ + + .(A.11) + 2s −4 rc2sc −4 ru2su −4 rv2sv −4 rρ ρ This estimate is similar to the estimate (5.31). After a proper selection of εci and the penalty parameter σc , we multiply (A.11) by 2, integrate with respect to time from 0 to t, apply Gronwall’s Lemma 2.6, and take the supremum over t. This results in an estimate, which is completely analogous to (5.33):  ZT  X r2 c c 2 c 2 + |||∇τ |||0,Ω + k[τ ]k0,e |e| e∈Γh 0   2 min(rρ +1,sρ )−2 h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 II h ,(A.12) + + + ≤C 2s −4 rc2sc −4 ru2su −4 rv2sv −4 rρ ρ

sup kτ c k20,Ω [0,T ]

where C II is a constant that depends on kρk(L∞ ([0,T ]);H 2 (Ω)) , kρt k(L∞ ([0,T ]);L2 (Ω)) , kck(L∞ ([0,T ]);H 2 (Ω)) , kct k(L∞ ([0,T ]);L2 (Ω)) , kuk(L∞ ([0,T ]);L2 (Ω)) , kvk(L∞ ([0,T ]);L2 (Ω)) , and T only. Hence, according to the definition on page 9, the estimate (A.12) implies that φcL ∈ S.

Next, we proceed with proving that φuL ∈ S. Once again, by the consistency Lemma (4.1), the exact solution satisfies the following equation (compare it with (3.5)): Z Ω

u

u ew +

=−

Z

XZ

E∈Eh E

e c(w )x +

ξ uwu −

Ω

u

XZ

X Z

u (−c)∗∗ u nx [w ]

e∈Γver h e

ξ c (w u )x +

E∈Eh E

X Z

e∈∂Ωver e

where

e∈∂Ωver e

ξ c nx w u − σu

u

e cnx w + σu

X

e∈Γh ∪∂Ωver

X

e∈Γh ∪∂Ωver

1

(−c)∗∗ u

−

X Z

ru2

|e|

Z

ru2 |e|

Z

[e u][w u ]

e

[ξ u ][w u ],

(A.13)

e

2

in in E aout aout cE L aL e − aL ce − [u]. := − L out in aL − ain aout L L − aL

Subtracting equation (A.13) from (5.6) and choosing w u = τ u , we obtain X

ru2 k[τ u ]k20,Ω |e| e∈Γh ∪∂Ωver Z Z X Z X Z X c u c ∗ ∗∗ u c u τ nx [τ ] + ξ u τ u ((−φL )u − (−c)u )nx [τ ] + τ (τ )x − =−

kτ u k20,Ω + σu

E∈Eh E

+

XZ

E∈Eh E

ξ c (τ u )x −

e∈Γver h e

X Z

e∈∂Ωver e

ξ c nx [τ u ] + σu

X

e∈Γh ∪∂Ωver

ru2

|e|

Z e

e∈∂Ωver e

Ω

[ξ u ][τ u ] =: T1u + ... + T7u , (A.14)

27

DG Methods for the Keller-Segel System and bound each term on the RHS of (A.14). To estimate the term T1u , we first integrate by parts and rewrite it as   XZ XZ X XZ c u c u u (τ c )x τ u − τ τ nx = (τ )x τ + T1 = − −

[τ c τ u ]nx .

e∈Γver h ∪∂Ωver e

E∈Eh E

E∈Eh e∈∂E e

E∈Eh E

Z

X

Then, using the formula for the jump and the average operators (2.1), we obtain X Z X Z X Z XZ u c c u c u u [τ ]{τ }nx − [τ ]{τ }nx − (τ )x τ − τ c [τ u ]nx . T1 = e∈Γver h e

e∈Γver h e

E∈Eh E

e∈∂Ωver e

Hence, using the Cauchy-Schwarz inequality, Young’s inequality, the inequality (2.5), and applying the assumption (3.2), we arrive at the following bound for T1u : |T1u | ≤

9 u 2 kτ k0,Ω + C1u 16

X

e∈Γh ∪∂Ωver

X r2 1 ru2 c k[τ u ]k20,e + |||(τ c )x |||20,Ω + C2u kτ c k20,Ω + C3u k[τ c ]k20,e . |e| 2 |e| e∈Γ h

(A.15) A bound for T2u can be obtained in a way similar to the one the bound on T3ρ has been established: Z     X  Z aout −ain L L u c E1 c E2 E1 u E2 u |T2 | ≤ aout − ain (φL )e − ce nx [τ ] + aout − ain (φL)e − ce nx [τ ] e∈Γver h

L

L

e  Z in out −a a L L u u [φL − u]nx [τ ] := I + II + III. + out in aL − aL

e

L

L

(A.16)

e

From (5.10) and (5.13), the first term on the RHS of (A.16) can be estimated by Z  X  Z 1 1 (τ c )E nx [τ u ] + (ξ c )E nx [τ u ] := eI. I≤ e e e∈Γh

e

e

Using then the Cauchy-Schwarz inequality, the trace inequality (2.3), the inequality (2.5), and the assumption (3.2), we estimate eI as follows: X r2 h2 min(rc +1,sc ) u u 2 eI ≤ 1 kτ c k2 + K . k[τ ]k0,e + C 0,Ω 2 |e| rc2sc e∈Γ h

A similar bound can be derived for the second term on the RHS of (A.16). The third term on the RHS of (A.16) is similar to the third term on the RHS of (5.19), hence it can be bounded by X 2  ru K1 h h2 min(ru +1,su ) u 2 . + K III ≤ k[τ ]k + C 2 0,e 2su ru2 |e| r u e∈Γ h

Combining the above bounds on I, II, and III, we arrive at   2 min(ru +1,su ) X r2 h2 min(rc +1,sc ) h u u c 2 u u 2 ∗ |T2 | ≤ kτ k0,Ω + C4 + . k[τ ]k0,e + C 2su 2sc |e| r r u c e∈Γ h

(A.17)

28

Y. Epshteyn and A. Kurganov

To bound the term T3u , we use the Cauchy-Schwarz inequality, Young’s inequality, and the inequality (2.5), which yield |T3u | ≤ C5u kτ c k20,Ω + C6u

X r2 u k[τ u ]k20,e . |e| e∈∂Ω

(A.18)

ver

The term T4u is bounded with the help of Cauchy-Schwarz inequality, Young’s inequality, and the approximation inequality (2.2): |T4u | ≤

h2 min(ru +1,su ) 1 u 2 . kτ k0,Ω + C ∗ 16 ru2su

(A.19)

Using the Cauchy-Schwarz inequality and the inverse inequality (2.7), we first bound T5u by X X kξ c k0,E h−1 ru2 kτ u k0,E := Te5u . (A.20) kξ ck0,E k(τ u )x k0,E ≤ |T5u | ≤ E∈Eh

E∈Eh

We then use Young’s inequality, the assumption (3.2), and the approximation inequality (2.2) to obtain 1 h2 min(rc +1,sc )−2 Te5u ≤ kτ u k20,Ω + C ∗ . (A.21) 16 rc2sc −4

The term T6u is bounded using the Cauchy-Schwarz inequality, the trace inequality (2.3), and the approximation inequality (2.2): |T6u |

≤

C7u

2 min(rc +1,sc ) X r2 u u 2 ∗h . k[τ ]k0,e + C |e| rc2sc e∈∂Ω

(A.22)

ver

ρ , estimated in (5.27). Hence, The last term T7u is similar to term T10

|T7u | ≤ C8u

X

e∈Γh ∪∂Ωver

ru2 h2 min(ru +1,su )−2 . k[τ u ]k20,e + C ∗ |e| ru2su −4

(A.23)

After obtaining the estimates (A.15) and (A.17)–(A.23), we plug them into (A.14) and use the assumption h < 1 to obtain X ru2 5 u 2 u u u u u kτ k0,Ω + (σu − C1 − C4 − C6 − C7 − C8 ) k[τ u ]k20,e ≤ (1 + C2u + C5u )kτ c k20,Ω 16 |e| e∈Γh ∪∂Ωver   X r2 1 h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 2 c u c c 2 ∗ + |||(τ )x |||0,Ω + C3 + . (A.24) k[τ ]k0,e + Cu 2sc −2 2su −4 2 |e| r r c u e∈Γ h

In the same way as we have derived the estimate (A.24), we can establish the following bound: X

rv2 k[τ v ]k20,e ≤ (1 + C2v + C5v )kτ c k20,Ω |e| e∈Γh ∪∂Ωhor   X r2 1 h2 min(rc +1,sc )−2 h2 min(rv +1,sv )−2 2 c c v c 2 ∗ + |||(τ )y |||0,Ω + C3 + . (A.25) k[τ ]k0,e + Cv 2sc −4 2sv −4 2 |e| r r c v e∈Γ

5 v 2 kτ k0,Ω + (σv − C1v − C4v − C6v − C7v − C8v ) 16

h

29

DG Methods for the Keller-Segel System

Next, we use Lemma 2.4 to bound kτ c k20,Ω on the RHS of (A.24) and (A.25). This results in X ru2 5 u 2 kτ k0,Ω + (σu − C1u − C4u − C6u − C7u − C8u ) k[τ u ]k20,e 16 |e| e∈Γh ∪∂Ωver 1  X r2 c ≤ + K(1 + C2u + C5u ) |||∇τ c |||20,Ω + (C3u + K(1 + C2u + C5u )) k[τ c ]k20,e 2 |e| e∈Γh   2 min(rc +1,sc )−2 h h2 min(ru +1,su )−2 , (A.26) + Cu∗ + rc2sc −2 ru2su −4 and X rv2 5 v 2 kτ k0,Ω + (σv − C1v − C4v − C6v − C7v − C8v ) k[τ v ]k20,e 16 |e| e∈Γh ∪∂Ωhor  1 X r2 c + K(1 + C2v + C5v ) |||∇τ c |||20,Ω + (C3v + K(1 + C2v + C5v )) k[τ c ]k20,e ≤ 2 |e| e∈Γh   2 min(rc +1,sc )−2 2 min(rv +1,sv )−2 h h + . (A.27) + Cv∗ 2s −2 rc c rv2sv −4 We then multiply both sides of (A.26) and (A.27) by 16/5, choose the appropriate penalty parameters σu and σv , integrate with respect to time from 0 to T , and use the estimate (A.12) to obtain ZT  0

kτ u k20,Ω

≤C

III

+



X

e∈Γh ∪∂Ωver

ru2 k[τ u ]k20,e |e|

h2 min(rρ +1,sρ )−2 2sρ −4

rρ



h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 + + + rc2sc −4 ru2su −4 rv2sv −4



(A.28)

and ZT  kτ v k20,Ω + 0

≤C

IV



X

e∈Γh ∪∂Ωhor

rv2 k[τ v ]k20,e |e|

h2 min(rρ +1,sρ )−2 2sρ −4

rρ



h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 + + + rc2sc −4 ru2su −4 rv2sv −4



where C III and C IV are constants that depend on kρk(L∞ ([0,T ]);H 2 (Ω)) , kρt k(L∞ ([0,T ]);L2 (Ω)) , kck(L∞ ([0,T ]);H 2 (Ω)) , kct k(L∞ ([0,T ]);L2 (Ω)) , kuk(L∞ ([0,T ]);L2 (Ω)) , kvk(L∞ ([0,T ]);L2 (Ω)) , and T only.

,(A.29)

We now estimate the RHS of (A.24) in a different way: we apply the inequality (2.5) and the inverse inequality (2.7), which yield X ru2 5 u 2 kτ k0,Ω + (σu − C1u − C4u − C6u − C7u − C7u ) k[τ u ]k20,e 16 |e| e∈Γh ∪∂Ωver   rc4 c 2 h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 ∗ + . ≤ Ku 2 kτ k0,Ω + Cu h rc2sc −2 ru2su −4

30

Y. Epshteyn and A. Kurganov

We then take the supremum over t, choose the appropriate penalty parameters σu and σv , and use the estimate (A.12) to obtain   X ru2 v 2 u 2 k[τ ]k0,e sup kτ k0,Ω + |e| [0,T ] e∈Γh ∪∂Ωver  2 min(rρ +1,sρ )−4  h h2 min(rc +1,sc )−4 h2 min(ru +1,su )−4 h2 min(rv +1,sv )−4 ≤C + + . + 2s −8 rc2sc −8 ru2su −8 rv2sv −8 rρ ρ Finally, using the assumptions on r, s, and h, we conclude that     X ru2 1 1 1 1 u 2 ∗ 2 v 2 sup kτ k0,Ω + + + + , k[τ ]k0,e ≤ Cu h |e| rρ4 rc4 ru4 rv4 [0,T ]

(A.30)

e∈Γh ∪∂Ωver

where the constant Cu∗ is independent of h and r. The bound on τ v is obtained similarly:     X 1 1 1 1 rv2 ∗ 2 v 2 v 2 + + + , k[τ ]k0,e ≤ Cv h sup kτ k0,Ω + |e| rρ4 rc4 ru4 rv4 [0,T ] e∈Γ ∪∂Ω h

(A.31)

hor

where Cv∗ is independent of h and r. According to the definition on page 9, the estimates (A.28)–(A.29) and (A.30)–(A.31) ensure that (φuL , φvL ) ∈ S. Acknowledgment: The research of Y.Epshteyn is based upon work supported by the Center for Nonlinear Analysis (CNA) under the National Science Foundation Grant # DMS-0635983. The research of A. Kurganov was supported in part by the NSF Grant # DMS-0610430.

References [1] J. Adler, Chemotaxis in bacteria, Ann. Rev. Biochem., 44 (1975), pp. 341–356. [2] S. Agmon, Lectures on Elliptic Boundary Value Problems, Van Nostrand, Princeton, NJ, 1965. [3] V. Aizinger, C. Dawson, B. Cockburn, and P. Castillo, Local discontinuous Galerkin methods for contaminant transport, Adv. in Water Res., 24 (2000), pp. 73–87. [4] D.N. Arnold, An interior penalty finite element method with discontinuous elements, SIAM J. Numer. Anal., 19 (1982), pp. 742–760. [5] I. Babuˇ ska and M. Suri, The h-p version of the finite element method with quasiuniform meshes, RAIRO Mod´el. Math. Anal. Num´er., 21 (1987), pp. 199–238. [6] I. Babuˇ ska and M. Suri, The optimal convergence rates of the p-version of the finite element method, SIAM J. Numer. Anal., 24 (1987), pp. 750–776. [7] G.A. Baker, W.N. Jureidini, and O.A. Karakashian, Piecewise solenoidal vector fields and the Stokes problem, SIAM J. Numer. Anal., 27 (1990), pp. 1466–1485.

DG Methods for the Keller-Segel System

31

[8] J.T. Bonner, The cellular slime molds, Princeton University Press, Princeton, New Jersey, 2nd ed., 1967. [9] S. Brenner, Poincar´e-Friedrichs inequalities for piecewise H 1 functions, SIAM J. Numer. Anal., 41 (2003), pp. 306–324. [10] E.O. Budrene and H.C. Berg, Complex patterns formed by motile cells of escherichia coli, Nature, 349 (1991), pp. 630–633. [11] E.O. Budrene and H.C. Berg, Dynamics of formation of symmetrical patterns by chemotactic bacteria, Nature, 376 (1995), pp. 49–53. [12] A. Chertock and A. Kurganov, A positivity preserving central-upwind scheme for chemotaxis and haptotaxis models, Numer. Math. submitted. [13] S. Childress and J.K. Percus, Nonlinear aspects of chemotaxis, Math. Biosc., 56 (1981), pp. 217–237. [14] B. Cockburn, G.E. Karniadakis, and C.-W. Shu, eds., First International Symposium on Discontinuous Galerkin Methods, vol. 11 of Lecture Notes in Computational Science and Engineering, Springer-Verlag, 2000. [15] B. Cockburn and C.-W. Shu., The local discontinuous Galerkin method for convectiondiffusion systems, SIAM J. Numer. Anal., 35 (1998), pp. 2440–2463. [16] M.H. Cohen and A. Robertson, Wave propagation in the early stages of aggregation of cellular slime molds, J. Theor. Biol., 31 (1971), pp. 101–118. [17] C. Dawson, E. Kubatko, and J. Westerink, hp Discontinuous Galerkin methods for advection-dominated problems in shallow water, Comp. Meth. Appl. Mech. Eng. to appear. [18] C. Dawson, S. Sun, and M.F. Wheeler, Compatible algorithms for coupled flow and transport, Comput. Meth. Appl. Mech. Engng, 193 (2004), pp. 2565–2580. [19] J. Douglas and T. Dupont, Lecture Notes in Physics, vol. 58, Springer-Verlag, 1976, ch. Interior penalty procedures for elliptic and parabolic Galerkin methods. `re, On the solution of incompressible two-phase flow by a p[20] Y. Epshteyn and B. Rivie version discontinuous Galerkin method, Comm. Numer. Methods Engrg., 22 (2006), pp. 741– 751. [21] F. Filbet, A finite volume scheme for the Patlak-Keller-Segel chemotaxis model, Numer. Math., 104 (2006), pp. 457–488. `re, and M.F. Wheeler, A discontinuous Galerkin method with [22] V. Girault, B. Rivie non-overlapping domain decomposition for the Stokes and Navier-Stokes problems, Math. Comp., 74 (2005), pp. 53–84. [23] S. Gottlieb, C.-W. Shu, and E. Tadmor, High order time discretization methods with the strong stability property, SIAM Review, 43 (2001), pp. 89–112.

32

Y. Epshteyn and A. Kurganov

[24] D.H. Griffel, Applied Functional Analysis, Dover, 2002. ´zquez, A blow-up mechanism for a chemotaxis model, [25] M.A. Herrero and J.J.L. Vela Ann. Scuola Normale Superiore, 24 (1997), pp. 633–683. [26] D. Horstmann, From 1970 until now: The Keller-Segel model in chemotaxis and its consequences I, Jahresber. DMV, 105 (2003), pp. 103–165. [27] D. Horstmann, From 1970 until now: The Keller-Segel model in chemotaxis and its consequences II, Jahresber. DMV, 106 (2004), pp. 51–69. [28] E.F. Keller and L.A. Segel, Initiation of slime mold aggregation viewed as an instability, J. Theor. Biol., 26 (1970), pp. 399–415. [29] E.F. Keller and L.A. Segel, Model for chemotaxis, J. Theor. Biol., 30 (1971), pp. 225– 234. [30] E.F. Keller and L.A. Segel, Traveling bands of chemotactic bacteria: A theoretical analysis, J. Theor. Biol., 30 (1971), pp. 235–248. [31] A. Kurganov and C.-T. Lin, On the reduction of numerical dissipation in central-upwind schemes, Commun. Comput. Phys., 2 (2007), pp. 141–163. [32] A. Kurganov, S. Noelle, and G. Petrova, Semi-discrete central-upwind schemes for hyperbolic conservation laws and Hamilton-Jacobi equations, SIAM J. Sci. Comput., 23 (2001), pp. 707–740. [33] A. Kurganov and G. Petrova, Central-upwind schemes on triangular grids for hyperbolic systems of conservation laws, Numer. Methods Partial Differential Equations, 21 (2005), pp. 536–552. [34] A. Marrocco, 2D simulation of chemotaxis bacteria aggregation, M2AN Math. Model. Numer. Anal., 37 (2003), pp. 617–630. [35] V. Nanjundiah, Chemotaxis, signal relaying and aggregation morphology, J. Theor. Biol., 42 (1973), pp. 63–105. ¨li, Discontinuous Galerkin finite element approximation of nonlinear [36] C. Ortner and E. Su second-order elliptic and hyperbolic systems., SIAM J. Numer. Anal., 45 (2007), pp. 1370– 1391. [37] C.S. Patlak, Random walk with persistence and external bias, Bull. Math: Biophys., 15 (1953), pp. 311–338. [38] L.M. Prescott, J.P. Harley, and D.A. Klein, Microbiology, Wm. C. Brown Publishers, Chicago, London, 3rd ed., 1996. [39] B. Rivi` ere, M.F. Wheeler, and V. Girault., A priori error estimates for finite element methods based on discontinuous approximation spaces for elliptic problems, SIAM J. Numer. Anal., 39 (2001), pp. 902–931.

DG Methods for the Keller-Segel System

33

[40] S. Sun and M.F. Wheeler, Symmetric and nonsymmetric discontinuous Galerkin methods for reactive transport in porous media, SIAM J. Numer. Anal., 43 (2005), pp. 195–219. [41] R. Tyson, S.R. Lubkin, and J.D. Murray, A minimal mechanism for bacterial pattern formation, Proc. Roy. Soc. Lond. B, 266 (1999), pp. 299–304. [42] R. Tyson, L.G. Stern, and R.J. LeVeque, Fractional step methods applied to a chemotaxis model, J. Math. Biol., 41 (2000), pp. 455–475.