You are currently viewing a beta version of our website. If you spot anything unusual, kindly let us know.

Preprint
Article

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

Altmetrics

Downloads

69

Views

28

Comments

0

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

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 ) 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 ( H ( c u r l ) 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 ) 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 Ω of N , N = 2 , 3 with boundary Γ , that we consider in this work is posed in the following setting.
  • First, we consider that Ω = Ω ¯ i n Ω 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 Ω ¯ i n with respect to Ω (the boundary of Γ o u t is Γ Γ i n ). The dielectric permittivity denoted by ε is assumed to belong to C 2 ( Ω ¯ ) and to fulfill the conditions:
ε 1   in   Ω out   and   ε 1   otherwise .
Let n be the unit outer normal vector on Γ . We denote by n ( · ) the outer normal derivative of a field on Γ . Taking into account conditions (2.1), we may prescribe wave-equation first order absorbing boundary conditions n e = t e on Γ × ( 0 , T ) [22,24], as seen hereafter.
Now, we are given e 0 { H 1 ( Ω ) } N and e 1 H ( d i v , Ω ) satisfying · ( ε e 0 ) = · ( ε e 1 ) 0 , together with f L 2 [ ( 0 , T ) ; H ( d i v , Ω ) ] satisfying · f 0 . Then, setting V : = { H 1 ( Ω ) } N , the problem to solve is:
Find   e V : = H 2 [ ( 0 , T ) ; { L 2 ( Ω ) } N ] L 2 [ ( 0 , T ) ; V ] such   that     ε t t e + × × e = f , · { ε e } = 0 . in   Ω i n × ( 0 , T ) , t t e Δ e = f in   Ω o u t × ( 0 , T ) , with   the   initial   conditions e ( · , 0 ) = e 0 ( · ) ,   and   t e ( · , 0 ) = e 1 ( · ) in   Ω and   the   boundary   conditions n e + t e = 0 on   Γ × ( 0 , T ) .
(2.2) tells us that the wave equation holds in Ω o u t × ( 0 , T ) . On the other hand, the equations in braces of (2.2) make up the Maxwell’s equations in Ω i n × ( 0 , T ) for the sole electric field. It is noticeable that, since f is solenoidal by assumption, the second equation in Ω i n × ( 0 , T ) is superfluous i.e. reduntant. Indeed, taking the divergence of both sides of the first equation in Ω i n × ( 0 , T ) and denoting by u the function · ( ε e ) , we have t t u = 0   in   Ω i n × ( 0 , T ) . Taking into account that u | t = 0 = t u | t = 0 = 0 by our assumption on e 0 and e 1 , it must hold · ( ε e ) = 0 in Ω i n × ( 0 , T ) . However, the same conclusion cannot be drawn for Ω o u t × ( 0 , T ) . This is because in this sub-domain u : = · 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 0 in Ω o u t × ( 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 ( Ω ¯ ) by | · | n , for n > 0 and the standard norm of C 0 ( Ω ¯ ) by · 0 , . A subset of Ω ¯ will be denoted by an upper case Latin letter with or without a subscript. For any D Ω we denote the standard inner product of { L 2 ( D ) } N by ( · , · ) D and the corresponding norm by { · } 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 ω L ( Ω ) we introduce the weighted L 2 ( Ω ) -semi-norm { · } ω : = Ω ω | { · } | 2 , which is actually a norm if ω 0 everywhere in Ω ¯ . The notation ( A , B ) ω expresses Ω ω A · B , A , B being two square-integrable fields in Ω . If ω is strictly positive this expression defines an inner product associated with the norm { · } ω .
  • We denote by < · , · > s , Γ i n the duality product of H s ( Γ i n ) H s ( Γ i n ) for s + .
In the sequel we shall repeatedly employ a well known operator identity applying to vector fields, namely,
Δ · + × × .
Let D be a bounded domain of N with boundary D . We recall that, according to the Trace Theorem (cf. [1]) and well known results, if a given function χ L 2 ( D ) has a well defined trace on D in the space H 1 / 2 ( D ) and Δ χ H 1 ( D ) , then necessarily χ 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 { H 1 ( Ω o u t ) } N , z H ( c u r l ; Ω i n ) H ( d i v ; Ω i n ) , Ω o u t u + Ω i n w = 0 } and let R ε be the subspace of R of the pairs ( v ; z ) satisfying · ( ε z ) = 0 in Ω i n . Noticing that the trace over Γ i n of a field in H ( c u r l ; Ω i n ) H ( d i v ; Ω i n ) } lies in { H 1 / 2 ( Γ i n ) } N (cf. [25]), recalling the space L 0 2 ( Ω ) : = { v | v L 2 ( Ω ) , Ω v = 0 } , we set the problem
Given   a   solenoidal   g { L 0 2 ( Ω ) } N   find   ( ( u ; w ) ; p ) R ε × { H 1 / 2 ( Γ i n ) } N such   that   ( ( v ; z ) ; q ) R ε × { H 1 / 2 ( Γ i n ) } N ( u , v ) Ω o u t + ( × w , × 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 .
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:
Given   a   solenoidal   g { L 0 2 ( Ω ) } N ,   find   ( ( u ; w ) ; r ; p ) R × L 2 ( Ω i n ) × { H 1 / 2 ( Γ i n ) } N such   that   ( ( v ; z ) ) ; s ; q ) R × L 2 ( Ω i n ) × { H 1 / 2 ( Γ i n ) } N ( u , v ) Ω o u t + ( × w , × z ) Ω i n + ( r , · ( ε z ) ) Ω i n + < p , v z > 1 / 2 , Γ i n = ( g | Ω o u t , v ) Ω o u t + ( g | Ω i n , z ) Ω i n , ( · ( ε w ) , s ) Ω i n = 0 , < q , u w > 1 / 2 , Γ i n = 0 . .
It is clear that the solution of (2.5) solves the following problem:
Given   g { L 0 2 ( Ω ) } N   fulfilling   · g = 0   in   Ω , find   ( u ; w ) R   such   that Δ u = g | Ω o u t in   Ω o u t     × × w = g | Ω i n , · ( ε w ) = 0 . in   Ω i n , u = w on   Γ i n , n u = 0 on   Γ , Ω o u t u + Ω i n w = 0 .
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 · w = log ( ε ) · w L 2 ( Ω i n ) from the assumed regularity of ε , we have · w { H 1 ( Ω i n ) } N . This implies that Δ w { H 1 ( Ω i n ) } N by the identity (2.3), since × × w { L 2 ( Ω i n } N . Therefore w { H 1 ( Ω i n ) } N owing to the coincidence of u and w on Γ 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 e ˜ | t = 0 = e 0 and { t e ˜ } | t = 0 = e 1 , we wish to,
Find   e V ˜ ,   where V ˜ : = { v | v H 2 [ ( 0 , T ) ; { L 2 ( Ω ) } N ] H 1 [ ( 0 , T ) ; V ] } such   that   v { H 1 ( Ω ) } N   and   t ( 0 , T ) t t e ˜ , v ε + ( e ˜ , v ) + ( · ε e ˜ , · v ) ( · e ˜ , · v ) + ( t e ˜ , v ) Γ = ( f , v ) .
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 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,
ε t t e Δ e + · e = f   in   Ω i n × ( 0 , T ) .
We know that the solution of (2.2) satisfies · ( ε e ) = 0 in Ω i n × ( 0 , T ) . If we subtract this equation from (2.8), we obtain:
ε t t e Δ e · ( ε e ) + · e = f   in   Ω i n × ( 0 , T ) .
Since ε 1 in Ω o u t it is readily seen that (2.9) also holds in the whole Ω × ( 0 , T ) .
  • Thus, taking an arbitrary v V and using Green’s first identity, together with the absorbing boundary conditions satisfied by e , we readily obtain v V ,
ε t t e , v + ( e , v ) ( · e , · v ) + ( · { ε e } , · v ) + ( t e , v ) Γ = ( f , v ) t ( 0 , T ) ,
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):
t ( 0 , T ] t e ˜ ε 2 + e ˜ 2 + · e ˜ ε 1 2 ( t ) E 0 T f ( · , t ) 2 d t + e 0 2 + e 0 2 + · e 0 ε 1 2 + e 1 ε 2 ,
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 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 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 V h such that, t ( 0 , T ]
t t e h , v ε + ( e h , v ) + ( · { ε e h } , · v ) ( · e h , · v ) + ( t e h , v ) Γ = 0 v V h , with e h ( · , 0 ) = e 0 h ( · )   and   t e h ( · , 0 ) = e 1 h ( · )   in   Ω .

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 τ ) by e h k for k = 1 , 2 , , M according to the following FE scheme for k = 1 , 2 , , M 1 :
e h k + 1 2 e h k + e h k 1 τ 2 , v ε + ( e h k , v ) + ( · ε e h k , · v ) ( · e h k , · v ) + e h k + 1 e h k 1 2 τ , v Γ = 0 v V h , with e h 0 = e 0 h   and   e h 1 = e h 0 + τ e 1 h   in   Ω .
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 K u | K · v | K d x , for every element K in T h , for ( u ; v ) { C 0 ( Ω ¯ ) } N × { C 0 ( Ω ¯ ) } 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 ) { C 0 ( Ω ¯ ) } N × { C 0 ( Ω ¯ ) } N .
( u , v ) ε h , h : = K T h ε ( G K ) m e a s ( K ) i = 1 N + 1 u ( S K , i ) · v ( S K , i ) N + 1 ,
where S K , i are the vertexes of K, i = 1 , , 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 Γ K . Then we set
( u , v ) Γ , h : = K S h F K Γ m e a s ( F K ) i = 1 N u ( R F K , i ) · v ( R F K , i ) N ,
where R F K , i are the vertexes of F K , i = 1 , , N .
For coherence ε h is defined to be the function whose value in each K T h is constant equal to ε ( G K ) . Furthermore we introduce the norms { · } ε h , h and { · } h of V h , given by ( { · } , { · } ) ε h , h 1 / 2 and ( { · } , { · } ) h 1 / 2 , respectively. Similarly we denote by { · } Γ , h the norm defined by ( { · } , { · } ) Γ , h 1 , 2 . Then still denoting the approximation of e h ( k τ ) by e h k , for k = 1 , 2 , . . . , M we determine e h k + 1 by,
e h k + 1 2 e h k + e h k 1 τ 2 , v ε h , h + ( e h k , v ) + ( · ε e h k , · v ) ( · e h k , · v ) + e h k + 1 e h k 1 2 τ , v Γ , h = 0 v V h , with e h 0 = e 0 h   and   e h 1 = e h 0 + τ e 1 h   in   Ω .
Now we adapt a result given in Lemma 3 of [12], which allows us to assert that the following upper bounds hold:
v ε h v ε h , h v V h .
Similarly we have,
v Γ v Γ , h v V h .

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,
e h k + 1 2 e h k + e h k 1 τ 2 , v ε h , h + ( e h k , v ) + ( · ε e h k , · v ) ( · e h k , · v ) + e h k + 1 e h k 1 2 τ , v Γ , h = F k ( v ) + ( d k , · v ) + G k ( v ) v V h , e h 0 = e 0 h   and   e h 1 = e h 0 + τ e 1 h   in   Ω .
where for every k { 1 , 2 , , 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 · h and · Γ , 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 { 1 , 2 , , M 1 } .
Taking v = e h k + 1 e h k 1 in (4.17) we obtain for k = 1 , 2 , , M 1 ,
e h k + 1 2 e h k + e h k 1 τ , e h k + 1 e h k 1 τ ε h , h + ( e h k , e h k + 1 e h k 1 ) + ( · { ε 1 } e h k , · e h k + 1 · e h k 1 ) + 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 , · { e h k + 1 e h k 1 } ) + G k ( e h k + 1 e h k 1 )
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):
e h k + 1 e h k τ ε h , h 2 e h k e h k 1 τ ε h , h 2 + I 1 + I 2 + 2 τ e h k + 1 e h k 1 2 τ Γ , h 2 | F k | h e h k + 1 e h k 1 h + d k · ( e h k + 1 e h k 1 ) + | G k | Γ , h e h k + 1 e h k 1 Γ , h where I 1 : = ( e h k , { e h k + 1 e h k 1 } ) ; and I 2 : = ( · { ε 1 } e h k , · { e h k + 1 e h k 1 } ) .
Next we estimate the terms I 1 and I 2 given by (4.19).
  • First of all it is easy to see that
