Applied Mathematics Letters 37 (2014) 77–81
Contents lists available at ScienceDirect
Applied Mathematics Letters journal homepage: www.elsevier.com/locate/aml
Uniqueness of numerical solutions to nonlinear parabolic equations by a fully implicit discontinuous Galerkin method Lunji Song a,b,∗ , Jinhu Ma a a
School of Mathematics and Statistics, Lanzhou University, Lanzhou 730000, PR China
b
Key Laboratory of Applied Mathematics and Complex Systems, Gansu Province, PR China
article
abstract
info
Article history: Received 15 April 2014 Accepted 31 May 2014 Available online 12 June 2014
We prove uniqueness of numerical solutions to nonlinear parabolic equations approximated by a fully implicit interior penalty discontinuous Galerkin (IPDG) method, with a mesh-independent constraint on time step. © 2014 Elsevier Ltd. All rights reserved.
Keywords: Discontinuous Galerkin method Uniqueness Nonlinear parabolic equation Fixed point
1. Introduction To describe the diffusion of a chemical species of the concentration u in a porous medium with a source term f (x, u), one can consider the following nonlinear parabolic equation ut − ∇ · (a(x, u)∇ u) = f (x, t ), a(x, u)∇ u · n = 0, u|t =0 = ψ(x),
in Ω × (0, T ),
in ∂ Ω × (0, T ),
on Ω × {0},
(1.1) (1.2) (1.3)
where Ω is an open interval in R1 , or a convex polygonal domain in R2 , n is the unit outward normal vector to the boundary ∂ Ω , a(x, u) is a bounded and positive function and T > 0 is arbitrary but fixed. We assume that f (x, t ) is uniformly Lipschitz ¯ × R), where C 2 (Ω ¯ × R) is the class of twice continuously continuous with respect to its variables, and a(x, v) ∈ C 2 (Ω ¯ × R. In this work, u is assumed to be a strong solution such that u ∈ C 2 (Ω × [0, T ]), a solution differentiable functions on Ω of the problem (1.1)–(1.3). Many contributions in the literature show that discontinuous Galerkin (DG) methods are a powerful simulation tool for solving partial differential equations (see e.g. [1–4]). However, in the literature, there are seldom uniqueness results on fully discrete implicit discontinuous Galerkin approximations of nonlinear parabolic problems. Because the problem is nonlinear, stability analysis of fully implicit interior penalty discontinuous Galerkin (IPDG) schemes cannot be directly applied to prove uniqueness of numerical solutions. So far there has been seldom theoretical analysis on uniqueness of numerical solutions to fully discrete implicit IPDG schemes for the nonlinear parabolic problems. Fortunately, Gudi and Pani [5] proved the existence and uniqueness of semi-discrete solutions for a class of elliptic equations by using Brouwer’s fixed point theorem. Hence, motivated by the work in [5], we would like to fill the gap on the uniqueness of numerical solutions for nonlinear parabolic equations by the IPDG approximations.
∗
Corresponding author at: School of Mathematics and Statistics, Lanzhou University, Lanzhou 730000, PR China. Tel.: +86 13609302099. E-mail address:
[email protected] (L. Song).
http://dx.doi.org/10.1016/j.aml.2014.05.016 0893-9659/© 2014 Elsevier Ltd. All rights reserved.
78
L. Song, J. Ma / Applied Mathematics Letters 37 (2014) 77–81
2. The discontinuous Galerkin method Let Th := {E1 , E2 , . . . , ENh } be a shape regular finite element subdivision of the domain Ω , where Ei is an interval in 1D or a triangle in 2D and Nh is the number of all elements. Here h denotes the maximal diameter of all elements. Note that ¯ = Ni=h1 E¯i , for each h > 0. For each element of Th , hi is the diameter of Ei , 1 ≤ i ≤ Nh . We introduce the set of all edges Ω M
P
by Fh including the set of interior edges {ei }i=h 1 and the set of boundary edges {ei }i=hPh +1 . On each edge ei (1 ≤ i ≤ Mh ) of Th , we fix a unit normal vector ni on ei by
ni =
the unit normal vector on ei , pointing from Ek to Ej , if ei = ∂ Ek ∩ ∂ Ej , the unit normal vector on ei , pointing outward of Ω , if ei ⊂ ∂ Ω ,
then we denote the average and jump operators as
{v} :=
v|Ej + v|Ek
[v] := v|E − v|E .
,
j k 2 We also set the average and jump on the boundary edge ek = ∂ Ei ∩ ∂ Ω as {v} := v|Ei ∩∂ Ω , [v] := v|Ei ∩∂ Ω . Along this article, the H m Sobolev norm on ω is defined by
∥ · ∥m,ω := | · |H m (ω) ,
∀ 0 ≤ m < ∞, ∀ ω ⊂ R1 (or R2 ).
(2.1)
Note that, by default, H (ω) denotes L (ω) and ∥ · ∥∞, ω is the standard L -norm on ω. Then using (2.1), we introduce the broken Sobolev space for any real number s: 0
∞
2
H s (Th ) = v ∈ L2 (Ω ) : v|Ei ∈ H s (Ei ), 1 ≤ i ≤ Nh ,
which is equipped with the broken space norm: Nh
|v|2H s (Th ) =
∥v∥2s,Ei .
i=1
The discontinuous finite element space is defined as
Dp (Th ) = v ∈ L2 (Ω ) : v|Ei ∈ Ppi (Ei ), 1 ≤ i ≤ Nh ,
(2.2)
where p = min{pi ≥ 1} for any 1 ≤ i ≤ Nh , and Ppi (Ei ) is the space of polynomials of (total) degree at most pi on Ei . An hp-quasi-uniformity assumption is imposed on the mesh Th , namely, there exists a positive constant CQ , independent of h and p, such that
max 1≤i≤Nh
pi
hi
max
hi
1≤i≤Nh
≤ CQ .
pi
We introduce the interior penalty term J0σ : Dp (Th ) × Dp (Th ) → R in the form: J0σ (v, w) =
Ph
σk
k=1
p2k
|ek |
[v][w]ds,
(2.3)
ek
which penalizes the jump of the functions across the edge ek , 1 ≤ k ≤ Ph . Here the penalty parameter σk is a positive real number and |ek | is the Lebesgue measure of the edge ek . We also define the energy norms on Dp (Th ): ∀ v ∈ Dp (Th ), 2
|||v||| =
Nh
∥∇v∥
2 0,Ek
σ
+ J0 (v, v),
k=1
2 Ph ∂v |ek | |||v|||+ = |||v||| + ds. 2 ∂ν p ek k k=1 2
2
Aiming to study the strong solution u ∈ C 2 (Ω × [0, T ]) of the problem (1.1)–(1.3), we proceed element by element and obtain the following consistent weak formulation of (1.1)–(1.3): Find u(t ) ∈ H s (Th ), s ≥ 2, such that
(ut , v) + Aϵ (u; u, v) = (f , v), ∀ v ∈ H s (Th ), (u(0), v) = (ψ, v), where Aϵ (ρ; v, w) is bilinear in the last two terms: Nh Ph Aϵ (ρ; v, w) = a(ρ)∇v · ∇w dx − {a(ρ)∇v · nk }[w]ds j =1
+ϵ
Ej
k=1
Ph k=1
(2.4) (2.5)
ek
{a(ρ)∇w · nk }[v]ds + J0σ (v, w),
∀ v, w ∈ H s (Th ),
(2.6)
ek
where the parameter ϵ in Aϵ may be taken the values −1, 0 or 1. Note that σk is chosen to satisfy the coercivity of Aϵ (see (3.11) of [6]).
L. Song, J. Ma / Applied Mathematics Letters 37 (2014) 77–81
79
To fully discrete (2.4), we use the discontinuous Galerkin method for the spatial variable and the Euler backward method for the time variable. Then a fully discrete IPDG formulation is to seek a sequence {unh }Nn≥0 of functions in Dp (Th ) such that ∀ n ≥ 0,
uhn+1 − unh
△t
, v + Aϵ (unh+1 ; unh+1 , v) = (f (tn+1 ), v),
∀ v ∈ Dp (Th ),
(2.7)
u0h = u0h ,
(2.8)
where u0h is taken as a L projection of ψ onto Dp (Th ). The main concern is about the existence and uniqueness of the numerical solution to (2.7)–(2.8). 2
3. Uniqueness of the numerical solutions in the DG schemes For a given unN at each time step, we will prove the existence and uniqueness of unN+1 of the fully implicit IPDG scheme (2.7)–(2.8). For simplification, we write a(x, w) := a(w). By Taylor series expansion, we obtain a(w) = a(u) + a˜ u (w)(w − u) = a(u) + au (w − u) + a˜ uu (w)(w − u)2 , where 1
a˜ u (w) =
a˜ uu (w) =
au (w + t (w − u))dt ,
1
0
(1 − t )auu (w + t (w − u))dt . 0
¯ × R), and auu ∈ C 0 (Ω ¯ × R), it is observed that a˜ u , a˜ uu ∈ L∞ (Ω ¯ × R). Thanks to the assumptions that au ∈ C 1 (Ω At time t = tn+1 , subtracting (2.4) and its fully discrete variational formula (2.7) with v = vh gives
∂ u unh+1 − unh − , vh + Aϵ (u; u, vh ) = Aϵ (uhn+1 ; unh+1 , vh ), ∂ t t =tn+1 △t
∀ vh ∈ Dp (Th ),
(3.1)
equivalently,
unh+1 − unh ∂ u − , vh + Aϵ (u; u − unh+1 , vh ) = Aϵ (unh+1 ; uhn+1 , vh ) − Aϵ (u; uhn+1 , vh ). ∂ t t =tn+1 △t
(3.2)
Since u is continuous on each edge ek , the right-hand side of the resultant equation above can be written as Aϵ (uhn+1 ; uhn+1 , vh ) − Aϵ (u; unh+1 , vh )
=
Nh Ej
j =1
+ϵ
(a(unh+1 ) − a(u))∇ unh+1 · ∇vh dx −
k=1
Ph k=1
=
Nh
{(a(unh+1 ) − a(u))∇vh · nk }[unh+1 ]ds
(a(unh+1 ) − a(u))∇(unh+1 − u) · ∇vh dx −
+ϵ
{(a(
unh+1
) − a(u))(∇
unh+1
{(a(unh+1 ) − a(u))∇ uhn+1 · nk }[vh ]ds ek
− ∇ u) · nk }[vh ]ds +
ek
Ph k=1
Ph k=1
Ph k=1
{(a(uhn+1 ) − a(u))∇ uhn+1 · nk }[vh ]ds
ek
ek
Ej
j =1
−
Ph
Nh j =1
(a(uhn+1 ) − a(u))∇ u · ∇vh dx Ej
{(a(unh+1 ) − a(u))∇vh · nk }[unh+1 − u]ds.
(3.3)
ek
Then, adding the following terms
−
Nh j =1
Ej
au (u)(unh+1 − u)∇ u · ∇vh dx +
Ph k=1
ek
au (u)(unh+1 − u)
∂u [vh ]ds, ∂ nk
to the both sides of (3.2) and splitting e := u − unh+1 = u − Ih u + Ih u − uhn+1 , we get another form of Eq. (3.2)
∂ u unh+1 − unh − , vh + Aˆ ϵ (u; u − unh+1 , vh ) = F (uhn+1 ; u − uhn+1 , vh ), ∂ t t =tn+1 △t
(3.4)
80
L. Song, J. Ma / Applied Mathematics Letters 37 (2014) 77–81
where Aˆ ϵ (ρ; w, v) = A(ρ; w, v) +
Nh
au (ρ)w∇ρ · ∇v dx − Ej
j =1
Ph k=1
au (ρ)
ek
∂ρ w [v]ds, ∂ nk
Nh F (unh+1 ; e, vh ) = (a(u) − a(unh+1 ))∇ e · ∇vh + au (u)e∇ u · ∇vh − (a(u) − a(uhn+1 ))∇ u · ∇vh dx Ej
j =1 Ph
∂e ∂u ∂u [vh ] + (a(u) − a(unh+1 )) [vh ] − au (u)e [vh ]ds ∂ nk ∂ nk ∂ nk k=1 ek Ph ∂v (a(u) − a(unh+1 )) +ϵ [e]ds. ∂ nk k=1 ek +
(a(unh+1 ) − a(u))
By Taylor series expansion, we write F (unh+1 ; e, vh ) =
Nh j =1
−
a˜ u (unh+1 )e∇ e · ∇vh + a˜ uu (unh+1 )e2 ∇ u · ∇vh dx
Ej
Ph
˜ (
au unh+1
k=1
ek
Ph ∂e ∂v n+1 2 ∂ u n +1 )e [vh ] + a˜ uu (uh )e [vh ]ds + ϵ [e]ds. a˜ u (uh )e ∂ nk ∂ nk ∂ nk k=1 ek
We see that e = u(tn+1 ) − unh+1 := η − χ , where η = u(tn+1 ) − Ih u(tn+1 ), χ = uhn+1 − Ih u(tn+1 ) and Ih u satisfies the local approximation properties; see Lemma 4.5 in [7]. Analogously to the proofs of Lemmas 4.2, 4.3 and Theorem 4.4 in [5], the following lemmas still hold. Lemma 3.1. Let uh , vh ∈ Dp (Th ). There exists a constant C > 0 which is independent of h and p such that
|F (uh ; u − uh , vh )| ≤ CCa
pi
max 1≤i≤Nh
1/2
1/2
2
|||χ||| + Cu h
hi
|||χ||| + |||η||| |||vh |||,
(3.5)
where Ca = max{∥˜au ∥L∞ (Ω ×R) , ∥˜auu ∥L∞ (Ω ×R) }, and Cu = max{∥u∥H 2 (Ω ) , ∥u∥H 1 (Ω ) |u|W 1,∞ (Ω ) }. Definition: For a given z ∈ Dp (Th ), let Sh : Dp (Th ) → Dp (Th ) be a map y = Sh z ∈ Dp (Th ) satisfying the stationary equation Aˆ ϵ (u; Ih u − y, vh ) = Aˆ ϵ (u; Ih u − u, vh ) + F (z ; u − z , vh ),
∀ vh ∈ Dp (Th ).
(3.6)
Lemma 3.2. Let β ≥ 1 and z ∈ Dp (Th ). Set y = Sh z. Then there exists a positive constant C which is independent of h and p such that
|||Ih u − y||| ≤ CCa
max 1≤i≤Nh
pi hi
1/2
2
|||Ih u − z ||| + Cu h
1/2
|||Ih u − z ||| + CCa 1 + Cu h1/2 |||Ih u − u|||+ .
Lemma 3.3. Set the ball Bδ (Ih u) = {v ∈ Dp (Th ) : |||Ih u − v||| ≤ δ} of radius δ . For some 0 < θ ≤ that
max 1≤i≤Nh
pi hi
1/2
1
|||Ih u − z |||2 ≤ CCQ Cu h 2 −2θ |||Ih u − u|||+ ,
1 4
(3.7)
and z ∈ Bδ (Ih u), it holds
(3.8)
where C is independent of h and p. Lemma 3.4. For sufficiently small h, there is a radius δ = map Sh maps Bδ (Ih u) into itself.
1 hθ
|||Ih u − u|||+ (0 < θ ≤ 41 ) and a positive constant C such that the
We also need to use the following Gårding-type inequality Aˆ ϵ (u; φ, vh ) ≥ C1 |||vh |||2 − C2 ∥φ∥2L2 (Ω ) , and the continuity of Aˆ ϵ Aˆ ϵ (u; φ, vh ) ≤ C |||φ||||||vh |||,
(3.9)
L. Song, J. Ma / Applied Mathematics Letters 37 (2014) 77–81
81
where C1 , C2 and C are positive constants independent of h and p. The proofs of Lemmas 3.1–3.4 are a parallel generalization as them appearing in [5]. Now we give a critical theorem for the existence and uniqueness of unh+1 in (3.4) as follows. Theorem 3.1. Assume that 0 < △t <
1 . C2
For sufficiently small h, at each time step, there is a radius δ =
(0 < θ ≤ ) and a positive constant C such that the following inequality holds for any given z1 , z2 ∈ Bδ (Ih u), 1 4
|||Sh z1 − Sh z2 ||| ≤ CCa CQ Cu hθ |||z1 − z2 |||.
1 hθ
|||Ih u − u|||+
(3.10)
Proof of Theorem 3.1. Setting y1 = Sh z1 and y2 = Sh z2 , for z1 , z2 ∈ Bδ (Ih u) and a fixed t, we have from (3.4)
∂ u y1 − unh − , vh + Aˆ ϵ (u; u − y1 , vh ) = F (z1 ; u − z1 , vh ), ∂ t t =tn+1 △t ∂u y2 − unh − , vh + Aˆ ϵ (u; u − y2 , vh ) = F (z2 ; u − z2 , vh ). ∂ t t =tn+1 △t
Subtracting these above identities gives
y1 − y2
△t
, vh + Aˆ ϵ (u; y2 − y1 , vh ) = F (z1 ; u − z1 , vh ) − F (z2 ; u − z2 , vh ).
(3.11)
Applying (3.9) and taking vh = y2 − y1 in (3.11), it is bounded below
−C 2 +
1
△t
∥y2 − y1 ∥2L2 (Ω ) + C1 |||y2 − y1 |||2 ≤ F (z1 ; u − z1 , y2 − y1 ) − F (z2 ; u − z2 , y2 − y1 ).
(3.12)
Setting Z = z1 − z2 , we bound the terms on the right hand side of (3.12) by Lemma 3.1
|F (z1 ; u − z1 , y2 − y1 ) − F (z2 ; u − z2 , y2 − y1 )| 1/2 pi ≤ CCa max |||Z |||2 + |||z1 − Ih u||||||Z ||| + |||Ih u − z2 ||||||Z ||| |||vh ||| 1≤i≤Nh
≤ CCa
max 1≤i≤Nh
hi pi
1/2
hi
(|||Z ||| + |||z1 − Ih u||| + |||Ih u − z2 |||) |||Z ||||||vh |||
1
≤ CCa CQ Cu h 2 −θ |||Z ||||||vh |||,
(3.13)
where we have used the inequality
|||Z ||| + |||z1 − Ih u||| + |||Ih u − z2 ||| ≤ 2(|||z1 − Ih u||| + |||Ih u − z2 |||) and Lemma 3.3 for the last step. As δ = h1θ |||η|||+ , 0 < θ < 41 and the time step satisfies that 0 < △t <
1 , C2
we conclude from (3.12) and (3.13) that
θ
|||y2 − y1 ||| ≤ CCa CQ Cu h |||z1 − z2 |||. This completes the proof.
For sufficiently small h, we deduce from Theorem 3.1 that there is a δ > 0 such that the map Sh : Bδ (Ih u) → Bδ (Ih u) is continuous. Hence, using Lemma 3.4 and Theorem 3.1, for small h, we conclude by Schauder’s fixed point theorem that there exists a uhn+1 ∈ Bδ (Ih u) such that Sh unh+1 = unh+1 . It is clear that uhn+1 is a unique fixed point of Sh . We conclude that at each time step there exists a unique solution to the fully discrete problem. Acknowledgment This work was supported by the National Natural Science Foundation of China (grant number 11101196). References [1] D.N. Arnold, An interior penalty finite element method with discontinuous elements, SIAM J. Numer. Anal. 19 (1982) 724–760. [2] B. Cockburn, J. Guzmán, H. Wang, Superconvergent discontinuous Galerkin methods for second-order elliptic problems, Math. Comp. 78 (265) (2009) 1–24. [3] S. Kaya, B. Riviére, A discontinuous subgrid eddy viscosity method for the time-dependent Navier–Stokes equations, SIAM J. Numer. Anal. 43 (4) (2005) 1572–1595. [4] S. Sun, M.F. Wheeler, Symmetric and nonsymmetric discontinuous Galerkin methods for reactive transport in porous media, SIAM J. Numer. Anal. 43 (1) (2005) 195–219. [5] T. Gudi, A.K. Pani, Discontinuous Galerkin methods for quasilinear elliptic problems of nonmonotone type, SIAM J. Numer. Anal. 45 (1) (2007) 163–192. [6] 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 Differential Equations 29 (4) (2013) 1341–1366. [7] I. Babuška, M. Suri, The h-p version of the finite element method with quasiunform meshes, RAIRO Modél. Math. Anal. Numeŕ. 21 (1987) 199–238.