Explicit P1 Finite Element Solution of the Maxwell-Wave Equation Coupling Problem with Absorbing b. c.

Preprint

Article

Explicit P1 Finite Element Solution of the Maxwell-Wave Equation Coupling Problem with Absorbing b. c.

Altmetrics

Downloads

Views

Comments

This version is not peer-reviewed

Submitted:

29 January 2024

Posted:

30 January 2024

You are already at the latest version

Alerts

Abstract

In this paper we address the approximation of the coupling problem for the wave equation and Maxwell's equations of electromagnetism in the time domain in terms of the electric field, by means of a nodal linear finite-element discretization in space, combined with a classical'explicit finite-difference scheme for the time discretization. Our study applies to the particular case where the dielectric permittivity has a constant value outside a sub-domain, whose closure does not intersect the boundary of the domain where the problem is defined. Inside this sub-domain Maxwell's equations hold. Outside this sub-domain the wave equation holds, which may correspond to the Maxwell's equations with a constant permittivity under certain conditions. We consider as a model the case of first order absorbing boundary conditions. Optimal error estimates that hold in natural norms under reasonable assumptions are given, among which lies a typical CFL condition for hyperbolic equations.

Keywords:

Subject: Computer Science and Mathematics - Applied Mathematics

MSC: 65M06; 65M12; 65M15; 65M55; 65M60

1. Introduction

The aim of this article is to show that an explicit

P_{1}

lumped-mass finite-element scheme can be a reliable method to solve the Maxwell’s equations of electromagnetism in the time domain and in a bounded domain in

ℜ^{N}

N = 2, 3

. This fact is illustrated here assuming a constant magnetic permeability. It is well known that in this case these equations can be expressed as a second order system in terms of the sole electric field. Moreover, the study is conducted in the particular case where the dielectric permittivity is constant in a region of the computational domain contiguous to its boundary. As a consequence, the wave equation holds therein and for this reason we address the case of a coupling problem for the Maxwell’s equations in the inner domain and the wave equation in the outer domain. However, if it happens that the Maxwell system also holds in the latter, then our study will directly apply to the Maxwell’s equations in the whole domain. Moreover, although it extends to other boundary conditions with minor modifications, as a model we consider in this work first-order absorbing boundaryboundary conditions.

Standard conforming linear finite elements are a priori an attractive tool to solve Maxwell’s equations, as an inexpensive method, especially in three-dimension space. However, this is not always a good choice for such a purpose. A primary explanation for this assertion is the fact that in geneoral the electric field cannot be searched for in the Sobolev space

{H^{1}}^{N}

, but rather in a subspace of

H (c u r l) \cap H (d i v)

consisting of vector fields satisfying certain boundary conditions. As a consequence, if the spacial domain in which the equations are defined has re-entrant corners, the subspace of

{H^{1}}^{N} \cap (H (c u r l) \cap H (d i v))

, incorporating for instance zero tangential boundary conditions, is a proper subspace of the corresponding subspace of

H (c u r l) \cap H (d i v)

owing to the so-called corner paradox (see e.g. [19,20]). This was one of the main motivations of different authors who looked into the design of finite element methods to cope with the issue. A celebrated contribution in this direction is due to Nédélec [29], namely, a family of

H (c u r l)

-conforming methods to solve these equations, commonly known as edge elements. The crucial point in the discussion on discretization methods to tackle the problem, is how to get rid of spurious solutions and other instabilities usually caused by nodal elements, such as the

P_{1}

FEM. A detailed study of these questions, together with methods especially designed to solve Maxwell’s equations is given for instance in [10,28]. A more recent approach to handle these equations using a Hermite variant of the lowest order Raviart-Thomas interpolation [30] was reported in [32].

In this work we rule out the aforementioned shortcoming, by taking advantage of the fact that the solution of the wave equation lies in

{H^{1}}^{N}

irrespective of boundary conditions. Hence in our case we can show that the

P_{1}

finite element method is indeed an accurate numerical solution tool, provided the coupled equations are written in a suitable VF (variational form). Actually, the VF employed in this work is similar but not equal to the AVF-Augmented Variational Formulation thoroughly studied by Ciarlet Jr. (cf. [16]) in the static and time-harmonic cases and in Jamelot [27] and Ciarlet Jr.-Jamelot [17]) in the time-dependent case. However, non negligible additional complexities must be dealt with, stemming from the fact that the dielectric permittivity varies in space. This is one of the main reasons that compelled the authors to carry out here a rigorous analysis of the

P_{1}

lumped-mass approximation of Maxwell’s equations. As a matter of fact, to the best of their knowledge, such results were lacking. Indeed, the case of a variable dielectric permittivity had been addressed in [18,23], but not for nodal finite elements, while in [2,15,17] nodal finite elements are dealt with, though for constant coefficients and formulations different from ours.

Conforming nodal finite elements were considered to handle the time-dependent case in the early nineties (cf. [2]) for convex domains. Later on specialists studied formulations of the static or the time-harmonic Maxwell’s equations suitable for a numerical solution with nodal elements, even for non convex domains. In this respect we refer to [3,17,21,27]. Such studies revealed the adequacy of nodal elements, at least in some relevant practical situations. Underlying this work lies precisely one such a case, characterized by the fact that the dielectric permittivity is constant in a neighborhood of the boundary of the domain of interest. On the other hand we emphasize again that there are not many studies of nodal elements as applied to the time-dependent Maxwell’s equations with variable coefficients. Hence this paper also gives a contribution in this direction. In short, by applying a lumped-mass explicit

P_{1}

finite-elements scheme to the Maxwell’s equations in terms of the electric field recast in a suitable VF, we establish that, at least in the above case, reliable numerical solutions are generated, as long as a classical stability condition is fulfilled.

Before pursuing, it should be pointed out that we had described and assessed our method in [7,8], without giving mathematical proofs of reliability. Here we present the main lines of its formal numerical analysis, though not for the same boundary conditions. For this reason the numerical examples shown in both articles are different from those presented in this work. We also observe that, akin to [8], we study here a Maxwell-wave equation coupling problem posed in two contiguous disjoint sub-domains, which happens to simplify into Maxwell’s equations in the union of both sub-domains under certain conditions. The mathematical grounds of our VF, in particular its equivalence with the strong form of the Maxwell’s equations, follows the main lines of [8].

As a matter of fact, this work somehow validates a numerical solution procedure described in [4] applied to a kind of CIP - Coefficient Inverse Problem - for the time-dependent Maxwell’s equations in the particular configuration underlined above. Let us briefly recall it below. For more details we refer to [9] and references therein.

Assume that in a vast nonconducting homogeneous medium with a given dielectric permittivity, an unknown object is searched for. The object’s material is also supposedly nonconducting, though with a different dielectric permittivity. Solenoidal electric waves of a regular pattern are emitted sufficiently far away from the search region during a certain time. In the absence of a hidden object such a pattern will be observed during all time in a certain location closer to the search zone. However, if there is indeed a hidden object, waves will be reflected and back-scattered onto the observation zone. Schematically such a process may be thought of as governed by the Maxwell system of equations of electromagnetism in an unbounded domain with a constant dielectric permittivity, except in a region surrounding the hidden object where the dielectric permittivity varies. The CIP consists of determining the spacial distribution of an unknown dielectric permittivity - i.e., a coefficient of Maxwell’s equations, - with the aim of locating the hidden object, on the basis of measurements of back-scattered waves at a small observation zone. In principle, the far electric field will still be as regular as the emitted one. However, we may conveniently solve the problem by taking a bounded computational domain consisting of two sub-domains, namely, an inner domain where the hidden object lies, and a surrounding outer domain on whose outer boundary - that is, the boundary of the computational domain, - absorbing boundary conditions are prescribed (see e.g. [22,24]). It is noticeable that, since the dielectric permittivity is constant in the outer domain, N wave equations hold therein. Moreover, since the far field is solenoidal, in practical terms it can be thought of as being also solenoidal on the boundary of the computational domain, as long as it is large enough. However, strictly speaking, such a condition cannot be prescribed, and hence it is mathematically inconsistent to consider that the full Maxwell system of equations hold in the outer domain, as it does in the inner domain. That is why we address in this article a Maxwell-wave coupling problem posed simultaneouly in the inner and the outer sub-domains, bearing in mind that, at least for the CIP in view, Maxwell’s equations are expected to hold in the union of both.

An outline of this paper is as follows: In Section 2 we describe the model problem being solved and study its equivalent VF employed in the sequel. In Section 3 we set up the discretizations of the model problem in both space and time. Section 4 is devoted to the formal reliability analysis of the explicit scheme considered in the previous section. A priori error estimates are given therein under the realistic assumption that the time step varies linearly with the mesh size as the mesh is refined. In Section 5 we present a numerical validation of our scheme. Finally we conclude in Section 6 with a few comments on the whole work.

2. The model problem

The Maxwell-wave equation coupling problem for a field

e

in a bounded Lipschitz domain

Ω

ℜ^{N}, N = 2, 3

with boundary

Γ

, that we consider in this work is posed in the following setting.

First, we consider that $Ω = {\bar{Ω}}_{i n} \cup Ω_{o u t}$ , where $Ω_{i n}$ is an interior open set whose boundary $Γ_{i n}$ does not intersect $Γ$ and $Ω_{o u t}$ is the complementary set of ${\bar{Ω}}_{i n}$ with respect to $Ω$ (the boundary of $Γ_{o u t}$ is $Γ \cup Γ_{i n}$ ). The dielectric permittivity denoted by $ε$ is assumed to belong to $C^{2} (\bar{Ω})$ and to fulfill the conditions:

ε \equiv 1 in Ω_{out} and ε \geq 1 otherwise .

(2.1)

Let

n

be the unit outer normal vector on

Γ

. We denote by

\partial_{n} (\cdot)

the outer normal derivative of a field on

Γ

. Taking into account conditions (2.1), we may prescribe wave-equation first order absorbing boundary conditions

\partial_{n} e = - \partial_{t} e

Γ \times (0, T)

[22,24], as seen hereafter.

Now, we are given

e_{0} \in {H^{1} (Ω)}^{N}

and

e_{1} \in H (d i v, Ω)

satisfying

\nabla \cdot (ε e_{0}) = \nabla \cdot (ε e_{1}) \equiv 0

, together with

f \in L^{2} [(0, T); H (d i v, Ω)]

satisfying

\nabla \cdot f \equiv 0

. Then, setting

V : = {H^{1} (Ω)}^{N}

, the problem to solve is:

\{\begin{matrix} Find e \in V : = H^{2} [(0, T); {L^{2} (Ω)}^{N}] \cap L^{2} [(0, T); V] \\ such that \\ \begin{matrix} ε \partial_{t t} e + \nabla \times \nabla \times e = f, \\ \nabla \cdot {ε e} = 0 . \end{matrix}\} & in Ω_{i n} \times (0, T), \\ \partial_{t t} e - Δ e = f & in Ω_{o u t} \times (0, T), \\ with the initial conditions \\ e (\cdot, 0) = e_{0} (\cdot), and \partial_{t} e (\cdot, 0) = e_{1} (\cdot) & in Ω \\ and the boundary conditions \\ \partial_{n} e + \partial_{t} e = 0 & on Γ \times (0, T) . \end{matrix}

(2.2)

(2.2) tells us that the wave equation holds in

Ω_{o u t} \times (0, T)

. On the other hand, the equations in braces of (2.2) make up the Maxwell’s equations in

Ω_{i n} \times (0, T)

for the sole electric field. It is noticeable that, since

f

is solenoidal by assumption, the second equation in

Ω_{i n} \times (0, T)

is superfluous i.e. reduntant. Indeed, taking the divergence of both sides of the first equation in

Ω_{i n} \times (0, T)

and denoting by u the function

\nabla \cdot (ε e)

, we have

\partial_{t t} u = 0 in Ω_{i n} \times (0, T)

. Taking into account that

u_{| t = 0} = {\partial_{t} u}_{| t = 0} = 0

by our assumption on

e_{0}

and

e_{1}

, it must hold

\nabla \cdot (ε e) = 0

Ω_{i n} \times (0, T)

. However, the same conclusion cannot be drawn for

Ω_{o u t} \times (0, T)

. This is because in this sub-domain

u : = \nabla \cdot e

solves a wave equation

u_{t t} - Δ u = 0

with zero initial conditions. Since zero boundary conditions do not necessarily hold for u, this is not sufficient to infer that

u \equiv 0

Ω_{o u t} \times (0, T)

. But, of course, nothing prevents the Maxwell’s equations from holding indeed in this domain too, as pointed out at the end of Section 1.

2.1. Notations and reminders

Before pursuing we introduce some notations and recall a couple of results to be used hereafter.

We denote the standard semi-norm of

C^{n} (\bar{Ω})

{| \cdot |}_{n, \infty}

for

n > 0

and the standard norm of

C^{0} (\bar{Ω})

{∥ \cdot ∥}_{0, \infty}

. A subset of

\bar{Ω}

will be denoted by an upper case Latin letter with or without a subscript. For any

D \subset Ω

we denote the standard inner product of

{L^{2} (D)}^{N}

{(\cdot, \cdot)}_{D}

and the corresponding norm by

{‖ {\cdot} ‖}_{D}

; if

D = Ω

we drop the subscript D.

m e a s (D)

represents the measure of D, which accounts for the length, the area or the volume of D, according to the case.

Any scalar function defined in $Ω$ will be denoted by a lower case Greek letter combined or not with other symbols different from upper case letters. For a given non negative function $ω \in L^{\infty} (Ω)$ we introduce the weighted $L^{2} (Ω)$ -semi-norm ${∥ {\cdot} ∥}_{ω} : = \sqrt{\int_{Ω} ω {| {\cdot} |}^{2}}$ , which is actually a norm if $ω \neq 0$ everywhere in $\bar{Ω}$ . The notation ${(A, B)}_{ω}$ expresses $\int_{Ω} ω A \cdot B$ , $A, B$ being two square-integrable fields in $Ω$ . If $ω$ is strictly positive this expression defines an inner product associated with the norm ${∥ {\cdot} ∥}_{ω}$ .
We denote by $< \cdot, \cdot >_{s, Γ_{i n}}$ the duality product of $H^{s} (Γ_{i n}) - H^{- s} (Γ_{i n})$ for $s \in ℜ^{+}$ .

In the sequel we shall repeatedly employ a well known operator identity applying to vector fields, namely,

- Δ \equiv - \nabla \nabla \cdot + \nabla \times \nabla \times .

(2.3)

Let D be a bounded domain of

ℜ^{N}

with boundary

\partial D

. We recall that, according to the Trace Theorem (cf. [1]) and well known results, if a given function

χ \in L^{2} (D)

has a well defined trace on

\partial D

in the space

H^{1 / 2} (\partial D)

and

Δ χ \in H^{- 1} (D)

, then necessarily

χ \in H^{1} (D)

2.2. Well-posedness considerations

The well-posedness of (2.2) will be taken for granted. However, it can be established by means of an argument similar to the one exploited in [8]. Since it is rather laborious, for the sake of brevity we skip details of such an analysis. Nevertheless we next sketch its main lines, which rely on the well-posedness of the following problem.

Let $R : = {(v, z) | v \in {H^{1} (Ω_{o u t})}^{N}, z \in H (c u r l; Ω_{i n}) \cap H (d i v; Ω_{i n}), \int_{Ω_{o u t}} u + \int_{Ω_{i n}} w = 0}$ and let $R_{ε}$ be the subspace of $R$ of the pairs $(v; z)$ satisfying $\nabla \cdot (ε z) = 0$ in $Ω_{i n}$ . Noticing that the trace over $Γ_{i n}$ of a field in $H (c u r l; Ω_{i n}) \cap H (d i v; Ω_{i n})}$ lies in ${H^{- 1 / 2} (Γ_{i n})}^{N}$ (cf. [25]), recalling the space $L_{0}^{2} (Ω) : = {v | v \in L^{2} (Ω), \int_{Ω} v = 0}$ , we set the problem