I 1 = 1 2 ( e h k + 1 2 + e h k 2 ( e h k + 1 e h k ) 2 ) 1 2 ( e h k 2 + e h k 1 2 ( e h k e h k 1 ) 2 ) .
Next we note that,
I 2 = J 1 + J 2 where J 1 : = ( · e h k , · { e h k + 1 e h k 1 } ) ε 1 and J 2 : = ( ε · e h k , · { e h k + 1 e h k 1 } ) .
Similarly to (4.20) we can write,
J 1 = 1 2 ( · e h k + 1 ε 1 2 + · e h k ε 1 2 · ( e h k + 1 e h k ) ε 1 2 ) 1 2 ( · e h k ε 1 2 + · e h k 1 ε 1 2 · ( e h k e h k 1 ) ε 1 2 )
Now observing that ε 0 on Γ , we integrate by parts J 2 given by (4.21), to get
J 2 = ( { ε · e h k } , e h k + 1 e h k 1 ) .
Let us rewrite J 2 as,
J 2 = M 1 + M 2 where M 1 : = ε e h k , e h k + 1 e h k 1 and M 2 : = { e h k } T ε , e h k + 1 e h k 1
M 1 in turn can be rewritten as follows:
M 1 = N 1 + N 2 where N 1 : = τ ε e h k , e h k + 1 e h k τ and N 2 : = τ ε e h k , e h k e h k 1 τ .
Then we further observe that
N 1 = τ 2 i = 1 k ε e h i e h i 1 τ , e h k + 1 e h k τ τ ε e h 0 , e h k + 1 e h k τ .
and hence,
N 1 τ 2 | ε | 2 , e h k + 1 e h k τ i = 1 k e h i e h i 1 τ + τ | ε | 2 , e h 0 e h k + 1 e h k τ
or yet,
N 1 τ | ε | 2 , e h k + 1 e h k τ τ k i = 1 k e h i e h i 1 τ 2 1 / 2 + e h 0 ,
and noting that k T / τ we get
N 1 τ | ε | 2 , e h k + 1 e h k τ τ T τ i = 1 k e h i e h i 1 τ 2 1 / 2 + e h 0 .
Applying to (4.27) Young’s inequality a b δ a 2 / 2 + b 2 / ( 2 δ ) a , b and δ > 0 with δ = 1 , we easily conclude that
N 1 τ 2 | ε | 2 , 2 e h k + 1 e h k τ 2 + τ T i = 1 k e h i e h i 1 τ 2 + e h 0 2 .
Similarly to (4.28),
N 2 τ 2 | ε | 2 , 2 e h k e h k 1 τ 2 + τ T i = 1 k e h i e h i 1 τ 2 + e h 0 2 .
Combining (4.28) and (4.29) we come up with
M 1 τ | ε | 2 , e h k + 1 e h k τ 2 + e h k e h k 1 τ 2 + L k , where L k : = τ T i = 1 k e h i e h i 1 τ 2 + e h 0 2 .
As for M 2 given by (4.24) we have:
M 2 ε 1 , e h k e h k + 1 e h k 1 τ | ε | 1 , e h k e h k + 1 e h k τ + e h k e h k 1 τ
or yet
M 2 τ 2 | ε | 1 , 2 e h k 2 + e h k + 1 e h k τ 2 + e h k e h k 1 τ 2 .
Now we recall (4.19) together with (3.16) and notice that for every square-integrable field A in Ω we have A A ε h , Then taking into account that · v N v v { H 1 ( Ω ) } N , and using Young’s inequality with δ = τ , δ = 1 / τ and δ = τ , respectively, we easily obtain the following estimates:
| F k | h e h k + 1 e h k 1 τ 2 | F k | h 2 + τ e h k + 1 e h k τ h 2 + e h k e h k 1 τ h 2 , d k · ( e h k + 1 e h k 1 ) N 2 τ d k 2 + τ ( e h k + 1 2 + e h k 1 ) 2 ) , | G k | Γ , h e h k + 1 e h k 1 Γ , h τ 2 | G k | Γ , h 2 + 2 τ e h k + 1 e h k 1 2 τ Γ , h 2 .
where in the first and the second inequality we also used the fact that A ± B 2 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 · · h · ε h · ε h , h we obtain for 1 k M 1 :
( A k B k ) + ( C k D k ) τ 2 | F k | h 2 + N 2 τ d k 2 + τ 2 G k Γ h 2 , where A k : = e h k + 1 e h k τ ε h , h 2 τ 2 2 + | ε | 1 , + 2 | ε | 2 , e h k + 1 e h k τ ε h , h 2 + 1 2 ( 1 2 τ ) e h k + 1 2 + ( 1 τ | ε | 1 , ) e h k 2 B k : = e h k e h k 1 τ ε h , h 2 + τ 2 2 + | ε | 1 , + 2 | ε | 2 , e h k e h k 1 τ ε h , h 2 + 1 2 1 + 2 τ e h k 1 2 + 1 + τ | ε | 1 , e h k 2 + τ | ε | 2 , L k C k : = 1 2 ( e h k + 1 e h k ) 2 + · ( e h k + 1 e h k ) ε 1 2 + 1 2 · e h k + 1 ε 1 2 + · e h k ε 1 2 D k : = 1 2 ( e h k e h k 1 ) 2 + · ( e h k e h k 1 ) ε 1 2 + 1 2 · e h k ε 1 2 + · e h k 1 ε 1 2 .
Setting
η : = 2 + | ε | 1 , + 2 | ε | 2 , ; θ : = | ε | 1 , ; ρ : = T 2 | ε | 2 , ,
we can rewrite A k and B k as follows:
A k : = e h k + 1 e h k τ ε h , h 2 τ 2 η e h k + 1 e h k τ ε h , h 2 + 1 2 e h k + 1 2 + e h k 2 τ e h k + 1 2 + θ 2 e h k 2 ; B k : = e h k e h k 1 τ ε h , h 2 + τ 2 η e h k e h k 1 τ ε h , h 2 + 1 2 e h k 1 2 + e h k 2 + τ e h k 1 2 + θ 2 e h k 2 + τ ρ T 2 L k
Then we note that for 1 k m 1 with 2 m M 1 ,
A k 1 B k = τ η e h k e h k 1 τ ε h , h 2 τ 1 + θ 2 e h k 2 + e h k 1 2 τ ρ T 2 L k .
It follows that
k = 1 m ( A k B k ) = A m B 1 + k = 2 m ( A k 1 B k ) = 1 2 2 τ η e h m + 1 e h m τ ε h , h 2 + 1 2 1 2 τ e h m + 1 2 + 1 2 1 τ θ e h m 2 1 2 2 + τ η e h 1 e h 0 τ ε h , h 2 1 2 1 + 2 τ e h 0 2 1 2 1 + τ θ e h 1 2 τ ρ T 2 L 1 k = 2 m τ η e h k e h k 1 τ ε h , h 2 + τ 2 2 + θ e h k 2 + e h k 1 2 + τ ρ T 2 L k .
Now we extend to k = 1 the summation range on the right hand side of (4.37) to obtain,
k = 1 m ( A k B k ) = 1 2 2 τ η e h m + 1 e h m τ ε h , h 2 + 1 2 1 2 τ e h m + 1 2 + 1 2 1 τ θ e h m 2 1 2 2 τ η e h 1 e h 0 τ ε h , h 2 1 2 1 τ θ e h 0 2 1 2 1 2 τ e h 1 2 k = 1 m τ η e h k e h k 1 τ ε h , h 2 + τ 2 2 + θ e h k 2 + e h k 1 2 + τ ρ T 2 L k .
Moreover since C k 1 = D k for all m k 2 , recalling (4.33) we easily come up with,
k = 1 m ( C k D k ) = C m D 1 = 1 2 ( e h m + 1 e h m ) 2 + · ( e h m + 1 e h m ) ε 1 2 + 1 2 · e h m + 1 ε 1 2 + · e h m ε 1 2 + 1 2 ( e h 1 e h 0 ) 2 + · ( e h 1 e h 0 ) ε 1 2 1 2 · e h 1 ε 1 2 + · e h 0 ε 1 2 .
Combining (4.38) and (4.39), we obtain for 2 m 1 :
k = 1 m ( A k + C k ) k = 1 m ( B k + D k ) = 1 2 2 τ η e h m + 1 e h m τ ε h , h 2 + 1 2 1 2 τ e h m + 1 2 + 1 2 1 τ θ e h m 2 1 2 2 τ η e h 1 e h 0 τ ε h , h 2 1 2 1 τ θ e h 0 2 1 2 1 2 τ e h 1 2 k = 1 m τ η e h k e h k 1 τ ε h , h 2 + τ 2 2 + θ e h k 2 + e h k 1 2 + τ ρ T 2 L k 1 2 ( e h m + 1 e h m ) 2 + · ( e h m + 1 e h m ) ε 1 2 + 1 2 · e h m + 1 ε 1 2 + · e h m ε 1 2 + 1 2 ( e h 1 e h 0 ) 2 + · ( e h 1 e h 0 ) ε 1 2 1 2 · e h 1 ε 1 2 + · e h 0 ε 1 2 .
On the other hand, recalling (4.30) we note that for 2 m M 1
k = 1 m L k k = 1 m τ T i = 1 m e h i e h i 1 τ 2 + e h 0 2 = T m τ k = 1 m e h k e h k 1 τ 2 + m e h 0 2
In short since m τ T , from (4.41) we easily obtain for 2 m M 1 :
τ ρ T 2 k = 1 m L k τ ρ k = 1 m e h k e h k 1 τ 2 + ρ T e h 0 2 .
Plugging (4.42) into (4.40) and summing up both sides of (4.33) from k = 1 through k = m for 2 m N 1 by using (4.35), (4.40) yields:
1 2 2 τ η e h m + 1 e h m τ ε h , h 2 + 1 2 1 2 τ e h m + 1 2 + 1 2 1 τ θ e h m 2 1 2 2 τ η e h 1 e h 0 τ ε h , h 2 1 2 1 τ θ e h 0 2 1 2 1 2 τ e h 1 2 k = 1 m τ η + ρ e h k e h k 1 τ ε h , h 2 + τ 2 2 + θ e h k 2 + e h k 1 2 1 2 ( e h m + 1 e h m ) 2 + · ( e h m + 1 e h m ) ε 1 2 + 1 2 · e h m + 1 ε 1 2 + · e h m ε 1 2 + 1 2 ( e h 1 e h 0 ) 2 + · ( e h 1 e h 0 ) ε 1 2 1 2 · e h 1 ε 1 2 + · e h 0 ε 1 2 k = 1 m τ 2 | F k | h 2 + N 2 τ d k 2 + τ 2 G k Γ h 2 + ρ T e h 0 2 .
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 e h j 2 for j = 0 , 1 and e 1 h ε h , h 2 , for 2 m M 1 we come up with:
1 2 2 τ η e h m + 1 e h m τ ε h , h 2 + 1 2 1 2 τ e h m + 1 2 + 1 2 1 τ θ e h m 2 1 2 ( e h m + 1 e h m ) 2 + · ( e h m + 1 e h m ) ε 1 2 + 1 2 · e h m + 1 ε 1 2 + · e h m ε 1 2 k = 1 m τ η + ρ e h k e h k 1 τ ε h , h 2 + τ 2 2 + θ e h k 2 + e h k 1 2 + τ 2 | F k | h 2 + N 2 τ d k 2 + τ 2 G k Γ , h 2 + e 1 h ε h , h 2 + τ 2 2 e 1 h 2 + · e 1 h ε 1 2 + 1 2 e h 1 2 + e h 0 2 + 1 2 · e h 1 ε 1 2 + · e h 0 ε 1 2 + ρ T e h 0 2 .
Now we recall a classical inverse inequality (cf. [14]), according to which,
v C h 1 v C h 1 v ε h , h   for   all   v V h ,
where C is a mesh-independent constant, and we apply the trivial upper bound · v ε 1 N ε 1 v for all v V h .
  • Let us assume that τ satisfies the following CFL-condition:
