An interior penalty discontinuous Galerkin finite element method for quasilinear parabolic problems

Finite Elements in Analysis and Design 95 (2015) 42–50

Contents lists available at ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

An interior penalty discontinuous Galerkin finite element method for quasilinear parabolic problems Ioannis Toulopoulos Johann Radon Institute for Computational and Applied Mathematics, Austrian Academy of Sciences, Linz Austria

art ic l e i nf o

a b s t r a c t

Article history: Received 27 June 2014 Received in revised form 7 November 2014 Accepted 15 November 2014

In this paper, an Interior Penalty Discontinuous Galerkin (IPDG) finite element method is analyzed for approximating quasilinear parabolic equations. The equations can be characterized as perturbed parabolic p-Laplacian equations. The fully discrete scheme is obtained by applying s-stage Diagonally Implicit Runge–Kutta (s-DIRK) methods for the time integration. The nonlinear systems of the algebraic equations appearing in s-DIRK cycles are solved by developing two low storage Picard iterative processes. A stability bound is shown for the semi-discrete IPDG solution in the broken ‖:‖DG;p
norm. Continuous in time a priori error estimates are proved in case of p 4 2, when linear approximation space is used. A numerical test is performed in order to compare the performance of the two Picard iterative processes. Also, the results presented in the theoretical analysis are confirmed by numerical examples. & 2014 Elsevier B.V. All rights reserved.

Keywords: Quasilinear parabolic problems Perturbed parabolic p-Laplace problem Interior penalty discontinuous Galerkin method Stability estimates A priori error estimates

1. Introduction In this paper, an Interior Penalty Discontinuous Galerkin (IPDG) method is studied for approximating solutions of quasilinear problems in Lp setting, which can be recognized as examples of the perturbed p-Laplace problem. The problems are described by nonlinear diffusion equations, where the diffusion coefficient has a standard p-exponent form, that is ðμþ j∇ujÞp
2 ; μ4 0, with most interesting case here μ¼1 [1]. Very often, these constitute the mathematical model in many practical applications, as in aerodynamic, non-Newtonian flows, plasticity and glaciology, see e.g. [2,3]. Over the last two decades, there has been an increasing interest on devising discontinuous Galerkin (DG) methods for the numerical solution of elliptic and parabolic problems. This interest comes from the advantages of the local approximation spaces without continuity requirements that DG methods offer. Finite element methods defined on discontinuous spaces with interior penalties for linear elliptic problems were first analyzed in [4,5]. These methods, for the construction of the penalty terms on the interfaces, use similar techniques as Nitsche's treatment of introducing penalty terms for imposing Dirichlet boundary conditions. These approaches are generalized by symmetric and non-symmetric IPDG methods, see [6–8], for a comprehensive analysis of IPDG methods for linear elliptic problems. Recently, DG methods have

E-mail address: [email protected] http://dx.doi.org/10.1016/j.finel.2014.11.001 0168-874X/& 2014 Elsevier B.V. All rights reserved.

been proposed and analyzed for applications to nonlinear elliptic problems formulated in W 1;2 ðΩÞ. For example, in [9], DG methods have been analyzed for second order elliptic and hyperbolic systems and in [10], DG symmetric/non-symmetric methods have been analyzed for non-Fickian diffusion problems. In [11], an incomplete IPDG is introduced for a class of second order monotone nonlinear elliptic problems and a priori error estimates are given under minimal regularity assumptions on the exact solution. We also refer [12], where an hp-DG method has been studied for monotone quasilinear elliptic problems. Using theory of monotone operators, the authors showed the uniqueness of the DG solution and derivea a priori error estimate in a mesh-dependent energy norm. In the literature, IPDG methods have been proposed for solving parabolic type problems, based on the theoretical results for steady problems. We refer, but not limited to, the following. In [13], the first analysis of a semi-discrete IPDG method was presented for linear problems and in [14], optimal error estimates for a semi-discrete symmetric IPDG method have obtained for nonlinear parabolic problems. We also mention [15,16], where error estimates are discussed for fully discrete IPDG methods, and furthermore, we refer [17] where three fully discrete IPDG methods are considered and analyzed. In contrast to the analysis of IPDG methods for elliptic problems with natural formulation in W 1;2 ðΩÞ, there are no contributions that are concerned with nonlinear problems formulated in W 1;p a 2 ðΩÞ, like the problem with p-exponent diffusion coefficient that is considered here. Maybe as one exception, we can refer the work presented in [18], where IPDG approximate solutions are

I. Toulopoulos / Finite Elements in Analysis and Design 95 (2015) 42–50

studied for the p-Laplace equation ðμ ¼ 0Þ. It is the purpose of this paper to make a first step in this direction. We point out that classical (continuous) finite element methods, for more general p-form problems, the so-called ðp; δÞ
structure problems, have been analyzed in the literature, see e.g. [19,20]. For parabolic ðp; δÞ
structure problems, we refer [21], where optimal convergence rates have been shown, in case of using linear finite element in space and implicit Euler scheme in time. The IPDG scheme proposed here, see (3.8), has the same form as the IPDG scheme in [12], but here the numerical flux is appropriately re-formulated in order to be compatible with the p-nature of the problem. As a first task, stability bounds are proved in ‖:‖DG;p
norm, for the case of μ¼ 0. Then, using the interpolation estimates presented in [20], a priori error estimates are given for the semi-discrete problem for p 4 2, assuming conventional regularity for the exact solution. The IPDG spatial discretization generates a nonlinear ODE system with respect to the degrees of freedom. We discretize in time this system by s-stage Diagonally Implicit Runge–Kutta (s-DIRK) methods. Every cycle of the Runge– Kutta method includes the solution of nonlinear algebraic systems. Two low computational cost Picard block-iterative methods are proposed for solving the nonlinear systems [22]. The Picard iterative methods are constructed based on the local (per element) approximation features of the IPDG method. The two different iterative methods are expected to have the same order of convergence (first order), but different performance speed, since the second one uses the latest available solution (and not the solution of the previous iteration) for updating the nonlinear parts of the system. The outline of the paper is as follows. It begins by presenting the model problem. Then inequalities for vectors a; b A R2 are shown, which are used later to derive the continuity–monotonicity properties of the scheme. In Section 3, the IPDG method is described. Section 4 includes the formulation of the s-DIRK method for the time discretization and the description of the two Picard methods. In Section 5, a stability bound for the discrete solution of the p-Laplace problem is presented. A priori error estimates for p 4 2 are shown in Section 6. The paper closes with the numerical tests in Section 7.

2. The model problem

W 1;p ðΩÞ. This fact motivates the need of further analysis and of developing numerical fluxes compatible with the p-exponent form of the diffusion coefficient. The goal of this paper is to make a first step in this direction. Assuming that f A Cð½0; T; L2 ðΩÞÞ; u0 A W 1;2 ðΩÞ \ W 1;p ðΩÞ, we call u weak solution of (2.1), if u A L1 ð0; T; W 1;p ðΩÞÞ \ W 1;2 ð0; T; W 1;2 ðΩÞÞ, ujΓD ≔uD satisfies the following formulation for any v A W 1;p 0 ðΩÞ: Z Z Z ut v dx þ Að∇uÞ  ∇v dx ¼ fv dx; ð2:2Þ 8 t 4 0; Ω

ut
div Að∇uÞ ¼ f u0 ðxÞ ¼ uðx; 0Þ u ¼ uD

in Ω  ð0; T

in Ω

ð2:1aÞ ð2:1bÞ

on Γ D  ð0; T;

ð2:1cÞ

Ω

 Z Lp ð0; T; VÞ ¼ v : ð0; TÞ-V :

ð2:3Þ