\{\begin{matrix} Given a solenoidal g \in {L_{0}^{2} (Ω)}^{N} find ((u; w); p) \in R_{ε} \times {H^{1 / 2} (Γ_{i n})}^{N} \\ such that \forall ((v; z); q) \in R_{ε} \times {H^{1 / 2} (Γ_{i n})}^{N} \\ {(\nabla u, \nabla v)}_{Ω_{o u t}} + {(\nabla \times w, \nabla \times z)}_{Ω_{i n}} + < p, v - z >_{1 / 2, Γ_{i n}} = {(g_{| Ω_{o u t}}, v)}_{Ω_{o u t}} + {(g_{| Ω_{i n}}, z)}_{Ω_{i n}}, \\ < q, u - w >_{1 / 2, Γ_{i n}} = 0 . \end{matrix}

(2.4)

Using the theory of saddle-point linear problems (cf. [11]) we have checked that (2.4) has a unique solution and moreover, by the same theory, it is equivalent to the following system:

\{\begin{matrix} Given a solenoidal g \in {L_{0}^{2} (Ω)}^{N}, find ((u; w); r; p) \in R \times L^{2} (Ω_{i n}) \times {H^{1 / 2} (Γ_{i n})}^{N} \\ such that \forall ((v; z)); s; q) \in R \times L^{2} (Ω_{i n}) \times {H^{1 / 2} (Γ_{i n})}^{N} \\ {(\nabla u, \nabla v)}_{Ω_{o u t}} + {(\nabla \times w, \nabla \times z)}_{Ω_{i n}} + {(r, \nabla \cdot (ε z))}_{Ω_{i n}} + < p, v - z >_{1 / 2, Γ_{i n}} \\ = {(g_{| Ω_{o u t}}, v)}_{Ω_{o u t}} + {(g_{| Ω_{i n}}, z)}_{Ω_{i n}}, \\ {(\nabla \cdot (ε w), s)}_{Ω_{i n}} = 0, \\ < q, u - w >_{1 / 2, Γ_{i n}} = 0 . . \end{matrix}

(2.5)

It is clear that the solution of (2.5) solves the following problem:

\{\begin{matrix} Given g \in {L_{0}^{2} (Ω)}^{N} fulfilling \nabla \cdot g = 0 in Ω, & find (u; w) \in R such that \\ - Δ u = g_{| Ω_{o u t}} & in Ω_{o u t} \\ \begin{matrix} \nabla \times \nabla \times w = g_{| Ω_{i n}}, \\ \nabla \cdot (ε w) = 0 . \end{matrix}\} & in Ω_{i n}, \\ u = w & on Γ_{i n}, \\ \partial_{n} u = 0 & on Γ, \\ \int_{Ω_{o u t}} u + \int_{Ω_{i n}} w = 0 . \end{matrix}

(2.6)

Thus, from classical results (cf. [13]), any linear second order hyperbolic counterpart of (2.6) assorted with proper initial conditions also has a unique solution. Let us consider a field

e

defined in

Ω

such that

e_{| Ω_{o u t}} = u

and

e_{| Ω_{i n}} = w

. Since

\nabla \cdot w = - \nabla log (ε) \cdot w \in L^{2} (Ω_{i n})

from the assumed regularity of

ε

, we have

\nabla \nabla \cdot w \in {H^{- 1} (Ω_{i n})}^{N}

. This implies that

Δ w \in {H^{- 1} (Ω_{i n})}^{N}

by the identity (2.3), since

\nabla \times \nabla \times w \in {L^{2} {(Ω_{i n}}}^{N}

. Therefore

w \in {H^{1} (Ω_{i n})}^{N}

owing to the coincidence of

u

and

w

Γ_{i n}

. Hence

e

lies in

{H^{1} (Ω)}^{N}

and moreover it solves a well-posed elliptic problem in

Ω

. Thus the well-posedness of (2.2) follows from the fact that it is a second order hyperbolic counterpart of (2.6).

Remark 2.1.

As already pointed out in Section 1, the study that follows also applies to several types of boundary conditions, for which such an

H^{1}

-regularity is known to hold. As pointeed out in Section 1, the choice of absorbing boundary conditions here was motivated by the fact that they correspond to practical situations addressed in [9] and references therein. ■

2.3. Variational form

Throughout this article we work with the variational problem (2.7) stated below, supposedly equivalent to (2.2).

Requiring that

{\tilde{e}}_{| t = 0} = e_{0}

and

{\partial_{t} \tilde{e}}_{| t = 0} = e_{1}

, we wish to,

\{\begin{matrix} Find e \in \tilde{V}, where \\ \tilde{V} : = {v | v \in H^{2} [(0, T); {L^{2} (Ω)}^{N}] \cap H^{1} [(0, T); V]} \\ such that \forall v \in {H^{1} (Ω)}^{N} and \forall t \in (0, T) \\ {(\partial_{t t} \tilde{e}, v)}_{ε} + (\nabla \tilde{e}, \nabla v) + (\nabla \cdot ε \tilde{e}, \nabla \cdot v) - (\nabla \cdot \tilde{e}, \nabla \cdot v) + {(\partial_{t} \tilde{e}, v)}_{Γ} = (f, v) . \end{matrix}

(2.7)

Under the conditions assumed in this work, problem (2.7) is equivalent to the coupling problem (2.2). Indeed, we have,

Proposition 2.1.

If the solution

e

to (2.2) belongs to

\tilde{V}

, the following assertions hold:

1.: $e$ is also a solution to (2.7).
2.: Any solution to (2.7) is unique, and thus it is the solution to equation (2.2).

Proof:

Using the operator identity (2.3) we rewrite the second equation of (2.2) as,

ε \partial_{t t} e - Δ e + \nabla \nabla \cdot e = f in Ω_{i n} \times (0, T) .

(2.8)

We know that the solution of (2.2) satisfies

\nabla \nabla \cdot (ε e) = 0

Ω_{i n} \times (0, T)

. If we subtract this equation from (2.8), we obtain:

ε \partial_{t t} e - Δ e - \nabla \nabla \cdot (ε e) + \nabla \nabla \cdot e = f in Ω_{i n} \times (0, T) .

(2.9)

Since

ε \equiv 1

Ω_{o u t}

it is readily seen that (2.9) also holds in the whole

Ω \times (0, T)

Thus, taking an arbitrary $v \in V$ and using Green’s first identity, together with the absorbing boundary conditions satisfied by $e$ , we readily obtain $\forall v \in V$ ,

(ε \partial_{t t} e, v) + (\nabla e, \nabla v) - (\nabla \cdot e, \nabla \cdot v) + (\nabla \cdot {ε e}, \nabla \cdot v) + {(\partial_{t} e, v)}_{Γ} = (f, v) \forall t \in (0, T),

(2.10)

which establishes that

e

solves (2.7).

2.: In order to prove the uniqueness of a solution to (2.7), we resort to the following energy estimate demonstrated in [5] for variational problem (2.7):

\begin{matrix} \forall t \in (0, T] \{‖ \partial_{t} \tilde{e} ‖_{ε}^{2} + ‖ \nabla \tilde{e} ‖^{2} + {‖ \nabla \cdot \tilde{e} ‖}_{ε - 1}^{2}\} (t) \\ \leq E \{\int_{0}^{T} {‖ f (\cdot, t) ‖}^{2} d t + ‖ e_{0} ‖^{2} + ‖ \nabla e_{0} ‖^{2} + ‖ \nabla \cdot e_{0} ‖_{ε - 1}^{2} + {‖ e_{1} ‖}_{ε}^{2}\}, \end{matrix}

(2.11)

where

E

is a constant independent of

e_{0}

and

e_{1}

The uniqueness follows from (2.11) owing to the linearity of problem (2.7). ■

3. Space-time discretization

We next describe our numerical scheme to solve (2.7). Henceforth, for the sake of simplicity we assume that both

Ω

and

Ω_{i n}

are polytopes and, without loss of generality, we take

f \equiv 0

3.1. Space semi-discretization

Let

T_{h}

be a mesh fitting

Ω

, consisting of N-simplices with maximum edge length h, belonging to a quasi-uniform family of meshes (cf. [14]). Each element

K \in T_{h}

is to be understood as a closed set. Practical calculations are certainly simplified in case

Ω_{i n}

is the union of N-simplices belonging to

T_{h}

, which we also assume, even though such an assumption is not essential. We denote by

V_{h}

the

P_{1}

FE-space of continuous functions related to

T_{h}

Setting $V_{h} : = {V_{h}}^{N}$ we define $e_{0 h}$ (resp. $e_{1 h}$ ) to be the standard $V_{h}$ -interpolate of $e_{0}$ (resp. $e_{1}$ ). Then the space semi-discretized problem to solve consists of finding $e_{h} \in V_{h}$ such that, $\forall t \in (0, T]$

\begin{matrix} {(\partial_{t t} e_{h}, v)}_{ε} + (\nabla e_{h}, \nabla v) + (\nabla \cdot {ε e_{h}}, \nabla \cdot v) - (\nabla \cdot e_{h}, \nabla \cdot v) + {(\partial_{t} e_{h}, v)}_{Γ} = 0 \\ \forall v \in V_{h}, \\ with \\ e_{h} (\cdot, 0) = e_{0 h} (\cdot) and \partial_{t} e_{h} (\cdot, 0) = e_{1 h} (\cdot) in Ω . \end{matrix}

(3.12)

3.2. Full discretization

To begin with we consider a natural centered time-discretization scheme to solve (3.12), namely: Given a number M of time steps we define the time increment

τ : = T / M

. Then we approximate

e_{h} (k τ)

e_{h}^{k}

for

k = 1, 2, \dots, M

according to the following FE scheme for

k = 1, 2, \dots, M - 1

\begin{matrix} {(\frac{e_{h}^{k + 1} - 2 e_{h}^{k} + e_{h}^{k - 1}}{τ^{2}}, v)}_{ε} + (\nabla e_{h}^{k}, \nabla v) + (\nabla \cdot ε e_{h}^{k}, \nabla \cdot v) - (\nabla \cdot e_{h}^{k}, \nabla \cdot v) \\ + {(\frac{e_{h}^{k + 1} - e_{h}^{k - 1}}{2 τ}, v)}_{Γ} = 0 \forall v \in V_{h}, \\ with \\ e_{h}^{0} = e_{0 h} and e_{h}^{1} = e_{h}^{0} + τ e_{1 h} in Ω . \end{matrix}

(3.13)

Owing to its coupling with

e_{h}^{k}

and

e_{h}^{k - 1}

on the left hand side of (3.13),

e_{h}^{k + 1}

cannot be determined explicitly by this scheme at every time step. In order to enable an explicit solution we resort to the classical mass-lumping technique. We recall that this consists of replacing on the left hand side the inner product

(u, v)

(resp.

{(u, v)}_{Γ}

) by an inner product

{(u, v)}_{h}

using the trapezoidal rule to approximate the integral of

\int_{K} u_{| K} \cdot v_{| K} d x

, for every element K in

T_{h}

, for

(u; v) \in {C^{0} (\bar{Ω})}^{N} \times {C^{0} (\bar{Ω})}^{N}

. In the case of (3.13),

u

stands for

ε (e_{h}^{k + 1} - 2 e_{h}^{k} + e_{h}^{k - 1})

and

v

is an arbitrary field in

V_{h}

. It is well-known that in this case the matrix associated with

{(ε e_{h}^{k + 1}, v)}_{h}

is a diagonal matrix. Similarly, a mass-lumping approximation

{(e_{h}^{k + 1} - e_{h}^{k - 1}, v)}_{Γ, h}

must be used for the inner product

{(e_{h}^{k + 1} - e_{h}^{k - 1}, v)}_{Γ}

The expression ${(ε u, v)}_{h}$ for continuous $u$ and $v$ gives rise to an approximation of the inner product ${(u, v)}_{ε}$ henceforth denoted by ${(u, v)}_{ε, h}$ . In order to simplify the calculations we further approximate the inner product ${(u, v)}_{ε, h}$ by the inner product ${(u, v)}_{ε_{h}, h}$ , whose definition is given below, followed by the expression ${(u, v)}_{Γ, h}$ approximating the inner product ${(u, v)}_{Γ}$ for every $(u; v) \in {C^{0} (\bar{Ω})}^{N} \times {C^{0} (\bar{Ω})}^{N}$ .

{(u, v)}_{ε_{h}, h} : = \sum_{K \in T_{h}} ε (G_{K}) m e a s (K) \sum_{i = 1}^{N + 1} \frac{u (S_{K, i}) \cdot v (S_{K, i})}{N + 1},

where

S_{K, i}

are the vertexes of K,

i = 1, \dots, N + 1

and

G_{K}

is the centroid of K.

Generically denoting by

F^{K}

the edges for

N = 2

or the faces for

N = 3

of an N -simplex K, let

S_{h}

be the subset of

T_{h}

consisting of K such that

Γ \cap K \neq \emptyset

. Then we set

{(u, v)}_{Γ, h} : = \sum_{K \in S_{h}} \sum_{F^{K} \cap Γ \neq \emptyset} m e a s (F^{K}) \sum_{i = 1}^{N} \frac{u (R_{F^{K}, i}) \cdot v (R_{F^{K}, i})}{N},

where

R_{F^{K}, i}

are the vertexes of

F^{K}

i = 1, \dots, N

For coherence

ε_{h}

is defined to be the function whose value in each

K \in T_{h}

is constant equal to

ε (G_{K})

. Furthermore we introduce the norms

{‖ {\cdot} ‖}_{ε_{h}, h}

and

{‖ {\cdot} ‖}_{h}

V_{h}

, given by

{({\cdot}, {\cdot})}_{ε_{h}, h}^{1 / 2}

and

{({\cdot}, {\cdot})}_{h}^{1 / 2}

, respectively. Similarly we denote by

{‖ {\cdot} ‖}_{Γ, h}

the norm defined by

{({\cdot}, {\cdot})}_{Γ, h}^{1, 2}

. Then still denoting the approximation of

e_{h} (k τ)

e_{h}^{k}

, for

k = 1, 2, . . ., M

we determine

e_{h}^{k + 1}

by,

\begin{matrix} {(\frac{e_{h}^{k + 1} - 2 e_{h}^{k} + e_{h}^{k - 1}}{τ^{2}}, v)}_{ε_{h}, h} + (\nabla e_{h}^{k}, \nabla v) + (\nabla \cdot ε e_{h}^{k}, \nabla \cdot v) - (\nabla \cdot e_{h}^{k}, \nabla \cdot v) \\ + {(\frac{e_{h}^{k + 1} - e_{h}^{k - 1}}{2 τ}, v)}_{Γ, h} = 0 \forall v \in V_{h}, \\ with \\ e_{h}^{0} = e_{0 h} and e_{h}^{1} = e_{h}^{0} + τ {e_{1}}_{h} in Ω . \end{matrix}

(3.14)

Now we adapt a result given in Lemma 3 of [12], which allows us to assert that the following upper bounds hold:

{‖ v ‖}_{ε_{h}} \leq {‖ v ‖}_{ε_{h}, h} \forall v \in V_{h} .

(3.15)

Similarly we have,

{‖ v ‖}_{Γ} \leq {‖ v ‖}_{Γ, h} \forall v \in V_{h} .

(3.16)

4. Reliability analysis

We next show that, under very reasonable conditions, optimal-order error estimates in a natural sense hold for approximations of the solution of problem (2.7) generated by the scheme (3.14).

4.1. Scheme stability

In order to conveniently prepare the subsequent steps of the reliability study of (3.14), following a technique thoroughly exploited in [31], we carry out the stability analysis of a more general form thereof, encompassing scheme (3.14) as a particular case, namely,

\begin{matrix} {(\frac{e_{h}^{k + 1} - 2 e_{h}^{k} + e_{h}^{k - 1}}{τ^{2}}, v)}_{ε_{h}, h} + (\nabla e_{h}^{k}, \nabla v) + (\nabla \cdot ε e_{h}^{k}, \nabla \cdot v) - (\nabla \cdot e_{h}^{k}, \nabla \cdot v) \\ + {(\frac{e_{h}^{k + 1} - e_{h}^{k - 1}}{2 τ}, v)}_{Γ, h} = F^{k} (v) + (d^{k}, \nabla \cdot v) + G^{k} (v) \forall v \in V_{h}, \\ e_{h}^{0} = e_{0 h} and e_{h}^{1} = e_{h}^{0} + τ e_{1 h} in Ω . \end{matrix}

(4.17)

where for every

k \in {1, 2, \dots, N - 1}

F^{k}

and

G^{k}

are given bounded linear functionals over

V_{h}

and the space of traces over

Γ

of fields in

V_{h}

equipped with the norms

{∥ \cdot ∥}_{h}

and

{∥ \cdot ∥}_{Γ, h}

respectively. We denote by

| F^{k} |_{h}

and

| G^{k} |_{Γ, h}

the underlying norms of both functionals.

d^{k}

in turn is a given function in

L^{2} (Ω)

for

k \in {1, 2, \dots, M - 1} .

Taking

v = e_{h}^{k + 1} - e_{h}^{k - 1}

in (4.17) we obtain for

k = 1, 2, \dots, M - 1

\begin{matrix} {(\frac{e_{h}^{k + 1} - 2 e_{h}^{k} + e_{h}^{k - 1}}{τ}, \frac{e_{h}^{k + 1} - e_{h}^{k - 1}}{τ})}_{ε_{h}, h} + (\nabla e_{h}^{k}, \nabla e_{h}^{k + 1} - \nabla e_{h}^{k - 1}) \\ + (\nabla \cdot {ε - 1} e_{h}^{k}, \nabla \cdot e_{h}^{k + 1} - \nabla \cdot e_{h}^{k - 1}) + {(\frac{e_{h}^{k + 1} - e_{h}^{k - 1}}{2 τ}, e_{h}^{k + 1} - e_{h}^{k - 1})}_{Γ, h} \\ = F^{k} (e_{h}^{k + 1} - e_{h}^{k - 1}) + (d^{k}, \nabla \cdot {e_{h}^{k + 1} - e_{h}^{k - 1}}) + G^{k} (e_{h}^{k + 1} - e_{h}^{k - 1}) \end{matrix}

(4.18)

Noting that

e_{h}^{k + 1} - 2 e_{h}^{k} + e_{h}^{k - 1} = (e_{h}^{k + 1} - e_{h}^{k}) - (e_{h}^{k} - e_{h}^{k - 1})

and that

e_{h}^{k + 1} - e_{h}^{k - 1} = (e_{h}^{k + 1} - e_{h}^{k}) + (e_{h}^{k} - e_{h}^{k - 1})

, the following estimate trivially holds for equation (4.17):

\begin{matrix} {∥\frac{e_{h}^{k + 1} - e_{h}^{k}}{τ}∥}_{ε_{h}, h}^{2} - {∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥}_{ε_{h}, h}^{2} + I_{1} + I_{2} + 2 τ {∥\frac{e_{h}^{k + 1} - e_{h}^{k - 1}}{2 τ}∥}_{Γ, h}^{2} \\ \leq | F^{k} |_{h} {∥e_{h}^{k + 1} - e_{h}^{k - 1}∥}_{h} + ∥ d^{k} ∥ ∥ \nabla \cdot (e_{h}^{k + 1} - e_{h}^{k - 1}) ∥ \\ + | G^{k} |_{Γ, h} {∥e_{h}^{k + 1} - e_{h}^{k - 1}∥}_{Γ, h} \\ where \\ I_{1} : = (\nabla e_{h}^{k}, \nabla {e_{h}^{k + 1} - e_{h}^{k - 1}}); \\ and \\ I_{2} : = (\nabla \cdot {ε - 1} e_{h}^{k}, \nabla \cdot {e_{h}^{k + 1} - e_{h}^{k - 1}}) . \end{matrix}

(4.19)

Next we estimate the terms

I_{1}

and

I_{2}

given by (4.19).

First of all it is easy to see that

\begin{matrix} I_{1} = \frac{1}{2} (∥ \nabla e_{h}^{k + 1} ∥^{2} + ∥ \nabla e_{h}^{k} ∥^{2} - ∥ \nabla (e_{h}^{k + 1} - e_{h}^{k}) ∥^{2}) \\ - \frac{1}{2} (∥ \nabla e_{h}^{k} ∥^{2} + ∥ \nabla e_{h}^{k - 1} ∥^{2} - ∥ \nabla (e_{h}^{k} - e_{h}^{k - 1}) ∥^{2}) . \end{matrix}

(4.20)

Next we note that,

\begin{matrix} I_{2} = J_{1} + J_{2} \\ where \\ J_{1} : = {(\nabla \cdot e_{h}^{k}, \nabla \cdot {e_{h}^{k + 1} - e_{h}^{k - 1}})}_{ε - 1} \\ and \\ J_{2} : = (\nabla ε \cdot e_{h}^{k}, \nabla \cdot {e_{h}^{k + 1} - e_{h}^{k - 1}}) . \end{matrix}

(4.21)

Similarly to (4.20) we can write,

\begin{matrix} J_{1} = \frac{1}{2} (∥ \nabla \cdot e_{h}^{k + 1} ∥_{ε - 1}^{2} + ∥ \nabla \cdot e_{h}^{k} ∥_{ε - 1}^{2} - ∥ \nabla \cdot (e_{h}^{k + 1} - e_{h}^{k}) ∥_{ε - 1}^{2}) \\ - \frac{1}{2} (∥ \nabla \cdot e_{h}^{k} ∥_{ε - 1}^{2} + ∥ \nabla \cdot e_{h}^{k - 1} ∥_{ε - 1}^{2} - ∥ \nabla \cdot (e_{h}^{k} - e_{h}^{k - 1}) ∥_{ε - 1}^{2}) \end{matrix}

(4.22)

Now observing that

\nabla ε \equiv 0

Γ

, we integrate by parts

J_{2}

given by (4.21), to get

\begin{matrix} J_{2} = - (\nabla {\nabla ε \cdot e_{h}^{k}}, e_{h}^{k + 1} - e_{h}^{k - 1}) . \end{matrix}

(4.23)

Let us rewrite

J_{2}

as,

\begin{matrix} J_{2} = M_{1} + M_{2} \\ where \\ M_{1} : = - (\nabla \nabla ε e_{h}^{k}, e_{h}^{k + 1} - e_{h}^{k - 1}) \\ and \\ M_{2} : = - ({\nabla e_{h}^{k}}^{T} \nabla ε, e_{h}^{k + 1} - e_{h}^{k - 1}) \end{matrix}

(4.24)

M_{1}

in turn can be rewritten as follows:

\begin{matrix} M_{1} = N_{1} + N_{2} \\ where \\ N_{1} : = - τ (\nabla \nabla ε e_{h}^{k}, \frac{e_{h}^{k + 1} - e_{h}^{k}}{τ}) \\ and \\ N_{2} : = - τ (\nabla \nabla ε e_{h}^{k}, \frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}) . \end{matrix}

(4.25)

Then we further observe that

N_{1} = - τ^{2} \sum_{i = 1}^{k} (\nabla \nabla ε \frac{e_{h}^{i} - e_{h}^{i - 1}}{τ}, \frac{e_{h}^{k + 1} - e_{h}^{k}}{τ}) - τ (\nabla \nabla ε e_{h}^{0}, \frac{e_{h}^{k + 1} - e_{h}^{k}}{τ}) .

(4.26)

and hence,

N_{1} \geq - τ^{2} {| ε |}_{2, \infty} ∥\frac{e_{h}^{k + 1} - e_{h}^{k}}{τ}∥ \sum_{i = 1}^{k} ∥\frac{e_{h}^{i} - e_{h}^{i - 1}}{τ}∥ + {τ | ε |}_{2, \infty} ∥ e_{h}^{0} ∥ ∥\frac{e_{h}^{k + 1} - e_{h}^{k}}{τ}∥

or yet,

N_{1} \geq - τ {| ε |}_{2, \infty} ∥\frac{e_{h}^{k + 1} - e_{h}^{k}}{τ}∥ \{τ \sqrt{k} {(\sum_{i = 1}^{k} {∥\frac{e_{h}^{i} - e_{h}^{i - 1}}{τ}∥}^{2})}^{1 / 2} + ∥ e_{h}^{0} ∥\},

and noting that

k \leq T / τ

we get

\begin{matrix} N_{1} \geq - τ {| ε |}_{2, \infty} ∥\frac{e_{h}^{k + 1} - e_{h}^{k}}{τ}∥ \{τ \sqrt{\frac{T}{τ}} {(\sum_{i = 1}^{k} {∥\frac{e_{h}^{i} - e_{h}^{i - 1}}{τ}∥}^{2})}^{1 / 2} + ∥ e_{h}^{0} ∥\} . \end{matrix}

(4.27)

Applying to (4.27) Young’s inequality

a b \leq δ a^{2} / 2 + b^{2} / (2 δ)

\forall a, b \in ℜ

and

δ > 0

with

δ = 1

, we easily conclude that

\begin{matrix} N_{1} \geq - \frac{τ}{2} {| ε |}_{2, \infty} (2 {∥\frac{e_{h}^{k + 1} - e_{h}^{k}}{τ}∥}^{2} + τ T \sum_{i = 1}^{k} {∥\frac{e_{h}^{i} - e_{h}^{i - 1}}{τ}∥}^{2} + {∥ e_{h}^{0} ∥}^{2}) . \end{matrix}

(4.28)

Similarly to (4.28),

\begin{matrix} N_{2} \geq - \frac{τ}{2} {| ε |}_{2, \infty} (2 {∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥}^{2} + τ T \sum_{i = 1}^{k} {∥\frac{e_{h}^{i} - e_{h}^{i - 1}}{τ}∥}^{2} + {∥ e_{h}^{0} ∥}^{2}) . \end{matrix}

(4.29)

Combining (4.28) and (4.29) we come up with

\begin{matrix} M_{1} \geq - τ {| ε |}_{2, \infty} ({∥\frac{e_{h}^{k + 1} - e_{h}^{k}}{τ}∥}^{2} + {∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥}^{2} + L_{k}), \\ where \\ L_{k} : = τ T \sum_{i = 1}^{k} {∥\frac{e_{h}^{i} - e_{h}^{i - 1}}{τ}∥}^{2} + {∥ e_{h}^{0} ∥}^{2} . \end{matrix}

(4.30)

As for

M_{2}

given by (4.24) we have:

M_{2} \geq - {|ε|}_{1, \infty} ∥\nabla e_{h}^{k}∥ ∥e_{h}^{k + 1} - e_{h}^{k - 1}∥ \geq - τ {| ε |}_{1, \infty} ∥\nabla e_{h}^{k}∥ (∥\frac{e_{h}^{k + 1} - e_{h}^{k}}{τ}∥ + ∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥)

or yet

M_{2} \geq - \frac{τ}{2} {| ε |}_{1, \infty} (2 ∥ \nabla e_{h}^{k} ∥^{2} + {∥\frac{e_{h}^{k + 1} - e_{h}^{k}}{τ}∥}^{2} + {∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥}^{2}) .

(4.31)

Now we recall (4.19) together with (3.16) and notice that for every square-integrable field

A

Ω

we have

∥ A ∥ \leq {∥ A ∥}_{ε_{h}}

, Then taking into account that

∥ \nabla \cdot v ∥ \leq \sqrt{N} ∥ \nabla v ∥

\forall v \in {H^{1} (Ω)}^{N}

, and using Young’s inequality with

δ = τ

δ = 1 / τ

and

δ = τ

, respectively, we easily obtain the following estimates:

\begin{matrix} | F^{k} |_{h} ∥ e_{h}^{k + 1} - e_{h}^{k - 1} ∥ \leq \frac{τ}{2} {| F^{k} |}_{h}^{2} + τ ({∥\frac{e_{h}^{k + 1} - e_{h}^{k}}{τ}∥}_{h}^{2} + {∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥}_{h}^{2}), \\ ∥ d^{k} ∥ ∥ \nabla \cdot (e_{h}^{k + 1} - e_{h}^{k - 1}) ∥ \leq \frac{N}{2 τ} ∥ d^{k} ∥^{2} + τ (∥ \nabla e_{h}^{k + 1} ∥^{2} + ∥ \nabla e_{h}^{k - 1} {) ∥}^{2}), \\ | G^{k} |_{Γ, h} {∥e_{h}^{k + 1} - e_{h}^{k - 1}∥}_{Γ, h} \leq \frac{τ}{2} {| G^{k} |}_{Γ, h}^{2} + 2 τ {∥\frac{e_{h}^{k + 1} - e_{h}^{k - 1}}{2 τ}∥}_{Γ, h}^{2} . \end{matrix}

(4.32)

where in the first and the second inequality we also used the fact that

{∥ A \pm B ∥}^{2} \leq {2 (∥ A ∥}^{2} + {∥ B ∥}^{2})

for all square-integrable fields

A

and

B

Now we collect (4.20), (4.21), (4.22), (4.23), (4.24) and (4.30), (4.31), (4.32) to plug them into (4.19). Using the fact that $∥ \cdot ∥ \leq {∥ \cdot ∥}_{h} \leq {∥ \cdot ∥}_{ε_{h}} \leq {∥ \cdot ∥}_{ε_{h}, h}$ we obtain for $1 \leq k \leq M - 1$ :

\begin{matrix} (A_{k} - B_{k}) + (C_{k} - D_{k}) \leq \frac{τ}{2} | F^{k} |_{h}^{2} + \frac{N}{2 τ} {∥ d^{k} ∥}^{2} + \frac{τ}{2} {|G^{k}|}_{Γ_{h}}^{2}, \\ where \\ A_{k} : = {∥\frac{e_{h}^{k + 1} - e_{h}^{k}}{τ}∥}_{ε_{h}, h}^{2} - \frac{τ}{2} (2 + {| ε |}_{1, \infty} + 2 {| ε |}_{2, \infty}) {∥\frac{e_{h}^{k + 1} - e_{h}^{k}}{τ}∥}_{ε_{h}, h}^{2} \\ + \frac{1}{2} \{(1 - 2 τ) ∥ \nabla e_{h}^{k + 1} ∥^{2} {+ (1 - τ | ε |}_{1, \infty}) ∥ \nabla e_{h}^{k} ∥^{2}\} \\ B_{k} : = {∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥}_{ε_{h}, h}^{2} + \frac{τ}{2} (2 + {| ε |}_{1, \infty} + 2 {| ε |}_{2, \infty}) {∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥}_{ε_{h}, h}^{2} \\ + \frac{1}{2} \{(1 + 2 τ) ∥ \nabla e_{h}^{k - 1} ∥^{2} + (1 + {τ | ε |}_{1, \infty}) {∥ \nabla e_{h}^{k} ∥}^{2}\} + τ {| ε |}_{2, \infty} L_{k} \\ C_{k} : = - \frac{1}{2} \{∥ \nabla (e_{h}^{k + 1} - e_{h}^{k}) ∥^{2} + {∥ \nabla \cdot (e_{h}^{k + 1} - e_{h}^{k}) ∥}_{ε - 1}^{2}\} + \frac{1}{2} (∥ \nabla \cdot e_{h}^{k + 1} ∥_{ε - 1}^{2} + {∥ \nabla \cdot e_{h}^{k} ∥}_{ε - 1}^{2}) \\ D_{k} : = - \frac{1}{2} \{∥ \nabla (e_{h}^{k} - e_{h}^{k - 1}) ∥^{2} + {∥ \nabla \cdot (e_{h}^{k} - e_{h}^{k - 1}) ∥}_{ε - 1}^{2}\} + \frac{1}{2} (∥ \nabla \cdot e_{h}^{k} ∥_{ε - 1}^{2} + {∥ \nabla \cdot e_{h}^{k - 1} ∥}_{ε - 1}^{2}) . \end{matrix}

(4.33)

Setting

\begin{matrix} η : = 2 + {| ε |}_{1, \infty} + 2 {| ε |}_{2, \infty}; \\ θ : = {| ε |}_{1, \infty}; \\ ρ : = T^{2} {| ε |}_{2, \infty}, \end{matrix}

(4.34)

we can rewrite

A_{k}

and

B_{k}

as follows:

\begin{matrix} A_{k} : = {∥\frac{e_{h}^{k + 1} - e_{h}^{k}}{τ}∥}_{ε_{h}, h}^{2} - \frac{τ}{2} η {∥\frac{e_{h}^{k + 1} - e_{h}^{k}}{τ}∥}_{ε_{h}, h}^{2} \\ + \frac{1}{2} (∥ \nabla e_{h}^{k + 1} ∥^{2} + {∥ \nabla e_{h}^{k} ∥}^{2}) - τ (∥ \nabla e_{h}^{k + 1} ∥^{2} + \frac{θ}{2} {∥ \nabla e_{h}^{k} ∥}^{2}); \\ B_{k} : = {∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥}_{ε_{h}, h}^{2} + \frac{τ}{2} η {∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥}_{ε_{h}, h}^{2} \\ + \frac{1}{2} (∥ \nabla e_{h}^{k - 1} ∥^{2} + {∥ \nabla e_{h}^{k} ∥}^{2}) + τ (∥ \nabla e_{h}^{k - 1} ∥^{2} + \frac{θ}{2} {∥ \nabla e_{h}^{k} ∥}^{2}) + \frac{τ ρ}{T^{2}} L_{k} \end{matrix}

(4.35)

Then we note that for

1 \leq k \leq m - 1

with

2 \leq m \leq M - 1

\begin{matrix} A_{k - 1} - B_{k} = - τ η {∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥}_{ε_{h}, h}^{2} - τ (1 + \frac{θ}{2}) (∥ \nabla e_{h}^{k} ∥^{2} + {∥ \nabla e_{h}^{k - 1} ∥}^{2}) - \frac{τ ρ}{T^{2}} L_{k} . \end{matrix}

(4.36)

It follows that

\begin{matrix} \sum_{k = 1}^{m} (A_{k} - B_{k}) = A_{m} - B_{1} + \sum_{k = 2}^{m} (A_{k - 1} - B_{k}) \\ = \frac{1}{2} (2 - τ η) {∥\frac{e_{h}^{m + 1} - e_{h}^{m}}{τ}∥}_{ε_{h}, h}^{2} + \frac{1}{2} (1 - 2 τ) ∥ \nabla e_{h}^{m + 1} ∥^{2} + \frac{1}{2} (1 - τ θ) {∥ \nabla e_{h}^{m} ∥}^{2} \\ - \frac{1}{2} (2 + τ η) {∥\frac{e_{h}^{1} - e_{h}^{0}}{τ}∥}_{ε_{h}, h}^{2} - \frac{1}{2} (1 + 2 τ) ∥ \nabla e_{h}^{0} ∥^{2} - \frac{1}{2} (1 + τ θ) {∥ \nabla e_{h}^{1} ∥}^{2} - \frac{τ ρ}{T^{2}} L_{1} \\ - \sum_{k = 2}^{m} \{τ η {∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥}_{ε_{h}, h}^{2} + \frac{τ}{2} (2 + θ) (∥ \nabla e_{h}^{k} ∥^{2} + {∥ \nabla e_{h}^{k - 1} ∥}^{2}) + \frac{τ ρ}{T^{2}} L_{k}\} . \end{matrix}

(4.37)

Now we extend to

k = 1

the summation range on the right hand side of (4.37) to obtain,

\begin{matrix} \sum_{k = 1}^{m} (A_{k} - B_{k}) \\ = \frac{1}{2} (2 - τ η) {∥\frac{e_{h}^{m + 1} - e_{h}^{m}}{τ}∥}_{ε_{h}, h}^{2} + \frac{1}{2} (1 - 2 τ) ∥ \nabla e_{h}^{m + 1} ∥^{2} + \frac{1}{2} (1 - τ θ) {∥ \nabla e_{h}^{m} ∥}^{2} \\ - \frac{1}{2} (2 - τ η) {∥\frac{e_{h}^{1} - e_{h}^{0}}{τ}∥}_{ε_{h}, h}^{2} - \frac{1}{2} (1 - τ θ) ∥ \nabla e_{h}^{0} ∥^{2} - \frac{1}{2} (1 - 2 τ) {∥ \nabla e_{h}^{1} ∥}^{2} \\ - \sum_{k = 1}^{m} \{τ η {∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥}_{ε_{h}, h}^{2} + \frac{τ}{2} (2 + θ) (∥ \nabla e_{h}^{k} ∥^{2} + {∥ \nabla e_{h}^{k - 1} ∥}^{2}) + \frac{τ ρ}{T^{2}} L_{k}\} . \end{matrix}

(4.38)

Moreover since

C_{k - 1} = D_{k}

for all

m \geq k \geq 2

, recalling (4.33) we easily come up with,

\begin{matrix} \sum_{k = 1}^{m} (C_{k} - D_{k}) = C_{m} - D_{1} \\ = - \frac{1}{2} \{∥ \nabla (e_{h}^{m + 1} - e_{h}^{m}) ∥^{2} + {∥ \nabla \cdot (e_{h}^{m + 1} - e_{h}^{m}) ∥}_{ε - 1}^{2}\} + \frac{1}{2} (∥ \nabla \cdot e_{h}^{m + 1} ∥_{ε - 1}^{2} + {∥ \nabla \cdot e_{h}^{m} ∥}_{ε - 1}^{2}) \\ + \frac{1}{2} \{∥ \nabla (e_{h}^{1} - e_{h}^{0}) ∥^{2} + {∥ \nabla \cdot (e_{h}^{1} - e_{h}^{0}) ∥}_{ε - 1}^{2}\} - \frac{1}{2} (∥ \nabla \cdot e_{h}^{1} ∥_{ε - 1}^{2} + {∥ \nabla \cdot e_{h}^{0} ∥}_{ε - 1}^{2}) . \end{matrix}

(4.39)

Combining (4.38) and (4.39), we obtain for

2 \leq m \leq - 1

\begin{matrix} \sum_{k = 1}^{m} (A_{k} + C_{k}) - \sum_{k = 1}^{m} (B_{k} + D_{k}) \\ = \frac{1}{2} (2 - τ η) {∥\frac{e_{h}^{m + 1} - e_{h}^{m}}{τ}∥}_{ε_{h}, h}^{2} + \frac{1}{2} (1 - 2 τ) ∥ \nabla e_{h}^{m + 1} ∥^{2} + \frac{1}{2} (1 - τ θ) {∥ \nabla e_{h}^{m} ∥}^{2} \\ - \frac{1}{2} (2 - τ η) {∥\frac{e_{h}^{1} - e_{h}^{0}}{τ}∥}_{ε_{h}, h}^{2} - \frac{1}{2} (1 - τ θ) ∥ \nabla e_{h}^{0} ∥^{2} - \frac{1}{2} (1 - 2 τ) {∥ \nabla e_{h}^{1} ∥}^{2} \\ - \sum_{k = 1}^{m} \{τ η {∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥}_{ε_{h}, h}^{2} + \frac{τ}{2} (2 + θ) (∥ \nabla e_{h}^{k} ∥^{2} + {∥ \nabla e_{h}^{k - 1} ∥}^{2}) + \frac{τ ρ}{T^{2}} L_{k}\} \\ - \frac{1}{2} \{∥ \nabla (e_{h}^{m + 1} - e_{h}^{m}) ∥^{2} + {∥ \nabla \cdot (e_{h}^{m + 1} - e_{h}^{m}) ∥}_{ε - 1}^{2}\} + \frac{1}{2} (∥ \nabla \cdot e_{h}^{m + 1} ∥_{ε - 1}^{2} + {∥ \nabla \cdot e_{h}^{m} ∥}_{ε - 1}^{2}) \\ + \frac{1}{2} \{∥ \nabla (e_{h}^{1} - e_{h}^{0}) ∥^{2} + {∥ \nabla \cdot (e_{h}^{1} - e_{h}^{0}) ∥}_{ε - 1}^{2}\} - \frac{1}{2} (∥ \nabla \cdot e_{h}^{1} ∥_{ε - 1}^{2} + {∥ \nabla \cdot e_{h}^{0} ∥}_{ε - 1}^{2}) . \end{matrix}

(4.40)

On the other hand, recalling (4.30) we note that for

2 \leq m \leq M - 1

\begin{matrix} \sum_{k = 1}^{m} L_{k} \leq \sum_{k = 1}^{m} (τ T \sum_{i = 1}^{m} {∥\frac{e_{h}^{i} - e_{h}^{i - 1}}{τ}∥}^{2} + {∥ e_{h}^{0} ∥}^{2}) = T m τ \sum_{k = 1}^{m} {∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥}^{2} + m {∥ e_{h}^{0} ∥}^{2} \end{matrix}

(4.41)

In short since

m τ \leq T

, from (4.41) we easily obtain for

2 \leq m \leq M - 1

\begin{matrix} \frac{τ ρ}{T^{2}} \sum_{k = 1}^{m} L_{k} \leq τ ρ \sum_{k = 1}^{m} {∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥}^{2} + \frac{ρ}{T} {∥ e_{h}^{0} ∥}^{2} . \end{matrix}

(4.42)

Plugging (4.42) into (4.40) and summing up both sides of (4.33) from

k = 1

through

k = m

for

2 \leq m \leq N - 1

by using (4.35), (4.40) yields:

\begin{matrix} \frac{1}{2} (2 - τ η) {∥\frac{e_{h}^{m + 1} - e_{h}^{m}}{τ}∥}_{ε_{h}, h}^{2} + \frac{1}{2} (1 - 2 τ) ∥ \nabla e_{h}^{m + 1} ∥^{2} + \frac{1}{2} (1 - τ θ) {∥ \nabla e_{h}^{m} ∥}^{2} \\ - \frac{1}{2} (2 - τ η) {∥\frac{e_{h}^{1} - e_{h}^{0}}{τ}∥}_{ε_{h}, h}^{2} - \frac{1}{2} (1 - τ θ) ∥ \nabla e_{h}^{0} ∥^{2} - \frac{1}{2} (1 - 2 τ) {∥ \nabla e_{h}^{1} ∥}^{2} \\ - \sum_{k = 1}^{m} \{τ (η + ρ) {∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥}_{ε_{h}, h}^{2} + \frac{τ}{2} (2 + θ) (∥ \nabla e_{h}^{k} ∥^{2} + {∥ \nabla e_{h}^{k - 1} ∥}^{2})\} \\ - \frac{1}{2} \{∥ \nabla (e_{h}^{m + 1} - e_{h}^{m}) ∥^{2} + {∥ \nabla \cdot (e_{h}^{m + 1} - e_{h}^{m}) ∥}_{ε - 1}^{2}\} + \frac{1}{2} (∥ \nabla \cdot e_{h}^{m + 1} ∥_{ε - 1}^{2} + {∥ \nabla \cdot e_{h}^{m} ∥}_{ε - 1}^{2}) \\ + \frac{1}{2} \{∥ \nabla (e_{h}^{1} - e_{h}^{0}) ∥^{2} + {∥ \nabla \cdot (e_{h}^{1} - e_{h}^{0}) ∥}_{ε - 1}^{2}\} - \frac{1}{2} (∥ \nabla \cdot e_{h}^{1} ∥_{ε - 1}^{2} + {∥ \nabla \cdot e_{h}^{0} ∥}_{ε - 1}^{2}) \\ \leq \sum_{k = 1}^{m} (\frac{τ}{2} | F^{k} |_{h}^{2} + \frac{N}{2 τ} {∥ d^{k} ∥}^{2} + \frac{τ}{2} {|G^{k}|}_{Γ_{h}}^{2}) + \frac{ρ}{T} {∥ e_{h}^{0} ∥}^{2} . \end{matrix}

(4.43)

Thus taking into account that

e_{h}^{1} - e_{h}^{0} = τ e_{1 h}

, leaving on the left hand side only the terms with superscripts

m + 1

and m, and increasing the coefficients of

∥ \nabla e_{h}^{j} ∥^{2}

for

j = 0, 1

and

∥ e_{1 h} ∥_{ε_{h}, h}^{2}

, for

2 \leq m \leq M - 1

we come up with:

\begin{matrix} \frac{1}{2} (2 - τ η) {∥\frac{e_{h}^{m + 1} - e_{h}^{m}}{τ}∥}_{ε_{h}, h}^{2} + \frac{1}{2} (1 - 2 τ) ∥ \nabla e_{h}^{m + 1} ∥^{2} + \frac{1}{2} (1 - τ θ) {∥ \nabla e_{h}^{m} ∥}^{2} \\ - \frac{1}{2} \{∥ \nabla (e_{h}^{m + 1} - e_{h}^{m}) ∥^{2} + {∥ \nabla \cdot (e_{h}^{m + 1} - e_{h}^{m}) ∥}_{ε - 1}^{2}\} + \frac{1}{2} (∥ \nabla \cdot e_{h}^{m + 1} ∥_{ε - 1}^{2} + {∥ \nabla \cdot e_{h}^{m} ∥}_{ε - 1}^{2}) \\ \leq \sum_{k = 1}^{m} \{τ (η + ρ) {∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥}_{ε_{h}, h}^{2} + \frac{τ}{2} (2 + θ) (∥ \nabla e_{h}^{k} ∥^{2} + {∥ \nabla e_{h}^{k - 1} ∥}^{2}) \\ + \frac{τ}{2} | F^{k} |_{h}^{2} + \frac{N}{2 τ} {∥ d^{k} ∥}^{2} + \frac{τ}{2} {|G^{k}|}_{Γ, h}^{2}\} + {∥e_{1 h}∥}_{ε_{h}, h}^{2} + \frac{τ^{2}}{2} (∥ \nabla e_{1 h} ∥^{2} + {∥ \nabla \cdot e_{1 h} ∥}_{ε - 1}^{2}) \\ + \frac{1}{2} (∥ \nabla e_{h}^{1} ∥^{2} + {∥ \nabla e_{h}^{0} ∥}^{2}) + \frac{1}{2} (∥ \nabla \cdot e_{h}^{1} ∥_{ε - 1}^{2} + {∥ \nabla \cdot e_{h}^{0} ∥}_{ε - 1}^{2}) + \frac{ρ}{T} {∥ e_{h}^{0} ∥}^{2} . \end{matrix}

(4.44)

Now we recall a classical inverse inequality (cf. [14]), according to which,

∥ \nabla v ∥ \leq C h^{- 1} ∥ v ∥ \leq C h^{- 1} {∥ v ∥}_{ε_{h}, h} for all v \in V_{h},

(4.45)

where C is a mesh-independent constant, and we apply the trivial upper bound

{∥ \nabla \cdot v ∥}_{ε - 1} \leq \sqrt{{N ∥ ε - 1 ∥}_{\infty}} ∥ \nabla v ∥

for all

v \in V_{h}

Let us assume that $τ$ satisfies the following CFL-condition:

τ \leq h / ν with ν = {C (1 + N ∥ ε - 1 ∥}_{\infty})^{1 / 2} .

(4.46)

Then we have,

\forall v \in V_{h}

{∥ \nabla v ∥}^{2} + {∥ \nabla \cdot v ∥}_{ε - 1}^{2} \leq ν^{2} \frac{τ^{2}}{h^{2}} \frac{{∥ v ∥}^{2}}{τ^{2}} \leq \frac{{∥ v ∥}_{ε_{h}, h}^{2}}{τ^{2}}

(4.47)

This means that

∥ \nabla (e_{h}^{m + 1} - e_{h}^{m}) ∥^{2} + {∥ \nabla \cdot (e_{h}^{m + 1} - e_{h}^{m}) ∥}_{ε - 1}^{2} \leq {∥\frac{e_{h}^{m + 1} - e_{h}^{m}}{τ}∥}_{ε_{h}, h}^{2}

(4.48)

Plugging (4.48) into (4.44) we obtain,

\begin{matrix} \frac{1}{2} (1 - τ η) {∥\frac{e_{h}^{m + 1} - e_{h}^{m}}{τ}∥}_{ε_{h}, h}^{2} + \frac{1}{2} (1 - 2 τ) ∥ \nabla e_{h}^{m + 1} ∥^{2} + \frac{1}{2} (1 - τ θ) {∥ \nabla e_{h}^{m} ∥}^{2} \\ \leq \sum_{k = 1}^{m} \{τ (η + ρ) {∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥}_{ε_{h}, h}^{2} + \frac{τ}{2} (2 + θ) (∥ \nabla e_{h}^{k} ∥^{2} + {∥ \nabla e_{h}^{k - 1} ∥}^{2}) \\ + \frac{τ}{2} | F^{k} |_{h}^{2} + \frac{N}{2 τ} {∥ d^{k} ∥}^{2} + \frac{τ}{2} {|G^{k}|}_{Γ, h}^{2}\} + {∥e_{1 h}∥}_{ε_{h}, h}^{2} + \frac{τ^{2}}{2} (∥ \nabla e_{1 h} ∥^{2} + {∥ \nabla \cdot e_{1 h} ∥}_{ε - 1}^{2}) \\ + \frac{1}{2} (∥ \nabla e_{h}^{1} ∥^{2} + {∥ \nabla e_{h}^{0} ∥}^{2}) + \frac{1}{2} (∥ \nabla \cdot e_{h}^{1} ∥_{ε - 1}^{2} + {∥ \nabla \cdot e_{h}^{0} ∥}_{ε - 1}^{2}) + \frac{ρ}{T} {∥ e_{h}^{0} ∥}^{2} . \end{matrix}

(4.49)

Next we note that both

1 - 2 τ

and

1 - τ θ

are bounded below by

1 - τ η

and moreover,

\frac{1}{2} (2 + θ)

is obviously bounded above by

η + ρ

. Therefore it is easy to see that (4.49) can be transformed into :

\begin{matrix} \frac{1}{2} (1 - τ η) ({∥\frac{e_{h}^{m + 1} - e_{h}^{m}}{τ}∥}_{ε_{h}, h}^{2} + ∥ \nabla e_{h}^{m + 1} ∥^{2} + {∥ \nabla e_{h}^{m} ∥}^{2}) \\ \leq S_{M} + E_{0} + \sum_{k = 1}^{m} τ (η + ρ) ({∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥}_{ε_{h}, h}^{2} + ∥ \nabla e_{h}^{k} ∥^{2} + {∥ \nabla e_{h}^{k - 1} ∥}^{2}), \end{matrix}

(4.50)

where

\begin{matrix} S_{M} : = \sum_{k = 1}^{M - 1} (\frac{τ}{2} | F^{k} |_{h}^{2} + \frac{N}{2 τ} {∥ d^{k} ∥}^{2} + \frac{τ}{2} {|G^{k}|}_{Γ, h}^{2}) \\ E_{0} : = {∥e_{1 h}∥}_{ε_{h}, h}^{2} + \frac{τ^{2}}{2} (∥ \nabla e_{1 h} ∥^{2} + {∥ \nabla \cdot e_{1 h} ∥}_{ε - 1}^{2}) \\ + \frac{1}{2} (∥ \nabla e_{h}^{1} ∥^{2} + {∥ \nabla e_{h}^{0} ∥}^{2}) + \frac{1}{2} (∥ \nabla \cdot e_{h}^{1} ∥_{ε - 1}^{2} + {∥ \nabla \cdot e_{h}^{0} ∥}_{ε - 1}^{2}) + \frac{ρ}{T} {∥ e_{h}^{0} ∥}^{2} . \end{matrix}

(4.51)

Let us assume that

τ \leq \frac{1}{2 η}

. Then from (4.51) we have

\begin{matrix} {∥\frac{e_{h}^{m + 1} - e_{h}^{m}}{τ}∥}_{ε_{h}, h}^{2} + ∥ \nabla e_{h}^{m + 1} ∥^{2} + {∥ \nabla e_{h}^{m} ∥}^{2} \\ \leq 4 \{S_{M} + E_{0} + \sum_{k = 1}^{m} τ (η + ρ) ({∥\frac{e_{h}^{k} - e_{h}^{k - 1}}{τ}∥}_{ε_{h}, h}^{2} + ∥ \nabla e_{h}^{k} ∥^{2} + {∥ \nabla e_{h}^{k - 1} ∥}^{2})\} . \end{matrix}

(4.52)

Now we set

β = 4 M τ (η + ρ) = 4 T \{2 + {| ε |}_{1, \infty} + (2 + T^{2}) {| ε |}_{2, \infty}\} .

(4.53)

From the discrete Grönwall’s Lemma and (3.15) from (4.52) we infer that

\forall m \leq M - 1

\begin{matrix} {∥\frac{e_{h}^{m + 1} - e_{h}^{m}}{τ}∥}^{2} + ∥ \nabla e_{h}^{m + 1} ∥^{2} + {∥ \nabla e_{h}^{m} ∥}^{2} \leq 4 (S_{M} + E_{0}) e^{β}, \end{matrix}

(4.54)

as long as

τ \leq min {h / ν, 1 / (2 η)}

, where

ν

η

ρ

and

β

are defined in (4.46), (4.34) and (4.53), respectively, and in the expression of

E_{0}

e_{h}^{1}

is to be replaced by

e_{h}^{0} + τ e_{1 h}

. ■

4.2. Scheme consistency

Before pursuing the reliability study of our scheme we need some approximation results related to the Maxwell’s equations. The arguments employed in this section found their inspiration in Thomée [33] and in Ruas [31].

4.2.1. Preliminaries

Henceforth we assume that

Ω

is a convex polygon for

N = 2

or a convex polyhedron for

N = 3

In this case one may reasonably assert that for every

g \in {L_{0}^{2} (Ω)}^{N}

the solution

v_{g} \in {H^{1} (Ω) \cap L_{0}^{2} (Ω)}^{N}

of the equation

\begin{matrix} - Δ v_{g} - \nabla {(ε - 1) \nabla \cdot v_{g}} = g & in Ω \\ \partial_{n} v_{g} = 0 & on Γ . \end{matrix}

(4.55)

belongs to

{H^{2} (Ω)}^{N}

Another result that we take for granted in this section is the existence of a constant $C_{ε}$ such that,

\forall g \in {L_{0}^{2} (Ω)}^{N} it holds ∥ H (v_{g}) ∥ \leq C_{ε} ∥ g ∥,

(4.56)

where

H (\cdot)

is the Hessian of a function or vector field.

(4.56) is a result whose grounds are found in analogous inequalities for convex polytopes applying to both the scalar Poisson problem and the linear elasticity system (or yet to the Stokes system) (cf. [26]). In fact (4.55) can be viewed as a problem half way between the vector Poisson equation with homogeneous Neumann boundary conditions whenever $ε \equiv 1$ , and a modified linear elasticity system with a smoothly varying Poisson ratio $ε - 1$ whenever $ε > 1$ . This fully justifies (4.56).

Now in order to establish the consistency of the explicit scheme (3.14) we first introduce an auxiliary field

{\hat{e}}_{h} (\cdot, t)

belonging to

V_{h}

for every

t \in [0, T]

, uniquely defined up to an additive field depending only on t as follows:

(\nabla {\hat{e}}_{h} {\cdot, t}, \nabla v) + {(\nabla \cdot {\hat{e}}_{h} {\cdot, t}, \nabla \cdot v)}_{ε - 1} = (\nabla e {\cdot, t}, \nabla v) + {(\nabla \cdot e {\cdot, t}, \nabla \cdot v)}_{ε - 1} \forall v \in V_{h} .

(4.57)

The time-dependent additive field up to which

{\hat{e}}_{h} {\cdot, t}

is defined can be determined by requiring that

\int_{Ω} {{\hat{e}}_{h} (\cdot, t) - e (\cdot, t)} d x = 0

\forall t \in [0, T]

Let us further assume that for every

t \in [0, T]

e (\cdot, t) \in {H^{2} (Ω)}^{N}

. In this case, from classical approximation results based on the interpolation error, we can assert that,

∥ \nabla {{\hat{e}}_{h} (\cdot, t) - e (\cdot, t)} ∥ + ∥ \nabla \cdot {{\hat{e}}_{h} (\cdot, t) - e (\cdot, t)} ∥_{ε - 1} \leq {\hat{C}}_{1} h ∥ H (e) (\cdot, t) ∥ \forall t \in [0, T],

(4.58)

where

{\hat{C}}_{1}

is a mesh-independent constant.

Let us show that there exists another mesh-independent constant ${\hat{C}}_{0}$ such that for every $t \in [0, T]$ it holds,

∥ {\hat{e}}_{h} (\cdot, t) - e (\cdot, t) ∥ \leq {\hat{C}}_{0} h^{2} ∥ H (e) (\cdot, t) ∥ \forall t \in [0, T] .

(4.59)

Since

{\hat{e}}_{h} (\cdot, t) - e (\cdot, t) \in {L_{0}^{2} (Ω)}^{N}

for every t, we may write:

∥ {\hat{e}}_{h} (\cdot, t) - e (\cdot, t) ∥ = sup_{g \in {L_{0}^{2} (Ω)}^{N} ∖ {0}} \frac{({\hat{e}}_{h} {\cdot, t} - e {\cdot, t}, g)}{∥ g ∥} \forall t \in [0, T] .

(4.60)

Defining

W : = {w | w \in {H^{2} (Ω) \cap L_{0}^{2} (Ω)}^{N} \partial_{n} w = 0 o n Γ}

, owing to (4.56) we have

\forall t \in [0, T]

∥ {\hat{e}}_{h} (\cdot, t) - e (\cdot, t) ∥ \leq C_{ε} sup_{w \in W ∖ {0}} \frac{({\hat{e}}_{h} {\cdot, t} - e {\cdot, t}, - Δ w - \nabla (ε - 1) \nabla \cdot w)}{∥ H (w) ∥} .

(4.61)

Since

\partial_{n} w = 0

and

ε = 1

Γ

we may integrate by parts the numerator in (4.61) to obtain for every

t \in [0, T]

∥ {\hat{e}}_{h} (\cdot, t) - e (\cdot, t) ∥ \leq C_{ε} sup_{w \in W ∖ {0}} \frac{(\nabla {{\hat{e}}_{h} {\cdot, t} - e {\cdot, t}}, \nabla w) + {(\nabla \cdot {{\hat{e}}_{h} (\cdot, t) - e (\cdot, t)}, \nabla \cdot w)}_{ε - 1}}{∥ H (w) ∥} .

(4.62)

Taking into account (4.57), the numerator of (4.62) can be rewritten as

(\nabla {{\hat{e}}_{h} {\cdot, t} - e {\cdot, t}}, \nabla {w - v}) + {(\nabla \cdot {{\hat{e}}_{h} {\cdot, t} - e {\cdot, t}}, \nabla \cdot {w - v})}_{ε - 1} \forall v \in V_{h} .

Then choosing

v

to be the

V_{h}

-interpolate of

w

, taking into account (4.58) we easily establish (4.59) with

{\hat{C}}_{0} = C_{ε} {\hat{C}}_{1}^{2}

To conclude these preliminary considerations, we refer to Chapter 5 of [31], to infer that the second order time-derivative $\partial_{t t} {\hat{e}}_{h} (\cdot, t)$ is well defined in ${H^{1} (Ω)}^{N}$ for every $t \in [0, T]$ , as long as $\partial_{t t} e (\cdot, t)$ lies in ${H^{1} (Ω)}^{N}$ for every $t \in [0, T]$ . Moreover, provided $\partial_{t t} e \in {H^{2} (Ω)}^{N}$ for every $t \in [0, T]$ , the following estimate holds:

∥ \partial_{t t} {{\hat{e}}_{h} (\cdot, t) - e (\cdot, t)} ∥ \leq {\hat{C}}_{0} h^{2} ∥ H (\partial_{t t} e) (\cdot, t) ∥ \forall t \in [0, T] .

(4.63)

In addition to the results given in this sub-section, we recall that, according to the Sobolev Embedding Theorem, there exists a constant

C_{T}

depending only on T such that it holds:

{∥ u ∥}_{L^{\infty} (0, T)} \leq C_{T} {{∥ u ∥}_{L^{2} (0, T)}^{2} + ∥ u_{t} {∥_{L^{2} (0, T)}^{2}}}^{1 / 2} \forall u \in H^{1} (0, T) .

(4.64)

In the remainder of this work we assume a certain regularity of

e

, namely,

Assumption* : The solution $e$ to equation (2.2) belongs to ${H^{4} (Ω \times (0, T))}^{N}$ .
Now taking $u = ∥ H (e) (\cdot, t) ∥$ we have $u_{t} (t) = \frac{\int_{Ω} {H (e) : H (\partial_{t} e)} (x, t) d x}{∥ {H (e)} (\cdot, t) ∥}$ , where : denotes the inner product of two constant tensors of order greater than or equal to three. Then by the Cauchy-Schwarz inequality and taking into account Assumption*, it trivially follows from (4.64) that the following upper bound holds:

∥ {H (e)} (\cdot, t) ∥ \leq C_{T} {{∥ ∥ H (e) ∥ ∥}_{L^{2} (0, T)}^{2} + ∥ ∥ H (\partial_{t} e) {∥ ∥}_{L^{2} (0, T)}^{2}}^{1 / 2} \forall t \in (0, T) .

(4.65)

In complement to the above ingredients we extend the inner products

{(\cdot, \cdot)}_{ε_{h}, h}

and

{(\cdot, \cdot)}_{Γ, h}

, and associated norms

{‖ \cdot ‖}_{ε_{h}, h}

and

{‖ \cdot ‖}_{Γ, h}

in a semi-definite manner, to fields in

{L^{2} (Ω)}^{N}

N = 2, 3

as follows:

First of all, $\forall K \in T_{h}$ let $Π_{K} : L^{2} (K) \to P_{1} (K)$ be the standard orthogonal projection operator onto the space $P_{1} (K)$ of linear functions in K. We set

\forall u and v \in {L^{2} (Ω)}^{N}, {(u, v)}_{ε_{h}, h} : = \sum_{K \in T_{h}} \frac{ε (G_{K}) m e a s (K)}{N + 1} \sum_{i = 1}^{N + 1} Π_{K} u (S_{K, i}) \cdot Π_{K} v (S_{K, i}) .

Let us generically denote by

F \subset Γ

an edge of a triangle or a face of a tetrahedron K such that

m e a s (K \cap Γ) > 0

. Moreover, we denote by

Π_{F} (v)

the standard orthogonal projection of a function

v \in L^{2} (F)

onto the space of linear functions on F. Similarly we define:

\forall u and v \in {L^{2} (Γ)}^{N}, {(u, v)}_{Γ, h} : = \sum_{F \subset Γ} \frac{m e a s (F)}{N} \sum_{i = 1}^{N} Π_{F} u (R_{F, i}) \cdot Π_{F} v (R_{F, i}) .

It is noteworthy that whenever

u

and

v

belong to

V_{h}

, both semi-definite inner products coincide with the inner products previously defined for such fields.

The following results hold in connection to the above inner products:

Lemma 4.1.

Let

Δ_{ε_{h}} (u, v) : = {(u, v)}_{ε_{h}, h} - {(u, v)}_{ε_{h}}

for

u, v \in {L^{2} (Ω)}^{N}

There exists a mesh independent constant

c_{Ω}

such that

\forall u \in {H^{1} (Ω)}^{N}

and

\forall v \in V_{h}

| Δ_{ε_{h}} (u, v) | \leq c_{Ω} {‖ ε ‖}_{0, \infty} {h ‖ \nabla u ‖ ‖ v ‖}_{h} ■ .

(4.66)

Lemma 4.2.

Let

Ξ_{h} (γ {u}, γ {v}) : = {(γ {u}, γ {v})}_{Γ, h} - {(γ {u}, γ {v})}_{Γ}

for

u, v \in {H^{1} (Ω)}^{N}

, where

γ {w}

represents the trace on Γ of a function

w \in H^{1} (Ω)

. Let also

\nabla_{Γ}

be the tangential gradient operator over Γ. There exists a mesh independent constant

{\tilde{c}}_{Γ}

such that

\forall u \in {H^{2} (Ω)}^{N}

and

\forall v \in V_{h}

| Ξ_{h} (γ {u}, γ {v}) | \leq {\tilde{c}}_{Γ} h ‖ \nabla_{Γ} {γ {u} ‖}_{Γ} {‖ γ {v} ‖}_{Γ, h} . ■ .

(4.67)

The proof of Lemma 4.1 is based on the Bramble-Hilbert Lemma, and we refer to [12] for more details. Lemma 4.2 in turn follows from the same arguments combined with the Trace Theorem, which ensures that

\nabla_{Γ} γ {u}

is well defined in

{L^{2} (Γ)}^{N}

u \in {H^{2} (Ω)}^{N}

. Incidentally the Trace Theorem allows us to bound above the right hand side of (4.67) in such a way that the following estimate also holds for another mesh independent constant

c_{Γ}

| Ξ_{h} (γ {u}, γ {v}) | \leq c_{Γ} {h {‖ \nabla u ‖}^{2} + {‖ H (u) ‖}^{2}}^{1 / 2} {‖ γ {v} ‖}_{Γ, h} . ■ .

(4.68)

To conclude we prove the validity of the following upper bounds:

Lemma 4.3.

\forall v \in {L^{2} (Ω)}^{N}

it holds

{∥ v ∥}_{ε_{h}, h} \leq \sqrt{N + 2} {‖ v ‖}_{ε_{h}} .

Proof: Denoting by

Π_{h} (v)

the function defined in

Ω

whose restriction to every

K \in T_{h}

Π_{K} (v)

for a given

v \in {L^{2} (Ω)}^{N}

, from an elementary property of the orthogonal projection we have

‖ Π_{h} {(v) ‖}_{ε_{h}} \leq {‖ v ‖}_{ε_{h}}, \forall v \in {L^{2} (Ω)} N .

(4.69)

Now taking v such that

v_{| K} \in P_{1} (K) \forall K \in T_{h}

, by a straightforward calculation using the expression of

v_{| K}

in terms of barycentric coordinates we have:

\int_{K} ε (G_{K}) v_{| K}^{2} d x = \frac{2 ε (G_{K}) m e a s (K)}{(N + 2) (N + 1)} \{\sum_{i = 1}^{N + 1} {| v (S_{K, i}) |}^{2} + \sum_{i = 1}^{N} \sum_{j = i + 1}^{N + 1} v (S_{K, i}) v (S_{K, j})\} .

It trivially follows that

\int_{K} ε (G_{K}) v_{| K}^{2} d x = \frac{ε (G_{K}) m e a s (K)}{(N + 2) (N + 1)} \{\sum_{i = 1}^{N + 1} {| v (S_{K, i}) |}^{2} + {(\sum_{i = 1}^{N + 1} v (S_{K, i}))}^{2}\} .

and finally

(N + 2) \int_{K} ε (G_{K}) v_{| K}^{2} d x \geq \frac{ε (G_{K}) m e a s (K)}{N + 1} \sum_{i = 1}^{N + 1} {| v (S_{K, i}) |}^{2} .

This immediately yields Lemma 4.3, taking into account (4.69). ■

Lemma 4.4.

\forall v \in {L^{2} (Γ)}^{N}

it holds

{∥ v ∥}_{Γ, h} \leq \sqrt{N + 1} {‖ v ‖}_{Γ}

. ■

The proof of this Lemma is based on arguments entirely analogous to Lemma 4.3.

4.2.2. Residual estimation

To begin with we define for

k = 0, 1, \dots, N

functions

{\hat{e}}_{h}^{k} \in V_{h}

{\hat{e}}_{h}^{k} (\cdot) : = {\hat{e}}_{h} (\cdot, k τ)

. In the sequel for any function or field

A

defined in

Ω \times (0, T)

A^{k} (\cdot)

denotes

A (\cdot, k τ)

, except for other quantities carrying the subscript h such as

e_{h}^{k}

Let us substitute

e_{h}^{k}

{\hat{e}}_{h}^{k}

for

k = 2, 3, \dots, N

on the left hand side of the first equation of (3.14) and take also as initial conditions

{\hat{e}}_{h}^{j}

instead of

e_{h}^{j}

j = 0, 1

The case of the initial conditions will be dealt with in the next section in the framework of the convergence analysis. As for the variational residual $E_{h}^{k} (v)$ resulting from the above substitution, where $E_{h}^{k}$ is a linear functional acting on $V_{h}$ , it can be expressed in the following manner:

\begin{matrix} E_{h}^{k} (v) = {({\partial_{t t} e}^{k}, v)}_{ε} + (\nabla e^{k}, \nabla v) + {(\nabla \cdot e^{k}, \nabla \cdot v)}_{ε - 1} + (\nabla ε \cdot e^{k}, \nabla \cdot v) + {({\partial_{t} e}^{k}, v)}_{Γ} \\ - {\bar{F}}^{k} (v) - ({\bar{d}}^{k}, \nabla \cdot v) - {\bar{G}}^{k} (v) \forall v \in V_{h}, \end{matrix}

(4.70)

where

\begin{matrix} {\bar{d}}^{k} = \nabla ε \cdot {{\hat{e}}_{h}^{k} - e^{k}}, \end{matrix}

(4.71)

and

{\bar{F}}^{k} (v)

and

{\bar{G}}^{k} (v)

can be written as follows:

\begin{matrix} {\bar{F}}^{k} (v) = F_{1}^{k} (v) + F_{2}^{k} (v) + F_{3}^{k} (v) + F_{4}^{k} (v), \\ with \\ F_{1}^{k} (v) = {(T^{k} ({\hat{e}}_{h} - e), v)}_{ε_{h}, h}; \\ F_{2}^{k} (v) = Δ_{ε_{h}} (T^{k} (e), v); \\ F_{3}^{k} (v) = ((ε_{h} - ε) T^{k} (e), v); \\ F_{4}^{k} (v) = {(T^{k} (e) - {\partial_{t t} e}^{k}, v)}_{ε}, \end{matrix}

(4.72)

T^{k}

being the finite-difference operator defined by,

T^{k} (\cdot) : = \frac{{\cdot}^{k + 1} - 2 {\cdot}^{k} + {\cdot}^{k - 1}}{τ^{2}},

(4.73)

and

\begin{matrix} {\bar{G}}^{k} (v) = G_{1}^{k} (v) + G_{2}^{k} (v) + G_{3}^{k} (v), \\ with \\ G_{1}^{k} (v) = {(Q^{k} ({\hat{e}}_{h} - e), v)}_{Γ, h}; \\ G_{2}^{k} (v) = Γ_{h} (Q^{k} (e), v); \\ G_{3}^{k} (v) = {(Q^{k} (e) - {\partial_{t} e}^{k}, v)}_{Γ}, \end{matrix}

(4.74)

Q^{k}

being the finite-difference operator defined by,

Q^{k} (\cdot) : = \frac{{\cdot}^{k + 1} - {\cdot}^{k - 1}}{2 τ} .

(4.75)

Notice that, under Assumption * both

\partial_{t} e (\cdot, t)

and

\partial_{t t} e (\cdot, t)

belong to

{H^{1} (Ω)}^{N}

for every

t \in [0, T]

. Hence we can define

\partial_{t} {\hat{e}}_{h}

from

\partial_{t} e

and

\partial_{t t} {\hat{e}}_{h}

from

\partial_{t t} e

in the same way as

{\hat{e}}_{h}

is defined from

e

. Moreover, straightforward calculations lead to,

T^{k} (\cdot) = \frac{1}{τ^{2}} \{\int_{(k - 1) τ}^{k τ} (\int_{t}^{k τ} \partial_{t t} {\cdot} d s) d t + \int_{k τ}^{(k + 1) τ} (\int_{k τ}^{t} \partial_{t t} {\cdot} d s) d t\} .

(4.76)

Furthermore another straightforward calculation allows us writing:

\begin{matrix} T^{k} (e) = {\partial_{t t} e}^{k} + R^{k} (e) \\ where \\ R^{k} (e) : = \frac{1}{2 τ^{2}} \{- \int_{(k - 1) τ}^{k τ} {t - (k - 1) τ}^{2} \partial_{t t t} e d t + \int_{k τ}^{(k + 1) τ} {(k + 1) τ - t}^{2} \partial_{t t t} e d t\} . \end{matrix}

(4.77)

Similarly,

Q^{k} (\cdot) = \frac{1}{2 τ} \int_{(k - 1) τ}^{(k + 1) τ} \partial_{t} {\cdot} d t,

(4.78)

and

\begin{matrix} Q^{k} (e) = {\partial_{t} e}^{k} + S^{k} (e), \\ where \\ S^{k} (e) = \frac{1}{2 τ} \{- \int_{(k - 1) τ}^{k τ} {t - (k - 1) τ} \partial_{t t} e d t + \int_{k τ}^{(k + 1) τ} {(k + 1) τ - t} \partial_{t t} e d t\} . \end{matrix}

(4.79)

Now we note that the sum of the terms on the first line of the expression of

E_{h}^{k} (v)

equals zero because they are just the left hand side of (2.7) at time

t = k τ

. Therefore the functions

{\hat{e}}_{h}^{k} \in V_{h}

are the solution of the following problem, for

k = 1, 2, \dots, M - 1

\begin{matrix} (ε \frac{{\hat{e}}_{h}^{k + 1} - 2 {\hat{e}}_{h}^{k} + {\hat{e}}_{h}^{k - 1}}{τ^{2}}, v) + (\nabla {\hat{e}}_{h}^{k}, \nabla v) + (\nabla \cdot {ε {\hat{e}}_{h}^{k}}, \nabla \cdot v) - (\nabla \cdot {\hat{e}}_{h}^{k}, \nabla \cdot v) \\ + {(\frac{{\hat{e}}_{h}^{k + 1} - {\hat{e}}_{h}^{k - 1}}{2 τ}, v)}_{Γ} = {\bar{F}}^{k} (v) + ({\bar{d}}^{k}, \nabla \cdot v) + {\bar{G}}^{k} (v) \forall v \in V_{h}, \\ {\hat{e}}_{h}^{0} (\cdot) = {\hat{e}}_{h} (\cdot, 0) and {\hat{e}}_{h}^{1} (\cdot) = {\hat{e}}_{h} (\cdot, τ) in Ω, \end{matrix}

(4.80)

{\bar{d}}^{k}

{\bar{F}}^{k}

and

{\bar{G}}^{k}

being given by (4.71), (4.72)-(4.73) and (4.74)-(4.75).

Estimating $∥ {\bar{d}}^{k} ∥$ is a trivial matter. Indeed, since $e^{k} \in {H^{2} (Ω)}^{N}$ , from (4.59) we immediately obtain,

∥ {\bar{d}}^{k} ∥ \leq {\hat{C}}_{0} h^{2} {| ε |}_{1, \infty} ∥ H (e^{k}) ∥ .

(4.81)

Therefore consistency of the scheme will result from suitable estimations of

| {\bar{F}}^{k} |_{h}

and

| {\bar{G}}^{k} |_{Γ, h}

in terms of

e

ε

τ

and h, which we next carry out.

First of all we search for upper bounds for the operators $T^{k} (\cdot)$ , $Q^{k} (\cdot)$ , $R^{k} (\cdot)$ and $S^{k} (\cdot)$ . With this aim we denote by $| \cdot |$ the euclidean norm of $ℜ^{n}$ , for $n \geq 1$ .
From (4.76) and the Cauchy-Schwarz inequality, we easily obtain for every $x \in Ω$ and $u$ such that ${\partial_{t t} u} (\cdot, t) \in {H^{2} (Ω)}^{N}$ $\forall t \in (0, T)$ ,

\begin{matrix} | T_{h}^{k} {u (x)} | \leq \frac{1}{τ^{2}} \{\int_{(k - 1) τ}^{k τ} (\int_{t}^{k τ} |{\partial_{t t} u} (x)| d s) d t + \int_{k τ}^{(k + 1) τ} (\int_{k τ}^{t} |{\partial_{t t} u} (x)| d s) d t\} \\ \leq \frac{1}{τ^{2}} \{\int_{(k - 1) τ}^{k τ} (\int_{(k - 1) τ}^{k τ} |{\partial_{t t} u} (x)| d s) d t + \int_{k τ}^{(k + 1) τ} (\int_{k τ}^{(k + 1) τ} |{\partial_{t t} u} (x)| d s) d t\} \\ = \frac{1}{τ} \int_{(k - 1) τ}^{(k + 1) τ} |{\partial_{t t} u} (x)| d t . \end{matrix}

It follows that, for every

u

such that

{\partial_{t t} u} (\cdot, t) \in {H^{2} (Ω)}^{N}

\forall t \in (0, T)

we have,

| T_{h}^{k} {u (x)} | \leq \sqrt{\frac{2}{τ}} {\{\int_{(k - 1) τ}^{(k + 1) τ} {|{\partial_{t t} u} (x)|}^{2} d t\}}^{1 / 2} .

(4.82)

Furthermore, from (4.77) and the inequality

a + b \leq {2 (a^{2} + b^{2})}^{1 / 2} \forall a, b \in ℜ

, for every

x \in Ω

we obtain:

\begin{matrix} | R^{k} {e (x)} | \leq \frac{1}{2} \{\int_{(k - 1) τ}^{k τ} |{\partial_{t t t} e} (x)| d t + \int_{k τ}^{(k + 1) τ} |{\partial_{t t t} e} (x)| d t\} \\ \leq \frac{\sqrt{τ}}{2} \{{(\int_{(k - 1) τ}^{k τ} {|{\partial_{t t t} e} (x)|}^{2} d t)}^{1 / 2} + {(\int_{k τ}^{(k + 1) τ} {|{\partial_{t t t} e} (x)|}^{2} d t)}^{1 / 2}\} . \end{matrix}

It follows that,

| R^{k} {e (x)} | \leq \sqrt{\frac{τ}{2}} {\{\int_{(k - 1) τ}^{(k + 1) τ} {|{\partial_{t t t} e} (x)|}^{2} d t\}}^{1 / 2} .

(4.83)

On the other hand from (4.78) and the Cauchy-Schwarz inequality we trivially have for every

x \in Γ

and

u

such that

γ {\partial_{t} u} (\cdot, t) \in {H^{3 / 2} (Γ)}^{N}

\forall t \in (0, T)

| Q^{k} {u (x)} | \leq \frac{1}{\sqrt{2 τ}} {\{\int_{(k - 1) τ}^{(k + 1) τ} {| {\partial_{t} u} (x) |}^{2} d t\}}^{1 / 2} .

(4.84)

Finally, similarly to (4.83), from (4.79) for every

x

Γ

we successively derive,

\begin{matrix} | S^{k} {e (x)} | \leq \frac{1}{2} \{\int_{(k - 1) τ}^{k τ} |{\partial_{t t} e} (x)| d t + \int_{k τ}^{(k + 1) τ} |{\partial_{t t} e} (x)| d t\} \\ \leq \frac{\sqrt{τ}}{2} \{{(\int_{(k - 1) τ}^{k τ} {|{\partial_{t t} e} (x)|}^{2} d t)}^{1 / 2} + {(\int_{k τ}^{(k + 1) τ} {|{\partial_{t t} e} (x)|}^{2} d t)}^{1 / 2}\} . \end{matrix}

Therefore it holds,

| S^{k} {e (x)} | \leq \sqrt{\frac{τ}{2}} {\{\int_{(k - 1) τ}^{(k + 1) τ} {|{\partial_{t t} e} (x)|}^{2} d t\}}^{1 / 2} .

(4.85)

Notice that bounds entirely analogous to (4.82) and (4.84) hold for

\nabla T^{k} (\cdot) = T^{k} (\nabla \cdot)

and

\nabla_{Γ} Q^{k} (\cdot) = Q^{k} (\nabla_{Γ} \cdot)

, that is,

\forall x \in Ω

\begin{matrix} | \nabla T_{h}^{k} (e) (x) | \leq \sqrt{\frac{2}{τ}} {\{\int_{(k - 1) τ}^{(k + 1) τ} {|{\partial_{t t} \nabla e} (x)|}^{2} d t\}}^{1 / 2}, \end{matrix}

(4.86)

and

\forall x \in Γ

| \nabla_{Γ} Q^{k} (e) (x) | \leq \frac{1}{\sqrt{2 τ}} {\{\int_{(k - 1) τ}^{(k + 1) τ} {| {\partial_{t} \nabla_{Γ} e} (x) |}^{2} d t\}}^{1 / 2} .

(4.87)

Next we estimate the four terms in the expression (4.72) of

F^{k} (v)

With the use of (4.82) and of Lemma 4.3 followed by a trivial manipulation, we successively have:

\begin{matrix} | F_{1}^{k} {(v) | \leq ‖ ε ‖}_{0, \infty} \sqrt{\frac{2}{τ}} {∥{\{\int_{(k - 1) τ}^{(k + 1) τ} {|{\partial_{t t} ({\hat{e}}_{h} - e)} (\cdot, t)|}^{2} d t\}}^{1 / 2}∥}_{h} {‖ v ‖}_{h} \\ \leq \sqrt{N + 2} {‖ ε ‖}_{0, \infty} \sqrt{\frac{2}{τ}} ∥{\{\int_{(k - 1) τ}^{(k + 1) τ} {|{\partial_{t t} ({\hat{e}}_{h} - e)} (\cdot, t)|}^{2} d t\}}^{1 / 2}∥ {‖ v ‖}_{h} \\ \leq \sqrt{N + 2} {‖ ε ‖}_{0, \infty} \sqrt{\frac{2}{τ}} {\{\int_{(k - 1) τ}^{(k + 1) τ} {∥{\partial_{t t} ({\hat{e}}_{h} - e)} (\cdot, t)∥}^{2} d t\}}^{1 / 2} {‖ v ‖}_{h} \end{matrix}

(4.88)

Recalling (4.72) and applying (4.63) to (4.88), we come up with,

| F_{1}^{k} (v) | \leq {\hat{C}}_{0} h^{2} \sqrt{\frac{2 (N + 2)}{τ}} {‖ ε ‖}_{0, \infty} {\{\int_{(k - 1) τ}^{(k + 1) τ} {‖ {H (\partial_{t t} e)} (\cdot, t) ‖}^{2} d t\}}^{1 / 2} {‖ v ‖}_{h} .

(4.89)

Next, combining (4.66) and (4.86) we immediately obtain.

| F_{2}^{k} (v) | \leq c_{Ω} {h ‖ ε ‖}_{0, \infty} \sqrt{\frac{2}{τ}} {\{\int_{(k - 1) τ}^{(k + 1) τ} {|{\partial_{t t} \nabla e} (\cdot, t)∥}^{2} d t\}}^{1 / 2} {‖ v ‖}_{h} .

(4.90)

Further, taking into account (3.15), (4.76), (4.86) and the standard estimate

‖ ε_{h} {- ε ‖}_{0, \infty} \leq C_{\infty} h {| ε |}_{1, \infty}

where

C_{\infty}

is a mesh-independent constant, it holds

| F_{3}^{k} (v) | \leq C_{\infty} {h | ε |}_{1, \infty} \sqrt{\frac{2}{τ}} {\{\int_{(k - 1) τ}^{(k + 1) τ} {∥{\partial_{t t} e} (\cdot, t)∥}^{2} d t\}}^{1 / 2} {‖ v ‖}_{h}

(4.91)

Finally by (3.15), (4.77) and (4.83), we have

| F_{4}^{k} {(v) | \leq ‖ ε ‖}_{0, \infty} \sqrt{\frac{τ}{2}} {\{\int_{(k - 1) τ}^{(k + 1) τ} {∥{\partial_{t t t} e} (\cdot, t)∥}^{2} d t\}}^{1 / 2} {‖ v ‖}_{h} .

(4.92)

Now we turn our attention to the three terms in the expression (4.74) of

G^{k} (v)

. First of all, owing to Assumption* and standard error estimates, we can write for a suitable mesh-independent constant

{\hat{C}}_{1}

‖ {\nabla \partial_{t} ({\hat{e}}_{h} - e)} (\cdot, t) ‖ \leq {\hat{C}}_{1} h ‖ {H (\partial_{t} e)} (\cdot, t)) ‖, \forall t \in [0, T] .

