Fully Discrete Analysis of a Discontinuous Finite Element ... - CMU Math

Comment

Report 1 Downloads 46 Views

c

2008 001

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

Fully Discrete Analysis of a Discontinuous Finite Element Method for the Keller-Segel Chemotaxis Model Yekaterina Epshteyn∗and Ahmet Izmirlioglu† Abstract This paper formulates and analyzes fully discrete schemes for the two-dimensional Keller-Segel chemotaxis model. The spatial discretization of the model is based on the discontinuous Galerkin methods and the temporal discretization is based either on Forward Euler or the second order explicit total variation diminishing (TVD) Runge-Kutta methods. We consider Cartesian grids and prove optimal fully discrete error estimates for the proposed methods. Our proof is valid for pre-blow-up times since we assume boundedness of the exact solution.

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

1

Introduction

The goal of this work is to formulate and analyze fully discrete discontinuous Galerkin (DG) methods for the solution of the two-dimensional (2-D) Keller-Segel chemotaxis model, [11, 26, 27, 28, 33, 34]. The underlying spatial discretization of the model is based on the methods proposed recently in [18] and the temporal discretization is based either on Forward Euler or the second order explicit total variation diminishing (TVD) Runge-Kutta methods. In this paper, we consider the classical formulation of the Keller-Segel system [11], 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 ∗

Department of Mathematical Sciences, Carnegie Mellon University, Pittsburgh, PA, [email protected] † Department of Mathematics, University of Pittsburgh, Pittsburgh, PA 15260; [email protected]

1

15213,

2

Y. Epshteyn and A. Izmirlioglu

vector. The system (1.1) is the basic step in the modeling of many biological processes. The Keller-Segel model (1.1) can be generalized to better describe the reality by taking into account some other factors such as growth and death of cells, presence of the food and other chemicals in the system, etc. It is well-known that solutions of the classical Keller-Segel system may blow up in finite time, see, e.g., [24, 25] 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, 6, 8, 9, 15, 35]. Capturing blowing up solutions numerically is a challenging problem. A finite-volume, [21], and a finite-element, [32], 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 [39] has been proposed in [40]. 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 [10], where the finitevolume Godunov-type central-upwind scheme was derived for (1.1) and extended to some other chemotaxis and haptotaxis models. The high-order discontinuous Galerkin method that is investigated here is based on the method proposed in [18]. The DG methods have recently become increasingly popular thanks to their flexibility for adaptive simulations, suitability for parallel computations, applicability to problems with discontinuous coefficients and/or solutions, and compatibility with other numerical methods. These methods have been successfully applied to a wide variety of problems, ranging from the solid mechanics to the fluid mechanics (see, e.g., [13, 14, 12, 19, 20, 22, 38] and references therein). Furthermore, the DG methods are among the methods that can be used for real biomedical problems, which are often considered in complex domains, have discontinuity in the coefficients, and incorporate PDEs of different mathematical nature. In order to develop high-order DG methods for (1.1) in [18], the Keller-Segel system is rewritten as a system of the nonlinear convection-diffusion-reaction equations by introducing new variables (u, v) := ∇c : kQt + F(Q)x + G(Q)y = k∆Q + R(Q), (1.2) 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.2), and k = 0 in the third and the fourth equations there. The methods proposed in [18] 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, [3, 16, 36]. The numerical fluxes in the proposed DG methods are the fluxes developed for the semidiscrete finite-volume central-upwind schemes in [30] (see also [29, 31]). 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

3

Analysis of a DG Method for the Chemotaxis Model

speeds of propagation—is incorporated into the central-upwind fluxes. In [18], Cartesian grids are considered and the continuous in time error estimates are proved for the proposed high-order DG methods under the assumption of boundedness of the exact solution. Some numerical tests that validate the methods are considered in [18] and [17]. In the following section, we introduce our notations, assumptions, and state some standard results. In §3-§4 we recall some results for the continuous in time scheme. In §5-§7 we formulate the explicit schemes and derive the error estimates under assumption of boundedness of the exact solution (some proof details are postponed to Appendix 9). The proof of the error estimates is based on the induction argument which simplifies the analysis significantly since we consider the coupled system of the nonlinear equations.

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 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

[w] := weE ,

2

1

2

[w] := weE − weE ,

1

2

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

{w} := weE ,

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 §6-§9. First, let us state some approximations properties and inequalities for the finite-element space.

4

Y. Epshteyn and A. Izmirlioglu

Lemma 2.1 (hp Approximation, [4, 5]) Let E ∈ Eh and ψ ∈ H s (E). 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

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

hµ− 2

h e

ψ − ψr ≤ C s− 1 kψks,E , 0,e r 2 ≤C

s ≥ 0,

(2.2)

1 s> , 2

(2.3)

where µ := min(r + 1, s) and e is the edge on ∂E.

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

where e is the edge on ∂E.

Lemma 2.3 ([36]) Let E be a mesh element with an edge e. Then there is a constant Ct independent of h and r such that 1

kwk0,e ≤ Ct h− 2 rkwk0,E , − 21

k∇w · ne k0,e ≤ Ct h

(2.6)

rk∇wk0,E , ,

∀w ∈ Pr (E)

(2.7)

Lemma 2.4 ([3, 7]) 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, [37]) 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.8) (2.9)

We also recall the following form of the discrete Gronwall’s lemma: Lemma 2.6 (discrete Gronwall) an , bn , cn , dP n (for integers n ≥ 0) be nonnegPl Let ∆t, H,Pand l ative numbers such that al + ∆t n=0 bn ≤ ∆t n=0 dn an + ∆t ln=0 cn + H for l ≥ 0. Suppose that ∆tdn < 1, ∀n. Then P P P l dn ∆t c + H for l ≥ 0. al + ∆t ln=0 bn ≤ exp ∆t ln=0 1−∆td n n=0 n

5

Analysis of a DG Method for the Chemotaxis Model

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 of 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).

3

Description of the Continuous in Time Numerical Scheme for the Keller-Segel Model

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

system (1.2). First, notice that the Jacobians of F and G are    0 χρ 0 χv 0 0 χρ      0 0 0  0 0 0 0 ∂G    =  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.2) 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, rmin

rmax := max{rρ , rc , ru , rv }, rmin := min{rρ , rc , ru , rv },

(3.2)

where a is a constant independent of rρ , rc , ru , and rv . 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 for the chemotaxis system (1.2): Z XZ XZ XZ DG ρ DG ρ DG ρ {∇w ρ · ne }[ρDG ] {∇ρ · ne }[w ] + ε ∇ρ ∇w − ρt w + Ω

+σρ

E∈Eh E X rρ2 Z

e∈Γh

−

XZ

E∈Eh E

|e|

e

e∈Γh e

e∈Γh e

DG

[ρ

ρ

][w ] −

XZ

DG DG

χρ

E∈Eh E

χρDG v DG (w ρ )y +

X Z

e∈Γhor h

e

u

ρ

(w )x +

X Z

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

e∈Γver h e

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

(3.3)

6

Y. Epshteyn and A. Izmirlioglu Z

c cDG t w

Ω

+σc Z

−

E∈Eh E X r2 Z c

uDG w u +

|e|

∇c

DG

c

∇w −

[cDG ][w c ] +

e XZ

Z Ω

cDG (w u )x +

X Z

c

DG

u

nx w + σu

X Z

XZ

cDG (w v )y +

DG

v

E∈Eh E

c

e∈Γh e

cDG w c − X Z

X Z

e∈Γh ∪∂Ωhor

|e|

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

DG

u

(·, 0)w =

Ω

Ω Z

u

u(·, 0)w ,

Ω

DG

Z

c

· ne }[w ] + ε

XZ

e∈Γh e

{∇w c · ne }[cDG ]

ρDG w c = 0,

(3.4)

Ω

(−cDG )∗u nx [w u ]

X

e∈Γhor e h 2 X rv

ny w + σv

e∈∂Ωhor e

Ω Z

{∇c

e∈Γh ∪∂Ωver

e∈∂Ωver e

v DG w v +

XZ

e∈Γver h e

E∈Eh E

−

Ω

+

e∈Γh

Ω

Z

XZ

ru2 |e|

Z

[uDG ][w u ] = 0,

(3.5)

e

(−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 ε corresponding 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 [30]: 1

2

in DG DG E aout (χρDG uDG )E u )e aout ain e − a (χρ (χρ u ) = − [ρDG ], 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 ain aout (cDG )E )e e − a (c DG ∗ − [uDG ], (−c )u = − aout − ain aout − ain 1 in DG E 2 bout (cDG )E )e bout bin e − b (c DG ∗ (−c )v = − − [v DG ]. out in out 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.2) is hyperbolic, these speeds can be estimated using the largest ∂F and ∂G (see (3.1)): and the smallest eigenvalues of the Jacobian ∂Q ∂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 DG E 2 in DG E 1 DG E 2 bout = max (χv DG )E , (χv ) , 0 , b = min (χv ) , (χv ) , 0 . e e e e

7

Analysis of a DG Method for the Chemotaxis Model 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 )u = − , (−c )v = − , 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)

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

4

Results for the Continuous in Time Scheme

Let us recall here the results that were obtained in [18] for the continuous in time scheme (3.3)(3.6). Lemma 4.1 (Consistency Lemma ) If the solution of the system (1.2) belongs to H 2 (Eh ), then it satisfies the formulation (3.3)–(3.6). Theorem 4.2 (L2 (H 1 ) and L∞ (L2 ) error estimates). Let the solution ρ, c, u and v of the Keller-Segel system (1.2) be sufficiently regular. Furthermore, we assume that penalty parameters σρ , σc , σu , σv are sufficiently large. Then there exists at least one DG solution to (3.3)-(3.6) and there exists constants Cρ and Cc , independent of h and r, such that

DG DG

