Advances in Numerical Analysis
Volume 2010 (2010), Article ID 923832, 28 pages
doi:10.1155/2010/923832
Research Article

Finite Element Approximation of Maxwell's Equations with Debye Memory

Simon Shaw

BICOM, Institute of Computational Mathematics, Brunel University, Uxbridge, Middlesex UB8 3PH, UK

Received 18 May 2010; Accepted 31 August 2010

Academic Editor: Dongwoo Sheen

Copyright © 2010 Simon Shaw. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Maxwell's equations in a bounded Debye medium are formulated in terms of the standard partial differential equations of electromagnetism with a Volterra-type history dependence of the polarization on the electric field intensity. This leads to Maxwell's equations with memory. We make a correspondence between this type of constitutive law and the hereditary integral constitutive laws from linear viscoelasticity, and we are then able to apply known results from viscoelasticity theory to this Maxwell system. In particular, we can show long-time stability by shunning Gronwall's lemma and estimating the history kernels more carefully by appeal to the underlying physical fading memory. We also give a fully discrete scheme for the electric field wave equation and derive stability bounds which are exactly analogous to those for the continuous problem, thus providing a foundation for long-time numerical integration. We finish by also providing error bounds for which the constant grows, at worst, linearly in time (excluding the time dependence in the norms of the exact solution). Although the first (mixed) finite element error analysis for the Debye problem was given by Li (2007), this seems to be the first time sharp constants have been given for this problem.

1. Introduction

The potential for noninvasive electromagnetic detection of biological anomalies in the human body (and of defects in structural artefacts) has been recently noted by Banks et al. in [1]. Such a diagnostic technology has the promise of being quick, painless, cheap, and readily deployable, [2]. Biological tissue is dispersive, and electromagnetic constitutive relationships represent this frequency dependence through “complex moduli” (see, e.g., the data for grey and white brain tissue in [3, 4]). These are manifested though hysteretic fading memory Volterra operators or as an equivalent set of evolution equations, [1, 5]. The basic idea behind such a diagnostic technology is to compare the result of the transmission of electromagnetic waves through a patient's tissue against a datum outcome from a healthy subject. The presence of unwelcome tumours or lesions would be signalled by a change in the lossy response of the biological dielectric.

A possibly more sinister side of the interaction of electromagnetic waves with biological tissue is the concern that wireless communication devices (e.g., cell phones and mobile internet) may have a long-term detrimental effect on health. As far back as 1977 one can find (at http://oai.dtic.mil/oai/oai?verb=getRecord&metadataPrefix=html&identifier=ADA051218 (accessed 26 August 2009)) the following abstract of a workshop dedicated to the concerns that increased use of microwave radiation might raise [6].

The wide application of industrial, commercial and military devices and systems which radiate frequencies in the radiofrequency and microwave portion of the electromagnetic spectrum plus numerous only partially understood indications of microwave effects upon living organisms have raised important questions of the physical basis of the interactions of electromagnetic fields with biological systems. These questions must be answered if the development of regulatory standards and of methods and techniques for controlling radiofrequency and microwave exposure is to be achieved. The same questions must be answered in connection with present and proposed therapeutic applications of these waves. The rapid increase in the use of these frequencies makes these questions matters of imperative concern, particularly in view of the possibilities of cumulative or delayed effects of exposure. The purpose of the Workshop on the Physical Basis of Electromagnetic Interactions with Biological Systems was to bring together the leading investigators in the field to present the results of recent research, to determine the present status of the field and the priority of significant problem areas, and to critically evaluate conflicting theoretical interpretations and experimental techniques.

We can also find more recent contextualising discussions of the potential effects, both good and bad, of the dramatic increase in the use of nonionizing electromagnetic radiation in our world: [7, 8].

The purpose of this paper is to formulate Maxwell's equations with a Debye (i.e., “lossy”) polarization model and give some stability estimates for two continuous models as well as a discrete model, and also some error bounds for the discrete model. We rely on a strong analogy with models of viscoelasticity to provide some high quality estimates in terms of long-time behaviour. To our knowledge this is the first time such an analogy has been used in this context. To elaborate further it is required that we first introduce the mathematical model.

Let Ω be a bounded polyhedral domain in ℝ 3 and 𝐼 ∶ = ( 0 , 𝑇 ] a finite time interval. We consider Maxwell's equations (e.g., Monk, [9]) in Ω × 𝐼 in the form, ̇ ̇ ∇ × 𝐄 + 𝐁 = 𝟎 , ∇ × 𝐇 − 𝐃 = 𝐉 , ∇ ⋅ 𝐃 = 𝜚 , ∇ ⋅ 𝐁 = 0 , ( 1 . 1 ) where here and throughout the overdot denotes partial time differentiation. In these 𝐄 is the electric field intensity, 𝐇 is the magnetic field intensity, 𝐃 is the electric displacement, and 𝐁 is the magnetic induction. The first equation is Faraday's law, the second is Ampere's modified circuital law, the last two are Gauss's law and Gauss's magnetic law, 𝜚 is the total charge density, and 𝐉 is the current density. Conservation of charge implies that ∇ ⋅ 𝐉 + ̇ 𝜚 = 0 and, combined with the above, means that if the Gauss laws hold for one time, 𝑡 , then they hold for all 𝑡 . We therefore build these conditions into the intial data, and take the opportunity to specify the boundary data. ̂ ̂ ∇ ⋅ 𝐃 = 𝜚 , ∇ ⋅ 𝐁 = 0 i n Ω a t 𝑡 = 0 , ( 1 . 2 ) 𝐧 × 𝐄 = 𝟎 , 𝐧 ⋅ 𝐇 = 0 o n 𝜕 Ω × 𝐼 , ( 1 . 3 ) where ̂ 𝐧 is the unit outward normal on the boundary.

We assume that 𝐁 = 𝜇 𝐇 where 0 < 𝜇 ∈ ℝ , which is standard, but are concerned here with the Debye model of a dielectric wherein 𝐃 is given as a linear functional of 𝐄 . This will introduce “history integrals” and, therefore, nonlocality in time (but the constitutive law will remain local in space). Our contribution lies in drawing an analogy between the Debye polarization model and the fading memory constitutive laws governing creep and relaxation in linear viscoelastic media. This will allow us to apply known estimates from viscoelasticity theory and derive long-time stability estimates for the solutions to Maxwell's equations (Gronwall's inequality is not used). We also demonstrate a fully discrete time stepping scheme, based spatially on finite elements and temporally on central finite differences, for which these stability estimates carry over exactly. This provides a foundation for long-time numerical integration of the system, and in this spirit we also present an outline error estimate to show that again Gronwall's inequality may be avoided. It is also worth noting that sharp stability estimates for the dual problem are often needed when deriving residual-based a posteriori error bounds for time dependent problems; see, for example, [10]. In such cases it is essential to avoid the use of Gronwall's inequality whenever possible.