(4.93)

On the other hand, by the Trace Theorem there exists a contant

C_{T r}

depending only on

Ω

such that,

‖ {\partial_{t} ({\hat{e}}_{h} - e)} (\cdot, t) ‖_{Γ} \leq C_{T r} {‖ {\partial_{t} ({\hat{e}}_{h} - e)} (\cdot, t) ‖^{2} + ‖ {\nabla \partial_{t} ({\hat{e}}_{h} - e)} (\cdot, t) {‖^{2}}}^{1 / 2} \forall t \in [0, T] .

(4.94)

Hence by (4.63), (4.93) and (4.94), we have, for a suitable mesh-independent constant

C_{B}

‖ {\partial_{t} ({\hat{e}}_{h} - e)} (\cdot, t) ‖_{Γ} \leq C_{B} h ‖ {H (\partial_{t} e)} (\cdot, t) ‖ \forall t \in [0, T] .

(4.95)

Now recalling (4.72) and taking into account (4.87) and Lemma 4.4, similarly to (4.88) we first obtain:

| G_{1}^{k} (v) | \leq \frac{2}{\sqrt{2 τ}} {\{\int_{(k - 1) τ}^{(k + 1) τ} {∥{\partial_{t} ({\hat{e}}_{h} - e)} (\cdot, t)∥}_{Γ}^{2} d t\}}^{1 / 2} {‖ v ‖}_{Γ, h} .