Through the paper C i ; i ¼ 1; …, will be generic constants with different values independent of crucial quantities. The explicit dependence on the problem data will be mentioned. 2.1. Helpful inequalities for vectors Working further on the results of Chapter I in [24,26], we prove special algebraic inequalities that are going to be used later. In the proofs, we use the function F : R2 -R2 FðaÞ ¼ ðμ þ jajÞðp
2Þ=2 a:

ð2:4Þ

We introduce the formula Z 1 d ðμþ ja þ tðb
aÞjÞp
2 ðaþ tðb
aÞÞ dt; AðbÞ
AðaÞ ¼ 0 dt

ð2:5Þ

and by an easy computation on the right-hand side of (2.5) we get AðbÞ
AðaÞ Z 1 ðμ þ ja þtðb
aÞjÞp
2 ðb
aÞ dt ¼ Z

1

þ ðp
2Þ 0


1 1  ðμ þjaþ tðb
aÞjÞp
3 a þ tðb
aÞ 2

2ðaþ tðb
aÞ; b
aÞðaþ tðb
aÞÞ dt:

ð2:6Þ

Multiplying (2.6) by b
a, we have ðAðbÞ
AðaÞ; b
aÞ Z 1 ðμ þ jaþ tðb
aÞjÞp
2 dt ¼ jb
aj2 0

Að∇uÞ ¼ ðμ þ j∇ujÞp
2 ∇u;

ða þtðb
aÞ; b
aÞ2 dt:

where j  j : R -R is the Euclidean measure and að∇uÞ ¼ ðμ þj∇ujÞp
2 is the diffusion coefficient. The nonlinear nature of the problem (2.1) comes by the appearance of j∇uj in the diffusive coefficient and this poses numerical challenges. The IPDG methods presented so far for nonlinear elliptic equations are referred to problems where the natural formulation is given in W 1;2 ðΩÞ, and either aðÞ is uniformly bounded, e.g. [23], or aðÞ satisfies a monotone condition, e.g. [11,12]. One cannot apply the same methodology for the problem (2.1) which is formulated in

0

 ‖vðtÞ‖pV dt o1 :

u A W 1;2 ð0; T; W 1;2 ðΩÞÞ \ Lp ð0; T; W s Z 2;p ðΩÞÞ:

þðp
2Þ

2

T

The existence–uniqueness of the solution of (2.2) (even with other assumptions on the data) are ensured by means of the monotone operators theory, see e.g. [1,24]. We refer [21,25], for regularity assumptions on the problem data for obtaining optimal rate of convergence for finite element solutions. For the analysis here, we assume the following conventional assumptions:

where ð0; T is the time interval, f : Ω  ð0; T-R, u0 : Ω-R, uD : Γ D  ð0; T-R are given smooth functions. The operator Að∇uÞ : R2 -R2 has the form p 4 1; μ Z 0;

Ω

where

0

Let Ω be a bounded domain in R2 , with smooth boundary Γ D ≔∂Ω. We consider the following scalar initial boundary value problem:

43

Z

1

ðμ þjaþ tðb
aÞjÞp
3 ja þ tðb
aÞj
1

0

ð2:7Þ

The last term on the right-hand side of (2.7) is positive and for p Z 2 we get Z 1 ðAðbÞ
AðaÞ; b
aÞ Z jb
aj2 ðμ þ ja þ tðb
aÞjÞp
2 dt ð2:8aÞ 0

by applying Cauchy–Schwarz inequalities, we further get Z 1 ðμþ ja þ tðb
aÞjÞp
2 dt: jAðbÞ
AðaÞj Z jb
aj 0

ð2:8bÞ

44

I. Toulopoulos / Finite Elements in Analysis and Design 95 (2015) 42–50

Also, applying Cauchy–Schwarz inequality on the last term on the right-hand side of (2.6), we have Z 1 ðμ þ ja þ tðb
aÞjÞp
3 ja þtðb
aÞj
1 ja þ tðb
aÞj2 jb
aj 0

Z r jb
aj

1

ð2:9Þ

0

jAðbÞ
AðaÞj r ðp
1Þjb
aj Z 1 ðμ þ ja þ tðb
aÞjÞp
2 dt:

ð2:10Þ

0

Recalling the forms of A and F and setting in (2.10) p≔ðp þ 2Þ=2 we obtain

r

ðp
2Þ=2

p2 2

fvg ¼ vjEin fvgD ¼

and

1 2 ðvjEin

The space ðμ þ ja þ tðb
aÞjÞp
2 dt:

Therefore, combining (2.9) and (2.6), we get

jðμ þjbjÞ

In case of e A E D , we define

ðp
2Þ=2

2

b
ðμ þ jajÞ aj !2 Z 1 ðp
2Þ=2 2 jb
aj ðμ þ ja þ tðb
aÞjÞ dt

ð2:11aÞ

0

and thus jFðbÞ
FðaÞj2 Z 1  p 2 r jb
aj2 ðμ þ jaþ tðb
aÞjÞp
2 dt: 2 0

jFðbÞ
FðaÞj2 rCðpÞðAðbÞ
AðaÞ; b
aÞ;

ð2:12aÞ

jFðbÞ
FðaÞj2 rjAðbÞ
AðaÞj2 :

ð2:12bÞ

Keeping p 42 and using that ðμ þ ja þ tðb
aÞjÞ2ðp
2Þ=2 r 2 max fðμ þ jajÞðp
2Þ=2 ; ðμ þ jbjÞðp
2Þ=2 gðμþ ja þ tðb
aÞjÞðp
2Þ=2 , we derive by (2.10) that Z 1 jAðbÞ
AðaÞj r ðp
1Þjb
aj ðμ þ ja þ tðb
aÞjÞ2ðp
2Þ=2 dt r CðpÞM ðjaj;jbjÞ

jb
ajðμþ ja þ tðb
aÞjÞ

ðp
2Þ=2

dt

0

r CðpÞM ðjaj;jbjÞ jFðbÞ
FðaÞj;

ð2:13Þ ðp
2Þ=2

where M ðjaj;jbjÞ ¼ 2 maxfðμþ jajÞ

; ðμ þ jbjÞ

ðp
2Þ=2

g.

3. The numerical scheme