τ h / ν   with   ν = C ( 1 + N ε 1 ) 1 / 2 .
Then we have, v V h :
v 2 + · v ε 1 2 ν 2 τ 2 h 2 v 2 τ 2 v ε h , h 2 τ 2
This means that
( e h m + 1 e h m ) 2 + · ( e h m + 1 e h m ) ε 1 2 e h m + 1 e h m τ ε h , h 2
Plugging (4.48) into (4.44) we obtain,
1 2 1 τ η e h m + 1 e h m τ ε h , h 2 + 1 2 1 2 τ e h m + 1 2 + 1 2 1 τ θ e h m 2 k = 1 m τ η + ρ e h k e h k 1 τ ε h , h 2 + τ 2 2 + θ e h k 2 + e h k 1 2 + τ 2 | F k | h 2 + N 2 τ d k 2 + τ 2 G k Γ , h 2 + e 1 h ε h , h 2 + τ 2 2 e 1 h 2 + · e 1 h ε 1 2 + 1 2 e h 1 2 + e h 0 2 + 1 2 · e h 1 ε 1 2 + · e h 0 ε 1 2 + ρ T e h 0 2 .
Next we note that both 1 2 τ and 1 τ θ are bounded below by 1 τ η and moreover, 1 2 ( 2 + θ ) is obviously bounded above by η + ρ . Therefore it is easy to see that (4.49) can be transformed into :
1 2 1 τ η e h m + 1 e h m τ ε h , h 2 + e h m + 1 2 + e h m 2 S M + E 0 + k = 1 m τ η + ρ e h k e h k 1 τ ε h , h 2 + e h k 2 + e h k 1 2 ,
where
S M : = k = 1 M 1 τ 2 | F k | h 2 + N 2 τ d k 2 + τ 2 G k Γ , h 2 E 0 : = e 1 h ε h , h 2 + τ 2 2 e 1 h 2 + · e 1 h ε 1 2 + 1 2 e h 1 2 + e h 0 2 + 1 2 · e h 1 ε 1 2 + · e h 0 ε 1 2 + ρ T e h 0 2 .
Let us assume that τ 1 2 η . Then from (4.51) we have
e h m + 1 e h m τ ε h , h 2 + e h m + 1 2 + e h m 2 4 S M + E 0 + k = 1 m τ η + ρ e h k e h k 1 τ ε h , h 2 + e h k 2 + e h k 1 2 .
Now we set
β = 4 M τ ( η + ρ ) = 4 T 2 + | ε | 1 , + ( 2 + T 2 ) | ε | 2 , .
From the discrete Grönwall’s Lemma and (3.15) from (4.52) we infer that m M 1 :
e h m + 1 e h m τ 2 + e h m + 1 2 + e h m 2 4 ( S M + E 0 ) e β ,
as long as τ 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 { L 0 2 ( Ω ) } N the solution v g { H 1 ( Ω ) L 0 2 ( Ω ) } N of the equation
Δ v g { ( ε 1 ) · v g } = g in   Ω n v g = 0 on   Γ .
belongs to { H 2 ( Ω ) } N .
  • Another result that we take for granted in this section is the existence of a constant C ε such that,