(4.96)

Then using (4.95) we immediately establish,

| G_{1}^{k} (v) | \leq C_{B} h \sqrt{\frac{2}{τ}} {\{\int_{(k - 1) τ}^{(k + 1) τ} {‖ {H (\partial_{t} e)} (\cdot, t) ‖}^{2} d t\}}^{1 / 2} {‖ v ‖}_{Γ, h} .

(4.97)

Next we switch to

G_{2}^{k}

. Similarly to (4.90), (4.68) and (4.87) yield

| G_{2}^{k} (v) | \leq c_{Γ} h \sqrt{\frac{2}{τ}} {\{\int_{(k - 1) τ}^{(k + 1) τ} (‖ {\nabla \partial_{t} e} (\cdot, t) ‖^{2} + {‖ {H (\partial_{t} e)} (\cdot, t) ‖}^{2}) d t\}}^{1 / 2} {‖ v ‖}_{Γ, h} .

(4.98)

As for

G_{3}^{k}

, taking into account (4.79) together with (4.85) and (3.16), we obtain:

| G_{3}^{k} (v) | \leq \sqrt{\frac{τ}{2}} {\{\int_{(k - 1) τ}^{(k + 1) τ} {∥{\partial_{t t} e} (\cdot, t)∥}_{Γ}^{2} d t\}}^{1 / 2} {‖ v ‖}_{Γ, h} .