The original motivations for this work are the papers by Banks et al. [1] and Young and Adams [5], and we will return to these below. We note also that Fabrizio and Morro in [11] have a wealth of information on dielectrics and conductors with memory, for example, general observations on memory effects in Chapter 4, dielectrics and conductors with memory in Chapter 5.2, as well as more general formulations where all three constitutive laws (including Ohm's law) have fading memory in Chapter 9.3.3. They also study the temporally-asymptotic behaviour of dielectrics with memory in Chapter 7.10 using the Fourier transform and in this sense the general flavour of the results given below will come as no surprise. However, we use “energy methods” and explicit representations of the history kernels which are far more revealing in terms of the damping mechanisms at work. These “energy methods” are also well suited to the application to numerical schemes in which the constants can be estimated more carefully (see, e.g., [12–15]).

Besides the temporal decay properties with which we are concerned below, we also would like to briefly remark on the spatial decay exhibited by Debye media. For example, in [16] Roberts and Petropoulos examine the amplitude asymptotics of short and long pulses with application to muscle mimicking materials. They are able to derive tight upper and lower bounds on the exponential spatial decay rates of the temporal 𝐿 𝑝 norm of the electric field, and they also discuss the nontrivial structure of the “wave hierarchy” induced in the dielectric. This structure was (seemingly) first uncovered in [17] where the 𝑀 + 1 Debye “modes” were combined to form a high-order PDE for the electric field.

To provide some more context for both the physical and mathematical aspects of this work, we note that Lakes in [18, Chapter 2] discusses constitutive equations for piezoelectric materials with memory. Here both strain and electric displacement are related to stress, electric field, and temperature through convolution integrals. On the other hand, in terms of numerical approximation schemes, it seems that the first time domain finite element methods (TDFEM, as in [19]) for Maxwell's equations in lossy media were proposed in [20]. They considered cold plasma, Debye and Lorentz media, with single or multiple poles, where the hysteretic effects were captured through a convolution, or “history”, integral. A discontinuous Galerkin TDFEM is given in [21] along with extensive numerics that illustrate also how the lossy medium can be coupled to a perfectly matched layer through the “auxiliary differential [relaxation] equations”, as opposed to a history integral. The first error estimates for Maxwell's equations with loss were developed by Li in [22] for a semidiscrete Galerkin approximation to the vector wave equation for the electric field in cold plasma. This work was extended to an implicit three-stage first-order time stepping scheme in [23], and then to the error analysis of a fully discrete mixed finite element method for cold plasma as well as one-pole Debye and two-pole Lorentz media in [24]. Furthermore, a posteriori error bounds for an interior penalty discontinuous Galerkin method are reported in [25].

Our presentation here is intended to complement Li's efforts in that we provide similar (Galerkin) formulations and error bounds but use the strong analogy to viscoelastic fading memory in order to avoid the use of Gronwall's inequality.

This paper is laid out as follows. In Section 2 we briefly review the aspects of linear viscoelasticity that are relevant and then in Section 3 discuss the Debye model and illustrate the similarities to viscoelasticity. This will allow us to make some general assumptions on the Debye kernels which will form the basis of our estimates. In Section 4 we formulate Maxwell's equations with memory as a problem for the electric displacement, 𝐃 , in Section 4.1 and also as a problem for the electric field, 𝐄 , in Section 4.2. In both cases, long-time stability estimates (Theorems 4.3 and 4.5) are derived without recourse to Gronwall's inequality. We give a fully discrete scheme for the second of these problems in Section 5 and show in Theorem 5.4 that the stability estimate continues to hold. A formal a priori error bound is stated in Theorem 5.6 which reveals that the constants are not time dependent. It follows from this that the only source of temporal error growth must come from the accumulation effects arising directly from the norms of the exact solution. By making standard assumptions on the approximation properties of the scheme this dependence is made explicit in Corollary 5.9, where we see that the convergence rate is optimal. Finally, Section 6 contains a few concluding remarks.

Everywhere below we will assume that Ohm's law holds so that 𝐉 = 𝜎 𝐄 + 𝐉 𝑎 , where 𝐉 𝑎 is known, but will not exclude the case where the conductivity, 𝜎 , is nonzero since the medium may not be a perfect dielectric and so may exhibit some limited conductivity. At the mathematical level this term is included for the sake of generality and to demonstrate that it can easily be dealt with by using the estimates given below. Also, for the sake of simplicity we assume that all physical parameters are constant in both space and time.

Specific notation will be explained at the point where it is introduced but otherwise our notation is standard. In particular, 𝐿 𝑝 ( Ω ) , 𝑊 𝑚 𝑝 ( Ω ) , and 𝐻 𝑚 ( Ω ) are the usual Lebesgue, Banach and Hilbert spaces and we set 𝐋 𝑝 ( Ω ) ∶ = 𝐿 𝑝 ( Ω ) 3 , 𝐇 𝑚 ( Ω ) ∶ = 𝐻 𝑚 ( Ω ) 3 , and so on. In general terms, if 𝑋 is a Banach space then its norm will be denoted by ‖ ⋅ ‖ 𝑋 . The only exception to this will be that ‖ ⋅ ‖ 0 will always mean the 𝐿 2 ( Ω ) norm as induced by ( ⋅ , ⋅ ) , the 𝐿 2 ( Ω ) inner product. We will not distinguish norms and inner products between the cases where they act on scalar- or vector-valued functions and time dependence is dealt with in the usual way. Indeed, given a map 𝑓 ∶ [ 0 , 𝑇 ] → 𝑋 , then we use standard notation such as ‖ 𝑓 ‖ 𝐿 𝑝 ( 0 , 𝑇 ; 𝑋 ) to signify the 𝐿 𝑝 ( 0 , 𝑇 ) norm of ‖ 𝑓 ( 𝑡 ) ‖ 𝑋 .

The natural space for Maxwell's equations as we will work with them is 𝐇 ( 𝐜 𝐮 𝐫 𝐥 ; Ω ) where 𝐇 ( 𝐜 𝐮 𝐫 𝐥 ; Ω ) ∶ = 𝐯 ∈ 𝐋 2 ( Ω ) ∶ ∇ × 𝐯 ∈ 𝐋 2 ( Ω ) ( 1 . 4 ) with graph norm, ‖ 𝐯 ‖ 𝐇 ( 𝐜 𝐮 𝐫 𝐥 ; Ω ) ∶ = ( ‖ 𝐯 ‖ 2 0 + ‖ ∇ × 𝐯 ‖ 2 0 ) 1 / 2 . This is a Hilbert space when equipped with the obvious norm-inducing inner product (see, e.g., [9]). We also define, 𝑉 ∶ = 𝐇 0 ̂ ( 𝐜 𝐮 𝐫 𝐥 ; Ω ) ∶ = { 𝐯 ∈ 𝐇 ( 𝐜 𝐮 𝐫 𝐥 ; Ω ) ∶ 𝐧 × 𝐯 = 𝟎 o n 𝜕 Ω } , ( 1 . 5 ) (also a Hilbert space) where ̂ 𝐧 is a.e. the unit outward normal on 𝜕 Ω . For Friedrich's-type properties of 𝑉 and associated spaces we refer to Monk in [9, Corollary 4.8].

Finally, in terms of preliminary notation, we recall Young's inequality for products, 𝑎 2 𝑎 𝑏 ⩽ 2 𝜖 + 𝜖 𝑏 2 , ∀ 𝑎 , 𝑏 ∈ ℝ , ∀ 𝜖 > 0 , ( 1 . 6 ) and Young's inequality for convolutions, ‖ 𝑤 ∗ 𝑣 ‖ 𝐿 𝑟 ( 0 , 𝑡 ) ⩽ ‖ 𝑤 ‖ 𝐿 𝑝 ( 0 , 𝑡 ) ‖ 𝑣 ‖ 𝐿 𝑞 ( 0 , 𝑡 ) w h e r e ( 𝑤 ∗ 𝑣 ) ( 𝑡 ) ∶ = 𝑡 0 𝑤 ( 𝑡 − 𝑠 ) 𝑣 ( 𝑠 ) 𝑑 𝑠 , ( 1 . 7 ) and 1 ⩽ 𝑝 , 𝑞 , 𝑟 ⩽ ∞ such that 𝑝 − 1 + 𝑞 − 1 = 1 + 𝑟 − 1 .

2. Elements of Viscoelasticity Theory

Before we recall the Debye model, it will be instructive to outline some basic facts from the constitutive theories of linear viscoelasticity theory as described in, for example, [26–28]. For a synchronous viscoelastic substance, the stress tensor, 𝝈 , is given as a linear functional of the strain tensor, 𝜺 , by the constitutive law in the following stress relaxation form: 𝝈 ( 𝑡 ) = 𝐂 𝜺 ( 𝑡 ) − 𝑡 0 𝜑 𝑠 ( 𝑡 − 𝑠 ) 𝐂 𝜺 ( 𝑠 ) 𝑑 𝑠 . ( 2 . 1 ) This can be inverted to the following creep form: 𝐂 𝜺 ( 𝑡 ) = 𝝈 ( 𝑡 ) − 𝑡 0 𝜓 𝑠 ( 𝑡 − 𝑠 ) 𝝈 ( 𝑠 ) 𝑑 𝑠 . ( 2 . 2 ) In these 𝜑 is a stress relaxation function and 𝜓 is a creep function. Also, the subscript denotes partial differentiation so that 𝜑 𝑠 ( 𝑡 − 𝑠 ) = − 𝜑 𝑡 ( 𝑡 − 𝑠 ) = − 𝜑 ′ ( 𝑡 − 𝑠 ) , for example, where the prime denotes partial differentiation with respect to the displayed time argument. Below we will use each of these notations with the aim of clarifying our reasoning in different places and contexts. For example, the 𝑠 -integral of an 𝑠 -derivative in the expressions above is useful later in Lemma 4.1, whereas the prime notation is useful in Assumptions 2.1 below where it neatly characterises the physical assumptions.

Causality and thermodynamics imply rather general conditions on 𝜑 and 𝜓 and, therefore, a good deal of freedom in ascribing specific forms. However, in many practical examples they typically have the Prony form, 𝜑 ( 𝑡 ) = 𝜑 0 + 𝑁 𝑖 = 1 𝜑 𝑖 𝑒 − 𝑡 / 𝜏 ′ 𝑖 , 𝜓 ( 𝑡 ) = 𝜓 0 + 𝑁 𝑖 = 1 𝜓 𝑖 𝑒 − 𝑡 / 𝜏 𝑖 ( 2 . 3 ) subject to normalization, 𝜑 ( 0 ) = 𝜓 ( 0 ) = 1 , and where 𝜑 𝑖 ⩾ 0 and 𝜏 𝑖 , 𝜏 𝑖 > 0 for each 𝑖 . When 𝜑 0 , 𝜓 0 > 0 we have a viscoelastic solid and when 𝜑 0 = 𝜓 0 = 0 we have a viscoelastic fluid (see, e.g., [14, 26] but note that the above form for 𝜓 will not apply in this case because 𝜓 ( 𝑡 ) ∼ 𝑡 as 𝑡 → ∞ ). The tensor 𝐂 is the Hooke tensor from linear elasticity theory and is subject to some standard symmetry and positivity requirements, but these are not important here.

If we denote the Laplace transformation as 𝔏 ∶ 𝑡 → 𝑝 and denote 𝔏 𝐹 by 𝐹 then it is readily shown that the relaxation and creep functions are related by 𝑝 2 𝜑 ( 𝑝 ) 𝜓 ( 𝑝 ) = 1 and it follows that, for a solid, 𝜑 ( ∞ ) 𝜓 ( ∞ ) = 𝜑 0 𝜓 0 = 1 . In particular, it is easy to derive the pair, 𝜑 ( 𝑡 ) = 𝜑 0 + 𝜑 1 𝑒 − 𝑡 / 𝜏 ⟺ 𝜓 ( 𝑡 ) = 𝜓 0 − 𝜓 1 𝑒 − 𝜑 0 𝑡 / 𝜏 w h e r e 𝜓 0 = 1 𝜑 0 , 𝜓 1 = 𝜑 1 𝜑 0 . ( 2 . 4 ) These simple observations illustrate the following more general properties (again, see [26]) which will be our assumptions for all of what follows.

Assumptions (creep-relaxation pairs). A normalized creep-relaxation pair subject to 𝜑 0 𝜓 0 = 1 satisfy: (i) 𝜑 ∈ 𝐶 3 ( 0 , ∞ ) with 𝜑 ( 0 ) = 1 , l i m 𝑡 ↑ ∞ 𝜑 ( 𝑡 ) = 𝜑 > 0 , 𝜑 ′ < 0 and 𝜑 > 0 ,(ii) 𝜓 ∈ 𝐶 3 ( 0 , ∞ ) with 𝜓 ( 0 ) = 1 , l i m 𝑡 ↑ ∞ 𝜓 ( 𝑡 ) = 𝜑 − 1 > 0 , 𝜓 ′ > 0 and 𝜓 < 0 , along with 𝑝 2 𝜓 ( 𝑝 ) 𝜑 ( 𝑝 ) = 1 .

We are now in a position to review the Debye model and draw the analogy to viscoelasticity.

3. The Debye Model

In the paper by Banks et al. [1], the electric displacement is written in terms of the electric field and polarization as 𝐃 = 𝜀 0 𝐏 𝐄 + where 𝜀 0 is the permittivity of free space and 𝐏 = 𝐏 𝐼 + 𝐏 . Here 𝐏 𝐼 = 𝜀 0 𝜒 𝐄 is the instantaneous polarization and 𝐏 is the “relaxation polarization”. With 𝜀 ∞ ∶ = 1 + 𝜒 , the constitutive model used in [1] is then 𝐃 = 𝜀 0 𝜀 ∞ ̇ 𝜀 𝐄 + 𝐏 w i t h 𝜏 𝐏 + 𝐏 = 𝑠 − 𝜀 ∞ 𝜀 0 ̆ 𝐄 s u b j e c t t o 𝐏 ( 0 ) = 𝐏 . ( 3 . 1 ) Note that here and throughout the dependence on the spatial variable 𝐱 ∈ Ω is suppressed. Integrating this ODE gives 𝐏 ( 𝑡 ) = 𝑒 − 𝑡 / 𝜏 ̆ 𝐏 + 𝑡 0 𝜀 𝑠 − 𝜀 ∞ 𝜀 0 𝜏 𝑒 − ( 𝑡 − 𝑠 ) / 𝜏 𝐄 ( 𝑠 ) 𝑑 𝑠 , ( 3 . 2 ) which, by henceforth assuming that ̆ 𝐏 = 𝟎 (since the instantaneous effects are already taken care of in 𝐏 𝐼 ), gives us 1 𝜀 0 𝜀 ∞ 𝜀 𝐃 ( 𝑡 ) = 𝐄 ( 𝑡 ) + 𝑠 − 𝜀 ∞ 𝜏 𝜖 ∞ 𝑡 0 𝑒 − ( 𝑡 − 𝑠 ) / 𝜏 𝐄 ( 𝑠 ) 𝑑 𝑠 . ( 3 . 3 ) The connection to viscoelasticity can now be made evident by defining the creep function 𝜓 ( 𝑡 ) ∶ = 1 − 𝜑 1 𝑒 − 𝜑 0 𝑡 / 𝜏 ′ 𝜑 0 = 𝜓 0 − 𝜓 1 𝑒 − 𝑡 / 𝜏 , ( 3 . 4 ) where 𝜏 ′ ∶ = 𝜑 0 𝜏 , 𝜓 0 = 1 / 𝜑 0 , 𝜓 1 = 𝜑 1 / 𝜑 0 , 𝜑 0 ∶ = 𝜀 ∞ / 𝜀 𝑠 , and 𝜑 1 ∶ = ( 𝜀 𝑠 − 𝜀 ∞ ) / 𝜀 𝑠 , and then for 𝐂 ∶ = ( 𝜀 0 𝜀 ∞ ) − 1 𝐈 , we have, 𝐂 𝐃 ( 𝑡 ) = 𝐄 ( 𝑡 ) − 𝑡 0 𝜓 𝑠 ( 𝑡 − 𝑠 ) 𝐄 ( 𝑠 ) 𝑑 𝑠 . ( 3 . 5 ) Immediately we can observe the analogy, or correspondence, between the electric displacement and field, 𝐃 , 𝐄 , and the viscoelastic strain and stress, 𝜺 , 𝝈 . We can also conclude that this relationship can be inverted to give, 𝐄 ( 𝑡 ) = 𝐂 𝐃 ( 𝑡 ) − 𝐂 𝑡 0 𝜑 𝑠 ( 𝑡 − 𝑠 ) 𝐃 ( 𝑠 ) 𝑑 𝑠 ( 3 . 6 ) for the relaxation function 𝜑 ( 𝑡 ) = 𝜑 0 + 𝜑 1 𝑒 − 𝑡 / 𝜏 ′ . It is clear that 𝜑 ( 0 ) = 𝜓 ( 0 ) = 1 and we also have 𝜑 0 = 𝜑 ( ∞ ) = 1 / 𝜓 ( ∞ ) = 1 / 𝜑 0 . On physical grounds, we expect that 𝜀 0 > 0 and 0 < 𝜀 ∞ ⩽ 𝜀 𝑠 . We then deduce that 𝜑 is monotone decreasing but bounded above zero by 𝜑 0 > 0 and that 𝜓 is monotone increasing but bounded above by 1 / 𝜑 0 . It therefore follows that this pair satisfy Assumptions 2.1.

Now, under the static conditions where 𝐄 is constant and 𝜀 𝑠 > 𝜀 ∞ we have 𝐃 ( 𝑡 ) = 1 − 𝑡 0 𝜓 𝑠 ( 𝜀 𝑡 − 𝑠 ) 𝑑 𝑠 0 𝜀 ∞ 𝐄 = 𝜓 ( 𝑡 ) 𝜀 0 𝜀 ∞ 𝐄 ( 3 . 7 ) and so, with ‖ ⋅ ‖ 𝔼 the Euclidean norm on ℝ 3 , we have ‖ 𝐃 ( 𝑡 1 ) ‖ 𝔼 < ‖ 𝐃 ( 𝑡 2 ) ‖ 𝔼 for all 0 ⩽ 𝑡 1 < 𝑡 2 < ∞ . We can see that this constitutive law represents a monotone creeping process for the electric displacement which is bounded by ‖ 𝐃 ( 0 ) ‖ 𝔼 = 𝜀 0 𝜀 ∞ ‖ 𝐄 ‖ 𝔼 and ‖ 𝐃 ( ∞ ) ‖ 𝔼 = 𝜑 0 − 1 𝜀 0 𝜀 ∞ ‖ 𝐄 ‖ 𝔼 .

We note that the Debye model presented in [11, Chapter 4, Section 2] agrees with that just presented if, in [11], we set 𝜀 0 ∶ = 𝜀 0 𝜀 ∞ and 𝜀 ′ ∶ = 𝜀 𝑠 𝜀 0 .

Just as in viscoelasticity theory the relaxation and creep functions can be represented as a Prony series of decaying exponentials rather than just the “single mode” model just given. This is the model presented by Young and Adams in [5]. They write 𝐃 = 𝜀 0 𝜀 ∞ 𝐄 + 𝑀 𝑚 = 1 𝐏 𝑚 w i t h 𝜏 𝑚 ̇ 𝐏 𝑚 + 𝐏 𝑚 = 𝜀 𝑠 𝑚 − 𝜀 ∞ 𝜀 0 𝐄 s u b j e c t t o 𝐏 𝑚 ̆ 𝐏 ( 0 ) = 𝑚 . ( 3 . 8 ) The connection between this model and the Banks et al. model is clear and we have, 𝐏 𝑚 ( 𝑡 ) = 𝑒 − 𝑡 / 𝜏 𝑚 ̆ 𝐏 𝑚 + 𝜀 𝑠 𝑚 − 𝜀 ∞ 𝜀 0 𝜏 𝑚 𝑡 0 𝑒 − ( 𝑡 − 𝑠 ) / 𝜏 𝑚 𝐄 ( 𝑠 ) 𝑑 𝑠 . ( 3 . 9 ) Assuming, as before, that each ̆ 𝐏 𝑚 = 𝟎 we can therefore write 1 𝜀 0 𝜀 ∞ 𝐃 ( 𝑡 ) = 𝐄 ( 𝑡 ) − 𝑡 0 𝜓 𝑠 ( 𝑡 − 𝑠 ) 𝐄 ( 𝑠 ) 𝑑 𝑠 , ( 3 . 1 0 ) where 𝜓 ( 𝑡 ) ∶ = 𝜓 0 − 𝑀 𝑚 = 1 𝜓 𝑚 𝑒 − 𝑡 / 𝜏 𝑚 w i t h 𝜓 𝑚 𝜀 ∶ = 𝑠 𝑚 − 𝜀 ∞ 𝜀 ∞ ( 3 . 1 1 ) and 𝜓 0 defined so that 𝜓 ( 0 ) = 1 . Drawing again on the viscoelasticity analogy we conclude that again this constitutive law can be inverted with the aid of a relaxation function, 𝜑 , to give 1 𝐄 ( 𝑡 ) = 𝜀 0 𝜀 ∞ 1 𝐃 ( 𝑡 ) − 𝜀 0 𝜀 ∞ 𝑡 0 𝜑 𝑠 ( 𝑡 − 𝑠 ) 𝐃 ( 𝑠 ) 𝑑 𝑠 . ( 3 . 1 2 ) Typically, for the creep function in use here we would expect that 𝜑 ( 𝑡 ) = 𝜑 0 + 𝑀 𝑚 = 1 𝜑 𝑚 𝑒 − 𝑡 / 𝜏 ′ 𝑚 ( 3 . 1 3 ) but this, of course, needs proof. We assume only that the pair satisfy Assumptions 2.1 and give specific examples later in Section 6.

4. Maxwell Systems with Memory

In this section, we introduce the Debye model in each of the forms (3.10) and (3.12) into the field equations, (1.1), the intial and boundary data, (1.2) and (1.3), and 𝐁 = 𝜇 𝐇 to derive partial differential Volterra equations. Firstly, in Section 4.1 we use (3.12) to obtain a PDE for the electric displacement, 𝐃 , and then in Section 4.2 we use (3.10) to derive a PDE for the electric field, 𝐄 . Stability estimates that exploit the “viscoelastic fading memory” in place of Gronwall's lemma are given in each case. We also recall that, for simplicity, all physical parameters are assumed to be constant in space-time.

4.1. A Formulation for the Electric Displacement

From (1.1) and Ohm's law, we obtain ̈ 𝐃 + ∇ × 𝜇 − 1 ̇ ̇ 𝐉 ∇ × 𝐄 + 𝜎 𝐄 = − 𝑎 ( 4 . 1 ) which, recalling (1.5), can be written in weak form as ̈ + 𝜇 𝐃 , 𝜙 − 1 + 𝜎 ̇ ̇ 𝐉 ∇ × 𝐄 , ∇ × 𝜙 𝐄 , 𝜙 = − 𝑎 , 𝜙 ∀ 𝜙 ∈ 𝑉 . ( 4 . 2 ) However, from (3.12) with 𝜀 ∶ = 𝜀 0 𝜀 ∞ we obtain after a partial integration 𝜀 ̇ ̇ 𝐄 ( 𝑡 ) = 𝐃 ( 𝑡 ) + 𝜑 ′ ( 𝑡 ) 𝐃 ( 0 ) − 𝑡 0 𝜑 𝑠 ( ̇ 𝑡 − 𝑠 ) 𝐃 ( 𝑠 ) 𝑑 𝑠 . ( 4 . 3 ) Note that in addition to using an overdot for partial time differentiation, we are using a prime to denote differentiation with respect to the displayed variable and subscripts to denote partial differentiation in the usual way.

By introducing these terms into the weak form, we can pose the problem as: find a twice differentiable map, 𝐃 ∶ 𝐼 → 𝑉 such that, ̈ + 𝜇 𝐃 ( 𝑡 ) , 𝜙 − 1 ∇ × 𝜀 − 1 − 𝐃 ( 𝑡 ) , ∇ × 𝜙 𝑡 0 𝜑 𝑠 ( 𝜇 𝑡 − 𝑠 ) − 1 ∇ × 𝜀 − 1 + 𝜎 𝐃 ( 𝑠 ) , ∇ × 𝜙 𝑑 𝑠 𝜀 ̇ 𝐃 ( 𝑡 ) − 𝑡 0 𝜑 𝑠 ̇ ( 𝑡 − 𝑠 ) 𝐃 ( 𝑠 ) 𝑑 𝑠 , 𝜙 = ( 𝐆 ( 𝑡 ) , 𝜙 ) ( 4 . 4 ) with ̇ 𝐉 𝐆 ( 𝑡 ) = − 𝑎 ( 𝑡 ) − 𝜎 𝜀 − 1 𝜑 ′ ( 𝑡 ) 𝐃 ( 0 ) known and given intial data.

Before we give the stability estimate, we need two lemmas which help us deal with the fading memory in a sharper way than Gronwall's inequality would admit. The first is a coercivity result as used in [29], and the second is a remarkable lemma given originally in [12] by Rivera and Menzala. (The author is very grateful to Professor Igor Bock (Slovak University of Technology, Bratislava) for pointing out Rivera and Menzala's result.) Outline proofs are given for the sake of completeness.