g { L 0 2 ( Ω ) } N   it   holds   H ( v g ) C ε g ,
where H ( · ) 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 ε 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 e ^ h ( · , t ) belonging to V h for every t [ 0 , T ] , uniquely defined up to an additive field depending only on t as follows:
( e ^ h { · , t } , v ) + ( · e ^ h { · , t } , · v ) ε 1 = ( e { · , t } , v ) + ( · e { · , t } , · v ) ε 1 v V h .
The time-dependent additive field up to which e ^ h { · , t } is defined can be determined by requiring that Ω { e ^ h ( · , t ) e ( · , t ) } d x = 0 t [ 0 , T ] .
Let us further assume that for every t [ 0 , T ] e ( · , t ) { H 2 ( Ω ) } N . In this case, from classical approximation results based on the interpolation error, we can assert that,
{ e ^ h ( · , t ) e ( · , t ) } + · { e ^ h ( · , t ) e ( · , t ) } ε 1 C ^ 1 h H ( e ) ( · , t ) t [ 0 , T ] ,
where C ^ 1 is a mesh-independent constant.
  • Let us show that there exists another mesh-independent constant C ^ 0 such that for every t [ 0 , T ] it holds,
e ^ h ( · , t ) e ( · , t ) C ^ 0 h 2 H ( e ) ( · , t ) t [ 0 , T ] .
Since e ^ h ( · , t ) e ( · , t ) { L 0 2 ( Ω ) } N for every t, we may write:
e ^ h ( · , t ) e ( · , t ) = sup g { L 0 2 ( Ω ) } N { 0 } ( e ^ h { · , t } e { · , t } , g ) g t [ 0 , T ] .
Defining W : = { w | w { H 2 ( Ω ) L 0 2 ( Ω ) } N n w = 0 o n Γ } , owing to (4.56) we have t [ 0 , T ] :
e ^ h ( · , t ) e ( · , t ) C ε sup w W { 0 } ( e ^ h { · , t } e { · , t } , Δ w ( ε 1 ) · w ) H ( w ) .
Since n w = 0 and ε = 1 on Γ we may integrate by parts the numerator in (4.61) to obtain for every t [ 0 , T ] ,
e ^ h ( · , t ) e ( · , t ) C ε sup w W { 0 } ( { e ^ h { · , t } e { · , t } } , w ) + ( · { e ^ h ( · , t ) e ( · , t ) } , · w ) ε 1 H ( w ) .
Taking into account (4.57), the numerator of (4.62) can be rewritten as
( { e ^ h { · , t } e { · , t } } , { w v } ) + ( · { e ^ h { · , t } e { · , t } } , · { w v } ) ε 1 v V h .
Then choosing v to be the V h -interpolate of w , taking into account (4.58) we easily establish (4.59) with C ^ 0 = C ε C ^ 1 2 .
  • To conclude these preliminary considerations, we refer to Chapter 5 of [31], to infer that the second order time-derivative t t e ^ h ( · , t ) is well defined in { H 1 ( Ω ) } N for every t [ 0 , T ] , as long as t t e ( · , t ) lies in { H 1 ( Ω ) } N for every t [ 0 , T ] . Moreover, provided t t e { H 2 ( Ω ) } N for every t [ 0 , T ] , the following estimate holds:
t t { e ^ h ( · , t ) e ( · , t ) } C ^ 0 h 2 H ( t t e ) ( · , t ) t [ 0 , T ] .
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 ( 0 , T ) C T { u L 2 ( 0 , T ) 2 + u t L 2 ( 0 , T ) 2 } 1 / 2 u H 1 ( 0 , T ) .
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 ( Ω × ( 0 , T ) ) } N .
  • Now taking u = H ( e ) ( · , t ) we have u t ( t ) = Ω { H ( e ) : H ( t e ) } ( x , t ) d x { H ( e ) } ( · , 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 ) } ( · , t ) C T { H ( e ) L 2 ( 0 , T ) 2 + H ( t e ) L 2 ( 0 , T ) 2 } 1 / 2 t ( 0 , T ) .
In complement to the above ingredients we extend the inner products ( · , · ) ε h , h and ( · , · ) Γ , h , and associated norms · ε h , h and · Γ , h in a semi-definite manner, to fields in { L 2 ( Ω ) } N , N = 2 , 3 as follows:
  • First of all, K T h let Π K : L 2 ( K ) P 1 ( K ) be the standard orthogonal projection operator onto the space P 1 ( K ) of linear functions in K. We set