(4.99)

Then using the Trace Theorem (cf. (4.94)), we finally establish,

| G_{3}^{k} (v) | \leq C_{T r} \sqrt{\frac{τ}{2}} {\{\int_{(k - 1) τ}^{(k + 1) τ} (‖ {\partial_{t t} e} (\cdot, t) ‖^{2} + {‖ {\nabla \partial_{t t} e} (\cdot, t) ‖}^{2}) d t\}}^{1 / 2} {‖ v ‖}_{Γ, h} .

(4.100)

Now collecting (4.89), (4.90), (4.91) and (4.92) we can write,

\begin{matrix} | {\bar{F}}^{k} |_{h} \leq {\hat{C}}_{0} h^{2} \sqrt{\frac{2 (N + 2)}{τ}} {‖ ε ‖}_{0, \infty} {\{\int_{(k - 1) τ}^{(k + 1) τ} {‖ {H (\partial_{t t} e)} (\cdot, t) ‖}^{2} d t\}}^{1 / 2} \\ + c_{Ω} h {‖ ε ‖}_{0, \infty} \sqrt{\frac{2}{τ}} {\{\int_{(k - 1) τ}^{(k + 1) τ} {|{\partial_{t t} \nabla e} (\cdot, t)∥}^{2} d t\}}^{1 / 2} \\ + C_{\infty} h {| ε |}_{1, \infty} \sqrt{\frac{2}{τ}} {\{\int_{(k - 1) τ}^{(k + 1) τ} {∥{\partial_{t t} \nabla e} (\cdot, t)∥}^{2} d t\}}^{1 / 2} \\ + {‖ ε ‖}_{0, \infty} \sqrt{\frac{τ}{2}} {\{\int_{(k - 1) τ}^{(k + 1) τ} {∥{\partial_{t t t} e} (\cdot, t)∥}^{2} d t\}}^{1 / 2} . \end{matrix}