ρ − ρ ∞ 2 ([0,T ];L2 (Ω)) + + |||∇(ρ − ρ)||| L 2 L ([0,T ];L (Ω)) ≤ Cρ (

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

rρρ

+

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

rρρ

where (rρ , rc , ru , rv ) ≥ 2.

+

T 0

X rρ2

1

[ρDG − ρ] 2 2 0,e |e| e∈Γ h

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

DG

DG

c − c ∞ 2 ([0,T ];L2 (Ω)) + + |||∇(c − c)||| L 2 L ([0,T ];L (Ω)) ≤ Cc (

Z

Z

0

T

X r2

2 12 c DG [c − c] 0,e |e| e∈Γ h

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

8

5

Y. Epshteyn and A. Izmirlioglu

Fully Discrete Schemes and Analysis

In the following sections, we formulate two explicit schemes and establish the convergence of the numerical solutions using an induction hypothesis. Existence of the discrete solution is trivial since the scheme is explicit in time. In the analysis below, we will assume that the exact solution of the system (1.2) 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 ≥ 5.5,

(5.2)

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 [24] and references therein. Let ∆t be a positive time step and let ti = i∆t denote the time step at the ith step. We denote by v i the function v evaluated at time ti .

6

Forward Euler Time Discretization

Find a discrete in time solution ρ i+1 i+1 i+1 c u v (ρi+1 DG , cDG , uDG , vDG ) ∈ Wrρ ,h × Wrc ,h × Wru ,h × Wrv ,h ,

which satisfies the following weak formulation for the chemotaxis system (1.2): Z Ω

Z

i XZ XZ XZ ρi+1 i ρ i ρ DG − ρDG ρ {∇w ρ · ne }[ρiDG ] {∇ρDG · ne }[w ] + ε ∇ρDG ∇w − w + ∆t e∈Γh e e∈Γh e E∈Eh E Z Z Z X X X rρ2 (χρiDG uiDG )∗ nx [w ρ] χρiDG uiDG (w ρ )x + [ρiDG ][w ρ] − +σρ |e| e∈Γver E∈Eh E e∈Γh h e e Z Z X X i i (χρiDG vDG )∗ ny [w ρ ] = 0, (6.3) χρiDG vDG (w ρ )y + − E∈Eh E i ci+1 DG − cDG

∆t

c

w +

XZ

E∈Eh E

Ω

e∈Γhor e h

∇ciDG ∇w c

−

XZ

e∈Γh e

{∇ciDG

· ne }[w ] + ε

Z Z X r2 Z c i c i c +σc [cDG ][w ] + cDG w − ρiDG w c = 0, |e| e∈Γh e Ω Ω Z X Z XZ ∗ u u u (−ci+1 ci+1 ui+1 DG )u nx [w ] DG (w )x + DG w + Ω

E∈Eh E

e∈Γver h e

c

XZ

e∈Γh e

{∇w c · ne }[ciDG ] (6.4)

9

Analysis of a DG Method for the Chemotaxis Model

− Z Ω

−

X Z

u ci+1 DG nx w

+ σu

i+1 v vDG w +

X Z

ru2 |e|

e∈Γh ∪∂Ωver

e∈∂Ωver e

XZ

X

v ci+1 DG (w )y +

E∈Eh E

v ci+1 DG ny w

X Z

e∈Γhor e h 2 X rv

+ σv

e∈Γh ∪∂Ωhor

e∈∂Ωhor e

|e|

Z

u [ui+1 DG ][w ] = 0,

(6.5)

e

∗ v (−ci+1 DG )v ny [w ]

Z

i+1 [vDG ][w v ] = 0.

(6.6)

e

To approximate the convective terms in (6.3), (6.5)-(6.6) we use the same central-upwind fluxes (3.8) as for the continuous in time scheme, with the one-sided local speeds given by: i E1 out i E2 i E1 in i E2 ai := max (χuDG )e , (χuDG)e , 0 , ai := min (χuDG )e , (χuDG )e , 0 , (6.7) i E1 i E2 in i E1 i E2 bout := max (χv ) , (χv ) , 0 , b := min (χv ) , (χv ) , 0 . i i DG e DG e DG e DG e

Notice that the inequalities similar to (3.10), aout i ≤ 1, aout − ain i i

−ain i ≤ 1, aout − ain i i

bout i ≤ 1, bout − bin i i

and

−bin i ≤ 1, bout − bin i i

(6.8)

which are needed in our convergence proof, are satisfied for the local speeds defined in (6.7) as in out in well (for simplicity, we assume that aout i −ai 6= 0 and bi −bi 6= 0 throughout the computational domain, see Remark in Section 3). Also, the initial conditions are: Z Z Z Z 0 ρ 0 ρ 0 c ρDG w = ρ w , cDG w = c0 w c , Ω Z

u0DG w u

Ω

=

Ω Z Ω

0

u

uw ,

ΩZ Ω

0 vDG wv

ΩZ

=

(6.9)

v0wv .

Ω

We denote by ρei , e ci , u ei, and vei the piecewise polynomial interpolants of the exact solution components ρi , ci , ui, and v i of the Keller-Segel system (1.2) and assume that these interpolants satisfy the approximation property (2.2). We then make the following induction hypothesis: assume that ∀ i : 0 ≤ i ≤ n, we have i ) ∈ Wruu ,h × Wrvv ,h : S = (uiDG , vDG n X i=0

2

∆tkuiDG − u eik0,Ω

≤ Cu

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

rρ

+Cut ∆t2 , n X 2 i ∆tkvDG − vei k0,Ω i=0

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

(6.10)

10

Y. Epshteyn and A. Izmirlioglu

≤ Cv

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

rρ

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

+Cvt ∆t2 , sup kuiDG − u ei k0,Ω ≤ Cu⋆

0≤i≤n

h 2 rmin

i sup kvDG − vei k0,Ω ≤ Cv⋆

,

0≤i≤n

h 2 rmin

(6.11)

,

(6.12)

where Cu , Cv ,, Cut , Cvt , Cu⋆ , and Cv⋆ are positive constants (which will be defined later) independent of h, the polynomial degrees (rρ , rc , ru , rv ), and n. The parameters sρ , sc , su , and sv denote the regularity of the corresponding components of the exact solution. Clearly, the induction hypothesis above holds true for i = 0. We need now to show that assuming that S holds true ∀ i : 0 ≤ i ≤ n, it will follow that it will be true for i + 1 = n + 1. i Let us first show that from the induction hypothesis S, it follows that functions uiDG , and vDG ,1 ≤ i ≤ n are bounded. i Lemma 6.1 For (uiDG , vDG ) ∈ S, there exist positive constants Mu , and Mv independent of h, ru , and rv , such that i sup kuiDGk∞,Ω ≤ Mu , sup kvDG k∞,Ω ≤ Mv . (6.13) 0≤i≤n

0≤i≤n

Proof: The result is derived from the definition of the subset S and the inverse inequality . Theorem 6.2 (l2 (H 1) and l∞ (L2 ) Forward Euler error estimates). Let the solution ρ, c, u and v of the Keller-Segel system (1.2) be sufficiently regular. Furthermore, we assume that penalty parameters σρ , σc , σu , σv are sufficiently large. Then the induction hypothesis holds true for n + 1. Furthermore, there exists constants Cρ and Cc , independent of h and r, such that 1 2

kρDG − ρkl∞ ([0,T ];L2(Ω)) + ∆t |||∇(ρDG − ρ)|||l2 ([0,T ];L2(Ω)) +

≤ Cρ (

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

+

rρ

hmin(rc +1,sc )−1 sc − 23

+

rc

su − 32

sρ − 32

rρ

+

hmin(rc +1,sc )−1 sc − 32

+

ru

kcDG − ckl∞ ([0,T ];L2 (Ω)) + ∆t |||∇(cDG − c)|||l2([0,T ];L2 (Ω)) + hmin(rρ +1,sρ )−1

i=0

hmin(ru +1,su )−1

1 2

≤ Cc (

N X

+

rc

su − 32

X rρ2

1

[ρi − ρi ] 2 2 DG 0,e |e| e∈Γ h

hmin(rv +1,sv )−1 sv − 23

+ ∆t)

rv N X i=0

hmin(ru +1,su )−1 ru

∆t

+

X r2

21 c i i 2 [cDG − c ] 0,e ∆t |e| e∈Γ h

hmin(rv +1,sv )−1 sv − 23

+ ∆t),

rv

where (rρ , rc , ru , rv ) ≥ 2. Proof: We introduce the following notation: τρi := ρiDG − ρei , ξρi := ρi − ρei , τui := uiDG − u ei , ξui = ui − u ei ,

τci := ciDG − e ci , ξci := ci − e ci , i v i. τvi := vDG −e v i , ξvi := v i − e

(6.14)

11

Analysis of a DG Method for the Chemotaxis Model

It follows from the consistency Lemma 4.1(see [18] for the details) that the exact solution of (1.2) satisfies the following weak formulation: 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 Ω Z Z Z Z X X X X (χρv)∗∗ ny [w ρ ] = 0,(6.15) χρv(w ρ )y + (χρu)∗∗ nx [w ρ] − χρu(w ρ)x + − e∈Γver h e

E∈Eh E

E∈Eh E

e∈Γhor e h

where 1

(χρu)∗∗ :=

2

E in E in aout aout i (χρu)e − ai (χρu)e i ai − [ρ], aout − ain aout − ain i i i i 1

(χρv)

∗∗

2

in E in bout (χρv)E bout e − bi (χρv)e i bi := i − [ρ], bout − bin bout − bin i i i i

in out in and the local speeds aout i , ai , bi , and bi are given by (6.7). Using (6.14), equation (6.15) can be rewritten as: Z X rρ2 Z XZ XZ XZ ρ i i ρ i ρ i ρ {∇w · ne }[e ρ ] + σρ {∇e ρ · ne }[w ] + ε ∇e ρ ∇w − ρet (t )w + [e ρi ][w ρ ] |e| e∈Γh e∈Γh e e∈Γh e E∈Eh E e Ω X Z XZ XZ XZ (χρi ui )∗∗ nx [w ρ ] χe ρi ξui (w ρ )x + χe ρi τui (w ρ )x − χe ρi uiDG (w ρ)x + −

−

XZ

i χe ρi vDG (w ρ )y

E∈Eh E

=−

Z

i

ρ

ξρ t (t )w −

Ω

− σρ

X rρ2 |e| e∈Γ h

Z

+

XZ

E∈Eh E

XZ

χe ρi τvi (w ρ )x

E∈Eh E

∇ξρi ∇w ρ

[ξρi ][w ρ ] +

XZ

+

XZ

e∈Γh e

−

XZ

χe ρi ξvi (w ρ )x

E∈Eh E

{∇ξρi

χξρi ui(w ρ )x +

+

X Z

(χρi v i )∗∗ ny [w ρ ]

e∈Γhor e h

ρ

· ne }[w ] − ε

XZ

XZ

e∈Γh e

{∇w ρ · ne }[ξρi ]

χξρi v i (w ρ)y .

(6.16)

E∈Eh E

E∈Eh E

e

e∈Γver h e

E∈Eh E

E∈Eh E

E∈Eh E

2

In order to obtain the estimate on the time step ∆t = O( hr4 ) we will proceed as follows, subtract equation (6.16) from (6.3) and choose w ρ = τρi , to obtain X rρ2 2

i+1 2

τ − τ i 2 + 2∆t|||∇τ i |||2 + 2∆tσρ

[τ i ] = τ i+1 − τ i 2 ρ 0,Ω ρ ρ ρ ρ ρ 0,Ω 0,Ω 0,Ω 0,e |e| e∈Γh XZ XZ XZ XZ i i i i i i i i χe ρi ξui (τρi )x χe ρ τu (τρ )x −2∆t {∇τρ ·ne }[τρ ]+2∆t χτρ uDG (τρ )x +2∆t +(1−ε)2∆t e∈Γh

−2∆t

X Z

e∈Γver h

e

e

E∈Eh E

(χρiDG uiDG )∗ − (χρi ui )∗∗

−2∆t

XZ

E∈Eh

E

χe ρi ξvi (τρi )y

nx [τρi ] + 2∆t

− 2∆t

XZ

E∈Eh E

X Z

e∈Γhor h

e

E∈Eh

E

i χτρi vDG (τρi )y

+ 2∆t

XZ

E∈Eh

i (χρiDG vDG )∗ − (χρi v i )∗∗ ny [τρi ]

E

E∈Eh

E

χe ρi τvi (τρi )y

12

Y. Epshteyn and A. Izmirlioglu

+2∆t

Z

i

Ω

+ε2∆t

ξρ t (t

)τρi + 2∆t

Z XZ XZ ρei+1 − ρei i i i i {∇ξρi · ne }[τρi ] ∇ξρ ∇τρ − 2∆t ρet (t ) − τρ + 2∆t ∆t Ω e∈Γ E∈E h

XZ

e∈Γh e

{∇τρi ·ne }[ξρi ]+2∆tσρ

X

e∈Γh

rρ2

|e|

=: R

Z

XZ

[ξρi ][τρi ]−2∆t

E∈Eh E

e

T1ρ

+

h

E

T2ρ

+ ... +

χξρi ui (τρi )x −2∆t

XZ

χξρi v i(τρi )y

E∈Eh E

(6.17)

ρ . T18

τρi+1 −τρi

τρi on the LHS of (6.17) was rewritten as:

i+1 2

i 2

i+1

Z i+1 i 2 τρ − τρi i τρ 0,Ω τρ 0,Ω τρ − τρ 0,Ω τρ = − − ∆t 2∆t 2∆t 2∆t Ω

Let us remark that term

e

Ω

∆t

(6.18)

Next, we bound each term on the RHS of (6.17) starting from term T2ρ using standard DG techniques, the term T1ρ will be bounded last. The quantities εi in the estimates below are positive real numbers, which will be defined later. Consider now, second term on the RHS T2ρ : X

{∇τ i } k[τ i ]k . |T ρ| ≤ 4∆t 2

ρ

ρ

0,e

0,e

e∈Γh

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

X 1 X

i i i E1 i E2 k{∇τρ }k0,e k[τρ ]k0,e ≤ 4∆t ∆t k[τρi ]k0,e

(∇τρ )e + (∇τρ )e 2 0,e 0,e e∈Γh e∈Γh X 2∆tCt rρ k∇τρi k0,E 1 + k∇τρi k0,E 2 k[τρi ]k0,e , ≤ √ h e∈Γver ∪Γhor h

and hence, using the fact that |e| ≤ the following bound on T2ρ :

√

h

h and using Cauchy-Schwarz inequality, we end up with

|T2ρ | ≤ ερ2 ∆t|||∇τρi |||20,Ω + R2 ∆t The term inequality:

T3ρ

X rρ2 2

[τ i ] ρ 0,e |e| e∈Γ

(6.19)

h

is bounded using Lemma 6.1 for uiDG and using Cauchy-Schwarz and Young’s

2 |T3ρ | ≤ ερ3 ∆t|||∇τρi |||20,Ω + R3 ∆t τρi 0,Ω

(6.20)

The term T4ρ is bounded using Cauchy-Schwarz and Young’s inequalities and assumption (3.2): X

2 K ∇τρi 0,E τui 0,E ≤ ερ4 ∆t|||∇τρi |||20,Ω + R4 ∆t τui 0,Ω |T4ρ | ≤ ∆t (6.21) E∈Eh

The term (2.2):

T5ρ

is bounded using Cauchy-Schwarz, Young’s inequality and approximation results |T5ρ | ≤ ερ5 ∆t|||∇τρi |||20,Ω + R5 ∆t

h2 min(ru +1,su ) ru2su

(6.22)

Analysis of a DG Method for the Chemotaxis Model Next, we bound T6ρ on the RHS of (6.17) as X Z aout ρ i i i E1 i i E1 i |T6 | ≤ 2∆t (χρ u ) − (χρ u ) n [τ ] x DG DG e e ρ aout − ain i i e∈Γver h e Z −ain i i i E2 i i E2 i + (χρDG uDG )e − (χρ u )e nx [τρ ] out in ai − ai e Z out −ain i ai i i i [ρ − ρ ]nx [τρ ] =: I + II + III. + out in DG ai − ai

13

(6.23)

e

Using (6.8) and (6.14), the first term on the RHS of (6.23) can be estimated by X Z 1 1 i i E i i E i (χρ u ) − (χρ u ) n [τ ] I ≤ 2∆tχ x DG DG e e ρ e∈Γver h

e Z X Z 1 1 (τ i ui )E nx [τ i ] + (ξ i ui )E nx [τ i ] ≤ 2∆tχ ρ DG e ρ DG e ρ ρ e∈Γver h

Z e Ze 1 1 i (ξ i ρi )E nx [τ i ] =: eI. + + (ρi τui )E n [τ ] x u e ρ e ρ e

e

We now use the Cauchy-Schwarz inequality, the trace inequality (2.4), the inequality (2.6), the assumption (3.2), the assumption that h < 1, r > 1, ∆t < 1, the approximation inequality (2.2), the bound on uiDG from Lemma 6.1: X rρ2 2

eI ≤ εe1 ∆t τ i 2 +(K1 +K2 )∆t

[τ i ] +K3 ∆tkτ i k2 +K4 ∆t ρ 0,Ω u 0,Ω ρ 0,e |e| e∈Γ h

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

rρ

h2 min(ru +1,su ) + ru2su

A similar bound can be derived for the second term II on the RHS of (6.23). To estimate the last term on the RHS of (6.23), we first use (6.14) and the definition of the one-sided local speeds (6.7) to obtain Z X 2 i i i f

[τρ ] 0,e + [ξρ ][τρ ] := III. III ≤ C∆t e∈Γver h

e

Then, using the Cauchy-Schwarz inequality, the assumption that h < 1 (and choosing h small f enough), r > 1, ∆t < 1, the approximation inequality (2.2), and inequality (2.4), we bound III as follows: 2 min(rρ +1,sρ )−1 X ∆trρ2 2 f

[τρi ] + ∆t h III ≤ . 2s 0,e |e| rρ ρ e∈Γver h

Combining the above bounds on I, II, and III, we arrive at |T6ρ |

≤

2 ερ6 ∆t τρi 0,Ω +R6 ∆t

2 min(rρ +1,sρ )−1 2 min(ru +1,su ) X rρ2 2 h h 2 i i

[τ ] +R7 ∆tkτ k +R8 ∆t . + u 0,Ω ρ 0,e 2s |e| ru2su rρ ρ e∈Γ h

(6.24)

.

14

Y. Epshteyn and A. Izmirlioglu

ρ The terms T7ρ - T10 are bounded in the same way as the terms T3ρ - T6ρ , respectively, and the bounds are:

2 (6.25) |T7ρ | ≤ ερ7 ∆t|||∇τρi |||20,Ω + R9 ∆t τρi 0,Ω

2 |T8ρ | ≤ ερ8 ∆t|||∇τρi |||20,Ω + R10 ∆t τvi 0,Ω (6.26)

|T9ρ | ≤ ερ9 ∆t|||∇τρi |||20,Ω + R11 ∆t

ρ | |T10

≤

h2 min(rv +1,sv ) rv2sv

X rρ2 2

[τρi ] + R13 ∆tkτvi k2 + R12 ∆t 0,Ω 0,e |e| e∈Γh 2 min(rρ +1,sρ ) h h2 min(rv +1,sv ) +R14 ∆t . + 2s rv2sv rρ ρ

(6.27)

2 ερ10 ∆t τρi 0,Ω

(6.28)

ρ we obtain the following bound: For the term T11

2 h2 min(rρ +1,sρ ) ρ |T11 | ≤ ερ11 ∆t τρi 0,Ω + R15 ∆t 2s rρ ρ

(6.29)

ρ Next, we bound term T12 . First using a Taylor expansion with integral remainder we can write

1 = ρe + ∆t∂t ρe(t ) + 2

i+1

ρe

1 2

Hence,

Z

i

ti+1

i

Z

ti+1

ti

1 (s − t )∂tt ρe(s)ds ≤ 2 i

ti

(s − ti )∂tt ρe(s)ds, and we can bound

Z

ti+1

ti

i 2

(s − t )

21 Z

ti+1

ti

21 (∂tt ρe(s))2 ds

Z ti+1 Z ti+1 i+1 i 2

1 ρ e − ρ e i 2 i

k∂tt ρe(s)k20,Ω ds, (s − t ) ≤ 4∆t ρet (t ) −

∆t ∆t ti ti 0,Ω

ρ Hence, using the above estimates we obtain the following bound for T12 : ρ |T12 |

≤

2 ερ12 ∆t τρi 0,Ω

2

+ R16 ∆t

Z

ti+1 ti

k∂tt ρe(s)k20,Ω ds.

(6.30)

ρ Next, we bound T13 using Cauchy-Schwarz and the approximation inequality (2.2): ρ |T13 | ≤ 2∆t

X

∇τ i ρ

E∈Eh

ρ Consider now term T14 :

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

i

∇ξ ≤ ερ ∆t|||∇τ i |||2 + R17 ∆t h 13 ρ 0,Ω 2s −2 0,E 0,E rρ ρ ρ |T14 |

≤ 2∆t

X

e∈Γh

||{∇ξρi

·

X Z i i ≤ 2∆t {∇ξρ · ne }[τρ ] e∈Γh

ne }||0,e ||[τρi ]||0,e

e

X |e| 21 rρ ||{∇ξρi · ne }||0,e 1 ||[τρi ]||0,e = 2∆t rρ |e| 2 e∈Γ h

(6.31)

Analysis of a DG Method for the Chemotaxis Model

15

Next, using Cauchy-Schwarz and trace inequality (2.5) we obtain: ≤ ερ14 ∆t

X rρ2 2 X

[τρi ] + C∆t |e| (h−1 ||∇ξρi ||20,E + h||∇2 ξρi ||20,E ) 2 0,e |e| r ρ E∈E e∈Γ h

h

Finally, using the approximation result (2.2) we derive: ρ | ≤ ερ14 ∆t |T14 ρ Next we bound term T15 :

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

[τρi ] + R18 ∆t h 2s −2 0,e |e| rρ ρ e∈Γ

(6.32)

h

ρ |T15 |

X Z ≤ ε2∆t {∇τρi · ne }[ξρi ] e∈Γh

e

Now using Lemma 2.3 and approximation result (2.3) we obtain: ≤ ε2∆t

X

E∈Eh

1

||∇τρi ||0,E

− 21

· Ch

· rρ

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

rρ

Next, applying Cauchy-Schwarz inequality we get the final estimate: ρ |T15 |

≤

ερ15 ∆t|||∇τρi |||20,Ω

+ R19 ∆t

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

rρ

(6.33)

ρ Now, we bound term T16 : ρ |T16 |

X rρ X rρ2 Z rρ i i ≤ 2∆tσρ [ξρi ][τρi ] ≤ 2∆tσρ 1 ||[τρ ]||0,e ||[ξρ ]||0,e 1 |e| e 2 |e| 2 e∈Γh |e| e∈Γh

ρ by using approximation result (2.3): Hence, we obtain the following bound for T16 ρ | ≤ ερ16 ∆t |T16

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

[τρi ] + R20 ∆t h 2s −3 0,e |e| rρ ρ e∈Γ

(6.34)

h

ρ ρ are bounded respectively using the approximation inequality (2.2) : and T18 Terms T17 ρ |T17 | ≤ ερ17 ∆t||τρi ||20,Ω + R21 ∆t

h2 min(rρ +1,sρ )

ρ |T18 | ≤ ερ18 ∆t|||∇τρi |||20,Ω + R22 ∆t

2sρ

rρ

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

rρ

(6.35)

(6.36)

The last term that needs to be bound is the term T1ρ . To bound this term we again subtract (6.16) from (6.3) and choose w ρ = τρi+1 − τρi . Using similar techniques as for the estimation of

16

Y. Epshteyn and A. Izmirlioglu

(6.17), except now use inverse inequality to estimate ∇(τρi+1 −τρi ) and inequality (2.6) to estimate ||[τρi+1 − τρi ]||0,e , we obtain the following estimate:

i+1

∆t2 rρ4 X rρ2 i 2 ∆t2 rρ4 i 2 k k i 2

[τ ]

τ |||∇τ ||| + C ≤ C − τ ρ 0,e ρ 0,Ω 2 1 ρ ρ 0,Ω h2 h2 e∈Γ |e| h

∆t2 rρ4 i 2 ∆t2 rρ4 i 2 ∆t2 rρ4 i 2

τρ + C4k

τu + C5k

τv 0,Ω 0,Ω 0,Ω h2 h2 h2 Z ti+1 ∆t2 rρ4 h2 min(rρ +1,sρ )−2 h2 min(ru +1,su ) h2 min(rv +1,sv ) k 3 k k∂tt ρe(s)k20,Ω ds + + C7 ∆t + +C6 2sρ −3 2su 2sv h2 r r i rρ t u v (6.37) Some details of the derivation of the estimate (6.37) can be found in Appendix 9. Next, we combine estimates (6.19)-(6.36) together with (6.37), plug it in (6.17), and using again the assumptions h < 1 and r > 1 to obtain: +C3k

i+1 2 ∆trρ4 ∆trρ4 X rρ2 i 2 i 2

τ − τ i 2 + ∆t(2 − C1

[τ ] )|||∇τρ |||0,Ω + ∆t(2σρ − C2 2 ) ρ 0,Ω ρ ρ 0,e 0,Ω h2 h |e| e∈Γ h

∆t2 rρ4 i 2 ∆t2 rρ4 i 2 ∆t2 rρ4 i 2

τρ + C4

τu + C5

τv 0,Ω 0,Ω 0,Ω h2 h2 h2 Z ti+1 ∆t2 rρ4 h2 min(rρ +1,sρ )−2 h2 min(ru +1,su ) h2 min(rv +1,sv ) 2 + + C7 ∆t k∂tt ρe(s)k20,Ω ds + +C6 2 2sρ −3 2su 2sv h r r rρ ti u v (6.38) Now sum (6.38) for i = 0, ..., n: ≤ C3

n n X rρ2 2

n+1 2 ∆trρ4 X ∆trρ4 X i 2

τρ − τρ0 2 + (2 − C1

[τρi ] ) ∆t|||∇τρ |||0,Ω + (2σρ − C2 2 ) ∆t 2 0,Ω 0,Ω 0,e h h |e| i=0 i=0 e∈Γ h

≤ C3

n n n

i 2

i 2

2 ∆trρ4 X ∆trρ4 X ∆trρ4 X

τ + C4

τ + C5 ∆t ∆t ∆t τvi 0,Ω ρ 0,Ω u 0,Ω 2 2 2 h i=0 h i=0 h i=0

2 min(rρ +1,sρ )−2 n Z ti+1 n X ∆trρ4 X h h2 min(ru +1,su ) h2 min(rv +1,sv ) 2 + +C7 ∆t ∆t k∂tt ρe(s)k20,Ω ds + +C6 2 2sρ −3 2su 2sv h i=0 r r rρ u v i=0 ti (6.39) Pn Pn i 2 i 2 Now using induction hypothesis S for i=0 ∆t kτu k0,Ω and i=0 ∆t kτv k0,Ω , (6.10)-(6.11), we obtain: n n X rρ2 2

n+1 2 ∆trρ4 X ∆trρ4 X i 2

[τ i ]

τ + (2 − C1 ) ∆t|||∇τ ||| + (2σ − C ) ∆t ρ 2 ρ 0,e ρ 0,Ω ρ 0,Ω h2 i=0 h2 i=0 |e| e∈Γ h

n

2

2 ∆trρ4 X ≤ C3 2 ∆t τρi 0,Ω + τρ0 0,Ω h i=0

17

Analysis of a DG Method for the Chemotaxis Model

∆trρ4 h2 min(rρ +1,sρ )−2 h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 + + + 2s −3 h2 rc2sc −3 ru2su −3 rv2sv −3 rρ ρ n Z ti+1 X ∆trρ4 2 2 t t k∂tt ρe(s)k20,Ω ds (6.40) +(C4 · Cu + C5 · Cv ) 2 ∆t + C7 ∆t h t i i=0

+(C4 ·Cu +C5 ·Cv +C6⋆ )

Let us now define Cmax := max(C3 , C4 · Cu + C5 · Cv + C6⋆ , C4 · Cut + C5 · Cvt ), where constant Cmax depends on the properties of the exact solutions ρ, u, v , T and domain Ω (it follows from the standard DG techniques and induction hypothesis for i = n). Hence, we obtain: n n X rρ2 2

n+1 2 ∆trρ4 X ∆trρ4 X i 2

τ + (2 − C1

[τ i ] ) ∆t|||∇τρ |||0,Ω + (2σρ − C2 2 ) ∆t ρ ρ 0,e 2 0,Ω h h |e| i=0 i=0 e∈Γ h

n

2

2 ∆trρ4 X ∆t τρi 0,Ω + τρ0 0,Ω ≤ Cmax 2 h i=0

∆trρ4 h2 min(rρ +1,sρ )−2 h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 + + + 2s −3 h2 rc2sc −3 ru2su −3 rv2sv −3 rρ ρ n Z ti+1 X ∆trρ4 2 2 +Cmax 2 ∆t + C7 ∆t (6.41) k∂tt ρe(s)k20,Ω ds h i=0 ti 2 h2 h2 h Next choosing ∆t ≤ min Cmax r4 , C1 r4 , C2 r4 , setting penalty parameter σρ to be large ρ ρ ρ ∆trρ4 ∆tr 4 enough, defining Cmin := min 1, 2−C1 h2 , 2σρ −C2 h2 ρ = 1, which can be done by choosR 1+C7 0T k∂tt ρe(s)k20,Ω ds 1 , = ing penalty parameter σρ large enough, and setting Cmm := max Cmin Cmin RT 1 + C7 0 k∂tt ρe(s)k20,Ω ds. Note that, now, the constant Cmm is independent of the induction hypothesis. Hence, we derive: +Cmax

n X X rρ2 2

n+1 2 i 2

τ +

[τ i ] ) ∆t(|||∇τ ||| + ρ ρ 0,Ω ρ 0,e 0,Ω |e| i=0 e∈Γ h

≤ Cmm

n X i=0

h2 min(rρ +1,sρ )−2 h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2

2 + ∆t τρi 0,Ω +Cmm + + 2s −3 rc2sc −3 ru2su −3 rv2sv −3 rρ ρ

2 +Cmm τ 0 + Cmm ∆t2 , (6.42) ρ

0,Ω

where Cmm depends only on the properties of the exact solutions ρ, T and domain Ω. Now, applying discrete Gronwall’s lemma to (6.42), we obtain the final estimate for τρn+1 : n 2 min(rρ +1,sρ )−2 X X rρ2 2

n+1 2 i ρ h i 2

τρ + [τ ] ) ≤ C ∆t(|||∇τ ||| + 0,Ω 2s −3 0,e 0,Ω |e| rρ ρ i=0 e∈Γ h

+

h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 + + + Cρt ∆t2 , rc2sc −3 ru2su −3 rv2sv −3

(6.43)

18

Y. Epshteyn and A. Izmirlioglu

where constants C ρ , Cρt depend only on the properties of the exact solutions ρ, domain Ω, T and independent of n. We now proceed in a similar way to the derivation of (6.43) for τci+1 = τcn+1 and we obtain: n 2 min(rρ +1,sρ )−2 h2 min(rc +1,sc )−2 X X r 2 2

n+1 2 c i i 2 c h

τ + ∆t(|||∇τc |||0,Ω + [τc ] 0,e ) ≤ C + c 2s −3 0,Ω |e| rc2sc −3 rρ ρ i=0 Γ h

+

h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 + + Cct ∆t2 ru2su −3 rv2sv −3

(6.44)

For the some details of the derivation of the above estimate (6.44) see Appendix 9. Next, we proceed by proving that the induction hypothesis S for un+1 DG . Once again, by the consistency Lemma (4.1), the exact solution satisfies the following equation (compare it with (3.5)): Z XZ X Z X Z i+1 u i+1 ∗∗ u i+1 u e ci+1 nx w u u e(t )w + (−c )u nx [w ] − e c (w )x + Ω

+σu

X

E∈Eh E Z ru2

e∈Γh ∪∂Ωver

+

X Z

e∈∂Ωver e

|e|

e∈Γver h e

[e ui+1 ][w u ] = −

e

Z Ω

ξci+1 nx w u − σu

X

e∈Γh ∪∂Ωver

where

ξui+1 w u −

ru2 |e|

Z

ξci+1(w u )x

E∈Eh E

[ξui+1 ][w u ],

(6.45)

e

1

(−ci+1 )∗∗ u := −

e∈∂Ωver e

XZ

2

i+1 E i+1 E in aout |e − ain |e aout i c i c i ai − [ui+1 ]. in out in aout − a a − a i i i i

Subtracting equation (6.45) from (6.5) and choosing w u = τui+1 , we obtain X

ru2 2 k[τui+1 ]k0,Ω |e| e∈Γh ∪∂Ωver Z Z X Z X Z X i+1 i+1 i+1 ∗ i+1 ∗∗ i+1 i+1 i+1 τc nx [τu ] + ξui+1 τui+1 ((−cDG )u − (−c )u )nx [τu ] + τc (τu )x − =− 2

kτui+1 k0,Ω + σu

e∈Γver h e

E∈Eh E

+

XZ

E∈Eh E

ξci+1 (τui+1 )x −

X Z

ξci+1 nx [τui+1 ] + σu

X

ru2

e∈Γh ∪∂Ωver

e∈∂Ωver e

|e|

Z

e∈∂Ωver e

Ω

[ξui+1 ][τui+1 ] =: T1u + ... + T7u ,

e

and bound each term on the RHS of (6.46). Let us now bound term by term on the right-hand side of (6.46). Consider the term T1u . Using integration by parts, we get: Z Z XZ i+1 i+1 i+1 i+1 τci+1 τui+1 nx ), − τc (τu )x = −(− τu (τc )x + E

E

e∈∂E

e

where nx denotes the x− component of the outward normal vector to element E.

(6.46)

19

Analysis of a DG Method for the Chemotaxis Model Summing over all the elements E we have: X XZ XZ i+1 i+1 τci+1 τui+1 nx ) (τc )x τu + −(− E

E∈Eh

=

XZ

E∈Eh

e

E∈Eh e∈∂E

(τci+1 )x τui+1 E

−

X

e∈Γver h ∪∂Ω

Z

[τci+1 τui+1 ]nx

e

Recalling the formula for the jump and average (2.1) we get: X Z X Z X Z XZ i+1 i+1 i+1 i+1 i+1 i+1 τci+1 [τui+1 ]nx [τu ]{τc }nx − [τc ]{τu }nx − (τc )x τu − = E

E∈Eh

e∈Γver h

e

e∈Γver h

e

e∈∂Ωver

e

Hence, using the inequality (2.6), Cauchy-Schwarz’s and Young’s inequalities, using lemma 2.4 for kτci+1 k and applying assumption (3.2), we have the following bound for T1u

2

|T1u | ≤ εu1 τui+1 0,Ω + U1

X

e∈Γh ∪∂Ωver

X r2

ru2 c i+1 2

[τui+1 ] 2 + U2 (|||∇τci+1 |||20,Ω + [τ ] ) (6.47) c 0,e 0,e |e| |e| e∈Γ h

A bound for T2u can be obtained in a way similar to the bound on T8ρ : Z out in X Z a −a 1 2 1 2 i i i+1 i+1 E E u i+1 E i+1 i+1 E i+1 |T2 | ≤ aout − ain (cDG )e − (c )e nx [τu ] + aout − ain (cDG )e − (c )e nx [τu ] i i i i e∈Γh e e Z in out −a a i i + [ui+1 − ui+1 ]nx [τui+1 ] := I + II + III. (6.48) out in DG ai − ai e

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

e

e

Using then the Cauchy-Schwarz inequality, the trace inequality (2.4), the inequality (2.6), and the assumption (3.2), we estimate eI as follows: eI ≤ KK1 kτ i+1 k2 + KK2 c 0,Ω

X r2 h2 min(rc +1,sc ) 2 u . k[τui+1 ]k0,e + KK3 2sc |e| r c e∈Γ h

A similar bound can be derived for the second term on the RHS of (6.48). The third term on the RHS of (6.48) is similar to the third term on the RHS of (6.23), hence it can be bounded by h2 min(ru +1,su )−1 KK4 · h 1 X ru2 i+1 2 . + k[τ ]k + KK · h III ≤ 6 0,e ru2 2 e∈Γ |e| u ru2su h

By choosing h small enough, we obtain: III ≤

X r2 h2 min(ru +1,su )−1 2 u . k[τui+1 ]k0,e + |e| ru2su e∈Γ h

20

Y. Epshteyn and A. Izmirlioglu Combining the above bounds on I, II, and III, and using lemma 2.4 we arrive at 2

|T2u | ≤ U3 (|||∇τci+1 |||0,Ω + +U5

X r2 X

c i+1 2 [τc ] 0,e ) + U4 |e| e∈Γ e∈Γ ∪∂Ω h

h

h2 min(ru +1,su )−1 h2 min(rc +1,sc ) + ru2su rc2sc

ver

ru2 2 k[τui+1 ]k0,e |e|

.

(6.49)

To bound the term T3u , we use the Cauchy-Schwarz inequality, Young’s inequality, the inequality (2.6), and lemma 2.4 which yield 2

|T3u | ≤ U6 (|||∇τci+1 |||0,Ω +

X r2 X r2

2 u c i+1 2 [τc ] 0,e ) + U7 k[τui+1 ]k0,e . |e| |e| e∈Γ e∈∂Ω

(6.50)

ver

h

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

|T4u | ≤ εu2 kτui+1 k0,Ω + U8

h2 min(ru +1,su ) . ru2su

(6.51)

Using similar idea as for the term T1u , we can rewrite T5u in the form X Z X Z XZ i+1 i+1 i+1 i+1 u [τui+1 ]{ξci+1}nx [ξc ]{τu }nx + (ξc )x τu + T5 = − E∈Eh

+

X Z

e∈∂Ωver

e

E

e∈Γver h

e

e∈Γver h

e

ξci+1 [τui+1 ]nx := Te5u .

(6.52)

We then use Lemma 2.1 and Lemma 2.3 to obtain X

2

f9 |Te5 | ≤ εu3 kτui+1 k0,Ω + U

e∈Γh ∪∂Ωver

ru2 h2 min(rc +1,sc )−2 ||[τui+1 ]||20,e + U9 . |e| rc2sc −3

(6.53)

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

X r2 h2 min(rc +1,sc ) 2 u . k[τui+1 ]k0,e + U11 2sc |e| r c e∈∂Ω

(6.54)

ver

The last term T7u is bounded using the approximation result (2.3). Hence, |T7u | ≤ U12

X

e∈Γh ∪∂Ωver

h2 min(ru +1,su )−2 ru2 2 . k[τui+1 ]k0,e + U13 |e| ru2su −3

(6.55)

After obtaining the estimates (6.47)–(6.55), we plug them into (6.46) and use the assumption h < 1, r > 1, assumption (3.2) , choose the penalty parameters large enough and an appropriate scaling to obtain

Analysis of a DG Method for the Chemotaxis Model

2

kτui+1 k0,Ω + ≤ Cu1 (|||∇τci+1 |||20,Ω +

X

e∈Γh ∪∂Ωver

21

ru2 2 k[τui+1 ]k0,e |e|

h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 X r2

c i+1 2 [τc ] 0,e ) + Cu2 + . |e| rc2sc −3 ru2su −3 e∈Γ

(6.56)

h

Next multiply both sides of (6.56) by ∆t and sum for i = −1, ..., n to get n X

i=−1

2 ∆tkτui+1 k0,Ω

+

n X

X

∆t

i=−1

e∈Γh ∪∂Ωver

≤ Cu3 (∆t|||∇τcn+1 |||20,Ω + ∆t +

n X i=0

Cu4 ∆t(|||∇τci |||20,Ω +

ru2 2 k[τui+1 ]k0,e |e|

X r2

c n+1 2 [τc ] 0,e ) |e| e∈Γ h

n h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 X r 2 2 X c i 5 Cu ∆t + . (6.57) [τc ] 0,e )+ |e| rc2sc −3 ru2su −3 i=−1 e∈Γ h

Next, using the inverse inequality, inequality (2.6) we obtain: n X

i=−1

2 ∆tkτui+1 k0,Ω

+

n X

X

∆t

e∈Γh ∪∂Ωver

i=−1

ru2 2 k[τui+1 ]k0,e |e|

4

∆trc4

τcn+1 2 + ∆trc τcn+1 2 ) 0,Ω 0,Ω 2 2 h h n n h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 X X X r 2 2 c i 8 7 i 2 [τc ] 0,e )+ + . (6.58) Cu ∆t + Cu ∆t(|||∇τc |||0,Ω + |e| rc2sc −3 ru2su −3 i=−1 i=0 e∈Γ

≤ Cu6 (

h

2

Final estimate is obtained by choosing ∆t ≤ C hr4 and by using the estimate (6.44): c

n X

i=−1

≤ Cu

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

rρ

+

2 ∆tkτui+1 k0,Ω

+

n X

∆t

X

e∈Γh ∪∂Ωver

i=−1

ru2 2 k[τui+1 ]k0,e |e|

h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 + + + Cut ∆t2 . (6.59) 2s −3 2s −3 2s −3 c u v rc ru rv

The estimate (6.59) proves the induction hypothesis (6.10) for i = n + 1: n+1 X i=0

2 ∆tkτui k0,Ω

+

≤ Cu

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

rρ

+

h2 min(rc +1,sc )−2 rc2sc −3

h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 + + Cut ∆t2 ru2su −3 rv2sv −3

22

Y. Epshteyn and A. Izmirlioglu

Induction hypothesis S, (6.12) can be shown using similar techniques as in (6.10), see details in Appendix 9. In the same way, we prove the induction hypothesis S for τv and show that: n+1 X i=0

≤ Cv

7

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

rρ

2 ∆tkτvi k0,Ω

n+1 X

+

i=0

∆t

X

e∈Γh ∪∂Ωhor

rv2 2 k[τvi ]k0,e |e|

h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 + + + Cvt ∆t2 . (6.60) + 2s −3 2s −3 2s −3 c u v rc ru rv

Runge-Kutta Second Order Time Discretization

Find a discrete in time solution ρ i+1 i+1 i+1 c u v (ρi+1 DG , cDG , uDG , vDG ) ∈ Wrρ ,h × Wrc ,h × Wru ,h × Wrv ,h ,

which satisfies the following weak formulation for the chemotaxis system (1.2): Z Z XZ XZ XZ i ρ i ρ i ρ i ρ {∇w ρ · ne }[ρiDG ] {∇ρDG · ne }[w ] + ε ∇ρDG ∇w − zDG,ρ w = ρDG w − ∆t Ω

Ω

+σρ

X

e∈Γh

− Z

XZ

E∈Eh E

i zDG,c wc =

Ω

rρ2

|e|

Z

[ρiDG ][w ρ ] −

e

XZ

X Z

e∈Γhor h

Ω

χρiDG uiDG (w ρ )x +

E∈Eh E

i χρiDG vDG (w ρ )y +

Z

e∈Γh e

E∈Eh E

e

X Z

e∈Γh e

(χρiDG uiDG )∗ nx [w ρ]

e∈Γver h e

i (χρiDG vDG )∗ ny [w ρ ] ,

(7.1)

XZ XZ XZ i c i c i c {∇w c · ne }[ciDG ] {∇cDG · ne }[w ] + ε ∇cDG ∇w − cDG w − ∆t e∈Γh e

e∈Γh e

E∈Eh E

Z Z X r2 Z c i c i c +σc [cDG ][w ] + cDG w − ρiDG w c , |e| e∈Γh e Ω Ω Z XZ 2 Z X Z X r u i i u i u i u (−zDG,c )∗u nx [w u ] zDG,c (w )x + zDG,u w + σu [zDG,u ][w ] = − |e| e∈Γh ∪∂Ωver E∈Eh E e∈Γhor e e Ω h Z X i zDG,c nx w u , −

Z

e∈∂Ωver e

i zDG,v wv

Ω

− Z Ω

+ σv

X Z

e∈∂Ωhor e

ρ ρi+1 DG w

X

1 = 2

Z

e∈Γh ∪∂Ωhor

rv2 |e|

i zDG,c ny w v , Ω

ρiDG w ρ

1 + 2

Z

i [zDG,v ][w v ]

(7.3)

XZ X Z i i v (−zDG,c )∗v ny [w v ] zDG,c (w )y + =− E∈Eh E

e

(7.2)

e∈Γhor e h

(7.4) Z

Ω

i zDG,ρ wρ

XZ XZ i i ρ {∇zDG,ρ · ne }[w ρ ] ∇zDG,ρ ∇w − − 0.5∆t E∈Eh E

e∈Γh e

23

Analysis of a DG Method for the Chemotaxis Model XZ

X rρ2 Z XZ i i i ρ {∇w · + σρ +ε χzDG,ρ zDG,u (w ρ )x [zDG,ρ ][w ] − |e| e∈Γh e∈Γh e E∈Eh E e X Z XZ X Z i i i i ρ i i ∗ ρ (χzDG,ρ zDG,v )∗ ny [w ρ] , χzDG,ρ zDG,v (w )y + (χzDG,ρ zDG,u ) nx [w ] − + Z

e∈Γver h

+ε Z Ω

Z

e

c ci+1 DG w

Ω

ρ

XZ

e∈Γh e

1 = 2

Z

i ne }[zDG,ρ ]

Ω

c

{∇w ·

ciDG w c

1 + 2

ru2

|e| X Z u ci+1 n w , − x DG e∈Γh ∪∂Ωver

e∈∂Ωver e

i+1 v vDG w

Ω

−

+ σv

X Z

X

e∈Γh ∪∂Ωhor

e∈∂Ωhor e

E∈Eh E

Ω

i zDG,c wc

rv2 |e|

v ci+1 , n w y DG

Z e

− 0.5

X r2 c + σc |e| e∈Γ

i ne }[zDG,c ]

X

u ui+1 DG w + σu

Z

h

XZ

Z

E∈Eh E

e∈Γhor h

i ∇zDG,c ∇w c

i [zDG,c ][w c ]

+

Z

e

−

XZ

e∈Γh e

i zDG,c wc

−

Ω

e

Z Ω

i {∇zDG,c · ne }[w c ]

i zDG,ρ wc ,

(7.6)

XZ X Z i+1 i+1 u ∗ u u [uDG ][w ] = − (−ci+1 cDG (w )x + DG )u nx [w ] e∈Γver h e

E∈Eh E

(7.7) Z

i+1 [vDG ][w v ]

XZ X Z ∗ v i+1 v (−ci+1 cDG (w )y + =− DG )v ny [w ] E∈Eh E

e

e∈Γhor e h

(7.8)

Again, we use (3.8) to approximate the convective terms in the scheme above, with the one-sided local speeds given by: i E1 i E2 in i E1 i E2 aout := max (χu ) , (χu ) , 0 , a := min (χu ) , (χu ) , 0 , i,u i,u DG e DG e DG e DG e i E1 i E2 in i E1 i E2 bout i,v := max (χvDG )e , (χvDG )e , 0 , bi,v := min (χvDG )e , (χvDG )e , 0 , (7.9) i E1 out i E2 i E1 in i E2 ai,zu := max (χzDG,u )e , (χzDG,u )e , 0 , ai,zu := min (χzDG,u )e , (χzDG,u )e , 0 , i E2 i E1 i E1 i E2 in := min (χz ) , (χz ) , 0 , , (χz ) , 0 , b := max (χz ) bout DG,v e DG,v e i,zv DG,v e DG,v e i,zv

Denote by m1 := u, zu , m2 := v, zv and notice that we make the same settings here as in the Remark in section 3. Also notice that the inequalities similar to (3.10), aout i,m1 ≤ 1, out ai,m1 − ain i,m1

−ain i,m1 ≤ 1, out ai,m1 − ain i,m1

bout i,m2 ≤ 1, out bi,m2 − bin i,m2

and

−bin i,m2 ≤ 1, out bi,m2 − bin i,m2

(7.10)

which are needed in our convergence proof, are satisfied for the local speeds defined in (7.9) in out in as well (for simplicity, we assume that aout i,m1 − ai,m1 6= 0 and bi,m2 − bi,m2 6= 0 throughout the computational domain) and the initial conditions: Z Z Z Z 0 ρ 0 ρ 0 c ρDG w = ρ w , cDG w = c0 w c , Ω Z Ω

(7.5)

u0DG w u

=

Ω Z Ω

0

u

uw ,

ΩZ Ω

0 vDG wv

ΩZ

=

Ω

v0wv .

(7.11)

24

Y. Epshteyn and A. Izmirlioglu Now, let us use a similar idea to [41] and define new variables zρ , zc , zu , zv : zρ (x, y, t) = ρ(x, y, t) + ρt (x, y, t)∆t

(7.12)

zc (x, y, t) = c(x, y, t) + ct (x, y, t)∆t

(7.13)

zu (x, y, t) = u(x, y, t) + ut (x, y, t)∆t

(7.14)

zv (x, y, t) = v(x, y, t) + vt (x, y, t)∆t

(7.15)

Next, using the regularity of the density ρ, the concentration c, and Taylor expansion, we obtain for ρ and c: ρ(x, y, t + ∆t) − ρ(x, y, t) − ρt (x, y, t)

∆t ∆t − ρt (x, y, t + ∆t) = O(∆t3 ) 2 2

c(x, y, t + ∆t) − c(x, y, t) − ct (x, y, t)

∆t ∆t − ct (x, y, t + ∆t) = O(∆t3 ) 2 2

Hence, we obtain: ρt (x, y, t+∆t) = −(χρ(x, y, t+∆t)u(x, y, t+∆t))x −(χρ(x, y, t+∆t)v(x, y, t+∆t))y +∆ρ(x, y, t+∆t) = (−χ(ρ(x, y, t) + ρt (x, y, t)∆t + O(∆t2 )))(u(x, y, t) + ut (x, y, t)∆t + O(∆t2 ))x +(−χ(ρ(x, y, t) + ρt (x, y, t)∆t + O(∆t2 )))(v(x, y, t) + vt (x, y, t)∆t + O(∆t2 ))y +∆(ρ(x, y, t) + ρt (x, y, t)∆t + O(∆t2 )) = −(χzρ zu )x − (χzρ zv )y + ∆zρ + O(∆t2 )

In the same way, it can be shown that

ct (x, y, t + ∆t) = ∆zc − zc + zρ + O(∆t2 ) From (7.12) and the Taylor expansion above for ρ, it follows that: zρi = ρi − ((χρi ui)x + (χρi v i )y − ∆ρi )∆t 1 1 ∆t ρi+1 = ρi + zρi − ((χzρi zui )x + (χzρi zvi )y − ∆zρi ) + O(∆t3 ), 2 2 2

(7.16)

Similarly, for c we get: zci = ci − (−∆ci + ci − ρi )∆t 1 1 ∆t ci+1 = ci + zci − (−∆zci + zci − zρi ) + O(∆t3 ), 2 2 2

(7.17)

Note that for u and v we will have: zu = u + ut ∆t = cx + (cx )t ∆t = (c + ct ∆t)x = (zc )x zv = (zc )y Hence, it follows that: zui = (zci )x ui+1 = (ci+1 )x ,

(7.18)

25

Analysis of a DG Method for the Chemotaxis Model zvi = (zci )y v i+1 = (ci+1 )y ,

(7.19)

Denoting by Err(x, y, i) = O(∆t3 ) and multiplying the above equation (7.16)- (7.19) by the test functions wρ , wc , wu and wv , integrating by parts, and using consistency of the DG scheme we obtain the scheme for zρi , ρi+1 , zci , ci+1 , zui , ui+1 and zvi , v i+1 : Z Z i ρ zρ w = ρi wρ + Aρ (ρi , ui, v i , w ρ)∆t, Ω Z Z Z ZΩ 1 1 i ρ i ρ i i i ρ ∆t i+1 ρ ρw + zρ w + Aρ (zρ , zu , zv , w ) + Err(x, y, i)w ρ (7.20) ρ w = 2 2 2 Ω Ω Ω Ω Z

zci w c

Z

ci w c + Ac (ci , ρi , w c )∆t, Ω Z ZΩ Z Z 1 1 i+1 c i c i c i i c ∆t c w = cw + zc w + Ac (zρ , zc , w ) + Err(x, y, i)w c 2 2 2 Ω Ω Ω Ω =

Z

zui w u

X

+ σu

Ω

e∈Γh ∪∂Ωver

Z

u

Z

zvi w v

i+1

X

u

w + σu

Ω

X

+ σv

e∈Γh ∪∂Ωhor

v

i+1

X

v

w + σv

Ω

Z

[zui ][w u ] = Au (zci , w u ), e

Z

ru2

e∈Γh ∪∂Ωver

Ω

Z

ru2 |e|

|e|

rv2 |e|

e∈Γh ∪∂Ωhor

Z

(7.21)

[ui+1 ][w u ] = Au (ci+1 , w u )

(7.22)

e

[zvi ][w v ] = Av (zci , w v ),

e

rv2 |e|

Z

[v i+1 ][w v ] = Av (ci+1 , w v ),

(7.23)

e

where we denoted by XZ XZ XZ ρ ρ {∇w ρ · ne }[x1 ] {∇x1 · ne }[w ] + ε ∇x1 ∇w − Aρ (x1 , x2 , x3 , w ) := − ρ

Z

e∈Γh e

E∈Eh E

XZ

X Z

X rρ2 χx1 x2 (w ρ )x + [x1 ][w ρ ] − |e| e∈Γh e∈Γver E∈Eh E h e Z Z X X (χx1 x3 )∗ ny [w ρ ] , χx1 x3 (w ρ )y + − +σρ

E∈Eh E

e∈Γhor h

e∈Γh e

(χx1 x2 )∗ nx [w ρ ]

e

(7.24)

e

XZ XZ XZ c c c {∇w c · ne }[x1 ] Ac (x1 , x2 , w ) := − {∇x1 · ne }[w ] + ε ∇x1 ∇w − +σc

X r2 c |e| e∈Γ h

Z e

E∈Eh E

[x1 ][w c ] +

Z Ω

e∈Γh e

x1 w c −

Z Ω

x2 w c ,

e∈Γh e

(7.25)

26

Y. Epshteyn and A. Izmirlioglu XZ X Z X Z u ∗ u u x1 nx w , (−x1 )u nx [w ] − x1 (w )x + Au (x1 , w ) := − u

E∈Eh E

e∈Γver h e

E∈Eh E

e∈Γver h

(7.26)

e∈∂Ωver e

XZ X Z X Z ∗ v v (−x1 )v ny [w ] − x1 (w )y + x1 ny w v , Av (x1 , w ) := − v

(7.27)

e∈∂Ωver e

e

Let us again introduce similar notations to (6.14)

and

τρi := ρiDG − ρei , ξρi := ρi − ρei , τui := uiDG − u ei , ξui = ui − u ei ,

i − zeρi , τziρ := zDG,ρ i τziu := zDG,u − zeui ,

ξzi ρ := zρi − zeρi , ξzi u = zui − zeui ,

τci := ciDG − e ci , ξci := ci − e ci , i τvi := vDG −e v i , ξvi := v i − e v i,

i − zeci , ξzi c := zci − zec i , τzic := zDG,c i τziv := zDG,v − zevi , ξzi v := zvi − zevi .

(7.28)

(7.29)

and subtract (7.20) from (7.1) and (7.5) respectively. We obtain the following error equations: Z Z i ρ τzρ w = τρi w ρ + Mρi (w ρ), (7.30) Ω Ω Z Z 1 1 1 τρi+1 w ρ = ( τρi + τziρ )w ρ + Nρi (w ρ) 2 2 2 Ω ZΩ 1 1 (7.31) = τρi w ρ + Mρi (w ρ ) + Nρi (w ρ ) 2 2 Ω where Mρi (w ρ)

=

Nρi (w ρ ) =

Z

ZΩ

i , w ρ) − Aρ (ρi , ui, v i , w ρ))∆t, (ξzi ρ − ξρi )w ρ + (Aρ (ρiDG , uiDG, vDG

(7.32)

(2ξρi+1 − ξρi − ξzi ρ − 2E(x, y, i))w ρ

Ω i i i +(Aρ (zDG,ρ , zDG,u , zDG,v , wρ)

− Aρ (zρi , zui , zvi , w ρ))∆t,

(7.33)

In the same way we obtain the error equations for concentration c, we obtain the following error equations: Z Z i c τzc w = τci w c + Mci (w c ), (7.34) Ω Ω Z Z 1 1 1 i+1 c τc w = ( τci + τzic )w c + Nci (w c ) 2 2 Ω Ω 2 Z 1 1 = τci w c + Mci (w c ) + Nci (w c ) (7.35) 2 2 Ω where Mci (w c )

=

Nci (w c ) =

Z

ZΩ Ω

(ξzi c − ξci )w c + (Ac (ciDG , ρiDG , w ρ) − Ac (ci , ρi , w c ))∆t,

(7.36)

i i (2ξci+1 − ξci − ξzi c − 2E(x, y, i))w c + (Ac (zDG,c , w c ) − Ac (zci , zρi , w c ))∆t, (7.37) , zDG,ρ

Analysis of a DG Method for the Chemotaxis Model The error equations for u and v are : Z τziu w u + σu Ω

Z

τui+1 w u

X

ru2 |e|

e∈Γh ∪∂Ωver

X

+ σu

Ω

Z

e

ru2 |e|

e∈Γh ∪∂Ωver

[τziu ][w u ] = Mui (w u ),

Z

[τui+1 ][w u ] = Nui (w u )

27

(7.38) (7.39)

e

where Mui (w u )

=

Z

Ω

Nui (w u ) =

Z

ξzi u w u

X

+ σu

e∈Γh ∪∂Ωver

ξui+1w u + σu Ω

X

ru2 |e|

Z

Ω

τziv w v

e

ru2 |e|

e∈Γh ∪∂Ωver

Z

Z

Z

e

X

+ σv

Ω

+ σv

u i+1 [ξui+1 ][w u ] + (Au (ci+1 , w u ))(7.41) DG , w ) − Au (c

rv2 |e|

e∈Γh ∪∂Ωhor

τvi+1 w v

i , w u ) − Au (zci , w u )), (7.40) [ξzi u ][w u ] + (Au (zDG,c

X

e∈Γh ∪∂Ωhor

Z

rv2 |e|

e

[τziv ][w v ] = Mvi (w v ),

Z

[τvi+1 ][w v ] = Nvi (w v )

(7.42) (7.43)

e

where Mvi (w v )

=

Z

Ω

Nui (w u ) =

Z

ξzi v w v

+ σv

X

e∈Γh ∪∂Ωhor

ξvi+1w v + σv

Ω

X

rv2 |e|

e∈Γh ∪∂Ωhor

Z

rv2 |e|

e

i , w v ) − Av (zci , w v )), (7.44) [ξzi v ][w v ] + (Av (zDG,c

Z

e

v i+1 [ξvi+1 ][w v ] + (Av (ci+1 , w v ))(7.45) DG , w ) − Av (c

Next, set w ρ = τρi in the equation of (7.30) and w ρ = τziρ in (7.31), add the two equations and after easy calculations, obtain the important equality that will be used to derive the error estimate for the density ρ:

2

i+1 2 i+1 i

τ − τ i 2 = (7.46)

τρ − τzρ + Mρi (τρi ) + Nρi (τziρ ) ρ 0,Ω ρ 0,Ω 0,Ω

Similarly, we obtain for the concentration c:

i+1 2

τc − τci 2 = τci+1 − τzi 2 + Mci (τci ) + Nci (τzi ) c c 0,Ω 0,Ω 0,Ω

As for Forward Euler Scheme, let us make the following induction hypothesis: i SR = (uiDG , vDG ) ∈ Wruu ,h × Wrvv ,h : n X i=0

2

∆tkuiDG − u ei k0,Ω

(7.47)

28

Y. Epshteyn and A. Izmirlioglu

≤ Cru

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

rρ

t +Cru ∆t4 , n X 2 i ∆tkvDG −e v i k0,Ω i=0

≤ Crv

(7.48)

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

rρ

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

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

t +Crv ∆t4 ,

(7.49) i ⋆ , sup kvDG −e v i k0,Ω ≤ Crv , 2 2 rmin rmin 0≤i≤n 2 2 ⋆ h i i ⋆ h ≤ Crρ 4 , sup kcDG − e c k0,Ω ≤ Crc 4 , rmin 0≤i≤n rmin

⋆ sup kuiDG − u ei k0,Ω ≤ Cru

0≤i≤n

sup kρiDG − ρei k0,Ω

0≤i≤n

h

h

(7.50) (7.51)

The induction hypothesis SR implies the following lemma.

i Lemma 7.1 For (ρiDG , ciDG , uiDG , vDG ) ∈ SR, there exist positive constants Mρ , Mc , Mu ,, Mv , Nρ , Nc , Nu , and Nv independent of h, rρ , rc , ru , and rv , such that

sup kρiDG k∞,Ω ≤ Mρ ,

sup kciDG k∞,Ω ≤ Mc ,

0≤i≤n

sup 0≤i≤n

i kzDG,ρ k∞,Ω

≤ Nρ ,

sup kuiDG k∞,Ω ≤ Mu ,

0≤i≤n

sup 0≤i≤n

i kzDG,c k∞,Ω

i sup kvDG k∞,Ω ≤ Mv ,

0≤i≤n

≤ Nc ,

sup 0≤i≤n

i kzDG,u k∞,Ω

0≤i≤n

≤ Nu ,

sup 0≤i≤n

(7.52) ≤ Nv .

i kzDG,v k∞,Ω

(7.53) Proof: For the details of the proof see Appendix 9.

Theorem 7.2 (l2 (H 1) and l∞ (L2 ) Runge-Kutta error estimates). Let the solution ρ, c, u and v of the Keller-Segel system (1.2) be sufficiently regular. Furthermore, we assume that penalty parameters σρ , σc , σu , σv are sufficiently large. Then the induction hypothesis holds true for n + 1. Furthermore, there exists constants Crρ and Crc , independent of h and r, such that 1 2

kρDG − ρkl∞ ([0,T ];L2(Ω)) + ∆t |||∇(ρDG − ρ)|||l2 ([0,T ];L2(Ω)) + ≤ Crρ (

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

+

hmin(rc +1,sc )−1 sc − 23

rρ

+

su − 32

kcDG − ckl∞ ([0,T ];L2 (Ω)) + ∆t |||∇(cDG − c)|||l2([0,T ];L2 (Ω)) + h

sρ − 32

min(rc +1,sc )−1

rρ

where (rρ , rc , ru , rv ) ≥ 2.

+

h

sc − 32

rc

+

N X i=0

min(ru +1,su )−1

+

h

su − 32

ru

+

X rρ2

1

[ρi − ρi ] 2 2 DG 0,e |e| e∈Γ h

hmin(rv +1,sv )−1 sv − 3 rv 2

ru

1 2

≤ Crc (

∆t

i=0

hmin(ru +1,su )−1

rc

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

N X

+ ∆t2 )

X r2

21 c i i 2 [cDG − c ] 0,e ∆t |e| e∈Γ h

min(rv +1,sv )−1

h

sv − 3 rv 2

+ ∆t2 ),

Analysis of a DG Method for the Chemotaxis Model

29

The proof of the above theorem is postponed to Appendix 9. Remark. The error estimates obtained in this paper are h-optimal and r-suboptimal (by 1/2), and the numerical tests reported in [18] confirm the theoretical error estimates.

8

Numerical Example

In this section, we consider the initial-boundary value problem for the Keller-Segel system and compare the solution obtained by the proposed Forward Euler DG method and by the RungeKutta DG method. We take the chemotactic sensitivity χ = 1 and the bell-shaped initial data ρ(x, y, 0) = 1300e−130(x

2 +y 2 )

c(x, y, 0) = 650e−65(x

,

2 +y 2 )

.

We consider the above problem in the square domain [− 21 , 12 ] × [− 21 , 21 ]. According to the results in [23], both components ρ and c of the solution are expected to blow up at the origin in finite time. In Figures 8.1–8.4, we plot the contours of the density ρ along the line L = [− 21 , 21 ] × 0, computed at different times before blow-up on uniform grid with h = 1/101 using quadratic polynomial approximation (Figures 8.1 and 8.2) and cubic polynomial approximation (Figures 8.3 and 8.4). The dashed line represents the solution obtained by the Forward Euler DG method and the solid line represents the solution obtained by the Runge-Kutta DG method. As one can see, the higher order time schemes are important for such type of problems with rapidly changing solution in time.

9

Appendix: Proof of Several Estimates

We collect in the Appendix details of the proofs of the several estimates.

9.1

Derivation of the Estimate (6.37)

To obtain the estimate (6.37), we again subtract (6.16) from (6.3) and choose w ρ = τρi+1 − τρi . XZ

i+1

2 i

τ ∇τρi ∇(τρi+1 − τρi ) − τρ 0,Ω = −∆t ρ E∈Eh

+∆t

XZ

e∈Γh

e

{∇τρi · ne }[τρi+1 − τρi ] − ε∆t +∆t

−∆t +∆t

XZ

E∈Eh

XZ

E∈Eh E

XZ

E∈Eh E

i χτρi vDG (τρi+1

−

e

e∈Γh

{∇(τρi+1 − τρi ) · ne }[τρi ] − ∆tσρ

χτρi uiDG (τρi+1 − τρi )x + ∆t

χe ρi ξui (τρi+1

E

XZ

E

−

τρi )x

τρi )y

XZ

E∈Eh

E

X rρ2 Z [τρi ][τρi+1 − τρi ] |e| e e∈Γ h

χe ρi τui (τρi+1 − τρi )x

X Z i i ∗ i i ∗∗ (χρDG uDG ) − (χρ u ) nx [τρi+1 − τρi ] − ∆t e∈Γver h

+ ∆t

XZ

E∈Eh

E

e

χe ρi τvi (τρi+1

−

τρi )y

− ∆t

XZ

E∈Eh

E

χe ρi ξvi (τρi+1 − τρi )y

30

Y. Epshteyn and A. Izmirlioglu Z X Z i i ∗ i i ∗∗ i+1 i (χρDG vDG ) − (χρ v ) ny [τρ − τρ ] + ∆t ξρt (ti )(τρi+1 − τρi ) −∆t e∈Γhor h

Ω

e

Z ρei+1 − ρei i+1 i +∆t ρet (t ) − (τρ − τρi ) ∆t Ω XZ XZ i i+1 i {∇ξρi · ne }[τρi+1 − τρi ] ∇ξρ ∇(τρ − τρ ) − ∆t +∆t e∈Γh e

E∈Eh E

+ε∆t

XZ

e∈Γh e

−∆t

{∇(τρi+1

XZ

E∈Eh E

−

τρi )

·

ne }[ξρi ]

X rρ2 Z + ∆tσρ [ξρi ][τρi+1 − τρi ] |e| e∈Γ

χξρi ui(τρi+1 − τρi )x − ∆t

h

XZ

E∈Eh E

(9.1)

e

χξρi v i (τρi+1 − τρi )y

ρ =: T T1ρ + T T2ρ + ... + T T20 .

Now let us bound each term on the RHS of (9.1). Using similar techniques as for the estimation of (6.17), except now we will use inverse inequality for the estimation of the ∇(τρi+1 − τρi ), and

inequality (2.6) to estimate [τρi+1 − τρi ] 0,e . Hence we obtain the following estimates: |T T1ρ | ≤

∆t2 rρ4 1

τρi+1 − τρi 2 + R1 |||∇τρi |||20,Ω 0,Ω 21 h2

∆t2 rρ4 1

τρi+1 − τρi 2 + R2 |||∇τρi |||20,Ω 0,Ω 21 h2

∆t2 rρ4 X rρ2 i 2 1 ρ 3 i+1 i 2

[τρ ] |T T3 | ≤ τ − τρ 0,Ω + R 0,e 21 ρ h2 e∈Γ |e| |T T2ρ | ≤

(9.2) (9.3)

(9.4)

h

|T T4ρ | ≤

∆t2 rρ4 X rρ2 i 2 1

τ i+1 − τ i 2 + R4

[τ ] ρ 0,Ω ρ 0,e 21 ρ h2 e∈Γ |e| h

∆t2 rρ4 i 2 1

τρi+1 − τρi 2 + R5

τρ 0,Ω 0,Ω 21 h2 2 4

1

τρi+1 − τρi 2 + R6 ∆t ru τui 2 |T T6ρ | ≤ 0,Ω 0,Ω 21 h2

∆t2 rρ4 h2 min(ru +1,su ) 1

τ i+1 − τ i 2 + R7 |T T7ρ | ≤ ρ 0,Ω 21 ρ h2 ru2su |T T5ρ | ≤

|T T8ρ| ≤

∆t2 rρ4 i 2 ∆t2 rρ4 i 2 1

τ i+1 − τ i 2 + R9

τ + R10

τ ρ 0,Ω ρ 0,Ω u 0,Ω 21 ρ h2 h2 ∆t2 rρ4 h2 min(rρ +1,sρ ) h2 min(ru +1,su ) 11 . + +R 2s h2 ru2su rρ ρ |T T9ρ|

(9.5)

∆t2 rρ4 i 2 1 12 i+1 i 2

τ τ − τρ 0,Ω + R ≤ ρ 0,Ω 21 ρ h2

(9.6) (9.7) (9.8)

(9.9) (9.10)

Analysis of a DG Method for the Chemotaxis Model 2 4

1

τ i+1 − τ i 2 + R13 ∆t rv τ i 2 ρ 0,Ω v 0,Ω 21 ρ h2

∆t2 rρ4 h2 min(rv +1,sv ) 1 ρ

τρi+1 − τρi 2 + R14 |≤ |T T11 0,Ω 21 h2 rv2sv

∆t2 rρ4 i 2 ∆t2 rρ4 i 2 1 ρ

τρi+1 − τρi 2 + R16

τρ + R17

τv |T T12 |≤ 0,Ω 0,Ω 0,Ω 21 h2 h2 ∆t2 rρ4 h2 min(rρ +1,sρ ) h2 min(rv +1,sv ) . +R18 2 + 2s h rv2sv rρ ρ 2 min(rρ +1,sρ )

1 ρ

τρi+1 − τρi 2 + R19 ∆t2 h |≤ |T T13 2s 0,Ω 21 rρ ρ Z ti+1

2 1 ρ 20 3 i+1 i

τ k∂tt ρe(s)k20,Ω ds. − τρ 0,Ω + R ∆t |T T14 | ≤ 21 ρ ti 2

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

τρi+1 − τρi 2 + R21 |≤ |T T15 2s −2 0,Ω 21 h2 rρ ρ

∆t2 rρ4 h2 min(rρ +1,sρ )−2 1 ρ 22 i+1 i 2

τ − τρ 0,Ω + R |T T16 | ≤ 2s −2 21 ρ h2 rρ ρ

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

τ i+1 − τ i 2 + R23 |T T17 |≤ ρ 0,Ω 2s −3 21 ρ h2 rρ ρ

∆t2 rρ4 h2 min(rρ +1,sρ )−2 1 ρ 24 i+1 i 2

|T T18 | ≤ τ − τρ 0,Ω + R 2s −3 21 ρ h2 rρ ρ

∆t2 rρ4 h2 min(rρ +1,sρ ) 1 ρ

τ i+1 − τ i 2 + R25 |T T19 |≤ ρ 0,Ω 2s 21 ρ h2 rρ ρ

∆t2 rρ4 h2 min(rρ +1,sρ ) 1 ρ

τρi+1 − τρi 2 + R26 |≤ |T T20 2s 0,Ω 21 h2 rρ ρ Now, combining all the bounds (9.2)-(9.21) and using the assumption that h < 1, r > estimate (6.37) follows. ρ |T T10 |≤

9.2

31

(9.11) (9.12)

(9.13) (9.14) (9.15) (9.16) (9.17) (9.18) (9.19) (9.20) (9.21) 1, the

Derivation of the Estimate for the Concentration (6.44)

First, from the consistency Lemma 4.1 we obtain that the exact solution of (1.2) satisfying the weak formulation of the form of the equation (3.4), which may be rewritten as Z X r2 Z XZ XZ XZ c c i i c i c i c {∇w · ne }[e c ] + σc {∇e c · ne }[w ] + ε ∇e c ∇w − e ct (t )w + [e ci ][w c ] |e| e∈Γh e∈Γh e e∈Γh e E∈Eh E e Ω Z Z Z XZ XZ + e ci w c − ρei w c = − ξc t (ti )w c − {∇ξci · ne }[w c ] ∇ξci ∇w c + Ω

−ε

XZ

e∈Γh e

Ω

E∈Eh E

Ω

c

{∇w ·

ne }[ξci ]

X r2 c − σc |e| e∈Γ h

Z e

[ξci ][w c ]

−

Z Ω

e∈Γh e

ξci w c

+

Z Ω

ξρi w c .

(9.22)

32

Y. Epshteyn and A. Izmirlioglu

We then subtract equation (9.22) from equation (6.4) and set w c = τci to obtain X r 2 2

i+1 2

i+1 c i

τc − τci 2 + 2∆t|||∇τci |||20,Ω + 2∆tσρ

τc − τci 2 [τ ] = c 0,Ω 0,Ω 0,e 0,Ω |e| e∈Γh Z Z Z Z

i 2 e ci+1 − e ci i i i i i i i

τc +2∆t ξc τc − 2∆t ξρ τc − 2∆t τρ τc − 2∆t τc 0,Ω + 2∆t e ct (ti ) − ∆t Ω Ω Ω Ω Z XZ XZ XZ i i i i i i {∇τc · ne }[τc ] + 2∆t ξc t (t )τc + 2∆t +2∆t(1 − ε) {∇ξci · ne }[τci ] ∇ξc ∇τc − 2∆t +2∆tε

XZ

e∈Γh e

e∈Γh e

Ω

{∇τci · ne }[ξci ] + 2∆tσc

Z

X r2 c |e| e∈Γ h

E∈Eh E

e∈Γh e

[ξci ][τci ] = T1c + Tc2 + ... + Tc12 .

(9.23)

e

Let us notice again that the term on the LHS of (9.23) is rewritten similar to the term on the LHS of (6.17): Z i+1 2 2 2 τc − τci i kτci+1 k0,Ω kτci k0,Ω kτci+1 − τci k0,Ω τc = − − (9.24) ∆t 2∆t 2∆t 2∆t Ω

The terms on the RHS of (9.23) are bounded using the same techniques as in (6.17). As before, we start with term T2c :

2 h2 min(rc +1,sc ) |T2c | ≤ εc2 ∆t τci 0,Ω + C2c ∆t (9.25) rc2sc

2 h2 min(rρ +1,sρ ) |T3c | ≤ εc3 ∆t τci 0,Ω + C3c ∆t 2s rρ ρ

2

2 |T4c | ≤ εc4 ∆t τci + C4c ∆t τρi 0,Ω

|T6c |

≤

0,Ω

2 |T5c | = 2∆t τci 0,Ω

2 εc6 ∆t τci 0,Ω

+

C6c ∆t2

Z

ti+1

ti

|T7c | ≤ εc7 ∆t|||∇τci |||20,Ω + C7c ∆t

k∂tt e c(s)k20,Ω ds.

rc2 X

[τci ] 2 0,e |e| e∈Γ

(9.26) (9.27) (9.28)

(9.29) (9.30)

h

2 h |T8c | ≤ εc8 ∆t τci 0,Ω + C8c ∆t

2 min(rc +1,sc )

|T9c | c |T10 |

≤

≤

εc9 ∆t|||∇τci |||20,Ω

εc10 ∆t

+

rc2sc

h2 min(rc +1,sc )−2 c C9 ∆t rc2sc −2

X r 2 2 h2 min(rc +1,sc )−2 c i c [τc ] 0,Ω + C10 ∆t |e| rc2sc −2 e∈Γ

(9.31) (9.32) (9.33)

h

c c |T11 | ≤ εc11 ∆t|||∇τci |||20,Ω + C11 ∆t

h2 min(rc +1,sc )−2 rc2sc −3

(9.34)

33

Analysis of a DG Method for the Chemotaxis Model

c |T12 | ≤ εc12 ∆t

X r 2 2 h2 min(rc +1,sc )−2 c i c [τc ] 0,e + C12 ∆t |e| rc2sc −3 e∈Γ

(9.35)

h

Next, we need to bound the term as it was done for the T1ρ . We subtract (9.22) from equation (6.4) and set w c = τci+1 − τci to obtain: Z Z Z Z

2

i+1 i i i+1 i i i+1 i i i+1 i

τc − τc = ∆t ξc (τc − τc ) − ∆t ξρ (τc − τc ) − ∆t τρ (τc − τc ) − ∆t τci (τci+1 − τci ) 0,Ω T1c

Ω

Ω

Ω

Ω

Z XZ XZ e ci+1 − e ci i+1 i i+1 i i i {∇τci · ne }[τci+1 − τci ] ∇τc ∇(τc − τc ) + ∆t (τc − τc ) − ∆t +∆t e ct (t ) − ∆t Ω e∈Γh e E∈Eh E Z Z Z 2 X r X c {∇(τci+1 − τci ) · ne }[τci ] − ∆tσc −ε∆t [τci ][τci+1 − τci ] + ∆t ξc t (ti )(τci+1 − τci ) |e| e e∈Γh e∈Γh e Ω Z Z XZ X X {∇(τci+1 − τci ) · ne }[ξci ] {∇ξci · ne }[τci+1 − τci ] + ∆tε ∇ξci ∇(τci+1 − τci ) − ∆t +∆t E∈Eh E

+∆tσc

X

e∈Γh

rc2

|e|

Z

e∈Γh e

e∈Γh e

[ξci ][τci+1 − τci ] = T T1c + T Tc2 + ... + T Tc14 .

(9.36)

e

Using similar techniques as in (9.1), we obtain the following estimates: |T T1c | ≤

2 min(rc +1,sc )

1

τ i+1 − τ i 2 + ∆t2 KC 1 h c 0,Ω 15 c rc2sc

2 min(rρ +1,sρ )

1

τci+1 − τci 2 + ∆t2 KC 2 h 2s 0,Ω 15 rρ ρ

1

τci+1 − τci 2 + ∆t2 KC 3 τρi 2 |T T3c| ≤ 0,Ω 0,Ω 15

1

τci+1 − τci 2 + ∆t2 KC 4 τci 2 |T T4c | ≤ 0,Ω 0,Ω 15 Z ti+1

1

τci+1 − τci 2 + KC 5 ∆t3 k∂tt e c(s)k20,Ω ds. |T T5c | ≤ 0,Ω 15 ti 2 4

1

τ i+1 − τ i 2 + KC 6 ∆t rc |||∇τ i |||2 |T T6c | ≤ c 0,Ω c 0,Ω 15 c h2 2 4

1

τ i+1 − τ i 2 + KC 7 ∆t rc |||∇τ i |||2 |T T7c | ≤ c 0,Ω c 0,Ω 15 c h2 2 4 X 2

rc i 1 2

τci+1 − τci 2 + KC 8 ∆t rc [τc ] 0,e |T T8c | ≤ 0,Ω 15 h2 e∈Γ |e|

|T T2c | ≤

(9.37) (9.38) (9.39) (9.40) (9.41) (9.42) (9.43) (9.44)

h

|T T9c | ≤

2 4 X 2

rc i 1 2

τci+1 − τci 2 + KC 9 ∆t rc [τc ] 0,e 0,Ω 2 15 h e∈Γ |e|

c |T T10 |≤

(9.45)

h

1

τci+1 − τci 2 + ∆t2 KC 10 h 0,Ω 15

2 min(rc +1,sc )

rc2sc

(9.46)

34

Y. Epshteyn and A. Izmirlioglu 2 4 2 min(rc +1,sc )−2

1

τ i+1 − τ i 2 + KC 11 ∆t rc h c 0,Ω 15 c h2 rc2sc −2 2 4 2 min(rc +1,sc )−2

1 c

τ i+1 − τ i 2 + KC 12 ∆t rc h |T T12 |≤ c 0,Ω 15 c h2 rc2sc −2 2 4 2 min(rc +1,sc )−2

1 c

τ i+1 − τ i 2 + KC 13 ∆t rc h |T T13 |≤ c 0,Ω 15 c h2 rc2sc −3 2 4 2 min(rc +1,sc )−2

1 i+1 i 2 14 ∆t rc h c

τ − τc 0,Ω + KC |T T14 | ≤ 15 c h2 rc2sc −3 c |T T11 |≤

(9.47) (9.48) (9.49) (9.50) 2

We combine all the bounds (9.25)-(9.35) and use bounds (9.37)-(9.50) to estimate kτci+1 − τci k0,Ω 2 (similar to (6.37)). Then we plug the estimates in (9.23). Choosing ∆t ≤ C hr4 , summing c

for i = 0, .., n, applying Lemma 2.4 for τρi 0,Ω , using estimate (6.43) and applying discrete Gronwall’s lemma we obtain the final estimate for the τci+1 : n 2 min(rρ +1,sρ )−2 h2 min(rc +1,sc )−2 X X r 2 2

n+1 2 c i i 2 c h

τc + ∆t(|||∇τ ||| + [τ ] ) ≤ C + c 0,Ω c 0,e 2sρ −3 0,Ω |e| rc2sc −3 r ρ i=0 e∈Γ h

+

9.3

h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 et 2 + + Cc ∆t ru2su −3 rv2sv −3

(9.51)

Proof of the Induction Hypothesis (6.12) for un+1 DG

Let us consider (6.46) and bound each term on the RHS of (6.46). Consider the term T1u . Using inverse inequality, Cauchy-Schwarz, assumption (3.2) and Young’s inequality, we obtain:

2

2 r4 (9.52) |T1u | ≤ εu1 τui+1 0,Ω + UU1 c2 τci+1 0,Ω h A bound for T2u can be obtained in a similar way to the bound on T8ρ : Z out in X Z a −a 1 1 2 2 i i u i+1 i+1 E i+1 E i+1 E i+1 E i+1 + |T2 | ≤ (c (c ) − (c ) n [τ ] ) − (c ) n [τ ] x x e e u e e u DG DG aout − ain aout − ain e∈Γver h

i

i

e Z in out −a a i i i+1 i+1 i+1 := I + II + III. [u − u ]n [τ ] + x u DG in aout − a i i

e

i

i

e

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

e

e

Then, using the Cauchy-Schwarz inequality, the trace inequality (2.4), the inequality (2.6), and the assumption (3.2), we estimate eI as follows: eI ≤ KU1 kτ i+1 k2 + KU2 c 0,Ω

X r2 h2 min(rc +1,sc ) 2 u . k[τui+1 ]k0,e + KU3 |e| rc2sc e∈Γ h

(9.53)

Analysis of a DG Method for the Chemotaxis Model

35

A similar bound can be derived for the second term on the RHS of (6.46). The third term on the RHS of (6.46) is similar to the third term on the RHS of (6.23), hence it can be bounded by KU1 · h 1 X ru2 h2 min(ru +1,su )−1 i+1 2 III ≤ . + k[τ ]k0,e + KU3 · h ru2 2 e∈Γ |e| u ru2su h

By choosing h small enough, we obtain: X r2 h2 min(ru +1,su )−1 2 u . k[τui+1 ]k0,e + III ≤ 2su |e| r u e∈Γ h

Combining the above bounds on I, II, and III, and using lemma 2.4 we arrive at 2 min(ru +1,su )−1 2 2 min(rc +1,sc ) X

i+1 2 r h h 2 u u i+1 |T2 | ≤ UU2 τc 0,Ω + UU3 + . k[τ ]k0,e + UU4 |e| u ru2su rc2sc e∈Γ ∪∂Ω h

ver

(9.54) To bound the term we use the Cauchy-Schwarz inequality, Young’s inequality and the inequality (2.6) which yield X r2

2

2 u k[τui+1 ]k0,e . (9.55) |T3u | ≤ UU5 τci+1 0,Ω + UU6 |e| e∈∂Ω T3u ,

ver

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

i+1 |T4u | ≤ εuu 2 kτu k0,Ω + UU7

h2 min(ru +1,su ) . ru2su

(9.56)

Using similar techniques as for (6.52) to obtain 2

i+1 g |T5u | ≤ εuu 3 kτu k0,Ω + UU8

X

e∈Γh ∪∂Ωver

ru2 h2 min(rc +1,sc )−2 . ||[τui+1 ]||20,e + UU8 |e| rc2sc −3

(9.57)

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

T7u

X r2 h2 min(rc +1,sc ) 2 u . k[τui+1 ]k0,e + UU10 |e| rc2sc e∈∂Ω

(9.58)

ver

is bounded using the approximation result (2.3). Hence,

|T7u |

≤ UU11

X

e∈Γh ∪∂Ωver

ru2 h2 min(ru +1,su )−2 i+1 2 . k[τ ]k0,e + UU12 |e| u ru2su −3

(9.59)

After obtaining the estimates (9.52)–(9.59), we plug them into (6.46) and use the assumption h < 1, r > 1, assumption (3.2) , choose the penalty parameters large enough, choose an appropriate scaling and use the estimate (6.44) for the kτcn+1 k to obtain 2 min(rρ +1,sρ )−4 h2 min(rc +1,sc )−4 h2 min(ru +1,su )−4 h2 min(rv +1,sv )−4 4

n+1 2 ft rc ∆t2 . fu h

τu ≤ C + + + + C u 2s −7 0,Ω rc2sc −7 ru2su −7 rv2sv −7 h2 rρ ρ (9.60)

36

Y. Epshteyn and A. Izmirlioglu

The estimate (9.60) proves the induction hypothesis on uDG (6.12) for n + 1 by making the 2 appropriate choice of ∆t = O( hr4 ): c

9.4

n+1

τ ≤ C ⋆ h u 2 u 0,Ω rmin

Proof of Lemma 7.1

(9.61)

Firstly, the proof of the statement (7.52) of Lemma 7.1 is the direct consequence of the induction hypothesis SR, and is the same as for the Lemma 6.1. Secondly, to prove the next statement (7.53) of Lemma 7.1, let us first consider (7.30), and set w ρ = τziρ , to obtain: Z

2

i (9.62) τρi τziρ + Mρi (τziρ ),

τzρ = 0,Ω

where

Mρi (w ρ)

=

Z

Ω

+ε

XZ

e∈Γh

−

e XZ

(ξzi ρ

−

ξρi )w ρ

XZ XZ i ρ {∇τρi · ne }[w ρ ] ∇τρ ∇w − − ∆t

{∇w ρ · ne }[τρi ] + σρ χτρi uiDG (wρi )x

E∈Eh E

−

Ω

E∈Eh E X rρ2 Z

e∈Γh

XZ

|e|

e∈Γh e

[τρi ][w ρ]

e

χe ρi τui (wρi )x

+

E

E∈Eh

XZ

E∈Eh

χe ρi ξui (wρi )x

E

XZ XZ X Z i i i i i i ∗ i i ∗∗ χe ρi τvi (wρi )y χτρ vDG (wρ )y − (χρDG uDG ) − (χρ u ) nx [wρ ] − + e

e∈Γver h

− −

XZ

χe ρi ξvi (wρi )y +

E

E∈Eh

XZ

E∈Eh E

e∈Γhor h

∇ξρi ∇w ρ

X rρ2 − σρ |e| e∈Γ h

=

T1ρ

+

X Z

T2ρ

Z

+

XZ

e∈Γh e

[ξρi ][w ρ ] ρ T19 .

E

i (χρiDG vDG )∗ − (χρi v i )∗∗ ny [wρi ]

{∇ξρi

· ne }[w ] − ε

XZ

χξρi ui(w ρ )x

E∈Eh E

e

+ ... +

+

e

E∈Eh

E∈Eh E

ρ

XZ

e∈Γh e

+

{∇w ρ · ne }[ξρi ]

XZ

E∈Eh E

χξρi v i (w ρ)y

(9.63)

Recalling the definitions (7.28), (7.29) of ξρi and ξzi ρ , we can then easily obtain the following bound, which will be used several times throughout the error analysis:

2

i+1

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

ξ − ξ i 2 + ≤ C ∆t (9.64) − ξ ξ

zρ ξρ ρ ρ ρ 0,Ω 2s 0,Ω rρ ρ

where positive constant Cρ depends on ∂t u and is independent of i, h and r (similar bound is valid for the concentration c which will be used in the derivation of the error estimate for c). ρ Using techniques similar to the estimation of the terms T T1ρ − T T20 , (9.2)-(9.21) and applying the already verified result (7.52), we obtain:

Analysis of a DG Method for the Chemotaxis Model

|T1ρ | ≤ ε1 kw ρ k20,Ω + ∆t2 R1 |T2ρ| |T3ρ|

≤

ε2 kw ρ k20,Ω

≤

ε3 kw ρ k20,Ω

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

rρ

∆t2 rρ8 i 2 + R2 4 τρ 0,Ω h

∆t2 rρ8 i 2 + R3 4 τρ 0,Ω h

∆t2 rρ8 i 2

τ ρ 0,Ω h4 ∆t2 rρ8 i 2 ρ ρ 2 |T5 | ≤ ε5 kw k0,Ω + R5 4 τρ 0,Ω h ∆t2 rρ4 2 |T6ρ| ≤ ε6 kw ρ k20,Ω + R6 2 τρi 0,Ω h 2 8 2 ∆t rρ h 2 |T7ρ | ≤ ε7 kw ρ k20,Ω + R7 4 4 τui 0,Ω h ru |T4ρ| ≤ ε4 kw ρ k20,Ω + R4

|T8ρ| ≤ ε8 kw ρ k20,Ω + R8

∆t2 rρ4 h2 min(ru +1,su ) h2 ru2su

∆t2 rρ4 i 2 ∆t2 rρ8 h2 i 2

τ

τρ + R10 0,Ω h2 h4 ru4 u 0,Ω ∆t2 rρ4 h2 min(rρ +1,sρ ) h2 min(ru +1,su ) . +R11 2 + 2s h ru2su rρ ρ

37

(9.65) (9.66) (9.67) (9.68) (9.69) (9.70) (9.71) (9.72)

|T9ρ | ≤ ε9 kw ρk20,Ω + R9

∆t2 rρ4 i 2 ≤ + R12 2 τρ 0,Ω h 2 8 2 ∆t rρ h ρ

τ i 2 |T11 | ≤ ε11 kw ρk20,Ω + R13 4 h rv4 v 0,Ω ρ | |T10

ρ |T12 |

≤

ε10 kw ρk20,Ω

ε12 kw ρk20,Ω

∆t2 rρ4 h2 min(rv +1,sv ) + R14 2 h rv2sv

∆t2 rρ8 h2 i 2 ∆t2 rρ4 i 2

τρ + R16

τ 0,Ω h2 h4 rv4 v 0,Ω ∆t2 rρ4 h2 min(rρ +1,sρ ) h2 min(rv +1,sv ) . + +R17 2 2s h rv2sv rρ ρ

(9.73) (9.74) (9.75) (9.76)

ρ |T13 | ≤ ε13 kw ρ k20,Ω + R15

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

(9.78)

∆t2 rρ4 h2 min(rρ +1,sρ )−2 + R19 2 2s −2 h rρ ρ

(9.79)

ρ |T14 | ≤ ε14 kw ρk20,Ω + R18 ρ |T15 |

≤

ε15 kw ρk20,Ω

(9.77)

38

Y. Epshteyn and A. Izmirlioglu ∆t2 rρ5 h2 min(rρ +1,sρ )−2 2s −2 h2 rρ ρ

(9.80)

∆t2 rρ5 h2 min(rρ +1,sρ )−2 + R20 2 2s −2 h rρ ρ

(9.81)

ρ |T16 | ≤ ε16 kw ρk20,Ω + R19 ρ |T17 |

≤

ε17 kw ρk20,Ω

ρ | ≤ ε18 kw ρ k20,Ω + R21 ∆t2 rρ2 |T18 ρ |T19 | ≤ ε19 kw ρ k20,Ω + R22 ∆t2 rρ2

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

(9.82)

2sρ −2

rρ

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

(9.83) 2s −2 rρ ρ Now combining the estimates (9.65)-(9.83), using the assumption that h < 1, r > 1, ∆t < 1, and plugging them into (9.63), we obtain the following estimate for Mρi (w ρ ): ∆t2 rρ8 i 2 ∆t2 rρ8 h2 i 2 ∆t2 rρ8 h2 i 2

τ + M3

τ

τ + M2 ρ 0,Ω h4 h4 ru4 u 0,Ω h4 rv4 v 0,Ω ∆t2 rρ5 h2 min(rρ +1,sρ )−2 h2 min(ru +1,su ) h2 min(rv +1,sv ) +M4 2 + (9.84) + 2s −2 h ru2su rv2sv rρ ρ

|Mρi (w ρ )| ≤ ε kw ρ k20,Ω + M1

2

Next, choosing w ρ = τziρ in (9.84), making ε small enough, taking ∆t ≤ C hr4 , and using the ρ induction hypothesis SR (7.50) and (7.51), we conclude that we obtain from (9.62)

2 4

i fρ h (9.85)

τzρ ≤ C 8 0,Ω rmin As with estimate (9.85), it can be shown that

i 2

τz c

0,Ω

4 fc h ≤C 8 rmin

(9.86)

Now, let us consider the equation of (7.38), it can be shown using the same techniques as in the derivation of (9.84) and (9.60): X

2

r 4 2 ru2

[τzi ] 2 |Mui (w u )| ≤ ε τziu 0,Ω + Mu1 c2 τzic 0,Ω + Mu2 u 0,e h |e| e∈Γh ∪∂Ωver 2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 h 3 + (9.87) +Mu rc2sc −3 ru2su −3 Considering w u = τziu , making ε small enough, taking penalty parameter σu large enough, and using result (9.86), we obtain from (7.38) that: 2

i 2 fu h

τz ≤ C u 0,Ω 4 rmin

In the same way, it can be shown that: 2

i 2

τz ≤ C fv h v 0,Ω 4 rmin

(9.88)

(9.89)

Applying these estimates (9.88) and (9.89), the inverse inequality, and the use of Lemma 2.1, concludes the proof of Lemma 7.1.

Analysis of a DG Method for the Chemotaxis Model

9.5

39

Proof of Theorem 7.2

First let us consider, (7.46), which can be rewritten in the following way:

X rρ2 X rρ2 2

i+1 2

i 2

[τ i ] ) + ∆t(|||∇τ i |||2 + σzρ

τ − τ i 2 + ∆t(|||∇τ i |||2 + σρ ] [τ

zρ ) zρ 0,Ω ρ 0,e ρ 0,Ω ρ 0,Ω ρ 0,Ω |e| |e| 0,e e∈Γh e∈Γh

2

fi (τρi ) + N e i (τzi ), = τρi+1 − τziρ + M (9.90) ρ 0,Ω

where

fi (τ i ) = M ρ ρ

e i (τ i ) = N ρ zρ

Z

Ω

Z

(ξzi ρ

Ω

−

ξρi )τρi

+

i (Aρ (ρiDG , uiDG, vDG , τρi )

i

i

− Aρ (ρ , u , v

i

, τρi )

−

(|||∇τρi |||20,Ω

X rρ2 2

[τ i ] ( + σρ ρ 0,e |e| e∈Γ h

i i i , zDG,u , zDG,v , τziρ ) − Aρ (zρi , zui , zvi , τziρ ) (2ξρi+1 − ξρi − ξzi ρ − 2E(x, y, i))τziρ + (Aρ (zDG,ρ

−(|||∇τziρ |||20,Ω + σzρ

X rρ2

i 2

[τzρ ] ))∆t, |e| 0,e e∈Γ

(

h

f and N. e To estimate these Let us estimate, the terms on the RHS of (9.90), starting with terms M terms, instead of techniques used to obtain estimates in (9.84), we will use techniques similar to ρ ρ ρ the ones used in estimating the terms T2ρ − T10 and T13 − T18 in (6.17) (including Lemma 6.1). fi (τ i ) : Thus, we will obtain the following bounds for the M ρ ρ

X rρ2 2

eM fi (τ i )| ≤ M f1 ∆t τ i 2 + εeM ∆t(|||∇τ i |||2 + C f2 ∆t τ i 2 + M

[τ i ] ) + M f3 ∆t τ i 2 |M ρ 0,Ω ρ ρ ρ 0,Ω u 0,Ω ρ 0,e v 0,Ω |e| e∈Γh 2 min(rρ +1,sρ )−2 h2 min(ru +1,su ) h2 min(rv +1,sv ) h f + + (9.93) +M4 ∆t 2s −3 ru2su rv2sv rρ ρ

eρi (τzi ) (applying now Lemma 7.1): Similarly, we can derive the following estimate N ρ

2 X rρ2

i 2 2 i i e3 ∆t τ i 2 + N e4 ∆t τ i 2 eN e1 ∆t e i (τ i )| ≤ N ) + N ||| + C + ε e ∆t(|||∇τ ] τ [τ |N

N zu 0,Ω zv 0, zρ 0,Ω zρ zρ ρ zρ 0,e 0,Ω |e| e∈Γh 2 min(rρ +1,sρ )−2 h h2 min(ru +1,su ) h2 min(rv +1,sv ) e +N5 ∆t + + O(∆t5 ) (9.94 + 2sρ −3 2s 2s u v ru rv rρ

Next, we need to bound the first term of the RHS (9.90). First, let us notice that by subtracting the equation (7.30) from the equation (7.31), and setting w ρ = τρi+1 − τziρ , we obtain that the first term on the RHS of (9.84) which can be expressed as :

2 1

i+1

(9.95)

τρ − τziρ = (Nρi (τρi+1 − τziρ ) − Mρi (τρi+1 − τziρ )) 2 0,Ω

2

Hence, in order to estimate τρi+1 − τziρ we need to estimate the RHS of (9.95). It can be 0,Ω

ρ shown, using ideas similar to the ideas used in the estimation of the terms T T1ρ −T T12 (9.2)-(9.13),

40

Y. Epshteyn and A. Izmirlioglu

ρ ρ and T T15 − T T20 (9.16)-(9.21), that the following estimate holds:

2

i+1

τρ − τziρ

0,Ω

+C2M N

≤ C1M N

X rρ2 2 ∆t2 rρ4 i 2

[τ i ] ) (|||∇τ ||| + ρ 0,e ρ 0,Ω h2 |e| e∈Γ h

∆t2 rρ4 i 2 ∆t2 rρ4 i 2 ∆t2 rρ4 i 2

τρ + C3M N

τu + C4M N

τv 0,Ω 0,Ω 0,Ω h2 h2 h2

X rρ2 ∆t2 rρ4

i 2 i 2 [τ ] (|||∇τ ||| + +C5M N

zρ ) zρ 0,Ω h2 |e| 0,e e∈Γ h

∆t2 rρ4 ∆t2 rρ4 i 2 ∆t2 rρ4 i 2

i 2 MN

τ + C M N

τ + C τ

7 8 zρ zu 0,Ω zv 0,Ω h2 h2 h2 0,Ω ∆t2 rρ4 h2 min(rρ +1,sρ )−2 h2 min(ru +1,su ) h2 min(rv +1,sv ) MN + + O(∆t5 ) + +C3 2sρ −3 2su 2sv h2 r r rρ u v +C6M N

(9.96)

At this point, we need to get estimates for τziu and τziv . Consider equation in (7.38). To obtain an estimate for τziu in terms of τzic , we apply similar techniques as in (6.46) to the terms T1u − T7u (6.47)-(6.55). Thus, we obtain the following estimate: 2

kτziu k0,Ω + ≤ Cz1u (|||∇τzic |||20,Ω + Similarly, for

τziv

e∈Γh ∪∂Ωver

ru2 2 k[τziu ]k0,e |e|

h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 X r2

c i 2 [τzc ] 0,e ) + Cz2u + . |e| rc2sc −3 ru2su −3 e∈Γ

(9.97)

h

we have:

2

kτziv k0,Ω + ≤ Cz1v (|||∇τzic |||20,Ω +

X

X

e∈Γh ∪∂Ωhor

rv2 2 k[τziv ]k0,e |e|

h2 min(rc +1,sc )−2 h2 min(rv +1,sv )−2 X r2

c i 2 . [τzc ] 0,e ) + Cz2v + 2sc −3 2sv −3 |e| r r c v e∈Γ

(9.98)

h

Next, let us plug the estimate (9.97) and (9.98) into (9.96), and after simple modifications we obtain:

2 X rρ2 2 ∆t2 rρ4

i+1

i 2 eM N

[τ i ] ) (|||∇τ ||| +

τρ − τziρ ≤ C 1 ρ 0,e ρ 0,Ω 0,Ω h2 |e| e∈Γ h

eM N +C 2

∆t2 rρ4 i 2 ∆t2 rρ4 i 2 ∆t2 rρ4 i 2 eM N eM N

τρ + C

τv τ + C 3 u 0,Ω 4 0,Ω 0,Ω h2 h2 h2

X rρ2 ∆t2 rρ4

i 2 i 2 MN e (|||∇τzρ |||0,Ω + C5

[τzρ ] ) h2 |e| 0,e e∈Γ h

∆t2 rρ4 MN e +C 6 ((|||∇τzic k||20,Ω 2 h

X r2

∆t2 rρ4

i 2 c i 2 MN e [τzc ] 0,e ) + C7 +

τzρ 0,Ω |e| h2 e∈Γ h

41

Analysis of a DG Method for the Chemotaxis Model

eM N +C 3

∆t2 rρ4 h2

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

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

+ O(∆t5 )

(9.99) Next, we combine estimate (9.94), (9.93) and (9.99), plug them into (9.90) and use the assumption that ∆t ≤ 1, h < 1, r > 1 to obtain (after some simplifications) the following:

X rρ2 X rρ2 2

i+1 2

i 2

[τρi ] ) + ∆t(|||∇τzi |||20,Ω + σzρ

τρ − τρi 2 + ∆t(|||∇τρi |||20,Ω + σρ

[τzρ ] ) ρ 0,e 0,Ω 0,Ω |e| |e| 0,e e∈Γ e∈Γ h

h

X rρ2 2 ∆t2 rρ4 i 2 ∆t2 rρ4 i 2

[τρi ] ) + C2r

τρ + C3r

τu ≤ + 0,Ω 2 2 0,Ω 0,e |e| h h e∈Γh

X rρ2 ∆t2 rρ4 2 ∆t2 rρ4

i 2 +C4r 2 τvi 0,Ω + C5r 2 (|||∇τziρ |||20,Ω +

[τzρ ] ) h h |e| 0,e e∈Γh

X r2

∆t2 rρ4 ∆t2 rρ4

2 c i 2 [τzc ] 0,e ) + C7r 2 τziρ +C6r 2 (|||∇τzic |||20,Ω + h |e| h 0,Ω e∈Γh ∆t2 rρ4 h2 min(rρ +1,sρ )−2 h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 r +C8 2 + + + 2s −3 h rc2sc −3 ru2su −3 rv2sv −3 rρ ρ ∆t2 rρ4 C1r 2 (|||∇τρi |||20,Ω h

+O(∆t5 )

(9.100)

Next, consider the error equations (7.47). Using similar techniques as for the estimates in the equation (7.46), we obtain the following estimate for τ c : X r2 X r 2 2

i+1 2

c i 2 c i 2 i

τ − τ i 2 + ∆t(|||∇τ i |||2 + σc ] [τ ] [τ ||| + σ ) + ∆t(|||∇τ z c c 0,e zc 0,e ) zc 0,Ω c c 0,Ω c 0,Ω 0,Ω |e| |e| e∈Γ e∈Γ h

h

X r 2 2

2 ∆t2 r 4 c i [τc ] 0,e ) + C2c ∆t2 τci 0,Ω ≤ C1c 2 c (|||∇τci k||20,Ω + h |e| e∈Γ h

2 X r2

c i 2 [τzc ] 0,e ) + C5c ∆t2 τziρ 0,Ω |e| e∈Γh

∆t2 rρ4 h2 min(rρ +1,sρ ) h2 min(rc +1,sc )−2 c 2 i 2 c + O(∆t5 ) (9.101) + +C6 ∆t τzc 0,Ω + C7 2 2sρ 2sc −3 h r rρ c

2 2

Now, choose ∆t ≤ C hr4 , apply lemma 2.4 to estimate τziρ , sum for i = 0, ..., n, and apply

2 +C3c ∆t2 τρi 0,Ω +

∆t2 rρ4 C4c 2 (|||∇τzic |||20,Ω h

+

c

0,Ω

discrete Gronwall’s lemma to (9.101) to obtain:

n n X X r 2 2 X X r2

n+1 2

c i c i 2 i 2 i 2

τc + ∆t(|||∇τ ||| + [τ ] [τ ] ) + ∆t(|||∇τ ||| + c 0,Ω c 0,e zc 0,e ) zc 0,Ω 0,Ω |e| |e| i=0 i=0 e∈Γ e∈Γ h

n

X X rρ2

i 2

i 2 c c

∆t(|||∇τziρ |||20,Ω + ∆t τρ 0,Ω + K2 ≤ K1

[τzρ ] ) |e| 0,e i=0 i=0 e∈Γh 2 min(rρ +1,sρ ) n X h2 min(rc +1,sc )−2 h c + +K3 ∆t 2sρ rc2sc −3 r ρ i=0 n X

h

42

Y. Epshteyn and A. Izmirlioglu +O(∆t4 )

(9.102)

Next, summing the equation (9.100) for i = 0, ...n, using the above estimate (9.102), con2 sidering ∆t ≤ C hr4 , choosing the penalty parameters large enough, and using the assumption ρ ∆t < 1, h < 1, r > 1, we obtain: n n

X X rρ2 2 X X rρ2

n+1 2

i 2 i 2 i i 2

τρ +

∆t(|||∇τ ||| + [τ ] ) + ∆t(|||∇τ [τ ||| + ]

zρ ) ρ 0,Ω ρ 0,e zρ 0,Ω 0,Ω 0,e |e| |e| i=0 i=0 e∈Γ e∈Γ h

≤ R1

+R4

n X i=0

n X i=0 4

2 ∆t τρi 0,Ω + R2

∆t

n X i=0

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

h

2 ∆t τui 0,Ω + R3

2 min(rc +1,sc )−2

+

2sρ −3

rρ

h

n X

h

rc2sc −3

i=0

2 ∆t τvi 0,Ω

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

+O(∆t )

(9.103)

Now, applying induction hypothesis SR (7.48) , (7.49) and using discrete Gronwall’s lemma we obtain: n n

X X rρ2 2 X rρ2 X

n+1 2

i 2 i 2 2 i i

τ +

∆t(|||∇τρ |||0,Ω + ∆t(|||∇τzρ |||0,Ω + [τρ ] 0,e ) +

[τzρ ] ) ρ 0,Ω |e| |e| 0,e i=0 i=0 e∈Γh e∈Γh 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 ≤C + + + 2s −3 rc2sc −3 ru2su −3 rv2sv −3 rρ ρ +O(∆t4 ) (9.104) The estimate (9.104) also confirms the induction hypothesis (7.51) for ρ for i + 1 = n + 1. Using the final estimate for τρn+1 (9.104) we derive the following bound for τcn+1 : n n X X r 2 2 X X r2

n+1 2

c i c i 2 i 2 i 2

τc + ∆t(|||∇τ ||| + [τ ] [τ ] ) + ∆t(|||∇τ ||| + ) c 0,Ω c z z 0,Ω c c 0,Ω 0,e 0,e |e| |e| i=0 i=0 e∈Γh e∈Γh 2 min(rρ +1,sρ )−2 2 min(rc +1,sc )−2 2 min(ru +1,su )−2 2 min(rv +1,sv )−2 h h h h + + + ≤C 2s −3 rc2sc −3 ru2su −3 rv2sv −3 rρ ρ +O(∆t4 ) (9.105)

Again, the above estimate (9.105) confirms the induction hypothesis (7.51) for c, for i + 1 = n. Next, to obtain the estimate for τun+1 we consider second error equation in (7.39). Employing similar techniques as in the case of the error equation (6.46) and multiplying by ∆t, it can be shown that : n n X X X ru2 2 2 k[τui+1 ]k0,e ∆tkτui+1 k0,Ω + C1u ∆t |e| i=−1 i=−1 e∈Γ ∪∂Ω h

≤

+C3u

n X i=0

∆t(|||∇τci |||20,Ω

∆tr 4 C2u ( 2 c h

n+1 2

τ + c 0,Ω

∆trc4 h2

ver

n+1 2

τ ) c 0,Ω

n h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 X X r 2 2 c i u ∆t [τc ] 0,e ) + C4 + . + |e| rc2sc −3 ru2su −3 i=0 e∈Γ h

(9.106)

Analysis of a DG Method for the Chemotaxis Model

43

2

Final estimate is obtained by choosing ∆t ≤ C hr4 and by using the estimate for τcn+1 (9.105): c

n X

2

i=−1

≤U

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

rρ

∆tkτui+1 k0,Ω +

n X

∆t

i=−1

X

e∈Γh ∪∂Ωver

ru2 2 k[τ i+1 ]k0,e |e| u

h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 + + + Ut ∆t4 . (9.107) + 2s −3 2s −3 2s −3 c u v rc ru rv

The estimate (9.107) proves the induction hypothesis (7.48) for i + 1 = n + 1: n+1 X i=0

2 ∆tkτui k0,Ω

+

≤U

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

rρ

h2 min(rc +1,sc )−2 + rc2sc −3

h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 + + Ut ∆t4 ru2su −3 rv2sv −3

Induction hypothesis SR (7.50) can be shown using similar techniques as in the proof of the induction hypothesis S (6.12). Similarly, we prove the induction hypothesis SR for τv and show that: n n X X X rv2 2 i+1 2 k[τvi+1 ]k0,e ∆tkτv k0,Ω + ∆t |e| i=0 i=0 e∈Γh ∪∂Ωhor

≤V

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

rρ

+

h2 min(rc +1,sc )−2 h2 min(ru +1,su )−2 h2 min(rv +1,sv )−2 + + + Vt ∆t4 . (9.108) rc2sc −3 ru2su −3 rv2sv −3

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.

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] D.N. Arnold, An interior penalty finite element method with discontinuous elements, SIAM J. Numer. Anal., 19 (1982), pp. 742–760. [4] 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. [5] 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. [6] J.T. Bonner, The cellular slime molds, Princeton University Press, Princeton, New Jersey, 2nd ed., 1967.

44

Y. Epshteyn and A. Izmirlioglu

[7] S. Brenner, Poincar´e-Friedrichs inequalities for piecewise H 1 functions, SIAM J. Numer. Anal., 41 (2003), pp. 306–324. [8] E.O. Budrene and H.C. Berg, Complex patterns formed by motile cells of escherichia coli, Nature, 349 (1991), pp. 630–633. [9] E.O. Budrene and H.C. Berg, Dynamics of formation of symmetrical patterns by chemotactic bacteria, Nature, 376 (1995), pp. 49–53. [10] A. Chertock and A. Kurganov, A positivity preserving central-upwind scheme for chemotaxis and haptotaxis models, Numer. Math. submitted. [11] S. Childress and J.K. Percus, Nonlinear aspects of chemotaxis, Math. Biosc., 56 (1981), pp. 217–237. [12] B. Cockburn, G. Kanschat, and D. Schotzau, A note on discontinuous Galerkin divergence-free solutions of the Navier-Stokes equations, Journal of Scientific Computing, 31 (2007), pp. 61–73. [13] 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. [14] B. Cockburn and C.-W. Shu., The local discontinuous Galerkin method for convectiondiffusion systems, SIAM J. Numer. Anal., 35 (1998), pp. 2440–2463. [15] 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. [16] J. Douglas and T. Dupont, Lecture Notes in Physics, vol. 58, Springer-Verlag, 1976, ch. Interior penalty procedures for elliptic and parabolic Galerkin methods. [17] Y. Epshteyn, Discontinuous Galerkin methods for the chemotaxis and haptotaxis models, Journal of Computational and Applied Mathematics, 224 (2008), pp. 168–181. [18] Y. Epshteyn and A. Kurganov, New Interior Penalty Discontinuous Galerkin Methods for the Keller-Segel Chemotaxis Model, Siam J. Numer. Anal., 47 (2008), pp. 386–408. also CNA report, http://www.math.cmu.edu/cna/pub2007.html. [19] Y. Epshteyn and B. Rivi` ere, On the solution of incompressible two-phase flow by a pversion discontinuous Galerkin method, Comm. Numer. Methods Engrg., 22 (2006), pp. 741– 751. [20] Y. Epshteyn and B. Rivi` ere, Convergence of High Order Methods for Miscible Displacement, International Journal of Numerical Analysis and Modeling, 5, Supp (2008), pp. 47–63. [21] F. Filbet, A finite volume scheme for the Patlak-Keller-Segel chemotaxis model, Numer. Math., 104 (2006), pp. 457–488.

Analysis of a DG Method for the Chemotaxis Model

45

[22] V. Girault, B. Rivi` ere, and M.F. Wheeler, A discontinuous Galerkin method with non-overlapping domain decomposition for the Stokes and Navier-Stokes problems, Math. Comp., 74 (2005), pp. 53–84. ´zquez, A blow-up mechanism for a chemotaxis model, [23] M.A. Herrero and J.J.L. Vela Ann. Scuola Normale Superiore, 24 (1997), pp. 633–683. [24] D. Horstmann, From 1970 until now: The Keller-Segel model in chemotaxis and its consequences I, Jahresber. DMV, 105 (2003), pp. 103–165. [25] D. Horstmann, From 1970 until now: The Keller-Segel model in chemotaxis and its consequences II, Jahresber. DMV, 106 (2004), pp. 51–69. [26] E.F. Keller and L.A. Segel, Initiation of slime mold aggregation viewed as an instability, J. Theor. Biol., 26 (1970), pp. 399–415. [27] E.F. Keller and L.A. Segel, Model for chemotaxis, J. Theor. Biol., 30 (1971), pp. 225– 234. [28] E.F. Keller and L.A. Segel, Traveling bands of chemotactic bacteria: A theoretical analysis, J. Theor. Biol., 30 (1971), pp. 235–248. [29] A. Kurganov and C.-T. Lin, On the reduction of numerical dissipation in central-upwind schemes, Commun. Comput. Phys., 2 (2007), pp. 141–163. [30] 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. [31] 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. [32] A. Marrocco, 2D simulation of chemotaxis bacteria aggregation, M2AN Math. Model. Numer. Anal., 37 (2003), pp. 617–630. [33] V. Nanjundiah, Chemotaxis, signal relaying and aggregation morphology, J. Theor. Biol., 42 (1973), pp. 63–105. [34] C.S. Patlak, Random walk with persistence and external bias, Bull. Math: Biophys., 15 (1953), pp. 311–338. [35] L.M. Prescott, J.P. Harley, and D.A. Klein, Microbiology, Wm. C. Brown Publishers, Chicago, London, 3rd ed., 1996. ere, M.F. Wheeler, and V. Girault., A priori error estimates for finite [36] B. Rivi` element methods based on discontinuous approximation spaces for elliptic problems, SIAM J. Numer. Anal., 39 (2001), pp. 902–931. [37] C. Schwab, p- and hp-finite element methods, Numerical Mathematics and Scientific Computation, (1998). The Clarendon Press, Oxford University Press, New York.

46

Y. Epshteyn and A. Izmirlioglu

[38] 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. [39] 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. [40] 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. [41] Q. Zhang and C-W. Shu, Error Estimates to Smooth Solutions of Runge–Kutta Discontinuous Galerkin Methods for Scalar Conservation Laws, SIAM J. Numer. Anal., 42 (2004), pp. 641–666.

47

Analysis of a DG Method for the Chemotaxis Model

3100

3500

3000

3000

2900 2500 2800 2000

2700

1500

2600

1000

2500 2400

500

2300 0 −0.5

0

0.5

−0.05

0

0.05

Figure 8.1: contours along the line L of the density ρ on 101 ×101 grid using quadratic polynomial approximation p = 2 at time t = 5 × 10−6 (left column) and zoom view of the same contours(right column). Forward Euler - dashed line, Runge - Kutta - solid line.

4

4

3.5

x 10

x 10

3.24 3

3.22 3.2

2.5

3.18

2

3.16 1.5 3.14 1

3.12

0.5

3.1

0 −0.5

0

0.5

3.08 −0.01

−0.005

0

0.005

Figure 8.2: contours along the line L of the density ρ on 101 ×101 grid using quadratic polynomial approximation p = 2 at later time t = 2×10−5 (left column) and zoom view of the same contours(right column). Forward Euler - dashed line, Runge - Kutta - solid line.

0.01

48

Y. Epshteyn and A. Izmirlioglu

3500

3100

3000

3000 2900

2500

2800 2000

2700

1500

2600

1000

2500 2400

500

2300 0 −0.5

0

0.5

−0.05

0

0.05

Figure 8.3: contours along the line L of the density ρ on 101 × 101 grid using cubic polynomial approximation p = 3 at time t = 5 × 10−6 (left column) and zoom view of the same contours(right column). Forward Euler - dashed line, Runge - Kutta - solid line.

4

4

3.5

x 10

x 10

3.24 3

3.22

2.5

3.2

2

3.18 3.16

1.5 3.14 1

3.12

0.5

3.1

0 −0.5

0

0.5

3.08 −0.01

−0.005

0

0.005

Figure 8.4: contours along the line L of the density ρ on 101 × 101 grid using cubic polynomial approximation p = 3 at later time t = 2×10−5 (left column) and zoom view of the same contours(right column). Forward Euler - dashed line, Runge - Kutta - solid line.

0.01

Recommend Documents

A Fully Discrete Evolving Surface Finite Element Method

A fully discrete stabilized finite element method for ... - Semantic Scholar

A NEW ERROR ANALYSIS FOR DISCONTINUOUS FINITE ELEMENT ...

ERROR ANALYSIS OF A MIXED FINITE ELEMENT ... - HKBU MATH

CONVERGENCE OF THE DISCONTINUOUS GALERKIN ... - CMU Math