u   and   v { L 2 ( Ω ) } N , ( u , v ) ε h , h : = K T h ε ( G K ) m e a s ( K ) N + 1 i = 1 N + 1 Π K u ( S K , i ) · Π K v ( S K , i ) .
Let us generically denote by F Γ an edge of a triangle or a face of a tetrahedron K such that m e a s ( K Γ ) > 0 . Moreover, we denote by Π F ( v ) the standard orthogonal projection of a function v L 2 ( F ) onto the space of linear functions on F. Similarly we define:
u   and   v { L 2 ( Γ ) } N , ( u , v ) Γ , h : = F Γ m e a s ( F ) N i = 1 N Π F u ( R F , i ) · Π 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 { L 2 ( Ω ) } N There exists a mesh independent constant c Ω such that u { H 1 ( Ω ) } N and v V h ,
| Δ ε h ( u , v ) | c Ω ε 0 , h u v h   .
Lemma 4.2.
Let Ξ h ( γ { u } , γ { v } ) : = ( γ { u } , γ { v } ) Γ , h ( γ { u } , γ { v } ) Γ for u , v { H 1 ( Ω ) } N , where γ { w } represents the trace on Γ of a function w H 1 ( Ω ) . Let also Γ be the tangential gradient operator over Γ. There exists a mesh independent constant c ˜ Γ such that u { H 2 ( Ω ) } N and v V h ,
| Ξ h ( γ { u } , γ { v } ) | c ˜ Γ h Γ γ { u } Γ γ { v } Γ , h .   .
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 Γ γ { u } is well defined in { L 2 ( Γ ) } N if u { 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 } ) | c Γ h { u 2 + H ( u ) 2 } 1 / 2 γ { v } Γ , h .   .
To conclude we prove the validity of the following upper bounds:
Lemma 4.3.
v { L 2 ( Ω ) } N it holds
v ε h , h N + 2 v ε h .
Proof: Denoting by Π h ( v ) the function defined in Ω whose restriction to every K T h is Π K ( v ) for a given v { L 2 ( Ω ) } N , from an elementary property of the orthogonal projection we have
Π h ( v ) ε h v ε h , v { L 2 ( Ω ) } N .
Now taking v such that v | K P 1 ( K ) K T h , by a straightforward calculation using the expression of v | K in terms of barycentric coordinates we have:
K ε ( G K ) v | K 2 d x = 2 ε ( G K ) m e a s ( K ) ( N + 2 ) ( N + 1 ) i = 1 N + 1 | v ( S K , i ) | 2 + i = 1 N j = i + 1 N + 1 v ( S K , i ) v ( S K , j ) .
It trivially follows that
K ε ( G K ) v | K 2 d x = ε ( G K ) m e a s ( K ) ( N + 2 ) ( N + 1 ) i = 1 N + 1 | v ( S K , i ) | 2 + i = 1 N + 1 v ( S K , i ) 2 .
and finally
( N + 2 ) K ε ( G K ) v | K 2 d x ε ( G K ) m e a s ( K ) N + 1 i = 1 N + 1 | v ( S K , i ) | 2 .
This immediately yields Lemma 4.3, taking into account (4.69). ■
Lemma 4.4.
v { L 2 ( Γ ) } N it holds v Γ , h 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 , , N functions e ^ h k V h by e ^ h k ( · ) : = e ^ h ( · , k τ ) . In the sequel for any function or field A defined in Ω × ( 0 , T ) , A k ( · ) denotes A ( · , k τ ) , except for other quantities carrying the subscript h such as e h k .
Let us substitute e h k by e ^ h k for k = 2 , 3 , , N on the left hand side of the first equation of (3.14) and take also as initial conditions 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:
E h k ( v ) = ( { t t e } k , v ) ε + ( e k , v ) + ( · e k , · v ) ε 1 + ( ε · e k , · v ) + ( { t e } k , v ) Γ F ¯ k ( v ) ( d ¯ k , · v ) G ¯ k ( v ) v V h ,
where
d ¯ k = ε · { e ^ h k e k } ,
and F ¯ k ( v ) and G ¯ k ( v ) can be written as follows:
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 ( 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 ) { t t e } k , v ) ε ,
T k being the finite-difference operator defined by,
T k ( · ) : = { · } k + 1 2 { · } k + { · } k 1 τ 2 ,
and
G ¯ k ( v ) = G 1 k ( v ) + G 2 k ( v ) + G 3 k ( v ) , with G 1 k ( v ) = ( Q k ( e ^ h e ) , v ) Γ , h ; G 2 k ( v ) = Γ h ( Q k ( e ) , v ) ; G 3 k ( v ) = ( Q k ( e ) { t e } k , v ) Γ ,
Q k being the finite-difference operator defined by,
Q k ( · ) : = { · } k + 1 { · } k 1 2 τ .
Notice that, under Assumption * both t e ( · , t ) and t t e ( · , t ) belong to { H 1 ( Ω ) } N for every t [ 0 , T ] . Hence we can define t e ^ h from t e and t t e ^ h from t t e in the same way as e ^ h is defined from e . Moreover, straightforward calculations lead to,
T k ( · ) = 1 τ 2 ( k 1 ) τ k τ t k τ t t { · } d s d t + k τ ( k + 1 ) τ k τ t t t { · } d s d t .
Furthermore another straightforward calculation allows us writing:
T k ( e ) = { t t e } k + R k ( e ) where R k ( e ) : = 1 2 τ 2 ( k 1 ) τ k τ { t ( k 1 ) τ } 2 t t t e d t + k τ ( k + 1 ) τ { ( k + 1 ) τ t } 2 t t t e d t .
Similarly,
Q k ( · ) = 1 2 τ ( k 1 ) τ ( k + 1 ) τ t { · } d t ,
and
Q k ( e ) = { t e } k + S k ( e ) , where S k ( e ) = 1 2 τ ( k 1 ) τ k τ { t ( k 1 ) τ } t t e d t + k τ ( k + 1 ) τ { ( k + 1 ) τ t } t t e d t .
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 e ^ h k V h are the solution of the following problem, for k = 1 , 2 , , M 1 :
ε e ^ h k + 1 2 e ^ h k + e ^ h k 1 τ 2 , v + ( e ^ h k , v ) + ( · { ε e ^ h k } , · v ) ( · e ^ h k , · v ) + e ^ h k + 1 e ^ h k 1 2 τ , v Γ = F ¯ k ( v ) + ( d ¯ k , · v ) + G ¯ k ( v ) v V h , e ^ h 0 ( · ) = e ^ h ( · , 0 )   and   e ^ h 1 ( · ) = e ^ h ( · , τ )   in   Ω ,
d ¯ k , F ¯ k and G ¯ k being given by (4.71), (4.72)-(4.73) and (4.74)-(4.75).
  • Estimating d ¯ k is a trivial matter. Indeed, since e k { H 2 ( Ω ) } N , from (4.59) we immediately obtain,
d ¯ k C ^ 0 h 2 | ε | 1 , H ( e k ) .
Therefore consistency of the scheme will result from suitable estimations of | F ¯ k | h and | 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 ( · ) , Q k ( · ) , R k ( · ) and S k ( · ) . With this aim we denote by | · | the euclidean norm of n , for n 1 .
  • From (4.76) and the Cauchy-Schwarz inequality, we easily obtain for every x Ω and u such that { t t u } ( · , t ) { H 2 ( Ω ) } N t ( 0 , T ) ,