Lemma 4.1 (coercivity of relaxation). If Assumptions 2.1, hold then 𝑡 0 𝐰 ( 𝑠 ) − 𝑠 0 𝜑 𝜉 ( 𝑠 − 𝜉 ) 𝐰 ( 𝜉 ) 𝑑 𝜉 , 𝐰 ( 𝑠 ) 𝑑 𝑠 ⩾ 𝜑 ‖ 𝐰 ‖ 2 𝐿 2 ( 0 , 𝑡 ; 𝐋 2 ( Ω ) ) ( 4 . 5 ) for all 𝑡 ∈ 𝐼 .

Proof. Since in general terms ( 𝑢 ± 𝑣 , 𝑢 ) ⩾ ‖ 𝑢 ‖ 2 − ‖ 𝑣 ‖ ‖ 𝑢 ‖ we need only to observe that ‖ ( 𝜑 ∗ 𝐰 , 𝐰 ) ‖ 𝐿 1 ( 0 , 𝑡 ) ⩽ ‖ ( − 𝜑 ∗ ‖ 𝐰 ‖ 0 ) ‖ 𝐰 ‖ 0 ‖ 𝐿 1 ( 0 , 𝑡 ) ⩽ ‖ − 𝜑 ‖ 𝐿 1 ( 0 , 𝑡 ) ‖ 𝐰 ‖ 2 𝐿 2 ( 0 , 𝑡 ; 𝐋 2 ( Ω ) ) from Young's inequality, (1.7). Assumptions 2.1 yield 1 − ‖ − 𝜑 ′ ‖ 𝐿 1 ( 0 , 𝑡 ) = 𝜑 ( 𝑡 ) ⩾ 𝜑 and this completes the proof.

Lemma 4.2 (Rivera and Menzala, [12]). Let 𝑎 ( 𝐰 , 𝐯 ) ∶ = ( ∇ × 𝐰 , ∇ × 𝐯 ) , then − 2 𝜏 0 𝑡 0 𝜑 𝑠 ( ̇ − 𝑡 − 𝑠 ) 𝑎 ( 𝐯 ( 𝑠 ) , 𝐯 ( 𝑡 ) ) 𝑑 𝑠 𝑑 𝑡 = 𝜑 ( 𝜏 ) 𝑎 ( 𝐯 ( 𝜏 ) , 𝐯 ( 𝜏 ) ) − 𝑎 ( 𝐯 ( 𝜏 ) , 𝐯 ( 𝜏 ) ) 𝜏 0 𝜑 ( − 𝜏 − 𝑠 ) 𝑎 ( 𝐯 ( 𝜏 ) − 𝐯 ( 𝑠 ) , 𝐯 ( 𝜏 ) − 𝐯 ( 𝑠 ) ) 𝑑 𝑠 𝜏 0 𝜑 ( + 𝑡 ) 𝑎 ( 𝐯 ( 𝑡 ) , 𝐯 ( 𝑡 ) ) 𝑑 𝑡 𝜏 0 𝑡 0 𝜑 ( 𝑡 − 𝑠 ) 𝑎 ( 𝐯 ( 𝑡 ) − 𝐯 ( 𝑠 ) , 𝐯 ( 𝑡 ) − 𝐯 ( 𝑠 ) ) 𝑑 𝑠 𝑑 𝑡 ( 4 . 6 ) for all 𝜏 ∈ 𝐼 and for all 𝐯 ∈ 𝐻 1 ( 𝐼 ; 𝑉 ) .

Proof. First we note that, 𝑑 𝑑 𝑡 𝑡 0 𝜑 𝑡 ( = 𝑡 − 𝑠 ) 𝑎 ( 𝐯 ( 𝑡 ) − 𝐯 ( 𝑠 ) , 𝐯 ( 𝑡 ) − 𝐯 ( 𝑠 ) ) 𝑑 𝑠 𝑡 0 𝜑 𝑡 𝑡 ( 𝑡 − 𝑠 ) 𝑎 ( 𝐯 ( 𝑡 ) − 𝐯 ( 𝑠 ) , 𝐯 ( 𝑡 ) − 𝐯 ( 𝑠 ) ) 𝑑 𝑠 − 2 𝑡 0 𝜑 𝑡 ̇ ( 𝑡 − 𝑠 ) 𝑎 ( 𝐯 ( 𝑠 ) , 𝐯 ( 𝑡 ) ) 𝑑 𝑠 + 2 𝑡 0 𝜑 𝑡 ( ̇ 𝑡 − 𝑠 ) 𝑎 ( 𝐯 ( 𝑡 ) , 𝐯 ( 𝑡 ) ) 𝑑 𝑠 , ( 4 . 7 ) and substitute for the last term on the right by using 𝑑 𝑑 𝑡 𝑡 0 𝜑 𝑡 ( 𝑡 − 𝑠 ) 𝑎 ( 𝐯 ( 𝑡 ) , 𝐯 ( 𝑡 ) ) 𝑑 𝑠 = 𝜑 ( 𝑡 ) 𝑎 ( 𝐯 ( 𝑡 ) , 𝐯 ( 𝑡 ) ) + 2 𝑡 0 𝜑 𝑡 ( ̇ 𝑡 − 𝑠 ) 𝑎 ( 𝐯 ( 𝑡 ) , 𝐯 ( 𝑡 ) ) 𝑑 𝑠 . ( 4 . 8 ) After a rearrangement of terms, this gives − 2 𝑡 0 𝜑 𝑠 ( ̇ = 𝑑 𝑡 − 𝑠 ) 𝑎 ( 𝐯 ( 𝑠 ) , 𝐯 ( 𝑡 ) ) 𝑑 𝑠 𝑑 𝑡 𝑡 0 𝜑 𝑡 ( 𝑡 − 𝑠 ) 𝑎 ( 𝐯 ( 𝑡 ) , 𝐯 ( 𝑡 ) ) 𝑑 𝑠 − 𝜑 − 𝑑 ( 𝑡 ) 𝑎 ( 𝐯 ( 𝑡 ) , 𝐯 ( 𝑡 ) ) 𝑑 𝑡 𝑡 0 𝜑 + ( 𝑡 − 𝑠 ) 𝑎 ( 𝐯 ( 𝑡 ) − 𝐯 ( 𝑠 ) , 𝐯 ( 𝑡 ) − 𝐯 ( 𝑠 ) ) 𝑑 𝑠 𝑡 0 𝜑 ( 𝑡 − 𝑠 ) 𝑎 ( 𝐯 ( 𝑡 ) − 𝐯 ( 𝑠 ) , 𝐯 ( 𝑡 ) − 𝐯 ( 𝑠 ) ) 𝑑 𝑠 . ( 4 . 9 ) The proof is then completed by integrating over ( 0 , 𝜏 ) , observing that 𝜏 0 𝜑 𝜏 ( 𝜏 − 𝑠 ) 𝑎 ( 𝐯 ( 𝜏 ) , 𝐯 ( 𝜏 ) ) 𝑑 𝑠 = ( 𝜑 ( 𝜏 ) − 𝜑 ( 0 ) ) 𝑎 ( 𝐯 ( 𝜏 ) , 𝐯 ( 𝜏 ) ) ( 4 . 1 0 ) and recalling that 𝜑 ( 0 ) = 1 .

Theorem 4.3 (stability for the 𝐃 equation). If Assumptions 2.1 hold and if 𝜇 > 0 , 𝜀 > 0 and 𝜎 ⩾ 0 are real numbers, then for every 𝑡 ∈ 𝐼 we have the following. Firstly, if 𝐆 = 𝟎 , ‖ ‖ ̇ ‖ ‖ 𝜇 𝜀 𝐃 ( 𝑡 ) 2 0 + 𝜑 ‖ ∇ × 𝐃 ( 𝑡 ) ‖ 2 0 ‖ ‖ ̇ 𝐃 ‖ ‖ + 2 𝜇 𝜎 𝜑 2 𝐿 2 ( 0 , 𝑡 ; 𝐋 2 ( Ω ) ) ‖ ‖ ̇ ‖ ‖ ⩽ 𝜇 𝜀 𝐃 ( 0 ) 2 0 + ‖ ∇ × 𝐃 ( 0 ) ‖ 2 0 . ( 4 . 1 1 ) Secondly, if 𝐆 ≠ 𝟎 and 𝜎 > 0 , ‖ ‖ ̇ ‖ ‖ 𝜇 𝜀 𝐃 ( 𝑡 ) 2 0 + 𝜑 ‖ ∇ × 𝐃 ( 𝑡 ) ‖ 2 0 ‖ ‖ ̇ 𝐃 ‖ ‖ + 𝜇 𝜎 𝜑 2 𝐿 2 ( 0 , 𝑡 ; 𝐋 2 ( Ω ) ) ‖ ‖ ̇ ‖ ‖ ⩽ 𝜇 𝜀 𝐃 ( 0 ) 2 0 + ‖ ∇ × 𝐃 ( 0 ) ‖ 2 0 + 𝜇 𝜀 2 𝜎 𝜑 ‖ 𝐆 ‖ 2 𝐿 2 ( 0 , 𝑡 ; 𝐋 2 ( Ω ) ) . ( 4 . 1 2 ) Thirdly, if 𝜎 = 0 , 𝜇 𝜀 2 ‖ ‖ ̇ 𝐃 ‖ ‖ 2 𝐿 ∞ ( 0 , 𝑡 ; 𝐋 2 ( Ω ) ) + 𝜑 ‖ ∇ × 𝐃 ( 𝑡 ) ‖ 2 0 ̇ 𝐃 ⩽ 2 𝜇 𝜀 ‖ ( 0 ) ‖ 2 0 + 2 ‖ ∇ × 𝐃 ( 0 ) ‖ 2 0 + 8 𝜇 𝜀 ‖ 𝐆 ‖ 2 𝐿 1 ( 0 , 𝑡 ; 𝐋 2 ( Ω ) ) . ( 4 . 1 3 )

Proof. We choose ̇ 𝜙 = 2 𝐃 ( 𝑡 ) in (4.4) and integrate over ( 0 , 𝜏 ) to get ̇ 𝜇 𝜀 ‖ 𝐃 ( 𝜏 ) ‖ 2 0 + ‖ ∇ × 𝐃 ( 𝜏 ) ‖ 2 0 − 2 𝜏 0 𝑡 0 𝜑 𝑠 ( ̇ 𝑡 − 𝑠 ) ∇ × 𝐃 ( 𝑠 ) , ∇ × 𝐃 ( 𝑡 ) 𝑑 𝑠 𝑑 𝑡 + 2 𝜇 𝜎 𝜏 0 ̇ 𝐃 ( 𝑡 ) − 𝑡 0 𝜑 𝑠 ̇ ̇ ̇ ( 𝑡 − 𝑠 ) 𝐃 ( 𝑠 ) 𝑑 𝑠 , 𝐃 ( 𝑡 ) 𝑑 𝑡 = 𝜇 𝜀 ‖ 𝐃 ( 0 ) ‖ 2 0 + ‖ ∇ × 𝐃 ( 0 ) ‖ 2 0 + 2 𝜇 𝜀 𝜏 0 ̇ 𝐆 ( 𝑡 ) , 𝐃 ( 𝑡 ) 𝑑 𝑡 . ( 4 . 1 4 ) We now use Lemmas 4.1 and 4.2 and note that the signs of the 𝜑 ′ and 𝜑 integrands mean that only nonnegative terms contribute to the left-hand side. This completes the proof for 𝐆 = 𝟎 . When 𝐆 ≠ 𝟎 and 𝜎 > 0 , Young's inequality, (1.6), gives, 2 𝜇 𝜀 𝜏 0 ̇ 𝐆 ( 𝑡 ) , 𝐃 ( 𝑡 ) 𝑑 𝑡 ⩽ 𝜇 𝜀 2 𝜎 𝜑 ‖ 𝐆 ‖ 2 𝐿 2 ( 0 , 𝑡 ; 𝐋 2 ( Ω ) ) ̇ + 𝜇 𝜎 𝜑 ‖ 𝐃 ‖ 2 𝐿 2 ( 0 , 𝑡 ; 𝐋 2 ( Ω ) ) ( 4 . 1 5 ) and we can complete the proof of the second claim. If 𝜎 = 0 and 𝐆 ≠ 𝟎 , we deduce that ̇ 𝜇 𝜀 ‖ 𝐃 ‖ 2 𝐿 ∞ ( 0 , 𝜏 ; 𝐋 2 ( Ω ) ) + ̆ 𝜑 ‖ ∇ × 𝐃 ( 𝜏 ) ‖ 2 0 ̇ ⩽ 2 𝜇 𝜀 ‖ 𝐃 ( 0 ) ‖ 2 0 + 2 ‖ ∇ × 𝐃 ( 0 ) ‖ 2 0 + 4 𝜇 𝜀 𝜏 0 | | ̇ | | ̇ 𝐆 ( 𝑡 ) , 𝐃 ( 𝑡 ) 𝑑 𝑡 , ⩽ 2 𝜇 𝜀 ‖ 𝐃 ( 0 ) ‖ 2 0 + 2 ‖ ∇ × 𝐃 ( 0 ) ‖ 2 0 + 8 𝜇 𝜀 ‖ 𝐆 ‖ 2 𝐿 1 ( 0 , 𝜏 ; 𝐋 2 ( Ω ) ) + 𝜇 𝜀 2 ‖ ̇ 𝐃 ‖ 2 𝐿 ∞ ( 0 , 𝜏 ; 𝐋 2 ( Ω ) ) , ( 4 . 1 6 ) and the third claim now follows.

4.2. A Formulation for the Electric Field