E Let T h ¼ fEi gN i ¼ 1 be a regular subdivision of Ω in triangular elements (without hanging nodes) with diameter hEi , where for simplicity we assume h≔minEi A T h hEi ¼ maxEi A T h hEi . We denote by E ¼ E I [ E D all the edges, where E I is the set of the interior edges of Th, that is E I ¼ fe : e ¼ ∂Ein \ ∂Eout ; for Ein ; Eout A T h g and E D is the set of the Dirichlet boundary edges E D ¼ fe : e ¼ ∂Ein \ Γ D ; Ein ; A T h g. For each of e A E I we associate a unit normal vector ne . For e A E D ; ne is considered to be the outward normal to ∂Ω. Define the following broken Sobolev spaces for s Z 2; p 4 1:

: vjE A W

s;p

ðEÞ; 8 E A T h g;

and the discontinuous finite element space V kh ðT h Þ≔fv A Lp ðΩÞ : vjE A Pk ðEÞ; 8 E A T h g;

ð3:1Þ V kh ðT h Þ 

W s;p ðT h Þ h ð3:2Þ

where Pk ðEÞ is the space of polynomials of degree less than or equal to k. Let e A E I , we define the average and the jump of v A W s;p ðT h Þ on h e by fvg ¼ 12 ðvjEin þ vjEout Þ

and

½vD ¼ vjEin
uD :

ð3:4Þ

is equipped with the broken DG norm [27,18], Z  p ½ϕ j∇ϕjp dx þ ∑ σh   ds J ϕ J pDG;p ¼ ∑ e h  e A EI E A Th E p Z   ½ϕD  ð3:5Þ þ ∑ σh   ds; e h  e A ED

where p 41 and σ 4 0 is a parameter. 3.2. Auxiliary results Next, we summarize some results from the literature, which are going to be frequently used. ðT h Þ Lemma 3.1 (Trace inequalities). For vh A V kh ðT h Þ, and v A W s;p h with p 4 1 there exist positive constants C 1 ðk; pÞ; C 2 ðk; pÞ; C 3 ðk; pÞ independent of the mesh size, such that (i) ‖h vh ‖pLp ðΓ Þ r C 1 ∑E A T h h‖vh ‖pLp ð∂EÞ , D
1 (ii) ‖vh ‖pLp ð∂EÞ r C 2 h ‖vh ‖pLp ðEÞ ,
1=2 (iii) ‖v‖L2 ð∂EÞ r C 3 h ð‖v‖L2 ðEÞ þ h‖∇v‖L2 ðEÞ Þ. Proof. The proofs of the above inequalities can be found in [7].

□

0

Hölder, Young and Poincare's inequalities: let 1 o p; p o 1 such 0 that 1=p þ 1=p0 ¼ 1, and ϵ 4 0, then for u A Lp ðΩÞ and v A Lp ðΩÞ we have Z juvj dx r ‖u‖Lp ðΩÞ ‖v‖Lp0 ðΩÞ ; ð3:6aÞ Ω

Z

0

ϵ ϵ
p =p juvj dx r ‖u‖pLp ðΩÞ þ ‖v‖p0p0 : L ðΩÞ p p0 Ω

ð3:6bÞ

ðT h Þ, see The generalized Poincare–Friedrichs inequality for v A W 1;2 h [28], !1=2 1 ‖½v‖2L2 ðeÞ : ð3:7Þ ‖v‖L2 ðΩÞ rC ∑ ‖∇v‖2L2 ðEÞ þ ∑ e A EI [ ED h E A Th 3.3. The IPDG discretization

3.1. Preliminaries – DG notation

ðT h Þ≔fv A Lp ðΩÞ W s;p h

and

W s;p ðT h Þ h Z

0

1

þ uD Þ

1=p

ð2:11bÞ

By (2.8a) and (2.8b), we have

Z

½v ¼ vjEin ;

½v ¼ vjEin
vjEout :

ð3:3Þ

For the simplification of the formulas below, we will often use R dx instead of ∑E E u dx. Inspired by the IPDG method in [12], we present the IPDG numerical scheme for discretizing the problem (2.1). We introduce the semi-linear form B : W s;p ðT h Þ  W s;p ðT h Þ-R, such h h s;p that for u; ϕ A W h ðT h Þ Z Bðu; ϕÞ ¼ ∑ að∇uÞ∇u∇ϕ dx EAT E Z   að∇uÞ∇u  ne ½ϕ ds
∑ R

Ωu

e A EI e

 Z 

½u ∇ϕ  ne ½u ds a h e A EI e Z

σ ½u ½u½ϕ ds þ ∑ a h e A EI h e

Z ½uD
∑ ∇ϕ  ne ½uD ds a h e A ED e

Z σ ½uD ½uD ϕ ds þ ∑ a h h e e A ED Z að∇uÞ∇u  ne ϕ ds;
∑
∑

e A ED e

ð3:8Þ

I. Toulopoulos / Finite Elements in Analysis and Design 95 (2015) 42–50

where σ≔σðk; pÞ is a positive parameter and will be specified in the error analysis. Giving an interpretation of the terms that appear in (3.8), we can say that the second integral in (3.8) gives an approximation of the trace of the nonlinear flux and ensures the consistency of the method. The third integral “symmetrizes” the flux form of Bð_; Þ_ which is important for the numerical computations. The fourth integral penalizes the jumps on the interfaces and helps for achieving _ The rest terms defined on the boundary the discrete coercivity of Bð_; Þ. edges have similar meaning with the formers. We also define the linear form Z LðϕÞ ¼ ∑ f ϕ dx: ð3:9Þ E A Th E

The semi-discrete problem is formulated as follows: find uh A W 1;2 ð0; T; V kh ðT h ÞÞ such that

Z Ω

∂uh ðtÞ ϕ dx þ Bðuh ; ϕÞ ¼ LðϕÞ; ∂t uh ð0Þ ¼ u0;DG ;

8 ϕ A V kh ðT h Þ ð3:10Þ

where u0;DG is the approximation of the initial condition to the V kh ðT h Þ space. Due to the assumed regularity (2.3) for the weak solution u (note that the jumps ½u ¼ 0 on the interfaces), it is easy to show that u satisfies the variational formulation (3.10), Z ∂uðtÞ ϕ dx þ Bðu; ϕÞ ¼ LðϕÞ; 8 ϕ A V kh ðT h Þ: ð3:11Þ Ω ∂t For every E A T h , the DG solution of (3.10) is expressed as uh ¼ ∑i U Ei ðtÞP i ðxÞ where UEi are the degrees of freedom and P i ðxÞ A Pk ðEÞ are the local polynomial basis functions. When this expression is substituted into (3.10), we obtain the following nonlinear ODE problem of finding the vector U ¼ ½…; U Ei ; … such that dUðtÞ þBðUðtÞÞ ¼ LðtÞ; dt Uð0Þ ¼ uh ð0Þ M

ð3:12Þ

where M is the block-diagonal mass matrix and the entries of B and L are specified by (3.8) and (3.9) respectively.

We discretize (3.12) with respect to time using s-stage Diagonally Implicit Runge–Kutta (s-DIRK) methods [29]. Hereafter, we denote by Δt the time step and with Un the approximation of Uðt n Þ at time t n ¼ nΔt; n ¼ 0; 1; 2; … . If τi ; i ¼ 1; …; s, are the quadrature points, bi are the weights and aij ; j ¼ 1; …; i, are the entries of Bucher's table, the s-DIRK method for the problem (3.12) is given by i

MΔUn;i ¼
Δt ∑ aij ðBðUn;j Þ
Lðt n;j ÞÞ;

i ¼ 1ð1Þs

ð4:1aÞ

j¼1

MU

nþ1

n

the Picard iterative process (4.2) are applied, (i) the element-Jacobi (PEJ) and (ii) the element Gauss–Seidel (PEGS). Both iterative approaches are simple applications of the Picard iterative method presented in [22]. In the element-Jacobi scheme, the full Picard matrix BP ðUn;l
1 Þ is approximated only by the block diagonal entries, neglecting in that way the contribution of the offdiagonal matrix blocks, which arise through the evaluation of the numerical fluxes on the interfaces. The numerical fluxes are computed using the previous solution vector Un;l
1 and are added to the right-hand residual R. The diagonal blocks of BP ðUn;l
1 Þ represent small dense matrices and are associated with each element E A T h . The solution of the resulting PEJ system of (4.2) is performed element by element using LU factorization method. The convergence of the previous proposed PEJ iterative method can be further accelerated by using Gauss–Seidel strategy, giving in this way, the second mentioned PEGS iterative method. PEGS method applies the same splitting of the matrix BP ðUn;l
1 Þ, but follows a passing over the interfaces by computing the numerical fluxes using the latest available solution Un;l or Un;l
1 where it is possible. In the numerical tests (see Section 7), the iterative process of (4.2) stops when J Un;l
Un;l
1 J o tol for a prescribed tolerance tol and then we set Un;i ≔Un;l . We point out that PEGS method is expected to have similar convergence rates per Runge– Kutta cycle as the PEJ method, but more improved performance behavior in terms of CPU (in fact the stopping criterion J Un;l
Un;l
1 J otol is achieved performing fewer iterations l than the PEJ method). Comparison between the two iterative methods will be shown in Section 7. Other higher-order iterative procedures (e.g. Newton) can be applied for computing the intermediate solutions of (4.1a). In many cases, the computation of the Jacobian matrix of BðUðtÞÞ may increase the CPU time of the whole ODE solver, see examples in [22], and more advanced numerical techniques must be applied, see e.g. [30,31]. Anyway, for the numerical tests presented in Section 7, the previous proposed Picard iterative methods have been found to be appropriate for solving (3.12).

5. A stability bound in ‖  ‖DG;p for the case of μ¼ 0

4. Fully discrete formulation

In this section, we give a stability estimate (a priori bound) for the DG solution uh in case of μ¼0 (p-Laplace problem) and note that Gronwall's lemma is not used. Stability bounds can also be obtained working in different directions using the monotonicity _ which are presented later. Here, the stability properties of Bð_; Þ, bound uses the ‖:‖DG;p
norm (3.5). Lemma 5.1. For the form (3.8) with μ¼0, there are constants κ 4 0; C D 4 0 such that Bðϕ; ϕÞ Z κ‖ϕ‖pDG;p


s

45

n;i

n;i

¼ MU
Δt ∑ bi ðBðU Þ
Lðt ÞÞ;

ð4:1bÞ

CD p
1

h

‖uD ‖pΓD ;

8 ϕ A V kh ðT h Þ:

ð5:1Þ

i¼1

where t ¼ t þ τi Δt and ΔUn;i ¼ Un;i
Un . The computation of the intermediate solutions Un;i in (4.1a) includes the solution of a nonlinear system, which is achieved by a Picard iterative process n;i

n

E A Th E

for l ¼ 1; …; lM ; compute Un;l by ðM þ aii ΔtBP ðU

n;l
1

set Un;i ¼ Un;lM ; n;l
1

ÞÞU

n;l

Proof. Choosing uh ¼ ϕ in (3.8) we obtain Z Z   Bðϕ; ϕÞ ¼ ∑ að∇ϕÞ∇ϕ  ∇ϕ dx
∑ að∇ϕÞ∇ϕ  ne ½ϕ ds

n

j

¼ RðU ; U Þ; ð4:2Þ

Þ is the iterative matrix produced by the Picard where BP ðU linearization and RðUn ; Uj Þ≔MUn
Δt∑ij
¼11 aij ðBðUn;j Þ
Lðt n;j ÞÞ is the residual computed using the previous solutions. In the present work, for computational efficiency, two low-storage variations of

e A EI e

 Z

Z 

½ϕ σ ½ϕ ∇ϕ  ne ½ϕ ds þ ∑ ½ϕ½ϕ ds a a
∑ h h e A EI e e A EI h e

Z ðϕ
uD Þ ∇ϕ  ne ðϕ
uD Þ ds
∑ a h e A ED e

Z σ ðϕ
uD Þ ðϕ
uD Þϕ ds a þ ∑ h e A ED h e

46

I. Toulopoulos / Finite Elements in Analysis and Design 95 (2015) 42–50

Z
∑

e A ED e

að∇ϕÞ∇ϕ  ne ϕ ds

¼ T 1
T 2
T 3 þ T 4
T 5 þ T 6
T 7 For the term T1, we have Z Z að∇ϕÞ∇ϕ  ∇ϕ dx ¼ ∑ j∇ϕjp dx: T1 ¼ ∑ E A Th E

E A Th E

For the rest of the terms, applying inequalities (3.6), Lemma 3.1 and introducing constants C i;ε ≔C i ðε; p; p0 Þ with ε 4 0, it follows that   Z     fað∇ϕÞ∇ϕ  ne g½ϕ ds T2 r ∑ e A E I e  Z    að∇ϕÞ∇ϕjg½ϕ ds r ∑ e A EI e

Z

r ∑

1=p0 

h

j∇ϕjp
1

e A EI e

1=p0

h

1=p0

h

0

fj∇ϕjp=p g

Similarly, adding uD
uD the term T7 can be bounded Z    T7 r ∑ j∇ϕjp
2 ∇ϕðϕ
uD þ uD Þ ds e A ED e

For T3, working in the same way as for T2 we have Z   ½ϕ ½ϕ 1=p0 þ 1=p   a ∇ϕ ds T3 r ∑  h h h e A EI e   0 Z  p=p ½ϕ   0   h1=p h1=p ∇ϕ ds r ∑ h e A EI e  0   11=p0 Z

1=p Z  p B ½ϕ 1=p r ∑ @   h dsA ðh fj∇ϕjgÞp ds e h  e e A EI   ! p Z   Z 1 ½ϕ r ∑   h ds þC 3;ε hfj∇ϕjgp ds e e A E I C 3;ε e  h  Z  p Z h ½ϕ ∑ j∇ϕjp ds: r   ds þ 3C 3;ε ∑ C 3;ε e A E I e  h  E A Th E A straightforward computation for the term T4 gives Z  p ½ϕ T 4 ¼ ∑ hσ   ds: e h  e A EI For the term T5 applying the same steps as for T3 yields Z Z h jϕ
uD jp T5 r ∑ dsþ 3C 5;ε ∑ j∇ϕjp ds; C 5;ε e A E D e h E D A T h ED where ED A T h are the boundary elements: fE A T h : ∂E \ Γ D a ∅g. Term T6 can be bounded as follows: Z jϕ
uD jp
2 ðϕ
uD Þ ðϕ
uD þ uD Þ ds T 6 ¼ ∑ σh h h h e e A ED

p Z jϕ
uD j ¼ ∑ σh ds h e e A ED

Z jϕ
uD j p
2 ðϕ
uD Þð
uD Þ
∑ σh ds 2 h e e A ED h !p Z   ϕ
uD  Z ∑ σh ds   h  e  e A ED

j∇ϕjp dx

ED A T h ED

ds

p0

Z

r C 7;ε ∑ þ

1=p0 Z

1=p

j½ϕj p ds ds 1=p0 e e h e A EI  p
2  0 Z

Z 1=p ½ϕ 0 ½ϕ2 1=p 0 p     ds r ∑ hfj∇ϕjp=p g ds h  h   h e e e A EI 0  p
2  2 1   Z Z   ½ϕ ½ϕ h B   dsC h  r ∑ @C 2;ε hj∇ϕjp dsþ A h C h 2;ε e  e e A EI      Z Z  p ½ϕ h   ds: r 3C 2;ε ∑ j∇ϕjp dx þ ∑ C 2;ε e A E I e  h  E A Th E r ∑

Z 

 j½ϕj

Z jϕ
uD j p
1 ð
uD Þ
∑ σh ds h h e e A ED   p Z   ϕ
uD  Z ð1
C 6;ε Þ ∑ σh   ds h  e e A ED Z 1 juD jp ∑ σh ds:
p
1 C 6;ε e A E D eh

 Z  ϕ
uD p juD jp 1  þ ∑ σh  ds: C 7 ; εe A E D h  hp
1 e

In the previous inequalities, choosing C i;ε such that 3C 2;ε þ 3C 3;ε þ 3C 5;ε þ C 7;ε r 12 ; and choosing the parameter σ to satisfy the following relations: σ 41

ð1
C 6;ε Þσ Z

1 1 þ ; C 5;ε C 7;ε

σ4

1 1 þ ; C 2;ε C 3;ε

while keeping h r1, we can find κ 4 0 and CD, in order (5.1) to be true. □ Now, choosing ϕ ¼ uh ðtÞ in (3.10) and using (5.1), we have d ‖u ‖22 þ κ‖uh ðtÞ‖pDG;p 2dt h L ðΩÞ    CD r Lðuh ðtÞÞ þ p
1 ‖uD ‖pΓD : h

ð5:2Þ

Applying inequality (3.6) on the right-hand side of (5.2), we get     Lðuh ðtÞÞ r 1 ‖f ðtÞ‖p0p0 þC 8;ε ‖uh ðtÞ‖pp : ð5:3Þ   C L ðΩÞ L ðΩÞ 8;ε Based on the discrete embeddings, see [27], Z Z   1 ‖ϕ‖pLp ðΩÞ r C p ∑ j∇ϕjp þ ∑ p
1 ½ϕjp ; e e A Eh E E 8 ϕ A V kh ðT h Þ;

ð5:4Þ p
1 ‖uh ‖pLp ðΩÞ r C p J uh J pDG;p þ ðC p =h Þ‖uD ‖pΓD .

we can easily show that Thus, inserting (5.3) into (5.2) and then applying inequality (5.4), we obtain for C 8;ε ¼ κ=2 that d κ ‖u ‖22 þ ‖u ‖p 2dt h L ðΩÞ 2 h DG;p 1 CD ‖f ‖p0p0 þ p
1 ‖uD ‖pΓD : r L ðΩÞ Cκ; p h

ð5:5Þ

Integrating from 0 to t, we get the following stability bound for uh: Z t ‖uh ðτÞ‖pDG;p dτ ‖uh ðtÞ‖2L2 ðΩÞ þ κ 0 Z t ‖f ðτÞ‖p0p0 r ‖u0h ‖2L2 ðΩÞ þ C 0

þ

CD h

p
1

L ðΩÞ

‖uD ðτÞ‖pΓD dτ:

ð5:6Þ

6. Continuous in time a priori error estimates Next, we give an error estimate on how close is the IPDG solution uh of (3.10) to u of (3.11), that is an estimate for



½u ½uh  2
F J L2 ðeÞ ; ‖u
uh ‖2F;DG ¼ ∑‖Fð∇uÞ
Fð∇uh Þ‖2L2 ðEÞ þ ∑ σh J F h h e A EI E

I. Toulopoulos / Finite Elements in Analysis and Design 95 (2015) 42–50

þ ∑ σh J F e A ED





½uD ½uh D 2
F J L2 ðeÞ ; h h

ð6:1Þ

where the function F has been defined in (2.4) and the jumps ½ in (3.3) and (3.4). We mention that similar error formula has been used in [32], where a LDG method studied for ðp; δÞ
structure problems. We consider the case where the solution u has the regularity (2.3), uh A V 1h ðT h Þ and I u A V 1h ðT h Þ is the Scott–Zhang interpolant of u [33]. For problem (2.1), we suppose that uD A P1 ðE D Þ and the parameter μ is such that (for example μ¼1) Z 1 ðμ þ ja þ sðb
aÞjÞp
2 ds Z 1; ð6:2Þ 0

where in (6.2), a represents u or uh, either their gradients and b takes the role of I u or its gradient. In the error analysis, we will make use of the following approximation result, which has been proved in [20]. Lemma 6.1. Let u A W s Z 2;p ðΩÞ with Fð∇uÞ A W 1;2 ðΩÞ, and I u A V 1h ðT h Þ its Scott–Zhang interpolant. Then there are constants C 1 ; C 2 4 0 independent of h such that 2

‖Fð∇uÞ
Fð∇I uÞ‖2L2 ðEÞ r C 1 h ‖∇Fð∇uÞ‖2L2 ðS Þ ; E

8 E A T h;

ð6:3Þ

where SE is a domain made of the neighboring elements of E in Th. Corollary 6.2. Under the assumptions of Lemma 6.1 and (2.3) the following estimate holds true for t 4 0: 2

‖u
I u‖2F;DG o Ch ∑ ‖∇Fð∇uÞ‖2L2 ðEÞ

ð6:4Þ

E A Th

for C 4 0 independent of h. Proof. We observe for u and the Scott–Zhang interpolant I u that ½u ¼ ½I u ¼ 0 on every e A E. The estimate (6.4) follows immediately by the (6.1) and the approximation result (6.3). □ Proposition 6.3. Under the assumptions (6.2), we can obtain the following estimates: ‖u
I u‖2L2 ðΩÞ r C Ω;p;u;I u ‖u
I u‖2F;DG ;

ð6:5aÞ

‖uh
I u‖2L2 ðΩÞ r C‖uh
I u‖2F;DG :

ð6:5bÞ

Proof. We recall inequality (3.7) for v≔u
I u and consequently we apply (2.8b) and (2.13) to obtain

Lemma 6.5. Let I u A V 1h ðT h Þ be the interpolant of u as in (6.3) and let ϕ ¼ uh
I u. For every edge e A E there are C 1;ε ; C 2;ε 4 0 such that Z       fað∇uh Þ∇uh
að∇I uÞ∇I ug  ne ½ϕ ds  e  



2 h  ½uh  ½I u 
F : r C 1;ε ‖Fð∇uh Þ
Fð∇I uÞ‖2L2 ðEin ⋃Eout Þ þ  F C 2;ε  h h  2 L ðeÞ

ð6:7Þ Proof. Let e ¼ ∂Ein ⋂∂Eout (or e A E D ). Applying sequentially the inequalities (3.6) and Lemma 3.1 on the left-hand side of (6.7), we have   Z  1   1=2  h  að∇uh Þ∇uh
að∇I uÞ∇I u  1=2 ½ϕ ds e h Z Z   1=2 1=2 1 j½ϕ2 ds hjfað∇uh Þ∇uh
að∇I uÞ∇I ug2 ds r e eh

1=2 h h r ‖að∇uh Þ∇uh
að∇I uÞ∇I u‖2L2 ðein Þ þ ‖…‖2L2 ðeout Þ 2 2

1=2  1  ‖½ϕ‖L2 ðeÞ rC trc ‖að∇uh Þ∇uh
að∇I uÞ∇I u‖2L2 ðE Þ in h

1=2 1=2 1 ‖½ϕ‖2L2 ðeÞ þ‖…‖2L2 ðE Þ out h   2     ½ϕ h   r C 1;ε J að∇uh Þ∇uh
að∇I uÞ∇I u‖2L2 ðE ⋃E Þ þ out in C 2;ε  h   L2 ðeÞ

h F ½uh  r ðby ð2:8bÞ;ð2:13ÞÞ C 1;ε ‖Fð∇uh Þ
Fð∇I uÞ‖2L2 ðE ⋃E Þ þ out in C 2;ε h

½I u 2
F L2 ðeÞ ; h where for simplicity, we used the notation: L2 ðein Þ≔L2 ðe  ∂Ein Þ.

□

I u A V 1h ðT h Þ

be the interpolant as in (6.3) of the Lemma 6.6. Let solution u and ϕ ¼ uh
I u. For every edge e ¼ ∂Ein ⋂∂Eout (or e A E D ) there are C 1;ε ; C 2;ε 4 0 such that Z 





    ½uh  ½I u   ½uh 
a ½I u ∇uh
∇I u  ne ds a   e  h h 



2 Z   ½uh  h ½I u  F
F r  dsþ C 1;ε ‖Fð∇uh Þ C 2;ε e  h h 
Fð∇I uÞ‖2L2 ðEin ⋃Eout Þ :

‖u
I u‖2L2 ðΩÞ r C Ω;p;u ‖u
I u‖2A;DG r C Ω;p;u ‖u
I u‖2F;DG :

47

ð6:8Þ

For the estimate (6.5b), let eA E D , then we have that juh
I uje r juh
uD jje þ juD
I ujje . Therefore, using the inequality ‖uh
I u‖L2 ðΩÞ r ‖uh
I u‖DG;2 (see (3.7)) and then applying (2.8b) and (2.13) for every term of ‖uh
I u‖DG;2 , we get

Proof. Following the same steps as in proof of Lemma 6.5, by applying Hölder's inequality, trace inequality, consequently Young's inequality and (2.8b) and (2.13), the relation (6.8) can be shown. □

‖uh
I u‖2L2 ðΩÞ r ‖uh
I u‖DG;2 r C Ω;p;uh ‖uh
I u‖2F;DG :

Theorem 6.7. Let u be the solution of (3.11) and let I u A V 1h ðT h Þ be its interpolant as in (6.3). Then for ϕ ¼ uh
I u and ε 4 0 there exist constants C 1;ε ; C 2;ε ; C 3;ε such that the form B of (3.8) satisfies

□

Lemma 6.4. Under the assumptions of Lemma 6.1, there exists a C 4 0 independent of h such that ∑

eAE

2 h‖Fð∇uÞ
Fð∇I uÞ‖2L2 ðeÞ r Ch

∑

E A Th

‖∇Fð∇uÞ‖2L2 ðEÞ :

ð6:6Þ

Proof. Using v≔Fð∇uÞ
Fð∇I uÞ in inequality (iii) of Lemma 3.1 and summing over all edges, we have that

∑ h‖Fð∇uÞ
Fð∇I uÞ‖2L2 ðeÞ r 3C ∑ ð‖Fð∇uÞ
Fð∇I uÞ‖2L2 ðEÞ

eAE

E A Th

2

þh

2 ‖∇ðFð∇uÞ
Fð∇I uÞÞ‖2L2 ðEÞ Þ r Ch

∑ ‖∇Fð∇uÞ‖L2 ðEÞ :

E A Th

□

1 ‖u
I u‖2F;DG þ C 2;ε ‖ϕ‖2F;DG jBðu; ϕÞ
BðI u; ϕÞj r C 1;ε Z 2 þ ∑ hjFð∇uÞ
Fð∇I uÞj2 ds: C 3;ε e A E e

ð6:9Þ

Proof. After a rearrangement of the terms of Bðu; ϕÞ
BðI u; ϕÞ, we have   Z      Bðu; ϕÞ
BðI u; ϕÞ r að∇uÞ∇u
að∇I uÞ∇I u∇ϕ dx   Ω Z    að∇uÞ∇u
að∇I uÞ∇I ujg½ϕ ds þ ∑ e A EI e

48

I. Toulopoulos / Finite Elements in Analysis and Design 95 (2015) 42–50

Z þ ∑

e A ED e

   að∇uÞ∇u
að∇I uÞ∇I uϕ ds

The last terms T6 and T7 are bounded by applying the same steps as before   

Z 

½I u ½I u 1=2  ½ϕ   ½u ½u
a T 6 rσ ∑ a h  1=2  ds h  h h h h  e A EI e  0  2 11=2 0Z   !1=2

Z 

 ½u ½u ½I u ½I u A @ 1  B 
a r σ ∑ @ ha  ½ϕ ds  h h h h  e  e h e A EI

( )

Z 

   ½I u  ½u    ½u
a ½I u ∇ϕ ds þ ∑ a    h h e A EI e   



Z    ½I uD  ½uD   ½uD
a ½I uD ∇ϕ ds þ ∑ a   h h e A ED e 



Z  ½u ½u ½I u ½I u  
a ½ϕ ds þ ∑ σ ðaÞ h h h h e e A EI





Z ½uD ½uD ½I uD ½I uD   þ ∑ σ
a ϕ ds a h h h h e e A ED







Z   ½u  ½uh  ½I u 2 
F hF hF  dsþ C 1;ε ∑ h h  h e  e A EI e  





2 Z   ½uD ½I u 2 1 ½I uD  
F ds:T 7 r ∑ hF
F  ds  h  C 2;ε e A E D e  h h 



2 Z   ½uh D ½I uD 
F hF þC 1;ε ∑  ds:  h h e A ED e 

1 r ∑ C 2;ε e A E I

¼ T 1 þ T 2 þT 3 þ T 4 þ T 5 þ T 6 þ T 7 : The first term T1 can be bounded by applying Hölder–Young's inequalities (3.6) and consequently (2.8b), (2.13), as follows: !1=2 !1=2 Z Z 2 2 T1 r ∑ jað∇uÞ∇u
að∇I uÞ∇I uj dx ∑ j∇ϕj dx E A Th E

1 r C 2;ε

E A Th E

Z

Choosing appropriate the constants C i;ε above and by gathering the bounds together, we can derive (6.9). □

Z 2

Ω

jFð∇uÞ
Fð∇I uÞj dx þ C 1;ε

2

Ω

jFð∇uh Þ
Fð∇I uÞj :

The term T2 can be bounded by applying the same steps as in Lemma 6.5,   Z  1   1=2  T2 r ∑ h  að∇uÞ∇u
að∇I uÞ∇I u  1=2 ½ϕ ds e A EI e h " Z

1=2 Z

1=2 # 1 j½ϕj2 ds hjfað∇uÞ∇u
að∇I uÞ∇I ugj2 ds r ∑ e eh e A EI 1 ∑ r C 2;ε e A E I

Z



2  ½u 

½I u   h 2
F hjFð∇uÞ
Fð∇I uÞj ds þ C 1;ε hF :   h h  2 e

Z

Theorem 6.8. Let I u A V 1h ðT h Þ be the interpolant of the solution u. The form Bð_; Þ_ is monotone with respect to second argument, in the sense that there is a κ0 40 such that Bðuh ; uh
I uÞ
BðI u; uh
I uÞ 4 κ 0 ‖uh
I u‖2F;DG :

Proof. Denoting ϕ ¼ uh
I u and after rearranging the terms, we obtain that Bðuh ; uh
I uÞ
BðI u; uh
I uÞ Z ðað∇uh Þ∇uh
að∇I uÞ∇I uÞ∇ϕ dx ¼ ∑ E A Th E

L ðeÞ

Z


∑

e A EI e

Analogously, for the term T3, we obtain that   Z  1   1=2  T3 r ∑ h að∇uÞ∇u
að∇I uÞ∇I u 1=2 ϕ ds e A ED e h  Z  ½u 

1  h D r ∑ hjFð∇uÞ
Fð∇I uÞj2 dsþ C 1;ε hF  C 2;ε e A E D e h

2 ½I uD 
F :   2 h

r

1 ∑ C 2;ε e A E I

Z





 Z  ½u ½I u 2 
F hF jFð∇uh Þ  ds þ 3C 1;ε ∑ h h  e  E E A Th


Fð∇I uÞj2 dx: Applying the same steps for T5, we get  0



2 Z   ½uD 1 ½I uD   @
F T5 r ∑ hF  ds  C 2;ε e A E D h h e  Z þ 3C 1;ε ∑ jFð∇uh Þ
Fð∇I uÞj2 dx: ED A T h

ED



 að∇uh Þ∇uh
að∇I uÞ∇I u ½ϕ ds

Z


∑

e A ED e

ðað∇uh Þ∇uh
að∇I uÞ∇I uÞϕ ds





Z   ½luh  ½I u ½uh 
a ½I u ∇ϕ ds a
∑ h h e A EI e





Z ½uh D ½I uD ½uh D
a ½I uD ∇ϕ ds
∑ a h h e A ED e





Z ½uh  ½uh  ½I u ½I u
a ½ϕ ds a þ ∑ σ h h h h e e A EI





Z ½uh D ½uh D ½I uD ½I uD
a ϕ ds: þ ∑ σ a h h h h e e A ED

L ðeÞ

The term T4 can be bounded working in the same way as in Lemma 6.5,  

( )

Z  ½u ½u ½I u ½I u 1=2   1=2 
a h a T4 r ∑ h ∇ϕ ds    h h h h  e A EI e 0  1 1=2 2

Z

1=2 Z 

 ½u ½u ½I u ½I u B
a hfj∇ϕjg2 ds r ∑ @ ha  dsA h h h h  e  e e A EI

ð6:10Þ

Using (2.12a), Lemmas 6.5 and 6.6, we have Z jFð∇uh Þ
Fð∇I uÞj2 dx Bðuh ; uh
I uÞ
BðI u; uh
I uÞ Z C p ∑ E A Th E


C 1;ε ∑ ‖Fð∇uh Þ
Fð∇I uÞ‖2L2 ðEÞ
E A Th





h ½uh  ½I u 2
F J L2 ðeÞ ∑ JF C 2;ε e A E I h h





h ½uh D ½I uÞD 2
F J L2 ðeÞ
C 1;ε ∑ ‖Fð∇uh Þ ∑ JF C 2;ε e A E D h h E A Th



½uh  ½I u 2
F J L2 ðeÞ
Fð∇I uÞ‖2L2 ðEÞ þ σ ∑ h J F h h e A EI



½uh D ½I uD 2
F J L2 ðeÞ : þ σ ∑ hJ F h h e A ED


Gathering the bounds and choosing appropriately the constants C i;ε and σ, for example 2C 1;ε o C p =2 and σ
1=C 2;ε 4 12, we can find κ0 such that the relation (6.10) to be true. □

I. Toulopoulos / Finite Elements in Analysis and Design 95 (2015) 42–50 2

1

1

0

0

y

y

2

-1

-2

49

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9

-1

-2

-1

0

1

2

-2

-2

-1

0

1

2

x

x

Fig. 1. Left: The domain Ω with the T h2 mesh. Right: The contours of uh computed for the p ¼ 2.3 test.

Applying (3.6) on the first term on the right-hand side of (6.14) and then using discrete embeddings (6.5b) yields

Table 1 CPUs for the two Picard iterative methods. – T hi

PEJ CPU for p ¼2.3

PEGS

i¼0 i¼1 i¼2 i¼3

21.6263 73.4690 195.4758 712.597

13.5540 56.5072 186.6697 625.4143

1∂ 1 κ0 ‖ϕ‖2L2 ðΩÞ þ κ0 ‖ϕ‖2F;DG r ‖∂t ðu
I uÞ‖2L2 ðΩÞ þ ‖ϕ‖2F;DG : 2 ∂t 4κ 0 4 1 2 2 2 ‖u
I u‖F;DG þ C 2;ε ‖ϕ‖F;DG þC 3;ε h ‖∇Fð∇uÞ‖L2 ðΩÞ : C 1;ε

Next, using (6.4) and choosing C 2;ε ¼ κ 0 =4 into (6.15), we have

Next, we give the estimate for the approximation error uh
u. Theorem 6.9. Under the assumptions (2.3) and (6.2), and choosing uh ð0Þ≔I u0 , there exist constants κ0 and C 4 0 such that for t A ð0; T Z κ0 t ‖uðτÞ
uh ðτÞ‖2F;DG dτ r ‖uðtÞ
I uðtÞ‖2L2 ðΩÞ ‖uðtÞ
uh ðtÞ‖2L2 ðΩÞ þ 2 0 Z t þC ‖∂t ðuðτÞ
I uðτÞÞ‖2L2 ðΩÞ dτ 0

þ Ch

2

Z 0

t

‖∇Fð∇uðτÞÞ‖L2 ðΩÞ dτ:

ð6:11Þ

Ω

Setting ϕ ¼ uh
I u and adding
sides of (6.12), we get

R

Ω ∂t I uϕ

dx
BðI u; ϕÞ on both

Z

Z Ω

∂t ðϕÞϕ dx þ Bðuh ; ϕÞ
BðI u; ϕÞ ¼

Ω

∂t ðu
I uÞϕ dx þ Bðu; ϕÞ
BðI u; ϕÞ:

ð6:13Þ Now, we use (6.6) in (6.9) and then we use the derived result on the right-hand side of (6.13). Consequently, we make use of (6.10) to the left-hand side of (6.13), we eventually end up with the following relation: Z 1∂ ‖ϕ‖2L2 ðΩÞ þ κ 0 ‖ϕ‖2F;DG r ∂t ðu
I uÞϕ dx 2 ∂t Ω þ

1∂ κ0 ‖u
I u‖2L2 ðΩÞ þ ‖uh
I u‖2F;DG 2 ∂t h 2 1 2 2 r ‖∂t ðu
I uÞ‖L2 ðΩÞ þ C ε h ‖∇Fð∇uÞ‖L2 ðΩÞ : 4κ0

ð6:16Þ

We integrate (6.16) from 0 to t: Z t ‖uh ðtÞ
I uðtÞ‖2L2 ðΩÞ þ κ0 ‖uh ðτÞ
I uðτÞ‖2F;DG dτ 0

Z r ‖u0;h
I u0 ‖2L2 ðΩÞ þ

t



0

1 ‖∂t ðuðτÞ
I uðτÞÞ‖2L2 ðΩÞ 2κ 0

 2 þ C ε h ‖∇Fð∇uðτÞÞ‖L2 ðΩÞ dτ:

ð6:17Þ

Observing that uh ð0Þ
I u0 ¼ 0 and applying the triangle inequality ‖uh ðtÞ
uðtÞ‖n r ‖uh ðtÞ
I uðtÞ‖n þ ‖uðtÞ
I uðtÞ‖n ;

Proof. We have by variational formulations (3.10) and (3.11) for t 4 0 that Z Z ∂t uh ϕ dx þ Bðuh ; ϕÞ ¼ ∂t uϕ dx þ Bðu; ϕÞ; 8 ϕ A V kh ðT h Þ: ð6:12Þ Ω

ð6:15Þ

1 2 ‖u
I u‖2F;DG þ C 2;ε ‖ϕ‖2F;DG þ C 3;ε h ‖∇Fð∇uÞ‖L2 ðΩÞ : C 1;ε

ð6:14Þ

in (6.17), we can deduce the estimate (6.11).

□

Using further the estimates (6.5a) in (6.11), we prove the following corollary. Corollary 6.10. Under the assumptions of Theorem 6.9, there is a C≔Cð‖∇Fð∇uðtÞÞ‖2L2 ð0;T;L2 ðΩÞÞ ; ‖∇Fð∇uðtÞÞ‖2L2 ðΩÞ ‖∇ut ‖2L2 ð0;T;L2 ðΩÞÞ Þ such that Z t 2 ‖uðτÞ
uh ðτÞ‖2F;DG dτ rCh : ð6:18Þ 0

Proof. The assertion follows by the application of the interpolation estimate of Lemma 6.1 and the estimate (6.5a) on the terms on the right-hand side of (6.11). □ 7. Numerical examples In this section, we present numerical results to illustrate the performance of the proposed IPDG method for solving problem (2.1) and to verify the theoretical results of the previous section. The numerical examples have been performed for p ¼ 2:3; p ¼ 2:5; p ¼ 3,

50

I. Toulopoulos / Finite Elements in Analysis and Design 95 (2015) 42–50

Table 2 Convergence rates for the three p-test cases. – T hi

p ¼2.3 Rates r

p¼ 2.5

p¼3

i¼ 0 i¼ 1 i¼ 2 i¼ 3

– 2.12 2.10 2.04

– 2.30 2.06 2.02

– 2.02 2.05 2.02

using μ¼1, σ¼2.5 (see (3.8)). The Picard iterative procedure was stopped until the tolerance value satisfied by tol r 1:E
07. The domain is Ω≔½
2; 2  ½
2; 2, where Γ D ¼ ∂Ω and the data f ; uD of (2.1) are specified so that the exact solution is uðx; y; tÞ ¼ BðtÞ sin ðx þ yÞ;

ð7:1Þ

where BðtÞ ¼ 1 þ expð
100tÞ. The initial unstructured mesh T h0 is generated by a triangular mesh generator with h0 ¼ 1 and the next finer meshes T hi are obtained by subdividing the triangles into four equal triangles, hi þ 1 ¼ hi =2. The problem has been solved up to final time T ¼0.5 using a second order, 1-stage DIRK method, [29]. In Fig. 1 left, the T h2 mesh of the domain Ω is presented and in Fig. 1 right, we plot the uh solution computed on T h2 mesh for the p ¼2.3 test case. In the first numerical test, the CPU time of the iterative methods PEJ and PEGS is compared. In Table 1, the CPU time for the p¼ 2.3 test case is given. As it was expected, for the same value of tol, PEGS performs faster and appears to be more efficient than the PEJ iterative method. Next, we give examples for the convergence rate of the error, Z t ehi ¼ ‖uðτÞ
uhi ðτÞ‖2F;DG dτ; ð7:2Þ 0

where uhi is the IPDG solution and u is the solution (7.1). All the numerical tests have been performed using the PEGS method with Δt o ðhi =10Þp . The numerical convergence rates r are computed by the formula r ¼ lnðehi =ehi þ 1 Þ=lnð2Þ. The results are shown in Table 2. We can observe that for all p-test cases the error (7.2) converges with the rate that has been predicted in Corollary 6.10.

8. Conclusions In this work, an IPDG method was presented for approximating the solution of a quasilinear parabolic problem formulated in Lpsetting. The resulting nonlinear ODE system was discretized in time by s-DIRK methods applying two low-storage Picard iterative schemes for solving the resulting nonlinear systems. A stability bound was shown in the broken ‖  ‖DG;p
norm for the IPDG solution. Optimal error estimates for the IPDG method were proved in the broken ‖  ‖F;DG
norm for the case of p 4 2. The theoretical results were validated by numerical tests for several values of p 4 2. Acknowledgments This work was supported by Austrian Science Fund (FWF) under the Grant NFN S117-03. References [1] T. Roubíček, Nonlinear partial differential equations with applications, ISNM, International Series of Numerical Mathematics, vol. 153, Birkhäuser Verlag, Basel, 2005.

[2] R. Glowinski, J. Rappaz, Approximation of a nonlinear elliptic problem arising in a non-Newtonian fluid flow model in glaciology, Math. Model. Numer. Anal. 37 (2003) 175–186. [3] M. Picasso, J. Rappaz, A. Reist, M. Funk, H. Blatter, Numerical simulation of the motion of a two dimensional glacier, Int. J. Numer. Methods Eng. 60 (2004) 995–1009. [4] D.N. Arnold, An interior penalty finite element method with discontinuous elements, SIAM J. Numer. Anal. 19 (1982) 742–760. [5] M.F. Wheeler, An elliptic collocation-finite element method with interior penalties, SIAM J. Numer. Anal. 15 (1978) 152–161. [6] B. Riviére, M. Wheeler, V. Girault, A priori error estimates for finite element methods based on discontinuous approximation spaces for elliptic problems, SIAM J. Numer. Anal. 39 (2) (2001) 902–931. [7] D.A. DiPietro, A. Ern, Mathematical Aspects of Discontinuous Galerkin Methods, Mathematiques et Applications, vol. 69, Springer-Verlag, Berlin, Heidelberg, 2012. [8] B. Riviére, Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations, SIAM, Society for Industrial and Applied Mathematics, Philadelphia, 2008. [9] C. Ortner, E. Süli, Discontinuous Galerkin finite element approximation of nonlinear second order elliptic and hyperbolic systems, SIAM J. Numer. Anal. 45 (2007) 1370–1397. [10] B. Riviére, S. Shaw, Discontinuous Galerkin finite element approximation of nonlinear non-Fickian diffusion in viscoelastic polymers, SIAM J. Numer. Anal. 44 (2006) 2650–2670. [11] C. Bi, Y. Lin, Discontinuous Galerkin method for monotone nonlinear elliptic problems, Int. J. Numer. Anal. Model. 4 (2012) 999–1024. [12] P. Huston, J. Robson, E. Suli, Discontinuous Galerkin finite element approximation of quasilinear elliptic boundary value problems I: the scalar case, IMA J. Numer. Anal. 25 (2005) 726–749. [13] D.N. Arnold, An interior penalty finite element method with discontinuous elements, SIAM J. Numer. Anal. 19 (4) (1982) 175–186. [14] M. Ohm, Error estimate for parabolic problem by backward Euler discontinuous Galerkin method, Int. J. Appl. Math. 22 (2) (2009) 117–128. [15] B. Riviére, M. Wheeler, A discontinuous Galerkin method applied to nonlinear parabolic equations, in: B. Cockburn, G. Karaniadakis, C.-W. Shu (Eds.), Discontinuous Galerkin Methods: Theory, Computation and Applications, Lecture Notes in Computer Science and Engineering, vol. 11, Springer, Berlin, 2000, pp. 231–244. [16] M. Ohm, H. Lee, J. Shin, Error estimates for discontinuous Galerkin method for nonlinear parabolic equations, J. Math. Anal. Appl. 315 (2006) 132–143. [17] L. Song, G.-M. Gie, M.-C. Shiue, Interior penalty discontinuous Galerkin methods with implicit time-integration techniques for nonlinear parabolic equations, Numer. Methods Partial Differ. Equ. 29 (4) (2013) 1341–1366. [18] E. Burman, D.A. DiPietro, Discontinuous Galerkin approximations with discrete variational principle for the nonlinear Laplacian, C. R. Acad. Sci. Paris Ser. I 346 (2008) 1013–1016. [19] L. Diening, C. Kreuzer, Linear convergence of an adaptive finite element for the p-Laplacian equation, SIAM J. Numer. Anal. 46 (2) (2008) 614–638. [20] L. Diening, M. Růžička, Interpolation operators in Orliz–Sobolev spaces, Numer. Math. 107 (2007) 107–129. [21] L. Diening, C. Ebmeyer, M. Růžička, Optimal convergence for the implicit space–time discretization of parabolic systems with p-structure, SIAM J. Numer. Anal. 45 (2) (2007) 457–472. [22] D. Kröner, M. Růžička, I. Toulopoulos, Numerical solutions of systems with (p,δ)structure using local discontinuous Galerkin finite element methods, Int. J. Numer. Methods Fluids, 76 (2014) 855–874, http://dx.doi.org/10.1002/fld.3955. [23] T. Gudi, N. Nataraj, A.K. Pani, An hp-local discontinuous Galerkin method for some quasilinear elliptic boundary value problems of nonmonotone type, Math. Comput. 77 (2008) 731–756. [24] E. DiBenedetto, Degenerate Parabolic Equations, Springer-Verlag, New York, 1993. [25] J.W. Barrett, W.B. Liu, Finite element approximation of the parabolic pLaplacian, SIAM J. Numer. Anal. 31 (2) (1994) 413–428. [26] G. Franzina, Existence, uniqueness, optimization and stability for low eigenvalues of some nonlinear operators (Ph.D. thesis), Athesina Stadiorum Universitas, 2012. [27] D.A. DiPietro, A. Ern, Discrete functional analysis tools for discontinuous Galerkin methods with application to the incompressible Navier–Stokes equations, Math. Comput. 79 (271) (2010) 1303–1330. [28] S.C. Brenner, Poincare–Friedrichs inequalities for piecewise H1 functions, SIAM J. Numer. Anal. 41 (2003) 306–324. [29] A. Roger, Diagonally implicit Runge–Kutta methods for stiff O.D.E.'S, SIAM J. Numer. Anal. 14 (6) (1977) 1006–1021. [30] I. Ly, An iterative method for solving cauchy problems for the p-Laplace operator, Complex Variables Elliptic Equ. 55 (11) (2010) 1079–1088. [31] Y.Q. Huang, L. Ruo, L. Wenbin, Preconditioned descent algorithms for pLaplacian, J. Sci. Comput. 32 (2) (2007) 343–371. [32] L. Diening, D. Kröner, M. Růžička, I. Toulopoulos, A local discontinuous Galerkin approximation for systems with p-structure, IMA J. Numer. Anal. doi:http://dx.doi.org/10.1093/imanum/drt040. [33] L.R. Scott, S. Zhang, Finite element interpolation of non-smooth functions satisfying boundary conditions, Math. Comput. 54 (190) (1990) 483–493.