| T h k { u ( x ) } | 1 τ 2 ( k 1 ) τ k τ t k τ { t t u } ( x ) d s d t + k τ ( k + 1 ) τ k τ t { t t u } ( x ) d s d t 1 τ 2 ( k 1 ) τ k τ ( k 1 ) τ k τ { t t u } ( x ) d s d t + k τ ( k + 1 ) τ k τ ( k + 1 ) τ { t t u } ( x ) d s d t = 1 τ ( k 1 ) τ ( k + 1 ) τ { t t u } ( x ) d t .
It follows that, for every u such that { t t u } ( · , t ) { H 2 ( Ω ) } N t ( 0 , T ) we have,
| T h k { u ( x ) } | 2 τ ( k 1 ) τ ( k + 1 ) τ { t t u } ( x ) 2 d t 1 / 2 .
Furthermore, from (4.77) and the inequality a + b { 2 ( a 2 + b 2 ) } 1 / 2 a , b , for every x Ω we obtain:
| R k { e ( x ) } | 1 2 ( k 1 ) τ k τ { t t t e } ( x ) d t + k τ ( k + 1 ) τ { t t t e } ( x ) d t τ 2 ( k 1 ) τ k τ { t t t e } ( x ) 2 d t 1 / 2 + k τ ( k + 1 ) τ { t t t e } ( x ) 2 d t 1 / 2 .
It follows that,
| R k { e ( x ) } | τ 2 ( k 1 ) τ ( k + 1 ) τ { t t t e } ( x ) 2 d t 1 / 2 .
On the other hand from (4.78) and the Cauchy-Schwarz inequality we trivially have for every x Γ and u such that γ { t u } ( · , t ) { H 3 / 2 ( Γ ) } N t ( 0 , T ) :
| Q k { u ( x ) } | 1 2 τ ( k 1 ) τ ( k + 1 ) τ | { t u } ( x ) | 2 d t 1 / 2 .
Finally, similarly to (4.83), from (4.79) for every x in Γ we successively derive,
| S k { e ( x ) } | 1 2 ( k 1 ) τ k τ { t t e } ( x ) d t + k τ ( k + 1 ) τ { t t e } ( x ) d t τ 2 ( k 1 ) τ k τ { t t e } ( x ) 2 d t 1 / 2 + k τ ( k + 1 ) τ { t t e } ( x ) 2 d t 1 / 2 .
Therefore it holds,
| S k { e ( x ) } | τ 2 ( k 1 ) τ ( k + 1 ) τ { t t e } ( x ) 2 d t 1 / 2 .
Notice that bounds entirely analogous to (4.82) and (4.84) hold for T k ( · ) = T k ( · ) and Γ Q k ( · ) = Q k ( Γ · ) , that is, x Ω ,
| T h k ( e ) ( x ) | 2 τ ( k 1 ) τ ( k + 1 ) τ { t t e } ( x ) 2 d t 1 / 2 ,
and x Γ ,
| Γ Q k ( e ) ( x ) | 1 2 τ ( k 1 ) τ ( k + 1 ) τ | { t Γ e } ( x ) | 2 d t 1 / 2 .
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:
| F 1 k ( v ) | ε 0 , 2 τ ( k 1 ) τ ( k + 1 ) τ { t t ( e ^ h e ) } ( · , t ) 2 d t 1 / 2 h v h N + 2 ε 0 , 2 τ ( k 1 ) τ ( k + 1 ) τ { t t ( e ^ h e ) } ( · , t ) 2 d t 1 / 2 v h N + 2 ε 0 , 2 τ ( k 1 ) τ ( k + 1 ) τ { t t ( e ^ h e ) } ( · , t ) 2 d t 1 / 2 v h
Recalling (4.72) and applying (4.63) to (4.88), we come up with,
| F 1 k ( v ) | C ^ 0 h 2 2 ( N + 2 ) τ ε 0 , ( k 1 ) τ ( k + 1 ) τ { H ( t t e ) } ( · , t ) 2 d t 1 / 2 v h .
Next, combining (4.66) and (4.86) we immediately obtain.
| F 2 k ( v ) | c Ω h ε 0 , 2 τ ( k 1 ) τ ( k + 1 ) τ { t t e } ( · , t ) 2 d t 1 / 2 v h .
Further, taking into account (3.15), (4.76), (4.86) and the standard estimate ε h ε 0 , C h | ε | 1 , where C is a mesh-independent constant, it holds
| F 3 k ( v ) | C h | ε | 1 , 2 τ ( k 1 ) τ ( k + 1 ) τ { t t e } ( · , t ) 2 d t 1 / 2 v h
Finally by (3.15), (4.77) and (4.83), we have
| F 4 k ( v ) | ε 0 , τ 2 ( k 1 ) τ ( k + 1 ) τ { t t t e } ( · , t ) 2 d t 1 / 2 v h .
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 C ^ 1 :
{ t ( e ^ h e ) } ( · , t ) C ^ 1 h { H ( t e ) } ( · , t ) ) , t [ 0 , T ] .
On the other hand, by the Trace Theorem there exists a contant C T r depending only on Ω such that,
{ t ( e ^ h e ) } ( · , t ) Γ C T r { { t ( e ^ h e ) } ( · , t ) 2 + { t ( e ^ h e ) } ( · , t ) 2 } 1 / 2 t [ 0 , T ] .
Hence by (4.63), (4.93) and (4.94), we have, for a suitable mesh-independent constant C B :
{ t ( e ^ h e ) } ( · , t ) Γ C B h { H ( t e ) } ( · , t ) t [ 0 , T ] .
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 ) | 2 2 τ ( k 1 ) τ ( k + 1 ) τ { t ( e ^ h e ) } ( · , t ) Γ 2 d t 1 / 2 v Γ , h .
Then using (4.95) we immediately establish,
| G 1 k ( v ) | C B h 2 τ ( k 1 ) τ ( k + 1 ) τ { H ( t e ) } ( · , t ) 2 d t 1 / 2 v Γ , h .
Next we switch to G 2 k . Similarly to (4.90), (4.68) and (4.87) yield
| G 2 k ( v ) | c Γ h 2 τ ( k 1 ) τ ( k + 1 ) τ { t e } ( · , t ) 2 + { H ( t e ) } ( · , t ) 2 d t 1 / 2 v Γ , h .
As for G 3 k , taking into account (4.79) together with (4.85) and (3.16), we obtain:
| G 3 k ( v ) | τ 2 ( k 1 ) τ ( k + 1 ) τ { t t e } ( · , t ) Γ 2 d t 1 / 2 v Γ , h .
Then using the Trace Theorem (cf. (4.94)), we finally establish,
| G 3 k ( v ) | C T r τ 2 ( k 1 ) τ ( k + 1 ) τ { t t e } ( · , t ) 2 + { t t e } ( · , t ) 2 d t 1 / 2 v Γ , h .
Now collecting (4.89), (4.90), (4.91) and (4.92) we can write,
| F ¯ k | h C ^ 0 h 2 2 ( N + 2 ) τ ε 0 , ( k 1 ) τ ( k + 1 ) τ { H ( t t e ) } ( · , t ) 2 d t 1 / 2 + c Ω h ε 0 , 2 τ ( k 1 ) τ ( k + 1 ) τ { t t e } ( · , t ) 2 d t 1 / 2 + C h | ε | 1 , 2 τ ( k 1 ) τ ( k + 1 ) τ { t t e } ( · , t ) 2 d t 1 / 2 + ε 0 , τ 2 ( k 1 ) τ ( k + 1 ) τ { t t t e } ( · , t ) 2 d t 1 / 2 .
On the other hand (4.97), (4.98) and (4.100) yield,
| G ¯ k | Γ , h C B h 2 τ ( k 1 ) τ ( k + 1 ) τ { H ( t e ) } ( · , t ) 2 d t 1 / 2 + c Γ h 2 τ ( k 1 ) τ ( k + 1 ) τ { t e } ( · , t ) 2 + { H ( t e ) } ( · , t ) 2 d t 1 / 2 + C T r τ 2 ( k 1 ) τ ( k + 1 ) τ { t t e } ( · , t ) 2 + { t t e } ( · , t ) 2 d t 1 / 2 .
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 e ¯ h k : = e ^ h k e h k for k = 0 , 1 , 2 , , 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 , , M , when the e h k s are replaced with the e ¯ h k s, and e h j is replaced with 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 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 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 e ¯ h k .
  • Let us first define,