(4.101)

On the other hand (4.97), (4.98) and (4.100) yield,

\begin{matrix} | {\bar{G}}^{k} |_{Γ, h} \leq C_{B} h \sqrt{\frac{2}{τ}} {\{\int_{(k - 1) τ}^{(k + 1) τ} {‖ {H (\partial_{t} e)} (\cdot, t) ‖}^{2} d t\}}^{1 / 2} \\ + c_{Γ} h \sqrt{\frac{2}{τ}} {\{\int_{(k - 1) τ}^{(k + 1) τ} (‖ {\nabla \partial_{t} e} (\cdot, t) ‖^{2} + {‖ {H (\partial_{t} e)} (\cdot, t) ‖}^{2}) d t\}}^{1 / 2} \\ + C_{T r} \sqrt{\frac{τ}{2}} {\{\int_{(k - 1) τ}^{(k + 1) τ} (‖ {\partial_{t t} e} (\cdot, t) ‖^{2} + {‖ {\nabla \partial_{t t} e} (\cdot, t) ‖}^{2}) d t\}}^{1 / 2} . \end{matrix}

(4.102)

Then, taking into account (4.80) and the stability condition (4.46), by inspection we can assert that the consistency of scheme (3.14) is an immediate consequence of (4.81), (4.101) and (4.102).

4.3. Convergence results

In order to establish the convergence of scheme (3.14) we combine the stability and consistence results obtained in the previous subsections. With this aim we define

{\bar{e}}_{h}^{k} : = {\hat{e}}_{h}^{k} - e_{h}^{k}

for

k = 0, 1, 2, \dots, M

. By linearity we can assert that the variational residual on the left hand side of the first equation of (3.14) for

k = 2, 3, \dots, M

, when the

e_{h}^{k}

s are replaced with the

{\bar{e}}_{h}^{k}

s, and

e_{h}^{j}

is replaced with

{\bar{e}}_{h}^{j}

for

j = 0, 1

, is exactly

E_{h}^{k} (v)

, since the residual corresponding to the

e_{h}^{k}

’s vanishes by definition. The initial conditions

{\bar{e}}_{h}^{j}

for

j = 0

and

j = 1

corresponding to the thus modified problem have to be estimated. This is the purpose of the next subsection.

4.3.1. Initial-condition deviations

Here we turn our attention to the estimate of

{\bar{E}}_{0}

, which accounts for the deviation in the initial conditions appearing in the stability inequality (4.54) that applies to the modification of (3.14) when

e_{h}^{k}

is replaced by

{\bar{e}}_{h}^{k}

Let us first define,

\begin{matrix} {\hat{e}}_{1 h} : = ({\hat{e}}_{h}^{1} - {\hat{e}}_{h}^{0}) / τ \\ {\bar{e}}_{1 h} : = {\hat{e}}_{1 h} - e_{1 h}, \\ {\bar{E}}_{0} : = {∥{\bar{e}}_{1 h}∥}_{ε_{h}, h}^{2} + \frac{τ^{2}}{2} (∥ \nabla {\bar{e}}_{1 h} ∥^{2} + {∥ \nabla \cdot {\bar{e}}_{1 h} ∥}_{ε - 1}^{2}) \\ + \frac{1}{2} (∥ \nabla {\bar{e}}_{h}^{1} ∥^{2} + {∥ \nabla {\bar{e}}_{h}^{0} ∥}^{2}) + \frac{1}{2} (∥ \nabla \cdot {\bar{e}}_{h}^{1} ∥_{ε - 1}^{2} + {∥ \nabla \cdot {\bar{e}}_{h}^{0} ∥}_{ε - 1}^{2}) + {T | ε |}_{2, \infty} {∥ {\bar{e}}_{h}^{0} ∥}^{2} \end{matrix}

(4.103)

Recalling that

e_{h}^{1} = e_{h}^{0} + τ e_{1 h}

we have

{\bar{e}}_{h}^{1} = {\bar{e}}_{h}^{0} + τ {\bar{e}}_{1 h}

. Thus, taking either

A = \nabla {\bar{e}}_{h}^{0}

A = \nabla \cdot {\bar{e}}_{h}^{0}

and either

B = τ \nabla {\bar{e}}_{1 h}

B = τ \nabla \cdot {\bar{e}}_{1 h}

, we apply twice the inequality

{∥ A + B ∥}^{2} / 2 \leq {∥ A ∥}^{2} + {∥ B ∥}^{2}

to (4.103) together with Lemma 4.3, to obtain,

\begin{matrix} {\bar{E}}_{0} \leq \frac{3}{2} \{∥ \nabla {\bar{e}}_{h}^{0} ∥^{2} + {∥ \nabla \cdot {\bar{e}}_{h}^{0} ∥}_{ε - 1}^{2} + τ^{2} (∥ \nabla {\bar{e}}_{1 h} ∥^{2} + {∥ \nabla \cdot {\bar{e}}_{1 h} ∥}_{ε - 1}^{2})\} \\ + {T | ε |}_{2, \infty} ∥ {\bar{e}}_{h}^{0} ∥^{2} + {(N + 2) ∥ ε ∥}_{0, \infty} {∥ {\bar{e}}_{1 h} ∥}^{2} . \end{matrix}

(4.104)

Finally using the inequality

{∥ \nabla \cdot v ∥}_{ε - 1}^{2} \leq {N ∥ ε - 1 ∥}_{0, \infty} {∥ \nabla v ∥}^{2}

\forall v \in {H^{1} (Ω)}^{N}

, after straightforward manipulations we easily infer from (4.104) that

\begin{matrix} {\bar{E}}_{0} \leq \frac{{3 + 3 N ∥ ε - 1 ∥}_{0, \infty}}{2} (∥ \nabla {\bar{e}}_{h}^{0} ∥^{2} + τ^{2} {∥ \nabla {\bar{e}}_{1 h} ∥}^{2}) + {T | ε |}_{2, \infty} ∥ {\bar{e}}_{h}^{0} ∥^{2} + {(N + 2) ∥ ε ∥}_{0, \infty} {∥ {\bar{e}}_{1 h} ∥}^{2} . \end{matrix}

(4.105)

We next use the obvious splitting

{\bar{e}}_{h}^{0} = ({\hat{e}}_{h}^{0} - e_{0}) + (e_{0} - e_{h}^{0})

, together with the trivial one,

{\bar{e}}_{1 h} = \frac{{\hat{e}}_{h}^{1} - {\hat{e}}_{h}^{0}}{τ} - e_{1 h} = ({\partial_{t} {\hat{e}}_{h}}_{| t = 0} - e_{1}) + \frac{1}{τ} \int_{0}^{τ} (τ - t) \partial_{t t} {\hat{e}}_{h} d t + (e_{1} - e_{1 h}) .

(4.106)

Then plugging (4.106) into (4.105), since

{\{\sum_{i = 1}^{p} ∥ A_{i} ∥\}}^{2} \leq p \sum_{i = 1}^{p} {∥ A_{i} ∥}^{2}

for any set of p functions or fields

A_{i}

, we obtain:

\begin{matrix} {\bar{E}}_{0} \leq \frac{{3 + 3 N ∥ ε - 1 ∥}_{0, \infty}}{2} \{2 \{∥ \nabla ({\hat{e}}_{h}^{0} - e_{0}) ∥^{2} + {∥ \nabla (e_{0} - e_{h}^{0}) ∥}^{2}\} \\ 3 τ^{2} \{∥ \nabla \partial_{t} ({\hat{e}}_{h} - e) (\cdot, 0) ∥^{2} + \frac{1}{τ^{2}} {∥\int_{0}^{τ} (τ - t) \partial_{t t} \nabla {\hat{e}}_{h} d t∥}^{2} + {∥ \nabla (e_{1} - e_{1 h}) ∥}^{2}\}\} \\ + {2 T | ε |}_{2, \infty} (∥ {\hat{e}}_{h}^{0} - e_{0} ∥^{2} + {∥ e_{0} - e_{h}^{0} ∥}^{2}) \\ + {3 (N + 2) ∥ ε ∥}_{0, \infty} \{∥ \partial_{t} ({\hat{e}}_{h} - e) (\cdot, 0) ∥^{2} + \frac{1}{τ^{2}} {∥\int_{0}^{τ} (τ - t) \partial_{t t} {\hat{e}}_{h} d t∥}^{2} + {∥ e_{1} - e_{1 h} ∥}^{2}\} . \end{matrix}

(4.107)

Owing to a trivial upper bound and to the Cauchy-Schwarz inequality we easily come up with

\begin{matrix} \frac{1}{τ^{2}} {∥\int_{0}^{τ} (τ - t) \nabla \partial_{t t} {\hat{e}}_{h} d t∥}^{2} \leq {∥\int_{0}^{τ} | \nabla \partial_{t t} {\hat{e}}_{h} | d t∥}^{2} \leq τ {∥{\{\int_{0}^{τ} {| \nabla \partial_{t t} {\hat{e}}_{h} |}^{2} d t\}}^{1 / 2}∥}^{2} . \end{matrix}

It follows that,

\frac{1}{τ^{2}} {∥\int_{0}^{τ} (τ - t) \nabla \partial_{t t} {\hat{e}}_{h} d t∥}^{2} \leq τ \int_{0}^{τ} {∥\nabla \partial_{t t} {\hat{e}}_{h}∥}^{2} d t \leq {(1 + N ∥ ε - 1 ∥}_{0, \infty}) τ \int_{0}^{τ} {∥\nabla \partial_{t t} e∥}^{2} d t .

