Journal of Applied Sciences

Year: 2017 | Volume: 17 | Issue: 2 | Page No.: 81-89
DOI: 10.3923/jas.2017.81.89

A Discontinuous Galerkin Method for the Wave Equation: A hp-a Priori Error Estimate

Abdelhamid Zaghdani, M. Ezzat Mohamed and A.I. El-Maghrabi

Abstract: Background: The discontinuous Galerkin method for the approximation of a partial differential equation solution has some advantages comparing to the classical finite element method. Objective: This study aimed to provide a numerical approximation of the wave equation solution derived from Maxwell’s equations. Methodology: This study applied the discontinuous Galerkin method for approximating the electric field which is solution of a wave equation that derives from Maxwell’s equations in a tridimensional domain. Results: Some discrete inequalities on discontinuous spaces for Maxwell’s equations were presented and a discontinuous Galerkin method for the numerical approximation of the solution of the wave equation was analyzed. Its hp-analysis was carried out and error estimates that were optimal in the mesh size and slightly suboptimal in the approximation degree were obtained. The DG spatial discretization was augmented with the second order Newmark scheme in time and some numerical results were obtained. Conclusion: The results of the study can be applied for approximating solutions of several partial differential equations.

Keywords: DG method, a priori error estimate, wave equation and Maxwell`s equations

INTRODUCTION

In this study, we present and prove some discrete inequalities on discontinuous space for Maxwell's equation and we analyze a hp-discontinuous Galerkin method for the wave equation in stable medium with perfect electric conductor boundary. Approximate continuity is imposed by including penalty terms in the form which defines the method. This method was analyzed for advection-diffusion-reaction problems^1,2 and for the approximation of second order elliptic equations³. An interior penalty finite element method with discontinuous element has been introduced and analysed⁴. This method is inspired from the DG method with the addition of penalty terms. For the time-harmonic Maxwell's equations, a hp-DG version has been presented and analysed⁵. For the plane-wave discontinuous Galerkin method we refer to Atcheson⁶.

The problem considered for the most of this study is the initial-boundary value problem derived from Maxwell's equations in stable medium with perfect electric conductor boundary:

where, Ω is a convex polyhedron included in R³, I = [0, T]⊂R. We suppose that the functions μ, ε are sufficiently smooth and satisfies 0<ε_min<∣ε(x)∣<ε_max and 0<μ_min<∣μ(x)∣<μ_max for all x in Ω with a constants ε_min, ε_max, μ_min, μ_max. The u₀ and u₁ are in H₀(∇×, Ω)∩H(∇⋅, Ω), g∈L² (I, L² (Ω)), f is defined on Ω×I and in L² (I, L² (Ω)³). Physically, u is the electric field, f and g are related to a current and charge density, respectively. Moreover, μ₀ε₀c² = 1, where μ₀≈4π10^–7 H⋅m^–1 and ε₀≈(36π10⁹)^–1 F⋅m^–1 are the magnetic permeability and the electric permittivity in vacuum, respectively. In fact, many of the results proved are valid so long as Ω is a bounded, convex domain with Lipchitz, connected and simply connected boundary. If we assume that Ω is a stable medium with perfect electric conductor boundary and if u is the exact solution of Maxwell's equations, then u belong to H¹(Ω)³ as described by Duvaut and Lions⁷.

Let ∏_h be a partition of Ω into tetrahedra. We denote by the set of all interior faces, the set of exterior faces and F_h the set of all faces of the partition. For e∈F_h we denote by 〈,〉_e; i.e., the scalar product in L²(e)³ or L²(e) furthermore we identify to 〈,〉 the scalar product in L² (∂Ω)³ or L² (∂Ω). The scalar product in L² (Ω)³ or L² (∂Ω) is denoted by (,) and is identified to wwith (,)_K is the scalar product in L² (K)³ or L² (K). For the other spaces and notations used in this study we refer to Zaghdani⁸.

In order to define the average of ∇×u in the formulation Eq. 5, we set for

Finite element spacesLet p = (p_K)K∈∏_h be a degree vector that assigns to each element K∈∏_h a polynomial approximation order p_K≥1. The generic hp-finite element space of piecewise polynomials is given by:

where, is the space of real polynomials of degree at most p_K in K.

Now, we set:

We define the local parameters h and p as h = min (h_K, h_K’), p = max (p_K, p_K’) in the case of interior faces and h = h_K, p = p_K in the case of boundary faces^8,9. In the following of this study σ denotes a stabilization parameter and will be chosen depending on the local mesh size and polynomial degree. We consider the same definition of this parameter as by Perugia and Schotzau⁹ and Zaghdani⁸ and define with a strictly positive constant κ.

MATERIALS AND METHODS

Formulation for the Maxwell problem: In the notations of Zaghdani⁸, we can state the basic integration by parts formulas:

We have¹⁰:

and:

In order to derive a weak formulation of Eq. 1, we note that Eq. 3 implies for any u with:

where, we have denoted by:

Now, we introduce the penalty term via the form:

Where:

and:

We also define:

Now, since J₀ (u, v) = (g, ∇⋅(εv)) for the exact solution u of Eq. 1, then u satisfies:

Properties of the bilinear form
Mesh-dependent norm: We now, introduce the discontinuous Galerkin norm and set for u∈H¹ (∇×, ∏_h):

We start by studying the continuity of the bilinear forms introduced above. We have:

Proposition 1: ∀v, u∈H¹ (∇×, ∏_h),

with a constant C independent of h and p.

Proof: The proof can be easily obtained from the definition of A, J^σ, ∥⋅∥ and the Cauchy-Schwarz inequality.

In order to study the coercivity of the bilinear form, we start by introducing the following inequality of Poincarré-Friedrichs type valid for u∈H¹ (∏_h)³.

Lemma: Let u∈H¹ (∏_h)³ and σ the stabilization parameter defined previously, then there exists C independent of h and p such that:

Proof: Since Γ is simply connected we have the following orthogonal decomposition¹¹:

Therefore, for u∈H¹ (∏_h)³, u∈L² (Ω)³ and we can write u = u₁+u₂ with:

Since u₁ = H₀ (∇×0, Ω) we can write u₁ = ∇q with ¹¹. Also u₂ = ∇×φ in Ω avec φ∈H(∇×, Ω)∩H₀(∇⋅0, Ω)². In particular, the traces of φ are well defined. Note that the inequality Eq. 6 implies that:

Using the integration by parts Eq. 3 and 4, we obtain:

Since,

We can get:

It is clear that:

Since, φ∈H (∇×, Ω)∩H₀ (∇⋅, Ω) and ∇⋅φ = 0. We obtain for the first inequality¹¹:

Using the trace inequality⁸, we get for every face e∈F_h:

By adding:

Similarly and using the fact that the embedding of H (∇×, Ω)∩ H₀ (∇⋅, Ω) in H¹ (Ω)³ is continuous, we can estimate:

and obtain:

Therefore, we get:

which is equivalent to:

Remark: If we use the assumptions on μ and ε, the definition of the bilinear forms A and J^σ we deduce in particular that there exists a constant C independent of h and p such that:

Now, the following coercivity result on the discrete space Σ_h holds.

Proposition 2: There exists two constants α>0 and independent of h and p such that:

Proof: Let us first recall the following inverse inequality:

with a constant C_inv>0, only depending on the shape regularity of the mesh. For the two dimensional elements, the proof of Eq. 7 can be found by Schwab¹². For the three dimension space, the proof is analogous².

Now, let α be an arbitrary real number and choose v∈Σ_h. Then:

Since, {∇×v} is the average of the flux at the face of two elements K_i and K_j the corresponding integral can be split into two integrals with integrands (∇×v)_i/σ and (∇×v)_j/σ each one associated with the elements K_i or K_j, respectively. Therefore, let e∈F_h and consider the integral associated with the element K. Using the inverse inequality, we have since ∇×Σ_h⊂Σ_h:

So that, selecting σ to be equal to κp²/h in Eq. 8, we obtain:

In particular:

Now, from the definition of a (v, v) we can get for all ∈>0:

From the previous inequality and the definition of J^σ we obtain:

It then follows that:

The previous inequality is true for all ∈>0, taking we obtain:

We can choose κ sufficiently large and:

and obtain:

Now, the following hp-approximation result to interpolate scalar function holds.

Proposition 3: Let K∈∏_h and suppose that Then there exists a sequence of polynomials satisfying:

Furthermore, if t_k≥1:

The constant C is independent of u, h_K and p_K but depends on the shape regularity of the mesh and on

Proof: The assertion in Eq. 9 has been proved by Babuska and Suri¹³ (Lemma) for two-dimensional domains. For three-dimensional domains, the proof is analogous². The assertion in Eq. 10 has been proved by Perugia and Schotzau⁹.

In order to interpolate the vector functions, we define the following.

Definition: For u = (u₁, u₂, u₃) we define by with is defined by where, is given in proposition 3.

Model problem: If u is the exact solution of the Maxwell problem, then u satisfies:

The interior penalty finite element approximation to u is to find u^h:I→Σ_h such that:

Upon choice of a basis for Σ_h and the data f and g, Eq. 11 determines uh as the only solution to an initial value problem for a linear system of ordinary differential equations. Note that if u is the exact solution of Eq. 1, then u satisfies the first line in Eq. 11 and thus the problem is consistent.

We now analyze the proposed procedure by the method of energy estimates.

A priori error estimate: In this study, u denote the exact solution of Eq. 1 and u^h the discrete solution of Eq. 11. The C is generic constant independent of h and p which takes different values at the different place sand depends of μ_min, μ_max, ε_min, ε_max, T, α, the coercivity constants of the form B and Ω.

Let ζ = u^h-u, then ζ satisfies:

Decompose ζ as η-υ where, and

Note that [η]_N = [η]_T on and [η]_T = 0 on thus:

Since, υ_t (t)∈Σ_h we can set v = υ_t (t), obtaining:

So:

Since, υ_t(0) = υ(0) = 0, integration over [0, t]⊂I, yields:

The final term may be integrated by parts in time. Hence:

Therefore, we can apply the coercivity of B and continuity of A to get:

In particular:

As this holds for all t∈I, Gronwall's Lemma implies that:

Since, ζ = η-υ and J^σ (η, η) = 0:

Then, error bounds for the finite element approximation to the true solution reduce to the error bounds for the piecewise polynomial interpolant. Thus, we start by estimating where is defined after proposition 3. By using proposition 3 and the definition of ∥⋅∥_h, we obtain the following proposition.

Proposition 4: Let u be the exact solution of Eq. 1 and suppose that for any t∈I with t_K≥2, then we have:

and

where, μ_K = min (p_K+1, t_k) and C is independent of h and p.

In order to obtain an estimation of we apply the previous proposition and get the following.

Proposition 5: Let u be the exact solution of Eq. 1 and suppose that with t_K≥2. Let u^h the discrete solution of Eq. 11, then the error ζ = u_h-u satisfies:

where, μ_K = min (p_K+1, t_K) and C is independent of h and p.

RESULTS AND DISCUSSION

We shall now present some numerical results which verify the sharpness of the theoretical error bounds stated in proposition 5. To obtain a full discretization of our wave equation, we choose to augment our DG spatial discretization with the second order Newmark scheme in time¹.

In our example, the DG stabilization parameter is set to κ = 10. The functions μ and ε in Eq. 1 are supposed constants and equal to 1.

Time discretizationThe discretization of Eq. 1 in space by the DG method Eq. 9 leads to the linear second order system of ordinary differential equations:

with initial conditions:

where, M is the mass matrix and A the stiffness matrix. To discretize Eq. 11 in time, we employ the Newmark time stepping scheme¹⁴. We let k denote the time step and set t_n = n.k. Then the Newmark method consists in finding approximation to u^h(t_n) such that:

and:

For n = 1, 2,...,N-1. Here f_n= f(t_n) while, β≥0 and are free parameters that still can be chosen. We recall that for the Newmark scheme is second order accurate in time, whereas, it is only first order accurate for For β = 0, the Newmark scheme Eq. 16 and 17 requires at each time step the solution of a linear system with the matrix M. However, because individual elements decouples, M is a bloc diagonal with a bloc size equal to the number of degrees of freedom per element. It can be inverted at very low computational cost and the scheme is essentially explicit. In fact, if the bases functions are chosen mutually orthogonal, M reduces to the identity¹⁵ and the references therein. Then, with the explicit Newmark method corresponds to the standard leap-frog scheme.

For β>0, the resulting scheme is implicit and involves the solution of a linear system with the symmetric positive definite stiffness matrix A at each time step. We finally note that the second order Newmark scheme with is unconditionally stable for whereas, for the time step k has to be restricted by a CFL condition. In the case β = 0 the condition is k²λ_max (A)≤4 (1-ε), ε∈(0, 1) where, λ_max (A) is the maximal eigen value of the DG stiffness matrix A.

In our test, we will employ the implicit second order Newmark scheme, setting and in Eq. 16 and 17.

Example: We consider the three dimensional wave Eq. 1 in Ω×I= (0.1)³×(0.1) and data f, g, u₀ and u₁ chosen such that the analytical solution is given by:

This solution is arbitrarily smooth so that our theoretical assumptions are satisfied. We discretize this problem using the polynomial spaces P^p (K)³, p = 1,2 on a sequence ∏_h of tetrahedral meshes. With decreasing mesh size h smaller time step k is not necessary, because the scheme is unconditionally stable.

We show the relative errors at time T = 1 in the energy norm, as we decrease h and we remark that the decrease of the energy norm as a function of the mesh size h is of order one for p = 1 and of order two for p = 2. Then the numerical results corroborate with the expected theoretical rates of O(h^p) as we decrease the mesh size (Table 1, Fig. 1, 2).

CONCLUSION

In this study, a discontinuous Galerkin method for the discretization of the wave has been proposed and its hp-error analysis has been carried out. The hp-error estimates obtained are optimal in the mesh size and suboptimal in the approximation degree. Some numerical results are given to confirm the convergence rates as a function of the mesh size.

ACKNOWLEDGMENT

The authors gratefully acknowledge the approval and the support of this study at from the Deanship of Scientific Research study by the grant No. 8-38-1436-5. K.S.A. Northern Border University, Arar.

REFERENCES