e ^ 1 h : = ( e ^ h 1 e ^ h 0 ) / τ e ¯ 1 h : = e ^ 1 h e 1 h , E ¯ 0 : = e ¯ 1 h ε h , h 2 + τ 2 2 e ¯ 1 h 2 + · e ¯ 1 h ε 1 2 + 1 2 e ¯ h 1 2 + e ¯ h 0 2 + 1 2 · e ¯ h 1 ε 1 2 + · e ¯ h 0 ε 1 2 + T | ε | 2 , e ¯ h 0 2
Recalling that e h 1 = e h 0 + τ e 1 h we have e ¯ h 1 = e ¯ h 0 + τ e ¯ 1 h . Thus, taking either A = e ¯ h 0 or A = · e ¯ h 0 and either B = τ e ¯ 1 h or B = τ · e ¯ 1 h , we apply twice the inequality A + B 2 / 2 A 2 + B 2 to (4.103) together with Lemma 4.3, to obtain,
E ¯ 0 3 2 e ¯ h 0 2 + · e ¯ h 0 ε 1 2 + τ 2 e ¯ 1 h 2 + · e ¯ 1 h ε 1 2 + T | ε | 2 , e ¯ h 0 2 + ( N + 2 ) ε 0 , e ¯ 1 h 2 .
Finally using the inequality · v ε 1 2 N ε 1 0 , v 2 v { H 1 ( Ω ) } N , after straightforward manipulations we easily infer from (4.104) that
E ¯ 0 3 + 3 N ε 1 0 , 2 e ¯ h 0 2 + τ 2 e ¯ 1 h 2 + T | ε | 2 , e ¯ h 0 2 + ( N + 2 ) ε 0 , e ¯ 1 h 2 .
We next use the obvious splitting e ¯ h 0 = ( e ^ h 0 e 0 ) + ( e 0 e h 0 ) , together with the trivial one,
e ¯ 1 h = e ^ h 1 e ^ h 0 τ e 1 h = ( { t e ^ h } | t = 0 e 1 ) + 1 τ 0 τ ( τ t ) t t e ^ h d t + ( e 1 e 1 h ) .
Then plugging (4.106) into (4.105), since i = 1 p A i 2 p i = 1 p A i 2 for any set of p functions or fields A i , we obtain:
E ¯ 0 3 + 3 N ε 1 0 , 2 2 ( e ^ h 0 e 0 ) 2 + ( e 0 e h 0 ) 2 3 τ 2 t ( e ^ h e ) ( · , 0 ) 2 + 1 τ 2 0 τ ( τ t ) t t e ^ h d t 2 + ( e 1 e 1 h ) 2 + 2 T | ε | 2 , e ^ h 0 e 0 2 + e 0 e h 0 2 + 3 ( N + 2 ) ε 0 , t ( e ^ h e ) ( · , 0 ) 2 + 1 τ 2 0 τ ( τ t ) t t e ^ h d t 2 + e 1 e 1 h 2 .
Owing to a trivial upper bound and to the Cauchy-Schwarz inequality we easily come up with
1 τ 2 0 τ ( τ t ) t t e ^ h d t 2 0 τ | t t e ^ h | d t 2 τ 0 τ | t t e ^ h | 2 d t 1 / 2 2 .
It follows that,
1 τ 2 0 τ ( τ t ) t t e ^ h d t 2 τ 0 τ t t e ^ h 2 d t ( 1 + N ε 1 0 , ) τ 0 τ t t e 2 d t .
The last inequality in (4.108) follows from the definition of e ^ h . Indeed we know that,
( t t e ^ h , v ) + ( · t t e ^ h , · v ) ε 1 = ( t t e , v ) + ( · t t e , · v ) ε 1 v V h .
Taking v = { t t e ^ h } ( · , t ) , by the Cauchy-Schwarz inequality we easily obtain:
{ t t e ^ h 2 + · t t e ^ h ε 1 2 } 1 / 2 { t t e 2 + · t t e ε 1 2 } 1 / 2 ,
or yet
t t e ^ h 1 + N ε 1 0 , t t e ,
which validates (4.108).
  • On the other hand according to (4.63) we have t t e ^ h t t e + C ^ 0 h 2 H ( t t e ) . This yields,
0 τ ( τ t ) t t e ^ h d t 2 τ 0 τ t t e ^ h 2 d t 2 τ 0 τ t t e 2 d t + C ^ 0 2 h 4 0 τ H ( t t e ) 2 d t .
Incidentally we point out that Assumption * and (4.64) allow us to assert that
  • t t e { L 2 ( Ω × ( 0 , T ) ) } N ;
  • t t e { L ( 0 , T ) } N ;
  • { t e } | t = 0 = e 1 { H 2 ( Ω ) } N ;
  • { e } | t = 0 = e 0 { 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,
t ( e ^ h e ) ( · , 0 ) C ^ 1 h H ( e 1 ) , t ( e ^ h e ) ( · , 0 ) C ^ 0 h 2 H ( e 1 ) .
Then taking into account that τ 1 / ( 2 η ) , from (4.107)-(4.108)-(4.112)-(4.113) and Assumption*, we conclude that there exists a constant C ¯ 0 depending on Ω , T and ε , but neither on h nor on τ , such that,
E ¯ 0 C ¯ 0 2 ( h 2 + τ h 2 + h 4 ) H ( e 0 ) 2 + ( h 2 + τ 2 h 2 ) H ( e 1 ) 2 + τ 3 t t e L 2 ( 0 , T ) 2 + τ 2 t t e L ( 0 , T ) 2 + τ h 4 0 τ H ( t t e ) 2 d t .
Notice that, starting from (4.64) with u = t t e , similarly to (4.65), we obtain
t t e L ( 0 , T ) C T { t t e L 2 ( 0 , T ) 2 + t t t e ) L 2 ( 0 , T ) 2 } 1 / 2 .
Thus noting that h d i a m ( Ω ) , using again the upper bound τ 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,
E ¯ 0 C 0 2 h 2 { H ( e 0 ) 2 + H ( e 1 ) 2 + H ( t t e ) L 2 ( Ω × ( 0 , T ) ) 2 } + τ 2 { t t e L 2 ( Ω × ( 0 , T ) ) 2 + t t t e L 2 ( Ω × ( 0 , T ) ) 2 + t t e L 2 ( Ω × ( 0 , T ) ) 2 } .

4.3.2. Error estimates