(4.108)

The last inequality in (4.108) follows from the definition of

{\hat{e}}_{h}

. Indeed we know that,

(\nabla \partial_{t t} {\hat{e}}_{h}, \nabla v) + {(\nabla \cdot \partial_{t t} {\hat{e}}_{h}, \nabla \cdot v)}_{ε - 1} = (\nabla \partial_{t t} e, \nabla v) + {(\nabla \cdot \partial_{t t} e, \nabla \cdot v)}_{ε - 1} \forall v \in V_{h} .

(4.109)

Taking

v = {\partial_{t t} {\hat{e}}_{h}} (\cdot, t)

, by the Cauchy-Schwarz inequality we easily obtain:

{∥ \nabla \partial_{t t} {\hat{e}}_{h} ∥^{2} + ∥ \nabla \cdot \partial_{t t} {\hat{e}}_{h} ∥_{ε - 1}^{2}}^{1 / 2} \leq {∥ \nabla \partial_{t t} {e ∥}^{2} + ∥ \nabla \cdot \partial_{t t} e {∥_{ε - 1}^{2}}}^{1 / 2},

(4.110)

or yet

∥ \nabla \partial_{t t} {\hat{e}}_{h} ∥ \leq \sqrt{{1 + N ∥ ε - 1 ∥}_{0, \infty}} ∥ \nabla \partial_{t t} e ∥,

(4.111)

which validates (4.108).

On the other hand according to (4.63) we have $∥ \partial_{t t} {\hat{e}}_{h} ∥ \leq ∥ \partial_{t t} e ∥ + {\hat{C}}_{0} h^{2} ∥ H (\partial_{t t} e) ∥$ . This yields,

{∥\int_{0}^{τ} (τ - t) \partial_{t t} {\hat{e}}_{h} d t∥}^{2} \leq τ \int_{0}^{τ} {∥\partial_{t t} {\hat{e}}_{h}∥}^{2} d t \leq 2 τ \{\int_{0}^{τ} {∥\partial_{t t} e∥}^{2} d t + {\hat{C}}_{0}^{2} h^{4} \int_{0}^{τ} {∥H (\partial_{t t} e)∥}^{2} d t\} .

(4.112)

Incidentally we point out that Assumption * and (4.64) allow us to assert that

$\partial_{t t} \nabla e \in {L^{2} (Ω \times (0, T))}^{N}$ ;
$∥ \partial_{t t} e ∥ \in {L^{\infty} (0, T)}^{N}$ ;
${\partial_{t} e}_{| t = 0} = e_{1} \in {H^{2} (Ω)}^{N}$ ;
${e}_{| t = 0} = e_{0} \in {H^{2} (Ω)}^{N}$ .

Clearly enough, besides (4.59) and (4.58), we will apply to (4.107) standard estimates based on the interpolation error in Sobolev norms (cf. [14]), together with the following obvious variants of (4.58) and (4.63), namely,

\begin{matrix} ∥\nabla \partial_{t} ({\hat{e}}_{h} - e) (\cdot, 0)∥ \leq {\hat{C}}_{1} h ∥H (e_{1})∥, \\ ∥\partial_{t} ({\hat{e}}_{h} - e) (\cdot, 0)∥ \leq {\hat{C}}_{0} h^{2} ∥H (e_{1})∥ . \end{matrix}

(4.113)

Then taking into account that

τ \leq 1 / (2 η)

, from (4.107)-(4.108)-(4.112)-(4.113) and Assumption*, we conclude that there exists a constant

{\bar{C}}_{0}

depending on

Ω

, T and

ε

, but neither on h nor on

τ

, such that,

\begin{matrix} {\bar{E}}_{0} \leq {\bar{C}}_{0}^{2} \{(h^{2} + τ h^{2} + h^{4}) ∥ H (e_{0}) ∥^{2} + (h^{2} + τ^{2} h^{2}) ∥ H (e_{1}) ∥^{2} + τ^{3} ∥ ∥ \nabla \partial_{t t} e {∥ ∥}_{L^{2} (0, T)}^{2} \\ + τ^{2} ∥ ∥ \partial_{t t} e {∥ ∥}_{L^{\infty} (0, T)}^{2} + τ h^{4} \int_{0}^{τ} {∥H (\partial_{t t} e)∥}^{2} d t\} . \end{matrix}

(4.114)

Notice that, starting from (4.64) with

u = ∥ \nabla \partial_{t t} e ∥

, similarly to (4.65), we obtain

∥ ∥ \partial_{t t} {e ∥ ∥}_{L^{\infty} (0, T)} \leq C_{T} {∥ ∥ \partial_{t t} {e ∥ ∥}_{L^{2} (0, T)}^{2} + ∥ ∥ \partial_{t t t} {e) ∥ ∥}_{L^{2} (0, T)}^{2}}^{1 / 2} .

Thus noting that

h \leq d i a m (Ω)

, using again the upper bound

τ \leq 1 / (2 η)

and extending the integral to the whole interval

(0, T)

in (4.114), from the latter inequality we infer the existence of another constant

C_{0}

independent of h and

τ

, such that,

\begin{matrix} {\bar{E}}_{0} \leq C_{0}^{2} \{h^{2} {∥ H (e_{0}) ∥^{2} + ∥ H (e_{1}) ∥^{2} + ∥ H (\partial_{t t} e) ∥_{L^{2} (Ω \times (0, T))}^{2}} \\ + τ^{2} {∥ \partial_{t t} {e ∥}_{L^{2} (Ω \times (0, T))}^{2} + ∥ \partial_{t t t} {e ∥}_{L^{2} (Ω \times (0, T))}^{2} + ∥ \nabla \partial_{t t} e ∥_{L^{2} (Ω \times (0, T))}^{2}}\} . \end{matrix}

(4.115)

4.3.2. Error estimates

In order to fully exploit the stability inequality (4.54) we further define,

{\bar{S}}_{M} : = \sum_{k = 1}^{M - 1} (\frac{τ}{2} | {\bar{F}}^{k} |_{h}^{2} + \frac{N}{2 τ} {∥ {\bar{d}}^{k} ∥}^{2} + \frac{τ}{2} {|{\bar{G}}^{k}|}_{Γ, h}^{2}) .

(4.116)

According to (4.80), in order to estimate

{\bar{S}}_{M}

under the regularity Assumption* on

e

, we resort to the estimates (4.81), (4.101) and (4.102). Using the inequality

| a \cdot b | \leq | a | | b |

for

a, b \in ℜ^{n}

, it is easy to see that there exists two constants

{\bar{C}}_{F}

and

{\bar{C}}_{G}

independent of h and

τ

such that

\sum_{k = 1}^{M - 1} \frac{τ}{2} | {\bar{F}}^{k} |_{h}^{2} \leq {\bar{C}}_{F} (h^{2} + τ^{2}) {| e |}_{H^{4} (Ω \times (0, T))}^{2},

(4.117)

\sum_{k = 1}^{M - 1} \frac{τ}{2} | {\bar{G}}^{k} |_{Γ, h}^{2} \leq {\bar{C}}_{G} (h^{2} + τ^{2}) {| e |}_{H^{3} (Ω \times (0, T))}^{2} .

(4.118)

On the other hand, recalling (4.65) we have

k = 1, 2, \dots, M - 1

∥ H (e^{k}) ∥^{2} \leq C_{T}^{2} \int_{0}^{T} \{{∥ H (e) ∥}^{2} + {∥ H (\partial_{t} e) ∥}^{2}\} d t .

(4.119)

Therefore, since

M = T / τ

we have:

\sum_{k = 1}^{M - 1} {∥ H (e^{k}) ∥}^{2} \leq \frac{C_{T}^{2} T}{τ} \int_{0}^{T} \{{∥ H (e) ∥}^{2} + {∥ H (\partial_{t} e) ∥}^{2}\} d t

(4.120)

It follows from (4.120) and (4.81) that,

\sum_{k = 1}^{M - 1} \frac{N}{2 τ} ∥ {\bar{d}}^{k} ∥^{2} \leq {\hat{C}}_{0}^{2} C_{T}^{2} T {| ε |}_{1, \infty}^{2} \frac{N h^{4}}{2 τ^{2}} ({| e |}_{H^{2} (Ω \times (0, T))}^{2} + {| e |}_{H^{3} (Ω \times (0, T))}^{2})

(4.121)

Plugging (4.117), (4.121) and (4.118) into (4.116) we can assert that there exists a constant

\tilde{C}

depending on

Ω

, T and

ε

, but neither on h nor on

τ

, such that,

\begin{matrix} {\bar{S}}_{M} \leq {\tilde{C}}^{2} (h^{2} + τ^{2} + \frac{h^{4}}{τ^{2}}) {∥ e ∥}_{H^{4} (Ω \times (0, T))}^{2} \end{matrix}

(4.122)

Now recalling (4.54)-(4.53), provided (4.46) holds, together with

τ \leq 1 / (2 η)

, we have:

{\{{∥\frac{{\bar{e}}_{h}^{m + 1} - {\bar{e}}_{h}^{m}}{τ}∥}^{2} + ∥ \nabla {\bar{e}}_{h}^{m + 1} ∥^{2} + {∥ \nabla {\bar{e}}_{h}^{m} ∥}^{2}\}}^{1 / 2} \leq 2 \sqrt{{\bar{S}}_{M} + {\bar{E}}_{0}} e^{β / 2} .

(4.123)

This implies that, for

m = 1, 2, \dots, M - 1

, it holds:

\begin{matrix} {\{{∥\frac{e_{h}^{m + 1} - e_{h}^{m}}{τ} - \frac{e^{m + 1} - e^{m}}{τ}∥}^{2} + ∥ \nabla (e_{h}^{m + 1} - e^{m + 1}) ∥^{2} + {∥ \nabla (e_{h}^{m} - e^{m})) ∥}^{2}\}}^{1 / 2} \\ \leq {\{{∥\frac{{\hat{e}}_{h}^{m + 1} - e^{m + 1}}{τ} - \frac{{\hat{e}}_{h}^{m} - e^{m}}{τ}∥}^{2} + ∥ \nabla ({\hat{e}}_{h}^{m + 1} - e^{m + 1}) ∥^{2} + {∥ \nabla ({\hat{e}}_{h}^{m} - e^{m}) ∥}^{2}\}}^{1 / 2} \\ + 2 \sqrt{{\bar{S}}_{M} + {\bar{E}}_{0}} e^{β / 2} . \end{matrix}

(4.124)

Let us define a function

e_{h}

\bar{Ω} \times [0, T]

whose value at

t = k τ

equals

e_{h}^{k}

for

k = 1, 2, \dots, M

and that varies linearly with t in each time interval

[(k - 1) τ, k τ]

, in such a way that

\partial_{t} e_{h} (x, t) = \frac{e_{h}^{k} (x) - e_{h}^{k - 1} (x)}{τ}

for every

x \in \bar{Ω}

and

t \in ((k - 1) τ, k τ)