We now derive an alternative formulation in which the electric field becomes the “unknown”. If the displacement formulation given above could be called the “relaxation problem”, then this could be called the “creep problem”. We again derive ̈ 1 𝐃 + ∇ × 𝜇 ̇ ̇ 𝐉 ∇ × 𝐄 + 𝜎 𝐄 = − 𝑎 , ( 4 . 1 7 ) but now use (3.10) to get ̇ ̇ 𝐃 ( 𝑡 ) = 𝜀 𝐄 ( 𝑡 ) + 𝜀 𝜓 ′ ( 𝑡 ) 𝐄 ( 0 ) − 𝜀 𝑡 0 𝜓 𝑠 ( ̇ 𝑡 − 𝑠 ) 𝐄 ( 𝑠 ) 𝑑 𝑠 , ( 4 . 1 8 ) from which, ̈ ̈ 𝐃 ( 𝑡 ) = 𝜀 𝐄 ( 𝑡 ) + 𝜀 𝜓 ( ̇ 𝑡 ) 𝐄 ( 0 ) + 𝜀 𝜓 ′ ( 0 ) 𝐄 ( 𝑡 ) + 𝜀 𝑡 0 𝜓 𝑠 𝑠 ( ̇ 𝑡 − 𝑠 ) 𝐄 ( 𝑠 ) 𝑑 𝑠 . ( 4 . 1 9 ) Introducing this results in the partial differential Volterra equation, 𝜀 ̈ ̇ 𝐄 ( 𝑡 ) + 𝜀 𝜓 ′ ( 0 ) 𝐄 ( 𝑡 ) + 𝜀 𝑡 0 𝜓 ( ̇ 1 𝑡 − 𝑠 ) 𝐄 ( 𝑠 ) 𝑑 𝑠 + ∇ × 𝜇 ̇ ̇ 𝐉 ∇ × 𝐄 ( 𝑡 ) + 𝜎 𝐄 ( 𝑡 ) = − 𝑎 ( 𝑡 ) − 𝜀 𝜓 ( 𝑡 ) 𝐄 ( 0 ) ( 4 . 2 0 ) or, more usefully, the weak problem: find a twice differentiable map 𝐄 ∶ 𝐼 → 𝑉 such that, 𝜀 ̈ + 𝜇 𝐄 ( 𝑡 ) , 𝜙 − 1 + 𝜎 ̇ + ̇ ∇ × 𝐄 ( 𝑡 ) , ∇ × 𝜙 𝐄 ( 𝑡 ) , 𝜙 𝜀 𝜓 ′ ( 0 ) 𝐄 ( 𝑡 ) + 𝜀 𝑡 0 𝜓 ̇ ( 𝑡 − 𝑠 ) 𝐄 ( 𝑠 ) 𝑑 𝑠 , 𝜙 = ( 𝐅 ( 𝑡 ) , 𝜙 ) ∀ 𝜙 ∈ 𝑉 , ( 4 . 2 1 ) with ̇ 𝐉 𝐅 ∶ = − 𝑎 ( 𝑡 ) − 𝜀 𝜓 ( 𝑡 ) 𝐄 ( 0 ) known and given intial data.

We require a slight adjustment to Lemma 4.1 before giving the stability estimate.

Lemma 4.4 (coercivity of creep). If Assumptions 2.1 hold, then 𝑡 0 𝜓 ′ ( 0 ) 𝐰 ( 𝑠 ) + 𝑠 0 𝜓 ( 𝑠 − 𝜉 ) 𝐰 ( 𝜉 ) 𝑑 𝜉 , 𝐰 ( 𝑠 ) 𝑑 𝑠 ⩾ 𝜓 ′ ( 𝑡 ) ‖ 𝐰 ‖ 2 𝐿 2 ( 0 , 𝑡 ; 𝐋 2 ( Ω ) ) ( 4 . 2 2 ) for all 𝑡 ∈ 𝐼 .

Proof. Since ‖ − 𝜓 ∗ ‖ 𝐰 ‖ 0 ‖ 𝐿 2 ( 0 , 𝑡 ) ⩽ ‖ 𝜓 ‖ 𝐿 1 ( 0 , 𝑡 ) ‖ 𝐰 ‖ 𝐿 2 ( 0 , 𝑡 ; 𝐋 2 ( Ω ) ) and 𝜓 ′ ( 0 ) − ‖ 𝜓 ‖ 𝐿 1 ( 0 , 𝑡 ) ∫ = 𝜓 ′ ( 0 ) + 𝑡 0 𝜓 ( 𝑠 ) 𝑑 𝑠 = 𝜓 ′ ( 𝑡 ) this is similar to Lemma 4.1.

Theorem 4.5 (stability for the 𝐄 equation). If Assumptions 2.1 hold and if 𝜇 > 0 , 𝜀 > 0 and 𝜎 ⩾ 0 are real numbers then for every 𝑡 ∈ 𝐼 we have first, ̇ 𝜇 𝜀 ‖ 𝐄 ( 𝑡 ) ‖ 2 0 + ‖ ∇ × 𝐄 ( 𝑡 ) ‖ 2 0 ̇ + 𝜇 ( 𝜎 + 𝜀 𝜓 ′ ( 𝑡 ) ) ‖ 𝐄 ‖ 2 𝐿 2 ( 0 , 𝑡 ; 𝐋 2 ( Ω ) ) ̇ ⩽ 𝜇 𝜀 ‖ 𝐄 ( 0 ) ‖ 2 0 + ‖ ∇ × 𝐄 ( 0 ) ‖ 2 0 + 𝜇 ‖ 𝐅 ‖ 2 𝐿 2 ( 0 , 𝑡 ; 𝐋 2 ( Ω ) ) , ( 𝜎 + 𝜀 𝜓 ′ ( 𝑡 ) ) ( 4 . 2 3 ) and second, 𝜇 𝜀 2 ‖ ̇ 𝐄 ‖ 2 𝐿 ∞ ( 0 , 𝑡 ; 𝐋 2 ( Ω ) ) + ‖ ∇ × 𝐄 ( 𝑡 ) ‖ 2 0 ̇ + 2 𝜇 ( 𝜎 + 𝜀 𝜓 ′ ( 𝑡 ) ) ‖ 𝐄 ‖ 2 𝐿 2 ( 0 , 𝑡 ; 𝐋 2 ( Ω ) ) ̇ ⩽ 2 𝜇 𝜀 ‖ 𝐄 ( 0 ) ‖ 2 0 + 2 ‖ ∇ × 𝐄 ( 0 ) ‖ 2 0 + 8 𝜇 𝜀 ‖ 𝐅 ‖ 2 𝐿 1 ( 0 , 𝑡 ; 𝐋 2 ( Ω ) ) . ( 4 . 2 4 )

Proof. We choose ̇ 𝜙 = 2 𝐄 ( 𝑡 ) in (4.21), integrate over ( 0 , 𝜏 ) and use Lemma 4.4 to get, ̇ 𝜇 𝜀 ‖ 𝐄 ( 𝜏 ) ‖ 2 0 + ‖ ∇ × 𝐄 ( 𝜏 ) ‖ 2 0 ̇ + 2 𝜇 ( 𝜎 + 𝜀 𝜓 ′ ( 𝜏 ) ) ‖ 𝐄 ‖ 2 𝐿 2 ( 0 , 𝜏 ; 𝐋 2 ( Ω ) ) ̇ ⩽ 𝜇 𝜀 ‖ 𝐄 ( 0 ) ‖ 2 0 + ‖ ∇ × 𝐄 ( 0 ) ‖ 2 0 + 2 𝜇 𝜏 0 ̇ 𝐅 ( 𝑡 ) , 𝐄 ( 𝑡 ) 𝑑 𝑡 . ( 4 . 2 5 ) Now by applying the Cauchy-Schwarz and Young inequalities in the following form: 2 𝜇 𝜏 0 ̇ 𝜇 𝐅 ( 𝑡 ) , 𝐄 ( 𝑡 ) 𝑑 𝑡 ⩽ ( 𝜎 + 𝜀 𝜓 ′ ( 𝜏 ) ) ‖ 𝐅 ‖ 2 𝐿 2 ( 0 , 𝜏 ; 𝐋 2 ( Ω ) ) ̇ + 𝜇 ( 𝜎 + 𝜀 𝜓 ′ ( 𝜏 ) ) ‖ 𝐄 ‖ 2 𝐿 2 ( 0 , 𝜏 ; 𝐋 2 ( Ω ) ) , ( 4 . 2 6 ) we complete the proof of the first claim. For the second, we easily see that ̇ 𝜇 𝜀 ‖ 𝐄 ‖ 2 𝐿 ∞ ( 0 , 𝜏 ; 𝐋 2 ( Ω ) ) + ‖ ∇ × 𝐄 ( 𝜏 ) ‖ 2 0 ̇ + 2 𝜇 ( 𝜎 + 𝜀 𝜓 ′ ( 𝜏 ) ) ‖ 𝐄 ‖ 2 𝐿 2 ( 0 , 𝜏 ; 𝐋 2 ( Ω ) ) ̇ ⩽ 2 𝜇 𝜀 ‖ 𝐄 ( 0 ) ‖ 2 0 + 2 ‖ ∇ × 𝐄 ( 0 ) ‖ 2 0 + 4 𝜇 𝜏 0 | | ̇ | | 𝐅 ( 𝑡 ) , 𝐄 ( 𝑡 ) 𝑑 𝑡 , ( 4 . 2 7 ) and with the Hölder and Young inequalities, 4 𝜇 𝜏 0 | | ̇ | | 𝐅 ( 𝑡 ) , 𝐄 ( 𝑡 ) 𝑑 𝑡 ⩽ 8 𝜇 𝜀 ‖ 𝐅 ‖ 2 𝐿 1 ( 0 , 𝜏 ; 𝐋 2 ( Ω ) ) + 𝜇 𝜀 2 ‖ ̇ 𝐄 ‖ 2 𝐿 ∞ ( 0 , 𝜏 ; 𝐋 2 ( Ω ) ) , ( 4 . 2 8 ) we arrive at the second claim.

Remark 4.6 (growth for 𝜎 = 0 ). If 𝜎 = 0 then the fact that 𝜓 ′ → 0 as 𝑡 → ∞ renders the first stability estimate in the above of limited interest.

We now move on to give a numerical scheme for the electric field equation and see that these stability estimates have an exact analogue.

5. A Numerical Scheme for the Electric Field

In this section, we suggest a fully discrete numerical scheme for the electric field problem, (4.21), provide a stability estimate that is essentially identical to that for the continuous problem, and also give an error bound showing that we still do not require Gronwall-type estimates. To prepare we first assume that the domain is “triangulated” with a quasiuniform mesh of hexahedra or tetrahedra, with meshwidth ℎ , and that a suitable finite element space, 𝑉 ℎ ⊂ 𝑉 , is specified with respect to this partition. (We would probably use Nédélec elements, [30, 31], but we do not yet need to be specific.) The time interval is then discretized into 𝑁 intervals of equal width 𝑘 so that 𝐼 = 𝐼 1 ∪ 𝐼 2 ∪ ⋯ ∪ 𝐼 𝑁 , where 𝐼 𝑖 ∶ = ( 𝑡 𝑖 − 1 , 𝑡 𝑖 ] and 𝑡 𝑖 = 𝑖 𝑘 . We continue to assume that 𝜇 > 0 , 𝜀 > 0 , and 𝜎 ⩾ 0 are real numbers but we recognise that some of the working and results given below could be made simpler by fixing on either of the more definite cases 𝜎 > 0 or 𝜎 = 0 .

Our numerical scheme employs a Galerkin finite element method for spatial discretization and a standard (for the wave equation) two-field centered finite difference approximation in time. For this we define the “velocity” field ̇ 𝐄 𝐖 ∶ = . Before we begin we just mention that numerical schemes for Maxwell's equations (see, e.g., [32, 33]) are often based on the first two of (1.1) and assume that 𝐄 ( 0 ) and 𝐇 ( 0 ) are given and (since 𝜀 and μ are assumed constant here) satisfy the divergence conditions in (1.2). We will need 𝐄 ( 0 ) in the scheme below and also ̇ 𝐄 ( 0 ) and we note that the latter quantity is available if 𝐇 ( 0 ) is because, from (4.3), 𝜀 ̇ ̇ 𝐄 ( 0 ) = 𝐃 ( 0 ) + 𝜑 ′ ( 0 ) 𝐃 ( 0 ) and from (3.10) and (1.1) we see that 𝐃 ( 0 ) = 𝜀 𝐄 ( 0 ) and ̇ 𝐃 ( 0 ) = ∇ × 𝐇 ( 0 ) − 𝜎 𝐄 ( 0 ) − 𝐉 𝑎 ( 0 ) . Based on this, from now on we just assume that we have physically admissible (in terms of (1.2)) intial data such that ̆ ̇ ̆ 𝐄 ( 0 ) ∶ = 𝐄 ∈ 𝑉 , 𝐄 ( 0 ) ∶ = 𝐖 ∈ 𝐋 2 ( Ω ) . ( 5 . 1 ) Our notation is 𝑤 𝑖 ∶ = 𝑤 ( 𝑡 𝑖 ) with 𝜕 𝑡 𝑤 𝑖 𝑤 ∶ = 𝑖 − 𝑤 𝑖 − 1 𝑘 , 𝜕 2 𝑡 𝑤 𝑖 = 𝜕 𝑡 𝜕 𝑡 𝑤 𝑖 = 𝑤 𝑖 − 2 𝑤 𝑖 − 1 + 𝑤 𝑖 − 2 𝑘 2 , 𝑤 𝑖 𝑤 ∶ = 𝑖 + 𝑤 𝑖 − 1 2 , Δ 𝑖 ̇ 𝑤 𝑤 ∶ = 𝑖 + ̇ 𝑤 𝑖 − 1 2 − 𝑤 𝑖 − 𝑤 𝑖 − 1 𝑘 . ( 5 . 2 ) Notice that Δ 𝑖 𝑤 = ̇ 𝑤 𝑖 − 𝜕 𝑡 𝑤 𝑖 and represents the trapezium rule.

The numerical scheme is introduced next in Section 5.1 where we establish stability and well-posedness. The error bound is then given in Section 5.2.

5.1. The Numerical Scheme

Our numerical approximation to the electric field wave equation, (4.21), is the problem: for each 𝑖 = 1 , 2 , … , 𝑁 in turn, find a pair ( 𝐄 ℎ 𝑖 , 𝐖 ℎ 𝑖 ) ∈ 𝑉 ℎ × 𝑉 ℎ such that, 𝜕 𝜇 𝜀 𝑡 𝐖 ℎ 𝑖 + , 𝜙 ∇ × 𝐄 ℎ 𝑖 𝜕 , ∇ × 𝜙 + 𝜇 𝜎 𝑡 𝐄 ℎ 𝑖 𝜓 , 𝜙 + 𝜇 𝜀 ( 0 ) 𝜕 𝑡 𝐄 ℎ 𝑖 + 𝜓 ⋆ 𝜕 𝑡 𝐄 ℎ 𝑖 − 1 / 2 , 𝜙 = 𝜇 𝐅 𝑖 , 𝜙 ∀ 𝜙 ∈ 𝑉 ℎ , 𝜕 ( 5 . 3 ) 𝑡 𝐄 ℎ 𝑖 = , 𝜻 𝐖 ℎ 𝑖 , 𝜻 ∀ 𝜻 ∈ 𝑉 ℎ . ( 5 . 4 ) These are supplemented with discrete intial data, 𝐄 ℎ 0 ≈ ̆ 𝐄 and 𝐖 ℎ 0 ≈ ̆ 𝐖 , about which we will say more later (see later below (5.28)), and the discrete convolution, denoted by “ ⋆ ”, is defined by, 𝜓 ⋆ 𝐗 𝑖 − 1 / 2 1 ∶ = 2 𝑡 𝑖 𝑡 𝑖 − 1 𝜓 𝑡 𝑖 − 𝑠 𝑑 𝑠 𝐗 𝑖 + 1 2 𝑖 − 1 𝑛 = 1 𝑡 𝑛 𝑡 𝑛 − 1 𝜓 𝑡 𝑖 − 𝑠 + 𝜓 𝑡 𝑖 − 1 − 𝑠 𝑑 𝑠 𝐗 𝑛 ( 5 . 5 ) with empty sums set to zero and for any subset { 𝐗 𝑛 } 𝑁 𝑛 = 1 ⊂ 𝐋 2 ( Ω ) . Of course this is merely ( 𝜓 ∗ 𝐗 ) 𝑖 , the average of ( 𝜓 ∗ 𝐗 ) at 𝑡 𝑖 and 𝑡 𝑖 − 1 , for a temporally piecewise constant 𝐗 . Also, for 𝑚 𝑛 ∶ = 𝜓 ′ ( 𝑡 𝑛 − 1 ) − 𝜓 ′ ( 𝑡 𝑛 + 1 ) > 0 (recall Assumptions 2.1), simple calculations give the alternative form 𝜓 ⋆ 𝐗 𝑖 − 1 / 2 1 = − 2 𝜓 ( 0 ) − 𝜓 𝐗 ( 𝑘 ) 𝑖 − 1 2 𝑖 − 1 𝑛 = 1 𝑚 𝑖 − 𝑛 𝐗 𝑛 ( 5 . 6 ) which will be useful below. We next establish a discrete version of (1.7).