In order to fully exploit the stability inequality (4.54) we further define,
S ¯ M : = k = 1 M 1 τ 2 | F ¯ k | h 2 + N 2 τ d ¯ k 2 + τ 2 G ¯ k Γ , h 2 .
According to (4.80), in order to estimate S ¯ M under the regularity Assumption* on e , we resort to the estimates (4.81), (4.101) and (4.102). Using the inequality | a · b | | a | | b | for a , b n , it is easy to see that there exists two constants C ¯ F and C ¯ G independent of h and τ such that
k = 1 M 1 τ 2 | F ¯ k | h 2 C ¯ F h 2 + τ 2 | e | H 4 ( Ω × ( 0 , T ) ) 2 ,
k = 1 M 1 τ 2 | G ¯ k | Γ , h 2 C ¯ G h 2 + τ 2 | e | H 3 ( Ω × ( 0 , T ) ) 2 .
On the other hand, recalling (4.65) we have k = 1 , 2 , , M 1 :
H ( e k ) 2 C T 2 0 T H ( e ) 2 + H ( t e ) 2 d t .
Therefore, since M = T / τ we have:
k = 1 M 1 H ( e k ) 2 C T 2 T τ 0 T H ( e ) 2 + H ( t e ) 2 d t
It follows from (4.120) and (4.81) that,
k = 1 M 1 N 2 τ d ¯ k 2 C ^ 0 2 C T 2 T | ε | 1 , 2 N h 4 2 τ 2 | e | H 2 ( Ω × ( 0 , T ) ) 2 + | e | H 3 ( Ω × ( 0 , T ) ) 2
Plugging (4.117), (4.121) and (4.118) into (4.116) we can assert that there exists a constant C ˜ depending on Ω , T and ε , but neither on h nor on τ , such that,
S ¯ M C ˜ 2 h 2 + τ 2 + h 4 τ 2 e H 4 ( Ω × ( 0 , T ) ) 2
Now recalling (4.54)-(4.53), provided (4.46) holds, together with τ 1 / ( 2 η ) , we have:
e ¯ h m + 1 e ¯ h m τ 2 + e ¯ h m + 1 2 + e ¯ h m 2 1 / 2 2 S ¯ M + E ¯ 0 e β / 2 .
This implies that, for m = 1 , 2 , , M 1 , it holds:
e h m + 1 e h m τ e m + 1 e m τ 2 + e h m + 1 e m + 1 2 + e h m e m ) 2 1 / 2 e ^ h m + 1 e m + 1 τ e ^ h m e m τ 2 + e ^ h m + 1 e m + 1 2 + e ^ h m e m 2 1 / 2 + 2 S ¯ M + E ¯ 0 e β / 2 .
Let us define a function e h in Ω ¯ × [ 0 , T ] whose value at t = k τ equals e h k for k = 1 , 2 , , M and that varies linearly with t in each time interval [ ( k 1 ) τ , k τ ] , in such a way that t e h ( x , t ) = e h k ( x ) e h k 1 ( x ) τ for every x Ω ¯ and t ( ( k 1 ) τ , k τ ) .
  • Now we define A ˜ m 1 / 2 ( · ) for any function or field A ( · , t ) to be the mean value of A ( · , t ) in ( m τ , ( m + 1 ) τ , that is A ˜ m + 1 / 2 = τ 1 m τ ( m + 1 ) τ A ( · , t ) d t . Clearly enough we have
e h m + 1 e h m τ e m + 1 e m τ 2 = { t e h ˜ } m + 1 / 2 { t e ˜ } m + 1 / 2 2 .
and also recalling (4.120) and (4.65)
e ^ h m + 1 e m + 1 τ e ^ h m e m τ 2 = m τ ( m + 1 ) τ t e ^ h t e τ d t 2 m τ ( m + 1 ) τ ( t e ^ h t e ) 2 d t 1 / 2 τ 2 = 1 τ m τ ( m + 1 ) τ t e ^ h t e 2 d t C ^ 0 2 h 4 τ m τ ( m + 1 ) τ H ( t e ) 2 d t C T 2 C ^ 0 2 h 4 0 T { H ( t e ) 2 + H ( t t e ) 2 } d t ,
which implies that
e ^ h m + 1 e m + 1 τ e ^ h m e m τ 2 C T 2 C ^ 0 2 h 2 | d i a m ( Ω ) | 2 e H 4 ( Ω × ( 0 , T ) ) 2 .
On the other hand from (4.58) and (4.65) we have for j = m or j = m + 1 :
( e ^ h j e j ) 2 C ^ 1 2 h 2 H ( e j ) 2 C ^ 1 2 C T 2 h 2 0 T H ( e ) 2 + H ( t e ) 2 d t .
which yields,
( e ^ h j e j ) 2 C T 2 C ^ 1 2 h 2 e H 3 ( Ω × ( 0 , T ) ) 2 .
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 C ˜ independent of h and τ it holds,
{ t e ˜ } m + 1 / 2 { t e } m + 1 / 2 C ˜ τ 2 e H 4 ( Ω × ( 0 , T ) ) ,
where for every function g defined in Ω ¯ × [ 0 , T ] and s + , g s is the function defined in Ω by g s ( · ) = g ( · , s τ ) .
  • Noting that { t e h ˜ } m + 1 / 2 ( · ) is nothing but { t e h } m + 1 / 2 ( · ) , 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 τ 1 / ( 2 η ) , under Assumption * on e , there exists a constant C depending only on Ω ,  ε and T such that
max 1 m M 1 { t ( e h e ) } m + 1 / 2 + max 2 m M ( e h m e m ) C ( τ + h + h 2 / τ ) e H 4 ( Ω × ( 0 , T ) ) + H ( e 0 ) + H ( e 1 ) .  
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 Ω by Ω 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 · is replaced by the norm · Ω h . Actually, unlike we did so far, in this section the notation · 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 : = 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 ) : = e r 2 t ε ( r ) .
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 = ( x 2 , x 1 ) e r ε ( r )   and   e 1 : = { t e } | t = 0 = 2 ( x 2 , x 1 ) e r ε ( r ) . Since · ( ε e ) 0 in Ω by construction, the right hand side f = ( f 1 ; f 2 ) is given by f : = ε t t e Δ e + ( · e ) , that is,
f = ( 4 x 2 e r 2 t Δ e 1 ; 4 x 1 e r 2 t Δ e 2 ) , where   ( Δ e 1 ; Δ e 2 ) = ( 3 x 2 r v θ x 2 v θ ; 3 x 1 r v θ + x 1 v θ ) ; with   v θ = ε ε ε 2 e r 2 t , v θ = ε 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 ) ) .
Figure 1. Function ε ( x , y ) in the domain Ω = ( 0 , 1 ) × ( 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 ) × ( 0 , 1 ) for different values of m in (5.130) on the mesh with h = 2 5 .
Preprints 97598 g001
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 × 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 ) × ( 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 if x 1 x 2 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 × 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,
e l 1 = max 1 k M e k e h k h max 1 k M e k h , e l 2 = max 1 k M ( e k e h k ) h max 1 k M e k h , e l 3 = max 1 k M 1 { t ( e e h ) } k + 1 / 2 h max 1 k M 1 { t e } k + 1 / 2 h .
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)
Preprints 97598 g002
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)
Preprints 97598 g003
Observation of these tables and figures clearly indicates that our scheme behaves like a first order method in the (semi-)norm of L [ ( 0 , T ) ; H 1 ( Ω ) ] for e and in the norm of L [ ( 0 , T ) ; L 2 ( Ω ) ] for 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 ( Ω × ( 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 [ ( 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 ( · ε e h k , · 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

  1. R.A. Adams, Sobolev Spaces, Academic Press, N.Y., 1975.
  2. 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]
  3. 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]
  4. 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]
  5. 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]
  6. L. Beilina, V. Ruas, An explicit P1 finite element scheme for Maxwell’s equations with constant permittivity in a boundary neighborhood, arXiv:1808.10720, 2020.
  7. 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.
  8. L. Beilina, V. Ruas, On the Maxwell-wave equation coupling problem and its explicit finite-element solution, Applications of Mathematics, open access (2022). [CrossRef]
  9. 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.
  10. A. Bossavit, Computational Electromagnetism, Variational Formulations, Complementary, Edge Elements, in: Electromagnetism, vol. 2, Academic Press, New York, 1998.
  11. 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.
  12. 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]
  13. 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]
  14. P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North Holland, 1978.
  15. P. Ciarlet Jr, J. Zou, Fully discrete finite element approaches for time-dependent Maxwell’s equations, Numerische Mathematik, 82-2 (1999), 193–219. [CrossRef]
  16. P. Ciarlet Jr, Augmented formulations for solving Maxwell equations, Computer Methods in Applied Mechanics and Engineering, 194-2–5 (2005), 559–586. [CrossRef]
  17. 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]
  18. G. Cohen, P. Monk, Gauss point mass lumping schemes for Maxwell’s equations, Numerical Methods for Partial Differential Equations, 14-1 (1998), 63–88.
  19. M. Costabel, A coercive bilinear form for Maxwell’s equations, J. Math. Anal. Appl., 157-2 (1991), 527–541. [CrossRef]
  20. 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.
  21. 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]
  22. A. Ditkowski and D. Gottlieb, On the Engquist–Majda absorbing boundary conditions for hyperbolic systems, Contemporary Mathematics, 330 (2003), 55–72.
  23. 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.
  24. B. Engquist, A. Majda, Absorbing boundary conditions for the numerical simulation of waves, Mathematics of Computation, 31 (1977), 629-651.
  25. V. Girault, P.A. Raviart, Finite Element Methods for Navier-Stokes Equations, Springer Series in Computational Mathematics, Vol. 5, Springer-Verlag, Berlin, 1986.
  26. P. Grisvard, Sigularities in Boundary Value Problems, Masson, Paris, 1992.
  27. E. Jamelot, Résolution des équations de Maxwell avec des éléments finis de Galerkin continus, thèse doctorale, École Polytechnique, Palaiseau, France, 2005.
  28. P. Monk. Finite Element Methods for Maxwell’s Equations, Clarendon Press, 2003.
  29. J.-C. Nédélec, Mixed finite elements in R3, Numerische Mathematik, 35 (1980), 315-341. [CrossRef]
  30. 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.
  31. V. Ruas, Numerical Methods for Partial Differential Equations; an Introduction, Wiley, 2016.
  32. V. Ruas and M.A. Silva Ramos, A Hermite Method for Maxwell’s Equations, Appl. Maths & Information Sci., 12-2 (2018), 271–283. [CrossRef]
  33. V. Thomée, Galerkin Finite Element Methods for Parabolic Problems, Springer Series in Computational Mathematics, Second edition, 1997.
  34. 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.
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.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

Subscribe

© 2024 MDPI (Basel, Switzerland) unless otherwise stated