Now we define ${\tilde{A}}^{m - 1 / 2} (\cdot)$ for any function or field $A (\cdot, t)$ to be the mean value of $A (\cdot, t)$ in $(m τ, (m + 1) τ$ , that is ${\tilde{A}}^{m + 1 / 2} = τ^{- 1} \int_{m τ}^{(m + 1) τ} A (\cdot, t) d t$ . Clearly enough we have

{∥\frac{e_{h}^{m + 1} - e_{h}^{m}}{τ} - \frac{e^{m + 1} - e^{m}}{τ}∥}^{2} = {∥{\tilde{\partial_{t} e_{h}}}^{m + 1 / 2} - {\tilde{\partial_{t} e}}^{m + 1 / 2}∥}^{2} .

(4.125)

and also recalling (4.120) and (4.65)

\begin{matrix} {∥\frac{{\hat{e}}_{h}^{m + 1} - e^{m + 1}}{τ} - \frac{{\hat{e}}_{h}^{m} - e^{m}}{τ}∥}^{2} = {∥\int_{m τ}^{(m + 1) τ} \frac{\partial_{t} {\hat{e}}_{h} - \partial_{t} e}{τ} d t∥}^{2} \\ \leq {∥\frac{{\{\int_{m τ}^{(m + 1) τ} {(\partial_{t} {\hat{e}}_{h} - \partial_{t} e)}^{2} d t\}}^{1 / 2}}{\sqrt{τ}}∥}^{2} = \frac{1}{τ} \int_{m τ}^{(m + 1) τ} {∥\partial_{t} {\hat{e}}_{h} - \partial_{t} e∥}^{2} d t \\ \leq \frac{{\hat{C}}_{0}^{2} h^{4}}{τ} \int_{m τ}^{(m + 1) τ} ∥ H (\partial_{t} e) ∥^{2} d t \leq C_{T}^{2} {\hat{C}}_{0}^{2} h^{4} \int_{0}^{T} {∥ H (\partial_{t} e) ∥^{2} + ∥ H (\partial_{t t} e) ∥^{2}} d t, \end{matrix}

which implies that

\begin{matrix} {∥\frac{{\hat{e}}_{h}^{m + 1} - e^{m + 1}}{τ} - \frac{{\hat{e}}_{h}^{m} - e^{m}}{τ}∥}^{2} \leq C_{T}^{2} {\hat{C}}_{0}^{2} h^{2} {| d i a m (Ω) |}^{2} {∥ e ∥}_{H^{4} (Ω \times (0, T))}^{2} . \end{matrix}

(4.126)

On the other hand from (4.58) and (4.65) we have for

j = m

j = m + 1

∥ \nabla ({\hat{e}}_{h}^{j} - e^{j}) ∥^{2} \leq {\hat{C}}_{1}^{2} h^{2} {∥ H (e^{j}) ∥}^{2} \leq {\hat{C}}_{1}^{2} C_{T}^{2} h^{2} \int_{0}^{T} \{{∥ H (e) ∥}^{2} + {∥ H (\partial_{t} e) ∥}^{2}\} d t .

which yields,

∥ \nabla ({\hat{e}}_{h}^{j} - e^{j}) ∥^{2} \leq C_{T}^{2} {\hat{C}}_{1}^{2} h^{2} {∥ e ∥}_{H^{3} (Ω \times (0, T))}^{2} .

(4.127)

Now using Taylor expansions about

t = (m + 1 / 2) τ

together with some arguments already exploited in this work, it is easy to establish that for a certain constant

\tilde{C}

independent of h and

τ

it holds,

∥{\tilde{\partial_{t} e}}^{m + 1 / 2} - {\partial_{t} e}^{m + 1 / 2}∥ \leq \tilde{C} τ^{2} {∥e∥}_{H^{4} (Ω \times (0, T))},

(4.128)

where for every function g defined in

\bar{Ω} \times [0, T]

and

s \in ℜ^{+}

g^{s}

is the function defined in

Ω

g^{s} (\cdot) = g (\cdot, s τ)

Noting that ${\tilde{\partial_{t} e_{h}}}^{m + 1 / 2} (\cdot)$ is nothing but ${\partial_{t} e_{h}}^{m + 1 / 2} (\cdot)$ , collecting (4.124), (4.125), (4.126), (4.127), (4.128), together with (4.115) and (4.122), we have thus proved the following convergence result for scheme (3.14):

Provided the CFL condition (4.46) is satisfied and

τ

also satisfies

τ \leq 1 / (2 η)

, under Assumption * on

e

, there exists a constant

C

depending only on

Ω

ε

and T such that

\begin{matrix} \begin{matrix} max_{1 \leq m \leq M - 1} ∥{\partial_{t} (e_{h} - e)}^{m + 1 / 2}∥ + max_{2 \leq m \leq M} ∥ \nabla (e_{h}^{m} - e^{m}) ∥ \\ \leq C (τ + h + h^{2} / τ) \{{∥ e ∥}_{H^{4} (Ω \times (0, T))} + ∥ H (e_{0}) ∥ + ∥ H (e_{1}) ∥\} . \end{matrix} \\ ■ \end{matrix}

(4.129)

In short, as long as

τ

varies linearly with h, first order convergence of scheme (3.14) in terms of either

τ

or h is thus established in the sense of the norms on the left hand side of (4.129).

5. Assessment of the scheme

The purpose of this section is to validate the theoretical results given in Section 4 by means of numerical experiments for

N = 2

. With this aim every partition

T_{h}

is assumed to fit both

Ω

and

Ω_{i n}

in the usual way. Let

Ω_{h}

be the union of all the triangles in

T_{h}

. If we replace

Ω

Ω_{h}

in the scheme (3.14), it is not difficult to see that, in the case of convex domains, the error estimates (4.129) extend to a curved domain

Ω

, as long as the norm

∥ \cdot ∥

is replaced by the norm

{∥ \cdot ∥}_{Ω_{h}}

. Actually, unlike we did so far, in this section the notation

{∥ \cdot ∥}_{h}

will rather stand for the later norm, for the sake of brevity.

We performed numerical tests for the model problem (2.2), taking

T = 0.5

and

Ω

to be the unit disk given by

Ω = {(x_{1}; x_{2}) | r < 1}

, where

r : = \sqrt{x_{1}^{2} + x_{2}^{2}}

H_{1 / 2}

being the Heaviside function, for an integer

m > 1

we take

ε (r) = 1 + {(1 - 4 r^{2})}^{m} (1 - H_{1 / 2} (r)) and set v_{θ} (r, t) : = \frac{e^{r - 2 t}}{ε (r)} .

(5.130)

We consider that the exact solution

e = (e_{1}, e_{2})

of (2.2) is given by

e_{1} = - x_{2} v_{θ} (r, t) and e_{2} = x_{1} v_{θ} (r, t) .

It is easy to check that

e

satisfies absorbing conditions on the boundary of the unit disk.

The initial data

e_{0}

and

e_{1}

are given by

e_{0} : = e_{| t = 0} = \frac{(- x_{2}, x_{1}) e^{r}}{ε (r)} and e_{1} : = {\partial_{t} e}_{| t = 0} = \frac{2 (x_{2}, - x_{1}) e^{r}}{ε (r)} .

Since

\nabla \cdot (ε e) \equiv 0

Ω

by construction, the right hand side

f = (f_{1}; f_{2})

is given by

f : = ε \partial_{t t} e - Δ e + \nabla (\nabla \cdot e)

, that is,

\begin{matrix} f = & (- 4 x_{2} e^{r - 2 t} - Δ e_{1}; 4 x_{1} e^{r - 2 t} - Δ e_{2}), \\ where & (Δ e_{1}; Δ e_{2}) = (- \frac{3 x_{2}}{r} v_{θ}^{'} - x_{2} v_{θ}^{''}; \frac{3 x_{1}}{r} v_{θ}^{'} + x_{1} v_{θ}^{''}); \\ with & v_{θ}^{'} = \frac{ε - ε^{'}}{ε^{2}} e^{r - 2 t}, v_{θ}^{''} = \frac{ε^{2} - 2 ε ε^{'} - ε ε^{''} + 2 {(ε^{'})}^{2}}{ε^{3}} e^{r - 2 t} \\ and & ε^{'} = - 8 m r {(1 - 4 r^{2})}^{m - 1} (1 - H_{1 / 2} (r)), \\ ε^{''} = 8 m (8 m^{2} r^{2} - 4 r^{2} - 1) {(1 - 4 r^{2})}^{m - 2} (1 - H_{1 / 2} (r)) . \end{matrix}

Figure 1. Function

ε (x, y)

in the domain

Ω = (0, 1) \times (0, 1)

for different values of m in (5.130) on the mesh with

h = 2^{- 5}

Figure 1. Function

ε (x, y)

in the domain

Ω = (0, 1) \times (0, 1)

for different values of m in (5.130) on the mesh with

h = 2^{- 5}

In our computations we used the software package Waves [34] for the finite element method applied to the solution of the model problem (2.2). The spatial domain

Ω

is discretized by a family of quasi-uniform meshes consisting of

2 \times 2^{2 l + 2}

triangles K for

l = 1, . . ., 6

, constructed as a certain mapping of the same number of triangles of the uniform mesh of the square

(- 1, 1) \times (- 1, 1)

, which is symmetric with respect to the cartesian axes and have their edges parallel to those axes and either to the line

x_{1} = x_{2}

x_{1} x_{2} \geq 0

or to the line

x_{1} = - x_{2}

otherwise. The above mapping is defined by means of a suitable transformation of cartesian into polar coordinates. For each value of l we define a reference mesh size

h_{l}

equal to

2^{- l}

We consider a partition of the time domain

(0, T)

into time intervals

J = (t_{k - 1}, t_{k}]

of equal length

τ_{l}

for a given number of time intervals

M, l = 1, . . ., 6

. We performed numerical tests taking

m = 2, 3, 4, 5

in (5.130). We choose the time step

τ_{l} = 0.025 \times 2^{- l}, l = 1, . . ., 6

, which provides numerical stability for all meshes. We computed the maximum value over the time steps of the relative errors measured in the

L^{2}

-norms of the function, its gradient and its time-derivative in the polygon

Ω_{h}

for the different meshes in use. Now

e

and

e_{h}

being the exact and approximate solution of (2.2), setting

M : = T / τ_{l}

, these quantities are represented by,

\begin{matrix} e_{l}^{1} & = \frac{max_{1 \leq k \leq M} {∥ e^{k} - e_{h}^{k} ∥}_{h}}{max_{1 \leq k \leq M} {∥ e^{k} ∥}_{h}}, \\ e_{l}^{2} & = \frac{max_{1 \leq k \leq M} {∥ \nabla (e^{k} - e_{h}^{k}) ∥}_{h}}{max_{1 \leq k \leq M} {∥ \nabla e^{k} ∥}_{h}}, \\ e_{l}^{3} & = \frac{max_{1 \leq k \leq M - 1} {∥ {\partial_{t} (e - e_{h})}^{k + 1 / 2} ∥}_{h}}{max_{1 \leq k \leq M - 1} {∥ {\partial_{t} e}^{k + 1 / 2} ∥}_{h}} . \end{matrix}

(5.131)

In Table 1, Table 2, Table 3 and Table 4 method’s convergence in these three senses is observed taking

m = 2, 3, 4, 5

in (5.130). Figure 2 shows convergence rates of our numerical scheme based on a

P_{1}

space discretization, taking the function

ε

defined by (5.130) with

m = 2

(on the left) and

m = 3

(on the right) for

ε (r)

. Similar convergence results are presented in Figure 3 taking

m = 4

(on the left) and

m = 5

(on the right) in (5.130).

Table 1. Computed maximum relative errors

e_{l}

in maximum energy, maximum

L^{2}

and maximum in broken time error on different meshes with mesh sizes

h_{l} = 2^{- l}, l = 1, . . ., 6

for the function

ε

with

m = 2

in (5.130).

Table 1. Computed maximum relative errors

e_{l}

in maximum energy, maximum

L^{2}

and maximum in broken time error on different meshes with mesh sizes

h_{l} = 2^{- l}, l = 1, . . ., 6

for the function

ε

with

m = 2

in (5.130).

l	$n e l$	$n n o$	$e_{l}^{(1)}$	$e_{l - 1}^{(1)} / e_{l}^{(1)}$	$e_{l}^{(2)}$	$e_{l - 1}^{(2)} / e_{l}^{(2)}$	$e_{l}^{(3)}$	$e_{l - 1}^{(3)} / e_{l}^{(3)}$
1	32	25	0.6057		2.2827		3.1375
2	128	81	0.1499	4.0418	1.0769	2.1198	1.1536	2.7196
3	512	289	0.0333	4.5007	0.4454	2.4178	0.5776	1.9972
4	2048	1089	0.0078	4.2466	0.2077	2.1449	0.2802	2.0617
5	8192	4225	0.0019	4.1288	0.1066	1.9483	0.1379	2.0313
6	32768	16641	0.0005	4.0653	0.0535	1.9905	0.0690	1.9981

Table 2. Computed maximum relative errors

e_{l}

in maximum energy, maximum

L^{2}

and maximum in broken time error on different meshes with mesh sizes

h_{l} = 2^{- l}, l = 1, . . ., 6

for the function

ε

with

m = 3

in (5.130).

Table 2. Computed maximum relative errors

e_{l}

in maximum energy, maximum

L^{2}

and maximum in broken time error on different meshes with mesh sizes

h_{l} = 2^{- l}, l = 1, . . ., 6

for the function

ε

with

m = 3

in (5.130).

l	$n e l$	$n n o$	$e_{l}^{(1)}$	$e_{l - 1}^{(1)} / e_{l}^{(1)}$	$e_{l}^{(2)}$	$e_{l - 1}^{(2)} / e_{l}^{(2)}$	$e_{l}^{(3)}$	$e_{l - 1}^{(3)} / e_{l}^{(3)}$
1	32	25	0.6144		1.8851		3.0462
2	128	81	0.1511	4.0666	1.0794	1.7464	1.1417	2.6682
3	512	289	0.0339	4.4553	0.4713	2.2904	0.5680	2.0099
4	2048	1089	0.0080	4.2216	0.2166	2.1753	0.2760	2.0583
5	8192	4225	0.0019	4.1207	0.1137	1.9049	0.1354	2.0381
6	32768	16641	0.0005	4.0615	0.0566	2.0092	0.0677	1.9997

Table 3. Computed maximum relative errors

e_{l}

in maximum energy, maximum

L^{2}

and maximum in broken time error on different meshes with mesh sizes

h_{l} = 2^{- l}, l = 1, . . ., 6

for the function

ε

with

m = 4

in (5.130).

Table 3. Computed maximum relative errors

e_{l}

in maximum energy, maximum

L^{2}

and maximum in broken time error on different meshes with mesh sizes

h_{l} = 2^{- l}, l = 1, . . ., 6

for the function

ε

with

m = 4

in (5.130).

l	$n e l$	$n n o$	$e_{l}^{(1)}$	$e_{l - 1}^{(1)} / e_{l}^{(1)}$	$e_{l}^{(2)}$	$e_{l - 1}^{(2)} / e_{l}^{(2)}$	$e_{l}^{(3)}$	$e_{l - 1}^{(3)} / e_{l}^{(3)}$
1	32	25	0.6122		1.9517		3.0545
2	128	81	0.1529	4.0027	1.0896	1.7912	1.1445	2.6689
3	512	289	0.0346	4.4266	0.4879	2.2331	0.5639	2.0296
4	2048	1089	0.0082	4.2069	0.2234	2.1839	0.2728	2.0667
5	8192	4225	0.0020	4.1151	0.1183	1.8879	0.1336	2.0418
6	32768	16641	0.0005	4.0585	0.0595	1.9890	0.0668	2.0008

Table 4. Computed maximum relative errors

e_{l}

in maximum energy, maximum

L^{2}

and maximum in broken time error on different meshes with mesh sizes

h_{l} = 2^{- l}, l = 1, . . ., 6

for the function

ε

with

m = 5

in (5.130).

Table 4. Computed maximum relative errors

e_{l}

in maximum energy, maximum

L^{2}

and maximum in broken time error on different meshes with mesh sizes

h_{l} = 2^{- l}, l = 1, . . ., 6

for the function

ε

with

m = 5

in (5.130).

l	$n e l$	$n n o$	$e_{l}^{(1)}$	$e_{l - 1}^{(1)} / e_{l}^{(1)}$	$e_{l}^{(2)}$	$e_{l - 1}^{(2)} / e_{l}^{(2)}$	$e_{l}^{(3)}$	$e_{l - 1}^{(3)} / e_{l}^{(3)}$
1	32	25	0.6107		1.9930		3.0603
2	128	81	0.1546	3.9505	1.1006	1.8108	1.1464	2.6696
3	512	289	0.0351	4.4031	0.4982	2.2090	0.5619	2.0403
4	2048	1089	0.0084	4.1954	0.2288	2.1777	0.2706	2.0765
5	8192	4225	0.0020	4.1106	0.1223	1.8715	0.1325	2.0417
6	32768	16641	0.0005	4.0561	0.0607	2.0139	0.0662	2.0011

Figure 2. Maximum in time of relative errors for

m = 2

(left) and

m = 3

(right)

Figure 2. Maximum in time of relative errors for

m = 2

(left) and

m = 3

(right)

Figure 3. Maximum in time of relative errors for

m = 4

(left) and

m = 5

(right)

Figure 3. Maximum in time of relative errors for

m = 4

(left) and

m = 5

(right)

Observation of these tables and figures clearly indicates that our scheme behaves like a first order method in the (semi-)norm of

L^{\infty} [(0, T); H^{1} (Ω)]

for

e

and in the norm of

L^{\infty} [(0, T); L^{2} (Ω)]

for

\partial_{t} e

for all the chosen values of m. As far as the values of m greater or equal to 4 are concerned this perfectly conforms to the a priori error estimates given in Section 4. However, those tables and figures also show that such theoretical predictions extend to cases not considered in our analysis such as

m = 2

and

m = 3

, in which the regularity of the exact solution is lower than assumed. Otherwise stated, some of our assumptions seem to be of academic interest only and a lower regularity of the solution such as

H^{2} (Ω \times (0, T))

should be sufficient to attain optimal first order convergence in both senses. On the other hand second-order convergence can be expected from our scheme in the norm of

L^{\infty} [(0, T); L^{2} (Ω)]

for

e

, according to Table 1, Table 2, Table 3 and Table 4 and Figure 2 and Figure 3.

6. Final remarks

As previously noted, the approach advocated in this work was extensively and successfully tested in the framework of the solution of CIPs governed by Maxwell’s equations. More specifically it was used with minor modifications to solve both the direct problem and the adjoint problem, as steps of an adaptive algorithm to determine the unknown dielectric permittivity. More details on this procedure can be found in [9] and references therein.

As a matter of fact, the method studied in this paper was designed to handle composite dielectrics structured in such a way that layers with higher permittivity are completely surrounded by layers with a (constant) lower permittivity, say

ε = 1

. It should be noted, however, that the assumption that

ε

attains a minimum in the outer layer was made here only to simplify things. Actually, under the same assumptions (4.129) also applies to the case where

ε

in inner layers is allowed to be smaller than in the outer layer, say

ε < 1

. We refer to [6] for further details.

Another issue that is worth a comment is the practical calculation of the term

(\nabla \cdot ε e_{h}^{k}, \nabla \cdot v)

in (3.14). Unless

ε

is a simple function such as a polynomial, it is not possible to compute this term exactly. That is why we recommend the use of the trapezoidal rule to carry out these computations. At the price of small adjustments in some terms involving norms of

ε

, the thus modified scheme is stable in the same sense as (4.54). Moreover, the qualitative convergence result (4.129) remains unchanged, provided a little more regularity is required from

ε

. We skip details for the sake of brevity.

Author Contributions

Conceptualization, L.B. and V.R.; methodology, L. B. and V.R.; software, L.B. and V.R.; validation, L.B. and V.R.; investigation, L.B. and V.R.; writing—original draft preparation, V.R.; writing—review and editing, L.B and V.R.; visualization, L.B. All authors have read and agreed to the published version of the manuscript.

Funding

The research of the first author is supported by the Swedish Research Council grant VR 2018-03661.

Institutional Review Board Statement

Not applicable

Informed Consent Statement

Not applicable

Data Availability Statement

Not applicable

Conflicts of Interest

The authors declare no conflict of interest.

References

R.A. Adams, Sobolev Spaces, Academic Press, N.Y., 1975.
F. Assous, P. Degond, E. Heintze, P.-A. Raviart, On a finite-element method for solving the three-dimensional Maxwell equations, J. Computational Physics, 109 (1993), 222–237. [CrossRef]
S. Badia, R. Codina, A Nodal-based Finite Element Approximation of the Maxwell Problem Suitable for Singular Solutions, SIAM J. Numer. Anal., 50-2 (2012), 398–417. [CrossRef]
L. Beilina, M. Grote, Adaptive Hybrid Finite Element/Difference Method for Maxwell’s equations, TWMS J. of Pure and Applied Mathematics, 1-2 (2010), 176-197. [CrossRef]
L. Beilina, Energy estimates and numerical verification of the stabilized Domain Decomposition Finite Element/Finite Difference approach for time-dependent Maxwell’s system, Cent. Eur. J. Math., 11-4 (2013) 702-733. [CrossRef]
L. Beilina, V. Ruas, An explicit P₁ finite element scheme for Maxwell’s equations with constant permittivity in a boundary neighborhood, arXiv:1808.10720, 2020.
L. Beilina, V. Ruas, Convergence of Explicit P1 Finite-Element Solutions to Maxwell’s Equations, Springer Proceedings in Mathematics and Statistics, vol 328, p.91-103, Springer, 2020.
L. Beilina, V. Ruas, On the Maxwell-wave equation coupling problem and its explicit finite-element solution, Applications of Mathematics, open access (2022). [CrossRef]
J. Bondestam-Malmberg, L. Beilina, An adaptive finite element method in quantitative reconstruction of small inclusions from limited observations, Appl. Maths & Information Sci., 12-1 (2018), 1–19.
A. Bossavit, Computational Electromagnetism, Variational Formulations, Complementary, Edge Elements, in: Electromagnetism, vol. 2, Academic Press, New York, 1998.
F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from Lagrange multipliers, RAIRO Analyse Numérique, 8-2 (1974), 129-151.
J.H. Carneiro de Araujo, P.D. Gomes, V. Ruas, Study of a finite element method for the time-dependent generalized Stokes system associated with viscoelastic flow, J. Computational Applied Mathematics, 234-8 (2010), 2562-2577. [CrossRef]
V.C. Chen and W.V. Wahl, Das Rand-Anfangswertproblem für quasilineare Wellen-gleichungen in Sobolevräumen niedriger Ordnung, Z. reine angew. Math., 337 (1982), 77-112. [CrossRef]
P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North Holland, 1978.
P. Ciarlet Jr, J. Zou, Fully discrete finite element approaches for time-dependent Maxwell’s equations, Numerische Mathematik, 82-2 (1999), 193–219. [CrossRef]
P. Ciarlet Jr, Augmented formulations for solving Maxwell equations, Computer Methods in Applied Mechanics and Engineering, 194-2–5 (2005), 559–586. [CrossRef]
P. Ciarlet Jr, E. Jamelot, Continuous Galerkin methods for solving the time-dependent Maxwell equations in 3D geometries, J. Computational Physics, 226-1 (2007), 1122-1135. [CrossRef]
G. Cohen, P. Monk, Gauss point mass lumping schemes for Maxwell’s equations, Numerical Methods for Partial Differential Equations, 14-1 (1998), 63–88.
M. Costabel, A coercive bilinear form for Maxwell’s equations, J. Math. Anal. Appl., 157-2 (1991), 527–541. [CrossRef]
M. Costabel, M. Dauge, Singularities of Maxwell’s equations on polyhedral domains, Analysis, Numerics and Applications of Differential and Integral Equations, Pitman Research Notes in Mathematics Series, Vol. 379, p.69–76, 1998.
M. Costabel, M. Dauge, Weighted regularization of Maxwell’s equations in polyhedral domains. A rehabilitation of nodal finite elements, Numerische Mathematik, 93-2 (2002), 239–277. [CrossRef]
A. Ditkowski and D. Gottlieb, On the Engquist–Majda absorbing boundary conditions for hyperbolic systems, Contemporary Mathematics, 330 (2003), 55–72.
A. Elmkies and P. Joly, Finite elements and mass lumping for Maxwell’s equations: the 2D case. Numerical Analysis, C. R. Acad.Sci.Paris, 324, pp. 1287–1293, 1997.
B. Engquist, A. Majda, Absorbing boundary conditions for the numerical simulation of waves, Mathematics of Computation, 31 (1977), 629-651.
V. Girault, P.A. Raviart, Finite Element Methods for Navier-Stokes Equations, Springer Series in Computational Mathematics, Vol. 5, Springer-Verlag, Berlin, 1986.
P. Grisvard, Sigularities in Boundary Value Problems, Masson, Paris, 1992.
E. Jamelot, Résolution des équations de Maxwell avec des éléments finis de Galerkin continus, thèse doctorale, École Polytechnique, Palaiseau, France, 2005.
P. Monk. Finite Element Methods for Maxwell’s Equations, Clarendon Press, 2003.
J.-C. Nédélec, Mixed finite elements in R³, Numerische Mathematik, 35 (1980), 315-341. [CrossRef]
P.-A. Raviart and J.-M. Thomas, Mixed Finite Element Methods for Second Order Elliptic Problems, Lecture Notes in Mathematics, Springer Verlag, 606: 292-315, 1977.
V. Ruas, Numerical Methods for Partial Differential Equations; an Introduction, Wiley, 2016.
V. Ruas and M.A. Silva Ramos, A Hermite Method for Maxwell’s Equations, Appl. Maths & Information Sci., 12-2 (2018), 271–283. [CrossRef]
V. Thomée, Galerkin Finite Element Methods for Parabolic Problems, Springer Series in Computational Mathematics, Second edition, 1997.
WavES, the software package, http://www.waves24.com.

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

Explicit P1 Finite Element Solution of the Maxwell-Wave Equation Coupling Problem with Absorbing b. c.

Abstract

1. Introduction

2. The model problem

2.1. Notations and reminders

2.2. Well-posedness considerations

2.3. Variational form

3. Space-time discretization

3.1. Space semi-discretization

3.2. Full discretization

4. Reliability analysis

4.1. Scheme stability

4.2. Scheme consistency

4.2.1. Preliminaries

4.2.2. Residual estimation

4.3. Convergence results

4.3.1. Initial-condition deviations

4.3.2. Error estimates

5. Assessment of the scheme

6. Final remarks

Author Contributions

Funding

Institutional Review Board Statement

Informed Consent Statement

Data Availability Statement

Conflicts of Interest

References

MDPI Initiatives

Important Links

Subscribe