Lemma 5.1 (Young's inequality for discrete convolutions). For sets { 𝑚 𝑖 } 𝑁 𝑖 = 1 and { 𝑤 𝑖 } 𝑁 𝑖 = 1 of nonnegative real numbers we have for every 𝑗 ∈ { 2 , 3 , … , 𝑁 } , ⎛ ⎜ ⎜ ⎝ 𝑗 𝑖 = 2 𝑖 − 1 𝑛 = 1 𝑚 𝑖 − 𝑛 𝑤 𝑛 𝑝 ⎞ ⎟ ⎟ ⎠ 1 / 𝑝 ⩽ 𝑗 − 1 𝑛 = 1 𝑚 𝑛 𝑗 − 1 𝑖 = 1 𝑤 𝑝 𝑖 1 / 𝑝 , ( 5 . 7 ) for 1 ⩽ 𝑝 ⩽ ∞ and the obvious interpretation if 𝑝 = ∞ .

Proof. Obviously, ∑ 𝑖 − 1 𝑛 = 1 𝑚 𝑖 − 𝑛 = ∑ 𝑖 − 1 𝑛 = 1 𝑚 𝑛 and so by two Hölder inequalities, 𝑗 𝑖 = 2 𝑖 − 1 𝑛 = 1 𝑚 𝑖 − 𝑛 𝑤 𝑛 𝑝 ⩽ 𝑗 − 1 𝑛 = 1 𝑚 𝑛 𝑗 𝑝 − 1 𝑖 = 2 𝑖 − 1 𝑛 = 1 𝑚 𝑖 − 𝑛 𝑤 𝑝 𝑛 . ( 5 . 8 ) We now interchange the order of summation with 𝑗 𝑖 = 2 𝑖 − 1 𝑛 = 1 𝑚 𝑖 − 𝑛 𝑤 𝑝 𝑛 = 𝑗 − 1 𝑖 = 1 𝑤 𝑝 𝑖 𝑗 − 𝑖 𝑛 = 1 𝑚 𝑛 , ( 5 . 9 ) apply Hölder's inequality again and take 𝑝 th roots.

Our first aim is to derive a stability estimate for the discrete problem and thereby establish its well-posedness. To this end, we need a discrete version of Lemma 4.4 but, before this, we need a bound on the discrete convolution.

Lemma 5.2. If Assumptions 2.1 hold, then | | | | | 𝑗 𝑖 = 1 𝜓 ⋆ 𝐰 𝑖 − 1 / 2 , 𝐰 𝑖 | | | | | ⩽ 𝜓 ( 0 ) − 𝜓 𝑗 𝑗 𝑖 = 1 ‖ ‖ 𝐰 𝑖 ‖ ‖ 2 0 ( 5 . 1 0 ) for every { 𝐰 𝑖 } 𝑁 𝑖 = 1 ⊂ 𝐋 2 ( Ω ) and every 𝑗 ∈ { 1 , 2 , … , 𝑁 } .

Proof. We have from the definition in (5.5) and (5.6) that 𝑗 𝑖 = 1 𝜓 ⋆ 𝐰 𝑖 − 1 / 2 , 𝐰 𝑖 1 = − 2 𝜓 ( 0 ) − 𝜓 ( 𝑘 ) 𝑗 𝑖 = 1 ‖ 𝐰 𝑖 ‖ 2 0 − 1 2 𝑗 𝑖 = 2 𝑖 − 1 𝑛 = 1 𝑚 𝑖 − 𝑛 𝐰 𝑛 , 𝐰 𝑖 , ( 5 . 1 1 ) and so applying Cauchy-Schwarz and then Lemma 5.1 gives | | | | | 𝑗 𝑖 = 1 𝜓 ⋆ 𝐰 𝑖 − 1 / 2 , 𝐰 𝑖 | | | | | ⩽ 1 2 𝜓 ′ ( 0 ) − 𝜓 ′ ( 𝑘 ) + 𝑗 − 1 𝑛 = 1 𝑚 𝑛 𝑗 𝑖 = 1 ‖ 𝐰 𝑖 ‖ 2 0 . ( 5 . 1 2 ) Now observe that ∑ 𝜓 ′ ( 0 ) − 𝜓 ′ ( 𝑘 ) + 𝑗 − 1 𝑛 = 1 𝑚 𝑛 = 2 𝜓 ′ ( 0 ) − 2 𝜓 𝑗 .

We can now state the discrete version of Lemma 4.4.

Lemma 5.3 (coercivity of discrete creep). If Assumptions 2.1 hold, then 𝑗 𝑖 = 1 𝜓 ( 0 ) 𝐰 𝑖 + 𝜓 ⋆ 𝐰 𝑖 − 1 / 2 , 𝐰 𝑖 ⩾ 𝜓 𝑗 𝑗 𝑖 = 1 ‖ 𝐰 𝑖 ‖ 2 0 ( 5 . 1 3 ) for every { 𝐰 𝑖 } 𝑁 𝑖 = 1 ⊂ 𝐋 2 ( Ω ) and every 𝑗 ∈ { 1 , 2 , … , 𝑁 } .

Proof. Noting that 𝑗 𝑖 = 1 𝜓 ( 0 ) 𝐰 𝑖 + 𝜓 ⋆ 𝐰 𝑖 − 1 / 2 , 𝐰 𝑖 ⩾ 𝜓 ′ ( 0 ) 𝑗 𝑖 = 1 ‖ 𝐰 𝑖 ‖ 2 0 − | | | | | 𝑗 𝑖 = 1 𝜓 ⋆ 𝐰 𝑖 − 1 / 2 , 𝐰 𝑖 | | | | | , ( 5 . 1 4 ) the result follows from Lemma 5.2.

Defining the norm ‖ 𝐰 ‖ ∞ , 𝑗 ∶ = m a x ‖ 𝐰 𝑖 ‖ 0 𝐰 ∶ 0 ⩽ 𝑖 ⩽ 𝑗 ∀ 𝑗 ⩽ 𝑁 , ∀ 𝑖 𝑁 𝑖 = 0 ⊂ 𝐋 2 ( Ω ) , ( 5 . 1 5 ) we can now state the main result of this subsection.

Theorem 5.4 (discrete well-posedness). If Assumptions 2.1 hold and we are given approximations ( 𝐄 ℎ 0 , 𝐖 ℎ 0 ) ∈ V ℎ × 𝑉 ℎ to ( ̆ ̆ 𝐄 , 𝐖 ) then there exists a unique solution { ( 𝐄 ℎ 𝑖 , 𝐖 ℎ 𝑖 ) } 𝑁 𝑖 = 1 ⊂ 𝑉 ℎ × 𝑉 ℎ to the discrete problem (5.3) and (5.4). Moreover, 𝜇 𝜀 2 ‖ 𝐖 ℎ ‖ 2 ∞ , 𝑗 + ‖ ∇ × 𝐄 ℎ 𝑗 ‖ 2 0 + 2 𝑘 𝜇 𝜎 + 𝜀 𝜓 𝑗 𝑗 𝑖 = 1 ‖ 𝜕 𝑡 𝐄 ℎ 𝑖 ‖ 2 0 ⩽ 2 𝜇 𝜀 ‖ 𝐖 ℎ 0 ‖ 2 0 + 2 ‖ ∇ × 𝐄 ℎ 0 ‖ 2 0 + 8 𝜇 𝜀 𝑘 𝑗 𝑖 = 1 ‖ 𝐅 𝑖 ‖ 0 2 ( 5 . 1 6 ) for each 𝑗 ∈ { 1 , 2 , … , 𝑁 } .

Proof. Since uniqueness implies existence for this finite dimensional problem, we need only prove the estimate and then take 𝐄 ℎ 0 = 𝐖 ℎ 0 = 𝐅 1 = ⋯ = 𝟎 . We select 𝜙 = 2 𝑘 𝜕 𝑡 𝐄 ℎ 𝑖 in (5.3) and 𝜻 = 2 𝑘 𝜕 𝑡 𝐖 ℎ 𝑖 in (5.4), merge them using the common term 2 𝑘 ( 𝜕 𝑡 𝐖 ℎ 𝑖 , 𝜕 𝑡 𝐄 ℎ 𝑖 ) , and then sum the result over 𝑖 = 1 , … , 𝑗 . Noting the usual telescoping sum that arises from products like 𝐰 𝑖 𝜕 𝑡 𝐰 𝑖 , we then use Lemma 5.3 to obtain 𝜇 𝜀 ‖ 𝐖 ℎ 𝑗 ‖ 2 0 + ‖ ∇ × 𝐄 ℎ 𝑗 ‖ 2 0 + 2 𝑘 𝜇 𝜎 + 𝜀 𝜓 𝑗 𝑗 𝑖 = 1 ‖ 𝜕 𝑡 𝐄 ℎ 𝑖 ‖ 2 0 ⩽ 𝜇 𝜀 ‖ 𝐖 ℎ 0 ‖ 2 0 + ‖ ∇ × 𝐄 ℎ 0 ‖ 2 0 + 2 𝑘 𝜇 𝑗 𝑖 = 1 𝐅 𝑖 , 𝜕 𝑡 𝐄 ℎ 𝑖 . ( 5 . 1 7 ) It is clear from (5.4) that 𝜕 𝑡 𝐄 ℎ 𝑖 = 𝐖 ℎ 𝑖 and so if 𝑖 ⩽ 𝑗 we easily conclude that ‖ 𝜕 𝑡 𝐄 ℎ 𝑖 ‖ 0 ⩽ ‖ 𝐖 ℎ ‖ ∞ , 𝑗 . Doubling the right-hand side to account for the ‖ ⋅ ‖ ∞ , 𝑗 norm a Young's inequality, | | | | | 4 𝑘 𝜇 𝑗 𝑖 = 1 𝐅 𝑖 , 𝜕 𝑡 𝐄 ℎ 𝑖 | | | | | ⩽ 𝜇 𝜀 2 ‖ 𝐖 ‖ 2 ∞ , 𝑗 + 8 𝜇 𝜀 𝑘 𝑗 𝑖 = 1 ‖ 𝐅 𝑖 ‖ 0 2 , ( 5 . 1 8 ) this completes the proof.

The stability estimate given above is an exact analogue of the second estimate given for the continuous problem in Theorem 4.5. We now examine the possibility of deriving an error estimate with a similar control over all of the “constants”.

5.2. An Error Bound

To arrive at an error bound for the discrete approximation, we introduce the 𝐇 ( 𝐜 𝐮 𝐫 𝐥 ; Ω ) projection as the linear continuous map, ℛ ∶ 𝑉 → 𝑉 ℎ such that, for every 𝐮 ∈ 𝑉 , ℛ 𝐮 ∈ 𝑉 ℎ is the unique solution to ( ℛ 𝐮 , 𝐯 ) 𝐇 ( 𝐜 𝐮 𝐫 𝐥 ; Ω ) = ( 𝐮 , 𝐯 ) 𝐇 ( 𝐜 𝐮 𝐫 𝐥 ; Ω ) ∀ 𝐯 ∈ 𝑉 ℎ . ( 5 . 1 9 ) It is clear that the best approximation property, ‖ 𝐮 − ℛ 𝐮 ‖ 𝐇 ( 𝐜 𝐮 𝐫 𝐥 ; Ω ) ⩽ ‖ 𝐮 − 𝐯 ‖ 𝐇 ( 𝐜 𝐮 𝐫 𝐥 ; Ω ) ∀ 𝐯 ∈ 𝑉 ℎ ( 5 . 2 0 ) is satisfied. Introducing also the 𝐋 2 ( Ω ) projection, 𝐱 𝒫 ∶ 𝐋 2 ( Ω ) → 𝑉 ℎ , through ( 𝒫 𝐮 , 𝐯 ) = ( 𝐮 , 𝐯 ) for all 𝐯 ∈ 𝑉 ℎ , we then (as usual) split the errors according to 𝜽 𝑖 = 𝐄 ℎ 𝑖 𝑡 − ℛ 𝐄 𝑖 ∈ 𝑉 ℎ , 𝝀 𝑖 = 𝐖 ℎ 𝑖 𝑡 − 𝒫 𝐖 𝑖 ∈ 𝑉 ℎ , 𝚯 ( 𝑡 ) = 𝐄 ( 𝑡 ) − ℛ 𝐄 ( 𝑡 ) ∈ 𝑉 , 𝚲 ( 𝑡 ) = 𝐖 ( 𝑡 ) − 𝒫 𝐖 ( 𝑡 ) ∈ 𝑉 , ( 5 . 2 1 ) so that 𝐄 ℎ 𝑖 − 𝐄 ( 𝑡 𝑖 ) = 𝜽 𝑖 − Θ 𝑖 and 𝐖 ℎ 𝑖 − 𝐖 ( 𝑡 𝑖 ) = 𝝀 𝑖 − Λ 𝑖 along with ( 𝚲 ( 𝑡 ) , 𝐰 ) = 0 , ( ∇ × 𝚯 ( 𝑡 ) , ∇ × 𝐰 ) = − ( 𝚯 ( 𝑡 ) , 𝐰 ) ∀ 𝐰 ∈ 𝑉 ℎ . ( 5 . 2 2 ) The first result is a representation formula for the spatially discrete components of the error and for this we define, 𝒢 𝑖 ∶ = 𝜇 𝜀 Δ 𝑖 𝐖 + 𝜇 𝜎 + 𝜀 𝜓 ( Δ 0 ) 𝑖 𝐄 + 𝜕 𝑡 𝚯 𝑖 − 𝚯 𝑖 𝜓 + 𝜇 𝜀 ⋆ 𝜕 𝑡 𝚯 𝑖 − 1 / 2 + 𝜇 𝜀 𝜓 ∗ ̇ 𝐄 − 𝜕 𝑡 𝐄 𝑖 , ℋ ( 5 . 2 3 ) 𝑖 ∶ = 𝜇 𝜀 Δ 𝑖 𝐄 + 𝜇 𝜀 𝜕 𝑡 𝚯 𝑖 . ( 5 . 2 4 )

Lemma 5.5. For each 𝑗 ∈ { 1 , 2 , … , 𝑁 } , we have 2 𝑘 𝜇 𝜀 𝑗 𝑖 = 1 𝜓 ′ ( 0 ) 𝜕 𝑡 𝜽 𝑖 + 𝜓 ⋆ 𝜕 𝑡 𝜽 𝑖 − 1 / 2 , 𝜕 𝑡 𝜽 𝑖 + 𝜇 𝜀 ‖ 𝝀 𝑗 ‖ 2 0 + ‖ ∇ × 𝜽 𝑗 ‖ 2 0 + 2 𝑘 𝜇 𝜎 𝑗 𝑖 = 1 ‖ 𝜕 𝑡 𝜽 𝑖 ‖ 2 0 = 𝜇 𝜀 ‖ 𝝀 0 ‖ 2 0 + ‖ ∇ × 𝜽 0 ‖ 2 0 + 2 𝑘 𝑗 𝑖 = 1 𝒢 𝑖 , 𝜕 𝑡 𝜽 𝑖 − ℋ 𝑖 , 𝜕 𝑡 𝝀 𝑖 . ( 5 . 2 5 )

Proof. We average the continuous problem, (4.21), at the times 𝑡 𝑖 and 𝑡 𝑖 − 1 and subtract the result from the discrete problem (5.3). Simultaneously adding and subtracting like terms in ℛ 𝐄 and 𝒫 𝐖 introduces the error components, 𝜽 , 𝝀 , 𝜽 , and Λ , and we collect the lower case quantities to the left of the equals sign and all of the other terms to the right. We also notice that the term ( 𝜕 𝑡 Λ 𝑖 , 𝜙 ) is zero while ( ∇ × Θ 𝑖 , ∇ × 𝜙 ) can be replaced by − ( Θ 𝑖 , 𝜙 ) . A simple manipulation provides 𝜓 ∗ ̇ 𝐄 𝑖 − 𝜓 ⋆ 𝜕 𝑡 ℛ 𝐄 𝑖 − 1 / 2 = 𝜓 ⋆ 𝜕 𝑡 𝚯 𝑖 − 1 / 2 + 𝜓 ∗ ̇ 𝐄 − 𝜕 𝑡 𝐄 𝑖 ( 5 . 2 6 ) and so choosing 𝜙 = 2 𝑘 𝜕 𝑡 𝜽 𝑖 generates the 𝒢 𝑖 terms on the right and with 𝜻 = 2 𝑘 𝜇 𝜀 𝜕 𝑡 𝝀 𝑖 we can use (5.4), with ( Λ 𝑖 , 𝜻 ) = 0 , to get 2 𝑘 𝜇 𝜀 ( 𝜕 𝑡 𝜽 𝑖 , 𝜕 𝑡 𝝀 𝑖 ) = 𝑘 𝜇 𝜀 𝜕 𝑡 ‖ 𝝀 𝑖 ‖ 2 0 + 2 𝑘 ( ℋ 𝑖 , 𝜕 𝑡 𝝀 𝑖 ) . Lastly we merge these and sum over 𝑖 = 1 , 2 , … , 𝑗 .

The first main effort in deriving the error bound is contained in the next lemma. For this we note first the “summation by parts” formula, 𝑘 𝑗 𝑖 = 1 𝑤 𝑖 , 𝜕 𝑡 𝑣 𝑖 = 𝑤 𝑗 , 𝑣 𝑗 − 𝑤 1 , 𝑣 0 − 𝑘 𝑗 𝑖 = 2 𝜕 𝑡 𝑤 𝑖 , 𝑣 𝑖 − 1 ( 5 . 2 7 ) and, secondly, that on taking 𝜻 = 𝜇 𝜀 𝜕 𝑡 𝜽 𝑖 we can use (5.4) to get, 𝜇 𝜀 ‖ 𝜕 𝑡 𝜽 𝑖 ‖ 0 ⩽ 𝜇 𝜀 ‖ 𝝀 ‖ ∞ , 𝑖 + ‖ ℋ 𝑖 ‖ 0 . ( 5 . 2 8 ) Thirdly, if the discrete intial data are defined by 𝐄 ℎ 0 ̆ 𝐄 ∶ = ℛ and 𝐖 ℎ 0 ̆ 𝐖 ∶ = 𝒫 , then it follows that 𝜽 0 = 𝝀 0 = 𝟎 .

We can now give a formal a priori error bound—formal because it deals only with the time-stepping dynamics of the error rather than with convergence orders. Of course this bound is meaningless until we know that the unresolved terms can be controlled but, in this way, by comparing with the stability estimates in Theorems 4.5 and 5.4 we can see that the error is governed (bounded) by essentially the same temporal dynamics as both the exact and discrete solutions. The “concrete form” of the error bound will follow later in Corollary 5.9.

Theorem 5.6 (“formal” a priori error bound). If Assumptions 2.1 hold and the discrete intial data are determined by 𝐄 ℎ 0 ̆ 𝐄 ∶ = ℛ and 𝐖 ℎ 0 ̆ 𝐖 ∶ = 𝒫 , we have 𝜇 𝜀 2 ‖ ‖ 𝐖 𝑗 − 𝐖 ℎ 𝑗 ‖ ‖ 2 0 + ‖ ‖ 𝐄 ∇ × 𝑗 − 𝐄 ℎ 𝑗 ‖ ‖ 2 0 + 2 𝑘 𝜇 𝜎 + 𝜀 𝜓 𝑗 𝑗 𝑖 = 1 ‖ ‖ 𝜕 𝑡 ( 𝐄 𝑖 − 𝐄 ℎ 𝑖 ) ‖ ‖ 2 0 ⩽ 𝜇 𝜀 ‖ 𝚲 𝑗 ‖ 2 0 + 2 ‖ 𝚯 𝑗 ‖ 2 𝐇 ( 𝐜 𝐮 𝐫 𝐥 ; Ω ) + 4 𝑘 𝜇 𝜎 + 𝜀 𝜓 𝑗 𝑗 𝑖 = 1 ‖ ‖ 𝜕 𝑡 𝚯 𝑖 ‖ ‖ 2 0 + 4 8 ‖ ‖ ℋ 𝜇 𝜀 𝑗 ‖ ‖ 2 0 + 4 8 𝑘 𝜇 𝜀 𝑗 𝑖 = 2 ‖ ‖ 𝜕 𝑡 ℋ 𝑖 ‖ ‖ 0 2 + 4 8 𝑘 𝜇 𝜀 𝑗 𝑖 = 1 ‖ ‖ 𝒢 𝑖 ‖ ‖ 0 2 + 8 𝑘 𝜇 𝜖 𝑗 𝑖 = 1 ‖ ‖ 𝒢 𝑖 ‖ ‖ 0 ‖ ‖ ℋ 𝑖 ‖ ‖ 0 ( 5 . 2 9 ) for every 𝑗 ∈ { 1 , 2 , … , 𝑁 } , where the 𝒢 𝑖 and ℋ 𝑖 are given in (5.23) and (5.24).

Proof. Lemmas 5.5 and 5.3 give 𝜇 𝜀 ‖ 𝝀 ‖ 2 ∞ , 𝑗 + ‖ ‖ ∇ × 𝜽 𝑗 ‖ ‖ 2 0 + 2 𝑘 𝜇 𝜎 + 𝜀 𝜓 𝑗 𝑗 𝑖 = 1 ‖ 𝜕 𝑡 𝜽 𝑖 ‖ 2 0 ⩽ | | | | | 4 𝑘 𝑗 𝑖 = 1 𝒢 𝑖 , 𝜕 𝑡 𝜽 𝑖 − ℋ 𝑖 , 𝜕 𝑡 𝝀 𝑖 | | | | | ( 5 . 3 0 ) and using (5.27) and Young's inequality then gives, | | | | | 4 𝑘 𝑗 𝑖 = 1 ℋ 𝑖 , 𝜕 𝑡 𝝀 𝑖 | | | | | ⩽ 4 𝜇 𝜀 𝜖 ‖ 𝝀 ‖ 2 ∞ , 𝑗 + 2 𝜖 𝜇 𝜀 ‖ ℋ 𝑗 ‖ 2 0 + 2 𝜖 𝑘 𝜇 𝜀 𝑗 𝑖 = 2 ‖ 𝜕 𝑡 ℋ 𝑖 ‖ 0 2 ( 5 . 3 1 ) for all 𝜖 > 0 . From (5.28), we find | | | | | 4 𝑘 𝑗 𝑖 = 1 𝒢 𝑖 , 𝜕 𝑡 𝜽 𝑖 | | | | | ⩽ 4 𝑘 𝜇 𝜀 𝑗 𝑖 = 1 ‖ 𝒢 𝑖 ‖ 0 ‖ ℋ 𝑖 ‖ 0 + 2 𝜇 𝜀 𝜖 ‖ 𝝀 ‖ 2 ∞ , 𝑗 + 2 𝜖 𝑘 𝜇 𝜀 𝑗 𝑖 = 1 ‖ 𝒢 𝑖 ‖ 0 2 ( 5 . 3 2 ) for all 𝜖 > 0 . In fact, choosing 𝜖 = 1 2 in the last two estimates gives 𝜇 𝜀 2 ‖ 𝝀 ‖ 2 ∞ , 𝑗 + ‖ ∇ × 𝜽 𝑗 ‖ 2 0 + 2 𝑘 𝜇 𝜎 + 𝜀 𝜓 𝑗 𝑗 𝑖 = 1 ‖ 𝜕 𝑡 𝜽 𝑖 ‖ 2 0 ⩽ 2 4 𝜇 𝜀 ‖ ℋ 𝑗 ‖ 2 0 + 2 4 𝑘 𝜇 𝜀 𝑗 𝑖 = 2 ‖ 𝜕 𝑡 ℋ 𝑖 ‖ 0 2 + 2 4 𝑘 𝜇 𝜀 𝑗 𝑖 = 1 ‖ 𝒢 𝑖 ‖ 0 2 + 4 𝑘 𝜇 𝜀 𝑗 𝑖 = 1 ‖ 𝒢 𝑖 ‖ 0 ‖ ℋ 𝑖 ‖ 0 , ( 5 . 3 3 ) and so by applying the triangle inequality in the form ‖ 𝑎 ± 𝑏 ‖ 2 ⩽ 2 ‖ 𝑎 ‖ 2 + 2 ‖ 𝑏 ‖ 2 , we complete the proof.

To turn this estimate into something more concrete, and to see what the convergence rate is, as well as what the potential growth in time of the error is, we need to estimate all of the terms on the right-hand side and for this we need some preparation. The following arises in a standard way and will be used without proof. ‖ ̇ 𝐯 ( 𝑠 ) − 𝜕 𝑡 𝐯 𝑖 ‖ 0 ̈ ⩽ 𝑘 ‖ 𝐯 ‖ 𝐿 1 ( 𝑡 𝑖 − 1 , 𝑡 𝑖 ; 𝐋 2 ( Ω ) ) 𝑡 , f o r 𝑠 ∈ 𝑖 − 1 , 𝑡 𝑖 , ( 5 . 3 4 ) ‖ Δ 𝑖 𝐯 ‖ 0 ⩽ 𝑘 2 ‖ 𝐯 𝑡 𝑡 𝑡 ‖ 𝐿 1 ( 𝑡 𝑖 − 1 , 𝑡 𝑖 ; 𝐋 2 ( Ω ) ) , ( 5 . 3 5 ) ‖ 𝜕 𝑡 𝐯 𝑖 ‖ 0 ⩽ 𝑘 − 1 ‖ ̇ 𝐯 ‖ 𝐿 1 ( 𝑡 𝑖 − 1 , 𝑡 𝑖 ; 𝐋 2 ( Ω ) ) , ( 5 . 3 6 ) ‖ 𝜕 𝑡 Δ 𝑖 𝐯 ‖ 0 ⩽ 𝑘 ‖ 𝐯 𝑡 𝑡 𝑡 𝑡 ‖ 𝐿 1 ( 𝑡 𝑖 − 2 , 𝑡 𝑖 ; 𝐋 2 ( Ω ) ) . ( 5 . 3 7 ) We will also need the following estimates for the continuous and discrete convolutions.

Lemma 5.7 (error in convolutions). If Assumptions 2.1 hold then we have 𝑘 𝑗 𝑖 = 1 ‖ ‖ 𝜓 ⋆ 𝜕 𝑡 𝐯 𝑖 − 1 / 2 ‖ ‖ 0 ⩽ 𝜓 ( 0 ) − 𝜓 𝑗 ‖ 𝐯 𝑡 ‖ 𝐿 1 ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) , 𝑘 𝑗 𝑖 = 1 ‖ ‖ 𝜓 ∗ 𝐯 𝑡 − 𝜕 𝑡 𝐯 𝑖 ‖ ‖ 0 ⩽ 𝜓 ( 0 ) − 𝜓 𝑗 𝑘 2 ‖ 𝐯 𝑡 𝑡 ‖ 𝐿 1 ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) , ( 5 . 3 8 ) for every 𝑖 ∈ { 1 , 2 , … , 𝑗 } and every 𝑗 ∈ { 1 , 2 , … , 𝑁 } .

Proof. First, from the definition in (5.6), we get after using (5.36) that 𝜓 𝑘 ‖ ⋆ 𝜕 𝑡 𝐯 𝑖 − 1 / 2 ‖ 0 ⩽ 1 2 ( 𝜓 ′ ( 0 ) − 𝜓 ′ ( 𝑘 ) ) ‖ 𝐯 𝑡 ‖ 𝐿 1 ( 𝑡 𝑖 − 1 , 𝑡 𝑖 ; 𝐋 2 ( Ω ) ) + 1 2 𝑖 − 1 𝑛 = 1 𝑚 𝑖 − 𝑛 ‖ 𝐯 𝑡 ‖ 𝐿 1 ( 𝑡 𝑛 − 1 , 𝑡 𝑛 ; 𝐋 2 ( Ω ) ) , ( 5 . 3 9 ) and then summing over 𝑖 = 1 , … , 𝑗 , we use Lemma 5.1 and arrive at 𝑘 𝑗 𝑖 = 1 ‖ 𝜓 ⋆ 𝜕 𝑡 𝐯 𝑖 − 1 / 2 ‖ 0 ⩽ 1 2 𝜓 ′ ( 0 ) − 𝜓 ′ ( 𝑘 ) + 𝑗 − 1 𝑛 = 1 𝑚 𝑛 ‖ 𝐯 𝑡 ‖ 𝐿 1 ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) ( 5 . 4 0 ) which is simplified to give the first result. The proof of the second is similar and begins with, 𝜓 ∗ 𝐯 𝑡 − 𝜕 𝑡 𝐯 𝑖 = 1 2 𝑡 𝑖 𝑡 𝑖 − 1 𝜓 𝑡 𝑖 ̇ − 𝑠 𝐯 ( 𝑠 ) − 𝜕 𝑡 𝐯 𝑖 + 1 𝑑 𝑠 2 𝑖 − 1 𝑛 = 1 𝑡 𝑛 𝑡 𝑛 − 1 𝜓 𝑡 𝑖 − 𝑠 + 𝜓 𝑡 𝑖 − 1 ̇ − 𝑠 𝐯 ( 𝑠 ) − 𝜕 𝑡 𝐯 𝑛 𝑑 𝑠 . ( 5 . 4 1 ) From (5.34) we define 𝑤 𝑛 ∶ = 𝑘 ‖ 𝐯 𝑡 𝑡 ‖ 𝐿 1 ( 𝑡 𝑛 − 1 , 𝑡 𝑛 ; 𝐋 2 ( Ω ) ) ̇ ⩾ ‖ 𝐯 ( 𝑠 ) − 𝜕 𝑡 𝐯 𝑛 ‖ 0 for 𝑠 ∈ [ 𝑡 𝑛 − 1 , 𝑡 𝑛 ] , and then after a little work, we get ‖ 𝜓 ∗ 𝐯 𝑡 − 𝜕 𝑡 𝐯 𝑖 ‖ 0 ⩽ 1 2 𝜓 ( 0 ) − 𝜓 𝑤 ( 𝑘 ) 𝑖 + 1 2 𝑖 − 1 𝑛 = 1 𝑚 𝑖 − 𝑛 𝑤 𝑛 . ( 5 . 4 2 ) Again we sum over 𝑖 and apply Lemma 5.1 to arrive at, 𝑘 𝑗 𝑖 = 1 ‖ 𝜓 ∗ 𝐯 𝑡 − 𝜕 𝑡 𝐯 𝑖 ‖ 0 ⩽ 𝑘 2 2 𝜓 ′ ( 0 ) − 𝜓 ′ ( 𝑘 ) + 𝑗 − 1 𝑛 = 1 𝑚 𝑛 ‖ 𝐯 𝑡 𝑡 ‖ 𝐿 1 ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) . ( 5 . 4 3 ) The same simplification completes the proof.

Corollary 5.9 contains the “concrete” form of the error bound but before we get to it we need to estimate the terms on the right of Theorem 5.6. From now on, we are not so concerned with the precise form of constants that will appear and so we adopt the common practice of writing “ 𝑎 ≲ 𝑏 ” to mean that “ 𝑎 ⩽ 𝐶 𝑏 ” where 𝐶 is a generic constant. In our working 𝐶 will depend only on 𝜇 , 𝜀 , 𝜓 ′ ( 0 ) > 0 and 𝜎 ⩾ 0 (but not on 𝑇 or 𝐄 ).

Lemma 5.8 (auxiliary bounds for 𝒢 𝑖 and ℋ 𝑖 ). If Assumptions 2.1 hold then for every 𝑗 ∈ { 1 , 2 , … , 𝑁 } , 4 8 𝜇 𝜀 ‖ ℋ 𝑗 ‖ 2 0 + 4 8 𝑘 𝜇 𝜀 𝑗 𝑖 = 2 ‖ 𝜕 𝑡 ℋ 𝑖 ‖ 0 2 + 4 8 𝑘 𝜇 𝜀 𝑗 𝑖 = 1 ‖ 𝒢 𝑖 ‖ 0 2 + 8 𝑘 𝜇 𝜀 𝑗 𝑖 = 1 ‖ 𝒢 𝑖 ‖ 0 ‖ ℋ 𝑖 ‖ 0 ≲ 𝑘 4 ‖ 𝐄 𝑡 𝑡 𝑡 ‖ 2 𝐿 ∞ ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) + ‖ 𝐄 𝑡 𝑡 ‖ 2 𝑊 2 1 ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) + 𝑗 𝑖 = 1 𝑘 ‖ 𝚯 𝑖 ‖ 0 2 + ‖ 𝚯 𝑡 ‖ 2 𝐿 ∞ ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) + ‖ 𝚯 𝑡 ‖ 2 𝑊 1 1 ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) . ( 5 . 4 4 )

Proof. First, from (5.24), (5.35), and (5.36) ‖ ℋ 𝑖 ‖ 0 ⩽ 𝜇 𝜀 𝑘 2 ‖ 𝐄 𝑡 𝑡 𝑡 ‖ 𝐿 1 ( 𝑡 𝑖 − 1 , 𝑡 𝑖 ; 𝐋 2 ( Ω ) ) + 𝜇 𝜀 𝑘 ‖ ̇ 𝚯 ‖ 𝐿 1 ( 𝑡 𝑖 − 1 , 𝑡 𝑖 ; 𝐋 2 ( Ω ) ) , ( 5 . 4 5 ) and with a Young's inequality, it follows from this that 4 8 𝜇 𝜀 ‖ ℋ 𝑗 ‖ 2 0 + 8 𝑘 𝜇 𝜀 𝑗 𝑖 = 1 ‖ 𝒢 𝑖 ‖ 0 ‖ ℋ 𝑖 ‖ 0 ≲ 𝜇 𝜀 𝑘 4 ‖ 𝐄 𝑡 𝑡 𝑡 ‖ 2 𝐿 ∞ ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) + 𝑘 2 ‖ 𝐄 𝑡 𝑡 𝑡 ‖ 𝐿 ∞ ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) + ‖ 𝚯 𝑡 ‖ 𝐿 ∞ ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) 2 + 𝑗 𝑖 = 1 𝑘 ‖ 𝒢 𝑖 ‖ 0 2 + 𝜇 𝜀 ‖ 𝚯 𝑡 ‖ 2 𝐿 ∞ ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) , ( 5 . 4 6 ) where the hidden constant is merely a positive real number. Now, with the observation 𝜕 2 𝑡 𝚯 𝑖 = 1 𝑘 𝜕 𝑡 𝑡 𝑖 𝑡 𝑖 − 1 ̇ 1 𝚯 ( 𝑠 ) 𝑑 𝑠 = 𝑘 2 𝑡 𝑖 𝑡 𝑖 − 1 ̇ ̇ 1 𝚯 ( 𝑠 ) − 𝚯 ( 𝑠 − 𝑘 ) 𝑑 𝑠 = 𝑘 2 𝑡 𝑖 𝑡 𝑖 − 1 𝑠 𝑠 − 𝑘 ̈ 𝚯 ( 𝜂 ) 𝑑 𝜂 𝑑 𝑠 ( 5 . 4 7 ) we see that ‖ 𝜕 2 𝑡 𝚯 𝑖 ‖ 0 ⩽ 1 𝑘 𝑡 𝑖 𝑡 𝑖 − 2 ‖ ̈ 𝚯 ( 𝜂 ) ‖ 0 𝑑 𝜂 . ( 5 . 4 8 ) Using this with (5.24) and (5.37) then gives ‖ 𝜕 𝑡 ℋ 𝑖 ‖ 0 ⩽ 𝜇 𝜀 𝑘 ‖ 𝐄 𝑡 𝑡 𝑡 𝑡 ‖ 𝐿 1 ( 𝑡 𝑖 − 2 , 𝑡 𝑖 ; 𝐋 2 ( Ω ) ) + 𝜇 𝜀 𝑘 ‖ 𝚯 𝑡 𝑡 ‖ 𝐿 1 ( 𝑡 𝑖 − 2 ; 𝑡 𝑖 ; 𝐋 2 ( Ω ) ) ( 5 . 4 9 ) which, in turn, produces 4 8 𝑘 𝜇 𝜀 𝑗 𝑖 = 2 ‖ 𝜕 𝑡 ℋ 𝑖 ‖ 0 2 𝑘 ⩽ 3 8 4 𝜇 𝜀 4 ‖ 𝐄 𝑡 𝑡 𝑡 𝑡 ‖ 2 𝐿 1 ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) + ‖ 𝚯 𝑡 𝑡 ‖ 2 𝐿 1 ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) . ( 5 . 5 0 ) The next task is to estimate the terms in ‖ 𝒢 𝑖 ‖ 0 . To begin we get from (5.23), (5.35), and (5.36) that ‖ 𝒢 𝑖 ‖ 0 ⩽ 𝜇 𝜀 ‖ Δ 𝑖 𝐖 ‖ 0 + 𝜇 𝜎 + 𝜀 𝜓 ( 0 ) ‖ Δ 𝑖 𝐄 ‖ 0 + ‖ 𝜕 𝑡 𝚯 𝑖 ‖ 0 + ‖ 𝚯 𝑖 ‖ 0 𝜓 + 𝜇 𝜀 ‖ ⋆ 𝜕 𝑡 𝚯 𝑖 − 1 / 2 ‖ 0 + 𝜇 𝜀 ‖ 𝜓 ∗ 𝐄 𝑡 − 𝜕 𝑡 𝐄 𝑖 ‖ 0 , ⩽ 𝜇 𝜀 𝑘 2 ‖ 𝐖 𝑡 𝑡 𝑡 ‖ 𝐿 1 ( 𝑡 𝑖 − 1 , 𝑡 𝑖 ; 𝐋 2 ( Ω ) ) + 𝜇 ( 𝜎 + 𝜀 𝜓 ′ ( 0 ) ) 𝑘 2 ‖ 𝐄 𝑡 𝑡 𝑡 ‖ 𝐿 1 ( 𝑡 𝑖 − 1 , 𝑡 𝑖 ; 𝐋 2 ( Ω ) ) + 𝜇 ( 𝜎 + 𝜀 𝜓 ′ ( 0 ) ) 𝑘 ‖ ̇ 𝚯 𝑖 ‖ 𝐿 1 ( 𝑡 𝑖 − 1 , 𝑡 𝑖 ; 𝐋 2 ( Ω ) ) + ‖ 𝚯 𝑖 ‖ 0 𝜓 + 𝜇 𝜀 ‖ ⋆ 𝜕 𝑡 𝚯 𝑖 − 1 / 2 ‖ 0 + 𝜇 𝜀 ‖ 𝜓 ∗ 𝐄 𝑡 − 𝜕 𝑡 𝐄 𝑖 ‖ 0 . ( 5 . 5 1 ) Introducing the results from Lemma 5.7 into this provides, 𝑘 𝑗 𝑖 = 1 ‖ 𝒢 𝑖 ‖ 0 ≲ 𝑘 2 ‖ 𝐄 𝑡 𝑡 ‖ 𝑊 2 1 ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) + ‖ 𝚯 𝑡 ‖ 𝐿 1 ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) + 𝑗 𝑖 = 1 𝑘 ‖ 𝚯 𝑖 ‖ 0 , ( 5 . 5 2 ) where 𝐶 depends only upon 𝜇 , 𝜀 , 𝜓 ′ ( 0 ) > 0 and 𝜎 ⩾ 0 but not on 𝑇 or 𝐄 . Squaring this yields, 𝑘 𝑗 𝑖 = 1 ‖ 𝒢 𝑖 ‖ 0 2 ≲ 𝑘 4 ‖ 𝐄 𝑡 𝑡 ‖ 2 𝑊 2 1 ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) + ‖ 𝚯 𝑡 ‖ 2 𝐿 1 ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) + 𝑗 𝑖 = 1 𝑘 ‖ 𝚯 𝑖 ‖ 0 2 . ( 5 . 5 3 ) The desired result then follows by merging (5.46), (5.50), (5.52) and (5.53).

The final step is to specify the finite element space, 𝑉 ℎ , and its approximation properties and for this we follow Monk in [34, Theorem 2.1 (part 1)] where, for the ℓ -th order first-type (since we are not aiming for optimal 𝐋 2 ( Ω ) bounds) curl-conforming Nédélec spaces in [30], we find the following approximation result. Let the mesh be tetrahedral and regular; then if 𝐄 ∈ 𝐇 ℓ + 1 ( Ω ) there is an interpolator Π ∶ 𝑉 → 𝑉 ℎ such that ‖ 𝐄 − Π 𝐄 ‖ 0 + ‖ ∇ × ( 𝐄 − Π 𝐄 ) ‖ 0 ⩽ 𝐶 ℎ ℓ ‖ 𝐄 ‖ 𝐇 ℓ + 1 ( Ω ) . ( 5 . 5 4 ) Hence, in view of (5.20) we have ‖ ‖ ‖ 𝜕 𝑛 𝚯 𝜕 𝑡 𝑛 ‖ ‖ ‖ 0 + ‖ ‖ ‖ 𝜕 ∇ × 𝑛 𝚯 𝜕 𝑡 𝑛 ‖ ‖ ‖ 0 ⩽ 𝐶 ℎ ℓ ‖ ‖ ‖ 𝜕 𝑛 𝐄 𝜕 𝑡 𝑛 ‖ ‖ ‖ 𝐇 ℓ + 1 ( Ω ) ( 5 . 5 5 ) and with this we can now state the a priori error bound.

Corollary 5.9 (“concrete” a priori error bound). If for the solution of (4.21) we have 𝐄 ∈ 𝑊 4 1 ( 0 , 𝑇 ; 𝐋 2 ( Ω ) ) ∩ 𝑊 2 1 ( 0 , 𝑇 ; 𝐇 ℓ + 1 ( Ω ) ) ∩ 𝑊 1 ∞ ( 0 , 𝑇 ; 𝐇 ℓ + 1 ( Ω ) ) with ̆ ̆ 𝐄 , 𝐖 ∈ 𝐇 ℓ + 1 ( Ω ) then, under the same conditions as in Theorem 5.6, we have ‖ 𝐖 𝑗 − 𝐖 ℎ 𝑗 ‖ 0 𝐄 + ‖ ∇ × 𝑗 − 𝐄 ℎ 𝑗 ‖ 0 + 𝜎 + 𝜀 𝜓 𝑗 𝑗 𝑖 = 1 𝑘 ‖ ‖ 𝜕 𝑡 𝐄 𝑖 − 𝐄 ℎ 𝑖 ‖ ‖ 2 0 1 / 2 ≲ ℎ ℓ 1 + 𝑡 𝑗 + 𝜎 + 𝜀 𝜓 𝑗 𝑡 𝑗 1 / 2 ‖ 𝐄 ‖ 𝑊 1 ∞ ( 0 , 𝑡 𝑗 ; 𝐇 ℓ + 1 ( Ω ) ) + ℎ ℓ ‖ 𝐄 𝑡 ‖ 𝑊 1 1 ( 0 , 𝑡 𝑗 ; 𝐇 ℓ + 1 ( Ω ) ) + 𝑘 2 ‖ 𝐄 𝑡 𝑡 𝑡 ‖ 𝐿 ∞ ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) + ‖ 𝐄 𝑡 𝑡 ‖ 𝑊 2 1 ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) ( 5 . 5 6 ) for every 𝑗 ∈ { 1 , 2 , … , 𝑁 } .

Proof. First use Lemma 5.8 to estimate the terms in 𝒢 𝑖 and ℋ 𝑖 in Theorem 5.6, and remove the squares by using norm equivalence. This gives ‖ 𝐖 𝑗 − 𝐖 ℎ 𝑗 ‖ 0 + ‖ ‖ 𝐄 ∇ × 𝑗 − 𝐄 ℎ 𝑗 ‖ ‖ 0 + 𝜎 + 𝜀 𝜓 𝑗 𝑗 𝑖 = 1 𝑘 ‖ ‖ 𝜕 𝑡 𝐄 𝑖 − 𝐄 ℎ 𝑖 ‖ ‖ 2 0 1 / 2 ≲ ‖ 𝚲 𝑗 ‖ 0 𝑡 + ‖ 𝚯 𝑗 ‖ 𝐇 ( 𝐜 𝐮 𝐫 𝐥 ; Ω ) + ‖ 𝚯 𝑡 ‖ 𝐿 ∞ ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) + ‖ 𝚯 𝑡 ‖ 𝑊 1 1 ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) + 𝑗 𝑖 = 1 𝑘 ‖ 𝚯 𝑖 ‖ 0 + ( 𝜎 + 𝜀 𝜓 ′ 𝑗 ) 𝑗 𝑖 = 1 𝑘 ‖ 𝜕 𝑡 𝚯 𝑖 ‖ 2 0 1 / 2 + 𝑘 2 ‖ 𝐄 𝑡 𝑡 𝑡 ‖ 𝐿 ∞ ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) + ‖ 𝐄 𝑡 𝑡 ‖ 𝑊 2 1 ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) . ( 5 . 5 7 ) Now ‖ Λ 𝑗 ‖ 0 ⩽ ‖ 𝐖 − Π 𝐖 ‖ 𝐿 ∞ ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) ≲ ℎ ℓ ‖ 𝐄 ‖ 𝑊 1 ∞ ( 0 , 𝑡 𝑗 ; 𝐇 ℓ + 1 ( Ω ) ) by the best approximation property of the 𝐋 2 ( Ω ) projection and with (5.20) and (5.36), we get 𝑡 ‖ 𝚯 𝑗 ‖ 𝐇 ( 𝐜 𝐮 𝐫 𝐥 ; Ω ) + ‖ 𝚯 𝑡 ‖ 𝐿 ∞ ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) + ‖ 𝚯 𝑡 ‖ 𝑊 1 1 ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) + 𝑗 𝑖 = 1 𝑘 ‖ 𝚯 𝑖 ‖ 0 + 𝜎 + 𝜀 𝜓 𝑗 𝑗 𝑖 = 1 𝑘 ‖ 𝜕 𝑡 𝚯 𝑖 ‖ 2 0 1 / 2 ≲ 1 + 𝑡 𝑗 ‖ 𝚯 ‖ 𝐿 ∞ ( 0 , 𝑡 𝑗 ; 𝐇 ( 𝐜 𝐮 𝐫 𝐥 ; Ω ) ) + ‖ 𝚯 𝑡 ‖ 𝑊 1 1 ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) + 1 + 𝜎 + 𝜀 𝜓 𝑗 𝑡 𝑗 1 / 2 ‖ 𝚯 ‖ 𝑊 1 ∞ ( 0 , 𝑡 𝑗 ; 𝐋 2 ( Ω ) ) . ( 5 . 5 8 ) Applying (5.55) and the triangle inequality then completes the proof.

6. Concluding Remarks

We have attempted to provide a foundation for the long-time numerical integration of Maxwell's equations with Debye polarization by analogy and appeal to known results in viscoelasticity theory. As a result the hidden constant in Corollary 5.9 does not depend upon 𝑇 in any way and the case 𝜎 ⩾ 0 is allowed. We require only that 𝜇 > 0 , 𝜖 > 0 and also that Assumptions 2.1 hold.

There is no suggestion that the formulations and scheme presented above are in any way “the best”, they have been chosen merely to exemplify the approach. Indeed there are many ways in which this material could be extended and improved upon. We close with some of the more obvious ones.

The application of Baker's technique in [35] for the scalar wave equation could be studied. This leads to lower regularity requirements on the intial data and, to our knowledge, has never been extended to deal with problems with Volterra-type memory. Other avenues worthy of investigation are: variable coefficients as in [33]; mixed FE schemes as in [24]; staggered time stepping as in [32]; the effect of, for example, nonlinear relaxation time; first-order formulations as in [32, 33]; and, the derivation of a discrete version of Rivera-Menzala's lemma.

The alternative approach where the hereditary constitutive law is replaced by evolution equations for a set of “internal variables”, as in [1, 5] and described earlier in Section 3, also has a well-established analgoue in viscoelasticity theory, see [36], and it seems that numerical schemes could be formulated and analysed in a similar way to that in [37–39].

Also interesting is whether the case of a viscoelastic fluid has a Debye analogy. In this case 𝜑 0 = 0 in (2.3) and 𝜓 alters so that 𝜓 ( 𝑡 ) ∼ 𝑡 as 𝑡 → ∞ (some sharp energy estimates are given for quasistatic viscoelasticity in [14]). We note though from [5, Table 1] that the theory presented here is good for water and methanol. In fact, in the notation developed around (3.10) and (3.12), and we have for water the (normalized) creep-relaxation pair 𝜓 ( 𝑡 ) = 2 3 . 0 4 5 9 − 2 1 . 3 3 9 3 𝑒 − 𝑡 / 1 7 . 6 7 − 0 . 7 0 6 6 𝑒 − 𝑡 / 0 . 9 , 𝜑 ( 𝑡 ) = 0 . 0 4 3 4 + 0 . 2 4 3 3 𝑒 − 𝑡 / 1 . 7 9 6 9 + 0 . 7 1 3 3 𝑒 − 𝑡 / 0 . 3 8 4 0 ; ( 6 . 1 ) and for methanol, 𝜓 ( 𝑡 ) = 1 3 . 5 2 3 3 − 1 0 . 6 4 8 7 𝑒 − 𝑡 / 5 1 . 5 − 1 . 1 1 8 3 𝑒 − 𝑡 / 7 . 0 9 − 0 . 7 5 6 3 𝑒 − 𝑡 / 1 . 1 2 , 𝜑 ( 𝑡 ) = 0 . 0 7 3 9 + 0 . 1 2 6 7 𝑒 − 𝑡 / 1 4 . 2 0 1 1 + 0 . 2 5 0 2 𝑒 − 𝑡 / 3 . 7 3 7 4 + 0 . 5 4 9 2 𝑒 − 𝑡 / 0 . 5 6 9 8 . ( 6 . 2 ) In these the creep functions were taken from [5, Table 1] (see also [40, 41] for data) and the relaxation functions were determined numerically (using matlab R2007b) with a technique similar to that described by Golden and Graham in [26]. All auxiliary calculations are quoted to four decimal places, and time is in picoseconds ( 1 0 − 1 2 s ). The symbolic capability of matlab can be used to check the requirement that ( 𝜑 ∗ 𝜓 ) ( 𝑡 ) = 𝑡 0 𝜑 ( 𝑡 − 𝑠 ) 𝜓 ( 𝑠 ) 𝑑 𝑠 = 𝑡 ( 6 . 3 ) and the values of 𝑡 − 𝜑 ∗ 𝜓 are shown in Figure 1. The quality could be improved by working to more than four decimal places, but we should bear in mind that the experimental errors are probably already swamping the numerical errors shown in the figures.

fig1
Figure 1: Graphical representation of 𝑡 − 𝜑 ∗ 𝜓 for the numerically determined water and methanol creep-relaxation pairs discussed in Section 6.

Another avenue of development is to generalise the frequency dependence and use the more general hereditary laws as given by [11, Chapter 9.3.3]. Here, for the constitutive laws we have 𝐃 ( 𝐱 , 𝑡 ) = 𝜀 0 ( 𝐱 ) 𝐄 ( 𝐱 , 𝑡 ) + 𝑡 0 𝜀 1 ( 𝐱 , 𝑡 , 𝑠 ) 𝐄 ( 𝐱 , 𝑠 ) 𝑑 𝑠 , 𝐁 ( 𝐱 , 𝑡 ) = 𝜇 0 ( 𝐱 ) 𝐇 ( 𝐱 , 𝑡 ) + 𝑡 0 𝜇 1 ( 𝐱 , 𝑡 , 𝑠 ) 𝐇 ( 𝐱 , 𝑠 ) 𝑑 𝑠 , 𝐉 ( 𝐱 , 𝑡 ) = 𝜎 0 ( 𝐱 ) 𝐄 ( 𝐱 , 𝑡 ) + 𝑡 0 𝜎 1 ( 𝐱 , 𝑡 , 𝑠 ) 𝐄 ( 𝐱 , 𝑠 ) 𝑑 𝑠 , ( 6 . 4 ) as well as the inverses of the first two (whenever they exist).

Lastly we mention that for kernels in the form (2.3), the numerical implementation of the Volterra operators is neither slow nor cumbersome. We refer to [42] for details.

References

  1. H. T. Banks, V. A. Bokil, and N. L. Gibson, “Analysis of stability and dispersion in a finite element method for Debye and Lorentz dispersive media,” Numerical Methods for Partial Differential Equations, vol. 25, no. 4, pp. 885–917, 2009. View at Publisher · View at Google Scholar
  2. C. De Cicco, L. Mariani, C. Vedruccio, et al., “Clinical application of spectral electromagnetic interaction in breast cancer: diagnosis results of a pilot study,” August 2009, http://www.piwe.it/elements_page1_files/BreastTumori.pdf.
  3. K. R. Foster, J. L. Schepps, R. D. Stoy, and H. P. Schwan, “Dielectric properties of brain tissue between 0.01 and 10 GHz,” Physics in Medicine and Biology, vol. 24, no. 6, pp. 1177–1187, 1979. View at Publisher · View at Google Scholar
  4. J.-Z. Bao, S.-T. Lu, and W. D. Hurt, “Complex dielectric measurements and analysis of brain tissues in the radio and microwave frequencies,” IEEE Transactions on Microwave Theory and Techniques, vol. 45, no. 10, pp. 1730–1741, 1997. View at Scopus
  5. J. L. Young and R. S. Adams, “On the time integration of Maxwell's equations associated with Debye relaxation processes,” IEEE Transactions on Antennas and Propagation, vol. 55, no. 8, pp. 2409–2412, 2007. View at Publisher · View at Google Scholar · View at MathSciNet
  6. L. S. Taylor and A. Y. Cheung, The Physical Basis of Electromagnetic Interactions with Biological Systems, University of Maryland, College Park, Md, USA, 1977.
  7. W. Grundler, F. Kaiser, F. Keilmann, and J. Walleczek, “Mechanisms of electromagnetic interaction with cellular systems,” Naturwissenschaften, vol. 79, no. 12, pp. 551–559, 1992. View at Publisher · View at Google Scholar
  8. A. H. Frey, “Electromagnetic field interactions with biological systems,” FASEB Journal, vol. 7, no. 2, pp. 272–281, 1993.
  9. P. Monk, Finite Element Methods for Maxwell’s Equations. Numerical Mathematics and Scientific Computation, Oxford University Press, Oxford, UK, 2003.
  10. C. Johnson, “Discontinuous Galerkin finite element methods for second order hyperbolic problems,” Computer Methods in Applied Mechanics and Engineering, vol. 107, no. 1-2, pp. 117–129, 1993.
  11. M. Fabrizio and A. Morro, Electromagnetism of Continuous Media: Mathematical Modelling and Applications, Oxford University Press, Oxford, UK, 2003.
  12. J. E. Muñoz Rivera and G. P. Menzala, “Decay rates of solutions to a von Kármán system for viscoelastic plates with memory,” Quarterly of Applied Mathematics, vol. 57, no. 1, pp. 181–200, 1999.
  13. V. Thomée and L. B. Wahlbin, “Long-time numerical solution of a parabolic equation with memory,” Mathematics of Computation, vol. 62, pp. 477–496, 1994.
  14. S. Shaw and J. R. Whiteman, “Optimal long-time Lp(0, T) stability and semidiscrete error estimates for the Volterra formulation of the linear quasistatic viscoelasticity problem,” Numerische Mathematik, vol. 88, no. 4, pp. 743–770, 2001, BICOM Tech. Rep. 98/7 http://www.brunel.ac.uk/bicom.
  15. S. Shaw, M. K. Warby, and J. R. Whiteman, “Error estimates with sharp constants for a fading memory Volterra problem in linear solid viscoelasticity,” SIAM Journal on Numerical Analysis, vol. 34, no. 3, pp. 1237–1254, 1997.
  16. T. M. Roberts and P. G. Petropoulos, “Asymptotics and energy estimates for electromagnetic pulses in dispersive media,” Journal of the Optical Society of America A, vol. 13, no. 6, pp. 1204–1217, 1996.
  17. P. G. Petropoulos, “The wave hierarchy for propagation in relaxing dielectrics,” Wave Motion, vol. 21, no. 3, pp. 253–262, 1995.
  18. R. Lakes, Viscoelastic Materials, Cambridge University Press, New York, NY, USA, 2009.
  19. J. Lee, R. Lee, and A. Cangellaris, “Time-domain finite-element methods,” IEEE Transactions on Antennas and Propagation, vol. 45, no. 3, pp. 430–442, 1997.
  20. D. Jiao and J.-M. Jin, “Time-domain finite-element modeling of dispersive media,” IEEE Microwave and Wireless Components Letters, vol. 11, no. 5, pp. 220–222, 2001. View at Publisher · View at Google Scholar
  21. T. Lu, P. Zhang, and W. Cai, “Discontinuous Galerkin methods for dispersive and lossy Maxwell's equations and PML boundary conditions,” Journal of Computational Physics, vol. 200, no. 2, pp. 549–580, 2004. View at Publisher · View at Google Scholar · View at MathSciNet
  22. J. Li, “Error analysis of finite element methods for 3-D Maxwell's equations in dispersive media,” Journal of Computational and Applied Mathematics, vol. 188, no. 1, pp. 107–120, 2006. View at Publisher · View at Google Scholar · View at MathSciNet
  23. J. Li and Y. Chen, “Analysis of a time-domain finite element method for 3-D Maxwell's equations in dispersive media,” Computer Methods in Applied Mechanics and Engineering, vol. 195, no. 33–36, pp. 4220–4229, 2006. View at Publisher · View at Google Scholar · View at MathSciNet
  24. J. Li, “Error analysis of fully discrete mixed finite element schemes for 3-D Maxwell's equations in dispersive media,” Computer Methods in Applied Mechanics and Engineering, vol. 196, no. 33-34, pp. 3081–3094, 2007. View at Publisher · View at Google Scholar · View at MathSciNet
  25. J. Li, “Posteriori error estimation for an interior penalty discontinuous Galerkin method for Maxwell's equations in cold plasma,” Advances in Applied Mathematicsand Mechanics, vol. 1, pp. 107–124, 2009.
  26. J. M. Golden and G. A. C. Graham, Boundary Value Problems in Linear Viscoelasticity, Springer, New York, NY, USA, 1988.
  27. J. D. Ferry, Viscoelastic Properties of Polymers, John Wiley and Sons, New York, NY, USA, 1970.
  28. W. N. Findley, J. S. Lai, and K. Onaran, Creep and Relaxation of Nonlinear Viscoelastic Materials (With an Introduction to Linear Viscoelasticity), Dover, New York, NY, USA, 1989.
  29. H. R. Hill, Adaptive finite elements for viscoelastic deformation problems, Ph.D. thesis, Brunel University, Uxbridge, UK, 2009, http://www.brunel.ac.uk/bicom.
  30. J. C. Nédélec, “Mixed finite elements in ℝ3,” Numerische Mathematik, vol. 35, no. 3, pp. 315–341, 1980. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  31. J. C. Nédélec, “A new family of mixed finite elements in ℝ3,” Numerische Mathematik, vol. 50, pp. 57–81, 1986.
  32. C. Ma, “Finite-element method for time-dependent Maxwell's equations based on an explicit-magnetic-field scheme,” Journal of Computational and Applied Mathematics, vol. 194, no. 2, pp. 409–424, 2006. View at Publisher · View at Google Scholar · View at MathSciNet
  33. J. Zhao, “Analysis of finite element approximation for time-dependent Maxwell problems,” Mathematics of Computation, vol. 73, no. 247, pp. 1089–1105, 2004. View at Publisher · View at Google Scholar · View at MathSciNet
  34. P. Monk, “A finite element method for approximating the time-harmonic Maxwell equations,” Numerische Mathematik, vol. 63, no. 1, pp. 243–261, 1992. View at Publisher · View at Google Scholar · View at MathSciNet
  35. G. A. Baker, “Error estimates for finite element methods for second order hyperbolic equations,” SIAM Journal on Numerical Analysis, vol. 13, no. 4, pp. 546–576, 1976.
  36. A. R. Johnson and A. Tessler, “A viscoelastic high order beam finite element,” in The Mathematics of Finite Elements and Applications (MAFELAP 1996), J. R. Whiteman, Ed., pp. 333–345, John Wiley and Sons, Chichester, UK, 1997.
  37. N. Bauermeister and S. Shaw, “Finite-element approximation of non-Fickian polymer diffusion,” IMA Journal of Numerical Analysis, vol. 30, pp. 702–730, 2010, BICOM Tech. Rep. 08/1, http://www.brunel.ac.uk/bicom. View at Publisher · View at Google Scholar
  38. B. Rivière, S. Shaw, and J. R. Whiteman, “Discontinuous Galerkin finite element methods for dynamic linear solid viscoelasticity problems,” Numerical Methods for Partial Differential Equations, vol. 23, no. 5, pp. 1149–1166, 2007. View at Publisher · View at Google Scholar · View at MathSciNet
  39. Report 05/7 at http://www.brunel.ac.uk/bicom.
  40. J. Barthel and R. Buchner, “High frequency permittivity and its use in the investigation of solution properties,” Pure and Applied Chemistry, vol. 63, no. 10, pp. 1473–1482, 1991.
  41. R. Buchner, J. Barthel, and J. Stauber, “The dielectric relaxation of water between 0°C and 35°C,” Chemical Physics Letters, vol. 306, no. 1-2, pp. 57–63, 1999.
  42. S. Shaw, M. K. Warby, J. R. Whiteman, C. Dawson, and M. F. Wheeler, “Numerical techniques for the treatment of quasistatic viscoelastic stress problems in linear isotropic solids,” Computer Methods in Applied Mechanics and Engineering, vol. 118, no. 3-4, pp. 211–237, 1994.