# Finite Element Approximation of Maxwell's Equations with Debye Memory

Authors: Simon Shaw
Journal: Advances in Numerical Analysis
 (2010)

---

## 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.


---

## Body

## 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 (athttp://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 I:=(0,T] a finite time interval. We consider Maxwell's equations (e.g., Monk, [9]) in Ω×I in the form,(1.1)∇×E+Ḃ=0,∇×H-Ḋ=J,∇⋅D=ϱ,∇⋅B=0,
where here and throughout the overdot denotes partial time differentiation. In these E is the electric field intensity, H is the magnetic field intensity, D is the electric displacement, and B 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 J is the current density. Conservation of charge implies that ∇·J+ϱ̇=0 and, combined with the above, means that if the Gauss laws hold for one time, t, then they hold for all t. We therefore build these conditions into the intial data, and take the opportunity to specify the boundary data.(1.2)∇⋅D=ϱ,∇⋅B=0inΩatt=0,(1.3)n̂×E=0,n̂⋅H=0on∂Ω×I,
where n̂ is the unit outward normal on the boundary.We assume thatB=μH where 0<μ∈ℝ, which is standard, but are concerned here with the Debye model of a dielectric wherein D is given as a linear functional of E. 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 Lp 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 M+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 Section2 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, D, in Section 4.1 and also as a problem for the electric field, E, 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 thatOhm's law holds so that J=σE+Ja, where Ja 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,Lp(Ω), Wpm(Ω), and Hm(Ω) are the usual Lebesgue, Banach and Hilbert spaces and we set Lp(Ω):=Lp(Ω)3, Hm(Ω):=Hm(Ω)3, and so on. In general terms, if X is a Banach space then its norm will be denoted by ∥·∥X. The only exception to this will be that ∥·∥0 will always mean the L2(Ω) norm as induced by (·,·), the L2(Ω) 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 f:[0,T]→X, then we use standard notation such as ∥f∥Lp(0,T;X) to signify the Lp(0,T) norm of ∥f(t)∥X.The natural space for Maxwell's equations as we will work with them isH(curl;Ω) where(1.4)H(curl;Ω):={v∈L2(Ω):∇×v∈L2(Ω)}
with graph norm, ∥v∥H(curl;Ω):=(∥v∥02+∥∇×v∥02)1/2. This is a Hilbert space when equipped with the obvious norm-inducing inner product (see, e.g., [9]). We also define,(1.5)V:=H0(curl;Ω):={v∈H(curl;Ω):n̂×v=0on∂Ω},
(also a Hilbert space) where n̂ is a.e. the unit outward normal on ∂Ω. For Friedrich's-type properties of V 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,(1.6)2ab⩽a2ϵ+ϵb2,∀a,b∈R,∀ϵ>0,
and Young's inequality for convolutions,(1.7)‖w*v‖Lr(0,t)⩽‖w‖Lp(0,t)‖v‖Lq(0,t)where(w*v)(t):=∫0tw(t-s)v(s)ds,
and 1⩽p,q,r⩽∞ such that p-1+q-1=1+r-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:(2.1)σ̲(t)=C̲ɛ̲(t)-∫0tφs(t-s)C̲ɛ̲(s)ds.
This can be inverted to the following creep form:(2.2)C̲ɛ̲(t)=σ̲(t)-∫0tψs(t-s)σ̲(s)ds.
In these φ is a stress relaxation function and ψ is a creep function. Also, the subscript denotes partial differentiation so that φs(t-s)=-φt(t-s)=-φ′(t-s), 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 s-integral of an s-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,(2.3)φ(t)=φ0+∑i=1Nφie-t/τi′,ψ(t)=ψ0+∑i=1Nψie-t/τi
subject to normalization, φ(0)=ψ(0)=1, and where φi⩾0 and τi′,τi>0 for each i. 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 ψ(t)~t as t→∞). The tensor C̲ 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𝔏:t→p and denote 𝔏F by F̂ then it is readily shown that the relaxation and creep functions are related by p2φ̂(p)ψ̂(p)=1 and it follows that, for a solid, φ(∞)ψ(∞)=φ0ψ0=1. In particular, it is easy to derive the pair, (2.4)φ(t)=φ0+φ1e-t/τ⟺ψ(t)=ψ0-ψ1e-φ0t/τwhereψ0=1φ0,ψ1=φ1φ0.
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)
φ∈C3(0,∞) with φ(0)=1, limt↑∞φ(t)=φ̌>0, φ′<0 and φ′′>0,(ii)
ψ∈C3(0,∞) with ψ(0)=1, limt↑∞ψ(t)=φ̌-1>0, ψ′>0 and ψ′′<0,along with p2ψ̂(p)φ̂(p)=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 D=ɛ0E+P̃ where ɛ0 is the permittivity of free space and P̃=PI+P. Here PI=ɛ0χE is the instantaneous polarization and P is the “relaxation polarization”. With ɛ∞:=1+χ, the constitutive model used in [1] is then(3.1)D=ɛ0ɛ∞E+PwithτṖ+P=(ɛs-ɛ∞)ɛ0EsubjecttoP(0)=P̆.
Note that here and throughout the dependence on the spatial variable x∈Ω is suppressed. Integrating this ODE gives(3.2)P(t)=e-t/τP̆+∫0t(ɛs-ɛ∞)ɛ0τe-(t-s)/τE(s)ds,
which, by henceforth assuming that P̆=0 (since the instantaneous effects are already taken care of in PI), gives us(3.3)1ɛ0ɛ∞D(t)=E(t)+ɛs-ɛ∞τϵ∞∫0te-(t-s)/τE(s)ds.
The connection to viscoelasticity can now be made evident by defining the creep function(3.4)ψ(t):=1-φ1e-φ0t/τ′φ0=ψ0-ψ1e-t/τ,
where τ′:=φ0τ, ψ0=1/φ0, ψ1=φ1/φ0, φ0:=ɛ∞/ɛs, and φ1:=(ɛs-ɛ∞)/ɛs, and then for C̲:=(ɛ0ɛ∞)-1I̲, we have,(3.5)C̲D(t)=E(t)-∫0tψs(t-s)E(s)ds.
Immediately we can observe the analogy, or correspondence, between the electric displacement and field, D, E, and the viscoelastic strain and stress, ɛ̲, σ̲. We can also conclude that this relationship can be inverted to give,(3.6)E(t)=C̲D(t)-C̲∫0tφs(t-s)D(s)ds
for the relaxation functionφ(t)=φ0+φ1e-t/τ′. 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<ɛ∞⩽ɛs. 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 whereE is constant and ɛs>ɛ∞ we have(3.7)D(t)=(1-∫0tψs(t-s)ds)ɛ0ɛ∞E=ψ(t)ɛ0ɛ∞E
and so, with ∥·∥𝔼 the Euclidean norm on ℝ3, we have ∥D(t1)∥𝔼<∥D(t2)∥𝔼 for all 0⩽t1<t2<∞. We can see that this constitutive law represents a monotone creeping process for the electric displacement which is bounded by ∥D(0)∥𝔼=ɛ0ɛ∞∥E∥𝔼 and ∥D(∞)∥𝔼=φ0-1ɛ0ɛ∞∥E∥𝔼.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 ɛ¯′:=ɛsɛ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(3.8)D=ɛ0ɛ∞E+∑m=1MPmwithτmṖm+Pm=(ɛsm-ɛ∞)ɛ0EsubjecttoPm(0)=P̆m.
The connection between this model and the Banks et al. model is clear and we have,(3.9)Pm(t)=e-t/τmP̆m+(ɛsm-ɛ∞)ɛ0τm∫0te-(t-s)/τmE(s)ds.
Assuming, as before, that each P̆m=0 we can therefore write(3.10)1ɛ0ɛ∞D(t)=E(t)-∫0tψs(t-s)E(s)ds,
where(3.11)ψ(t):=ψ0-∑m=1Mψme-t/τmwithψm:=ɛsm-ɛ∞ɛ∞
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(3.12)E(t)=1ɛ0ɛ∞D(t)-1ɛ0ɛ∞∫0tφs(t-s)D(s)ds.
Typically, for the creep function in use here we would expect that(3.13)φ(t)=φ0+∑m=1Mφme-t/τm′
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 B=μH to derive partial differential Volterra equations. Firstly, in Section 4.1 we use (3.12) to obtain a PDE for the electric displacement, D, and then in Section 4.2 we use (3.10) to derive a PDE for the electric field, E. 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(4.1)D̈+∇×μ-1∇×E+σĖ=-J̇a
which, recalling (1.5), can be written in weak form as(4.2)(D̈,ϕ)+(μ-1∇×E,∇×ϕ)+(σĖ,ϕ)=-(J̇a,ϕ)∀ϕ∈V.
However, from (3.12) with ɛ:=ɛ0ɛ∞ we obtain after a partial integration(4.3)ɛĖ(t)=Ḋ(t)+φ′(t)D(0)-∫0tφs(t-s)Ḋ(s)ds.
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,D:I→V such that,(4.4)(D̈(t),ϕ)+(μ-1∇×ɛ-1D(t),∇×ϕ)-∫0tφs(t-s)(μ-1∇×ɛ-1D(s),∇×ϕ)ds+(σɛ(Ḋ(t)-∫0tφs(t-s)Ḋ(s)ds),ϕ)=(G(t),ϕ)
with G(t)=-J̇a(t)-σɛ-1φ′(t)D(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 Assumptions2.1, hold then
(4.5)∫0t(w(s)-∫0sφξ(s-ξ)w(ξ)dξ,w(s))ds⩾φ̌‖w‖L2(0,t;L2(Ω))2
for all t∈I.Proof.
Since in general terms(u±v,u)⩾∥u∥2-∥v∥∥u∥ we need only to observe that ∥(φ′*w,w)∥L1(0,t)⩽∥(-φ′*∥w∥0)∥w∥0∥L1(0,t)⩽∥-φ′∥L1(0,t)∥w∥L2(0,t;L2(Ω))2 from Young's inequality, (1.7). Assumptions 2.1 yield 1-∥-φ′∥L1(0,t)=φ(t)⩾φ̌ and this completes the proof.Lemma 4.2 (Rivera and Menzala, [12]).
Leta(w,v):=(∇×w,∇×v), then
(4.6)-2∫0τ∫0tφs(t-s)a(v(s),v̇(t))dsdt=φ(τ)a(v(τ),v(τ))-a(v(τ),v(τ))-∫0τφ′(τ-s)a(v(τ)-v(s),v(τ)-v(s))ds-∫0τφ′(t)a(v(t),v(t))dt+∫0τ∫0tφ′′(t-s)a(v(t)-v(s),v(t)-v(s))dsdt
for all τ∈I and for all v∈H1(I;V).Proof.
First we note that,(4.7)ddt∫0tφt(t-s)a(v(t)-v(s),v(t)-v(s))ds=∫0tφtt(t-s)a(v(t)-v(s),v(t)-v(s))ds-2∫0tφt(t-s)a(v(s),v̇(t))ds+2∫0tφt(t-s)a(v(t),v̇(t))ds,
and substitute for the last term on the right by using
(4.8)ddt∫0tφt(t-s)a(v(t),v(t))ds=φ′(t)a(v(t),v(t))+2∫0tφt(t-s)a(v(t),v̇(t))ds.
After a rearrangement of terms, this gives
(4.9)-2∫0tφs(t-s)a(v(s),v̇(t))ds=ddt∫0tφt(t-s)a(v(t),v(t))ds-φ′(t)a(v(t),v(t))-ddt∫0tφ′(t-s)a(v(t)-v(s),v(t)-v(s))ds+∫0tφ′′(t-s)a(v(t)-v(s),v(t)-v(s))ds.
The proof is then completed by integrating over (0,τ), observing that
(4.10)∫0τφτ(τ-s)a(v(τ),v(τ))ds=(φ(τ)-φ(0))a(v(τ),v(τ))
and recalling that φ(0)=1.Theorem 4.3 (stability for theD equation).
If Assumptions2.1 hold and if μ>0, ɛ>0 and σ⩾0 are real numbers, then for every t∈I we have the following. Firstly, if G=0,
(4.11)μɛ‖Ḋ(t)‖02+φ̌‖∇×D(t)‖02+2μσφ̌‖Ḋ‖L2(0,t;L2(Ω))2⩽μɛ‖Ḋ(0)‖02+‖∇×D(0)‖02.
Secondly, if G≠0 and σ>0,
(4.12)μɛ‖Ḋ(t)‖02+φ̌‖∇×D(t)‖02+μσφ̌‖Ḋ‖L2(0,t;L2(Ω))2⩽μɛ‖Ḋ(0)‖02+‖∇×D(0)‖02+μɛ2σφ̌‖G‖L2(0,t;L2(Ω))2.
Thirdly, if σ=0,
(4.13)μɛ2‖Ḋ‖L∞(0,t;L2(Ω))2+φ̌∥∇×D(t)∥02⩽2μɛ∥Ḋ(0)∥02+2∥∇×D(0)∥02+8μɛ∥G∥L1(0,t;L2(Ω))2.Proof.
We chooseϕ=2Ḋ(t) in (4.4) and integrate over (0,τ) to get
(4.14)μɛ∥Ḋ(τ)∥02+∥∇×D(τ)∥02-2∫0τ∫0tφs(t-s)(∇×D(s),∇×Ḋ(t))dsdt+2μσ∫0τ(Ḋ(t)-∫0tφs(t-s)Ḋ(s)ds,Ḋ(t))dt=μɛ∥Ḋ(0)∥02+∥∇×D(0)∥02+2μɛ∫0τ(G(t),Ḋ(t))dt.
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 G=0. When G≠0 and σ>0, Young's inequality, (1.6), gives,
(4.15)2μɛ∫0τ(G(t),Ḋ(t))dt⩽μɛ2σφ̌∥G∥L2(0,t;L2(Ω))2+μσφ̌∥Ḋ∥L2(0,t;L2(Ω))2
and we can complete the proof of the second claim. If σ=0 and G≠0, we deduce that
(4.16)μɛ∥Ḋ∥L∞(0,τ;L2(Ω))2+φ̆∥∇×D(τ)∥02⩽2μɛ∥Ḋ(0)∥02+2∥∇×D(0)∥02+4μɛ∫0τ|(G(t),Ḋ(t))|dt,⩽2μɛ∥Ḋ(0)∥02+2∥∇×D(0)∥02+8μɛ∥G∥L1(0,τ;L2(Ω))2+μɛ2∥Ḋ∥L∞(0,τ;L2(Ω))2,
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(4.17)D̈+∇×1μ∇×E+σĖ=-J̇a,
but now use (3.10) to get(4.18)Ḋ(t)=ɛĖ(t)+ɛψ′(t)E(0)-ɛ∫0tψs(t-s)Ė(s)ds,
from which,(4.19)D̈(t)=ɛË(t)+ɛψ′′(t)E(0)+ɛψ′(0)Ė(t)+ɛ∫0tψss(t-s)Ė(s)ds.
Introducing this results in the partial differential Volterra equation, (4.20)ɛË(t)+ɛψ′(0)Ė(t)+ɛ∫0tψ′′(t-s)Ė(s)ds+∇×1μ∇×E(t)+σĖ(t)=-J̇a(t)-ɛψ′′(t)E(0)
or, more usefully, the weak problem: find a twice differentiable map E:I→V such that,(4.21)(ɛË(t),ϕ)+(μ-1∇×E(t),∇×ϕ)+(σĖ(t),ϕ)+(ɛψ′(0)Ė(t)+ɛ∫0tψ′′(t-s)Ė(s)ds,ϕ)=(F(t),ϕ)∀ϕ∈V,
with F:=-J̇a(t)-ɛψ′′(t)E(0) known and given intial data.We require a slight adjustment to Lemma4.1 before giving the stability estimate.Lemma 4.4 (coercivity of creep).
If Assumptions2.1 hold, then
(4.22)∫0t(ψ′(0)w(s)+∫0sψ′′(s-ξ)w(ξ)dξ,w(s))ds⩾ψ′(t)∥w∥L2(0,t;L2(Ω))2
for all t∈I.Proof.
Since∥-ψ′′*∥w∥0∥L2(0,t)⩽∥ψ′′∥L1(0,t)∥w∥L2(0,t;L2(Ω)) and ψ′(0)-∥ψ′′∥L1(0,t)=ψ′(0)+∫0tψ′′(s)ds=ψ′(t) this is similar to Lemma 4.1.Theorem 4.5 (stability for theE equation).
If Assumptions2.1 hold and if μ>0, ɛ>0 and σ⩾0 are real numbers then for every t∈I we have first,
(4.23)μɛ∥Ė(t)∥02+∥∇×E(t)∥02+μ(σ+ɛψ′(t))∥Ė∥L2(0,t;L2(Ω))2⩽μɛ∥Ė(0)∥02+∥∇×E(0)∥02+μ∥F∥L2(0,t;L2(Ω))2(σ+ɛψ′(t)),
and second,
(4.24)μɛ2∥Ė∥L∞(0,t;L2(Ω))2+∥∇×E(t)∥02+2μ(σ+ɛψ′(t))∥Ė∥L2(0,t;L2(Ω))2⩽2μɛ∥Ė(0)∥02+2∥∇×E(0)∥02+8μɛ∥F∥L1(0,t;L2(Ω))2.Proof.
We chooseϕ=2Ė(t) in (4.21), integrate over (0,τ) and use Lemma 4.4 to get,
(4.25)μɛ∥Ė(τ)∥02+∥∇×E(τ)∥02+2μ(σ+ɛψ′(τ))∥Ė∥L2(0,τ;L2(Ω))2⩽μɛ∥Ė(0)∥02+∥∇×E(0)∥02+2μ∫0τ(F(t),Ė(t))dt.
Now by applying the Cauchy-Schwarz and Young inequalities in the following form:
(4.26)2μ∫0τ(F(t),Ė(t))dt⩽μ(σ+ɛψ′(τ))∥F∥L2(0,τ;L2(Ω))2+μ(σ+ɛψ′(τ))∥Ė∥L2(0,τ;L2(Ω))2,
we complete the proof of the first claim. For the second, we easily see that
(4.27)μɛ∥Ė∥L∞(0,τ;L2(Ω))2+∥∇×E(τ)∥02+2μ(σ+ɛψ′(τ))∥Ė∥L2(0,τ;L2(Ω))2⩽2μɛ∥Ė(0)∥02+2∥∇×E(0)∥02+4μ∫0τ|(F(t),Ė(t))|dt,
and with the Hölder and Young inequalities,
(4.28)4μ∫0τ|(F(t),Ė(t))|dt⩽8μɛ∥F∥L1(0,τ;L2(Ω))2+μɛ2∥Ė∥L∞(0,τ;L2(Ω))2,
we arrive at the second claim.Remark 4.6 (growth forσ=0).
Ifσ=0 then the fact that ψ′→0 as t→∞ 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.


## 4.1. A Formulation for the Electric Displacement

From (1.1) and Ohm's law, we obtain(4.1)D̈+∇×μ-1∇×E+σĖ=-J̇a
which, recalling (1.5), can be written in weak form as(4.2)(D̈,ϕ)+(μ-1∇×E,∇×ϕ)+(σĖ,ϕ)=-(J̇a,ϕ)∀ϕ∈V.
However, from (3.12) with ɛ:=ɛ0ɛ∞ we obtain after a partial integration(4.3)ɛĖ(t)=Ḋ(t)+φ′(t)D(0)-∫0tφs(t-s)Ḋ(s)ds.
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,D:I→V such that,(4.4)(D̈(t),ϕ)+(μ-1∇×ɛ-1D(t),∇×ϕ)-∫0tφs(t-s)(μ-1∇×ɛ-1D(s),∇×ϕ)ds+(σɛ(Ḋ(t)-∫0tφs(t-s)Ḋ(s)ds),ϕ)=(G(t),ϕ)
with G(t)=-J̇a(t)-σɛ-1φ′(t)D(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 Assumptions2.1, hold then
(4.5)∫0t(w(s)-∫0sφξ(s-ξ)w(ξ)dξ,w(s))ds⩾φ̌‖w‖L2(0,t;L2(Ω))2
for all t∈I.Proof.
Since in general terms(u±v,u)⩾∥u∥2-∥v∥∥u∥ we need only to observe that ∥(φ′*w,w)∥L1(0,t)⩽∥(-φ′*∥w∥0)∥w∥0∥L1(0,t)⩽∥-φ′∥L1(0,t)∥w∥L2(0,t;L2(Ω))2 from Young's inequality, (1.7). Assumptions 2.1 yield 1-∥-φ′∥L1(0,t)=φ(t)⩾φ̌ and this completes the proof.Lemma 4.2 (Rivera and Menzala, [12]).
Leta(w,v):=(∇×w,∇×v), then
(4.6)-2∫0τ∫0tφs(t-s)a(v(s),v̇(t))dsdt=φ(τ)a(v(τ),v(τ))-a(v(τ),v(τ))-∫0τφ′(τ-s)a(v(τ)-v(s),v(τ)-v(s))ds-∫0τφ′(t)a(v(t),v(t))dt+∫0τ∫0tφ′′(t-s)a(v(t)-v(s),v(t)-v(s))dsdt
for all τ∈I and for all v∈H1(I;V).Proof.
First we note that,(4.7)ddt∫0tφt(t-s)a(v(t)-v(s),v(t)-v(s))ds=∫0tφtt(t-s)a(v(t)-v(s),v(t)-v(s))ds-2∫0tφt(t-s)a(v(s),v̇(t))ds+2∫0tφt(t-s)a(v(t),v̇(t))ds,
and substitute for the last term on the right by using
(4.8)ddt∫0tφt(t-s)a(v(t),v(t))ds=φ′(t)a(v(t),v(t))+2∫0tφt(t-s)a(v(t),v̇(t))ds.
After a rearrangement of terms, this gives
(4.9)-2∫0tφs(t-s)a(v(s),v̇(t))ds=ddt∫0tφt(t-s)a(v(t),v(t))ds-φ′(t)a(v(t),v(t))-ddt∫0tφ′(t-s)a(v(t)-v(s),v(t)-v(s))ds+∫0tφ′′(t-s)a(v(t)-v(s),v(t)-v(s))ds.
The proof is then completed by integrating over (0,τ), observing that
(4.10)∫0τφτ(τ-s)a(v(τ),v(τ))ds=(φ(τ)-φ(0))a(v(τ),v(τ))
and recalling that φ(0)=1.Theorem 4.3 (stability for theD equation).
If Assumptions2.1 hold and if μ>0, ɛ>0 and σ⩾0 are real numbers, then for every t∈I we have the following. Firstly, if G=0,
(4.11)μɛ‖Ḋ(t)‖02+φ̌‖∇×D(t)‖02+2μσφ̌‖Ḋ‖L2(0,t;L2(Ω))2⩽μɛ‖Ḋ(0)‖02+‖∇×D(0)‖02.
Secondly, if G≠0 and σ>0,
(4.12)μɛ‖Ḋ(t)‖02+φ̌‖∇×D(t)‖02+μσφ̌‖Ḋ‖L2(0,t;L2(Ω))2⩽μɛ‖Ḋ(0)‖02+‖∇×D(0)‖02+μɛ2σφ̌‖G‖L2(0,t;L2(Ω))2.
Thirdly, if σ=0,
(4.13)μɛ2‖Ḋ‖L∞(0,t;L2(Ω))2+φ̌∥∇×D(t)∥02⩽2μɛ∥Ḋ(0)∥02+2∥∇×D(0)∥02+8μɛ∥G∥L1(0,t;L2(Ω))2.Proof.
We chooseϕ=2Ḋ(t) in (4.4) and integrate over (0,τ) to get
(4.14)μɛ∥Ḋ(τ)∥02+∥∇×D(τ)∥02-2∫0τ∫0tφs(t-s)(∇×D(s),∇×Ḋ(t))dsdt+2μσ∫0τ(Ḋ(t)-∫0tφs(t-s)Ḋ(s)ds,Ḋ(t))dt=μɛ∥Ḋ(0)∥02+∥∇×D(0)∥02+2μɛ∫0τ(G(t),Ḋ(t))dt.
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 G=0. When G≠0 and σ>0, Young's inequality, (1.6), gives,
(4.15)2μɛ∫0τ(G(t),Ḋ(t))dt⩽μɛ2σφ̌∥G∥L2(0,t;L2(Ω))2+μσφ̌∥Ḋ∥L2(0,t;L2(Ω))2
and we can complete the proof of the second claim. If σ=0 and G≠0, we deduce that
(4.16)μɛ∥Ḋ∥L∞(0,τ;L2(Ω))2+φ̆∥∇×D(τ)∥02⩽2μɛ∥Ḋ(0)∥02+2∥∇×D(0)∥02+4μɛ∫0τ|(G(t),Ḋ(t))|dt,⩽2μɛ∥Ḋ(0)∥02+2∥∇×D(0)∥02+8μɛ∥G∥L1(0,τ;L2(Ω))2+μɛ2∥Ḋ∥L∞(0,τ;L2(Ω))2,
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(4.17)D̈+∇×1μ∇×E+σĖ=-J̇a,
but now use (3.10) to get(4.18)Ḋ(t)=ɛĖ(t)+ɛψ′(t)E(0)-ɛ∫0tψs(t-s)Ė(s)ds,
from which,(4.19)D̈(t)=ɛË(t)+ɛψ′′(t)E(0)+ɛψ′(0)Ė(t)+ɛ∫0tψss(t-s)Ė(s)ds.
Introducing this results in the partial differential Volterra equation, (4.20)ɛË(t)+ɛψ′(0)Ė(t)+ɛ∫0tψ′′(t-s)Ė(s)ds+∇×1μ∇×E(t)+σĖ(t)=-J̇a(t)-ɛψ′′(t)E(0)
or, more usefully, the weak problem: find a twice differentiable map E:I→V such that,(4.21)(ɛË(t),ϕ)+(μ-1∇×E(t),∇×ϕ)+(σĖ(t),ϕ)+(ɛψ′(0)Ė(t)+ɛ∫0tψ′′(t-s)Ė(s)ds,ϕ)=(F(t),ϕ)∀ϕ∈V,
with F:=-J̇a(t)-ɛψ′′(t)E(0) known and given intial data.We require a slight adjustment to Lemma4.1 before giving the stability estimate.Lemma 4.4 (coercivity of creep).
If Assumptions2.1 hold, then
(4.22)∫0t(ψ′(0)w(s)+∫0sψ′′(s-ξ)w(ξ)dξ,w(s))ds⩾ψ′(t)∥w∥L2(0,t;L2(Ω))2
for all t∈I.Proof.
Since∥-ψ′′*∥w∥0∥L2(0,t)⩽∥ψ′′∥L1(0,t)∥w∥L2(0,t;L2(Ω)) and ψ′(0)-∥ψ′′∥L1(0,t)=ψ′(0)+∫0tψ′′(s)ds=ψ′(t) this is similar to Lemma 4.1.Theorem 4.5 (stability for theE equation).
If Assumptions2.1 hold and if μ>0, ɛ>0 and σ⩾0 are real numbers then for every t∈I we have first,
(4.23)μɛ∥Ė(t)∥02+∥∇×E(t)∥02+μ(σ+ɛψ′(t))∥Ė∥L2(0,t;L2(Ω))2⩽μɛ∥Ė(0)∥02+∥∇×E(0)∥02+μ∥F∥L2(0,t;L2(Ω))2(σ+ɛψ′(t)),
and second,
(4.24)μɛ2∥Ė∥L∞(0,t;L2(Ω))2+∥∇×E(t)∥02+2μ(σ+ɛψ′(t))∥Ė∥L2(0,t;L2(Ω))2⩽2μɛ∥Ė(0)∥02+2∥∇×E(0)∥02+8μɛ∥F∥L1(0,t;L2(Ω))2.Proof.
We chooseϕ=2Ė(t) in (4.21), integrate over (0,τ) and use Lemma 4.4 to get,
(4.25)μɛ∥Ė(τ)∥02+∥∇×E(τ)∥02+2μ(σ+ɛψ′(τ))∥Ė∥L2(0,τ;L2(Ω))2⩽μɛ∥Ė(0)∥02+∥∇×E(0)∥02+2μ∫0τ(F(t),Ė(t))dt.
Now by applying the Cauchy-Schwarz and Young inequalities in the following form:
(4.26)2μ∫0τ(F(t),Ė(t))dt⩽μ(σ+ɛψ′(τ))∥F∥L2(0,τ;L2(Ω))2+μ(σ+ɛψ′(τ))∥Ė∥L2(0,τ;L2(Ω))2,
we complete the proof of the first claim. For the second, we easily see that
(4.27)μɛ∥Ė∥L∞(0,τ;L2(Ω))2+∥∇×E(τ)∥02+2μ(σ+ɛψ′(τ))∥Ė∥L2(0,τ;L2(Ω))2⩽2μɛ∥Ė(0)∥02+2∥∇×E(0)∥02+4μ∫0τ|(F(t),Ė(t))|dt,
and with the Hölder and Young inequalities,
(4.28)4μ∫0τ|(F(t),Ė(t))|dt⩽8μɛ∥F∥L1(0,τ;L2(Ω))2+μɛ2∥Ė∥L∞(0,τ;L2(Ω))2,
we arrive at the second claim.Remark 4.6 (growth forσ=0).
Ifσ=0 then the fact that ψ′→0 as t→∞ 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 h, and that a suitable finite element space, Vh⊂V, 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 N intervals of equal width k so that I=I1∪I2∪⋯∪IN, where Ii:=(ti-1,ti] and ti=ik. 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” fieldW:=Ė. 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 E(0) and H(0) are given and (since ɛ and μ are assumed constant here) satisfy the divergence conditions in (1.2). We will need E(0) in the scheme below and also Ė(0) and we note that the latter quantity is available if H(0) is because, from (4.3), ɛĖ(0)=Ḋ(0)+φ′(0)D(0) and from (3.10) and (1.1) we see that D(0)=ɛE(0) and Ḋ(0)=∇×H(0)-σE(0)-Ja(0). Based on this, from now on we just assume that we have physically admissible (in terms of (1.2)) intial data such that(5.1)E(0):=Ĕ∈V,Ė(0):=W̆∈L2(Ω).
Our notation is wi:=w(ti) with(5.2)∂twi:=wi-wi-1k,∂t2wi=∂t∂twi=wi-2wi-1+wi-2k2,w¯i:=wi+wi-12,Δiw:=ẇi+ẇi-12-wi-wi-1k.
Notice that Δiw=ẇ¯i-∂twi and represents the trapezium rule.The numerical scheme is introduced next in Section5.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 i=1,2,…,N in turn, find a pair (Eih,Wih)∈Vh×Vh such that,(5.3)μɛ(∂tWih,ϕ)+(∇×E¯ih,∇×ϕ)+μσ(∂tEih,ϕ)+μɛ(ψ′(0)∂tEih+(ψ′′⋆∂tEh)i-1/2,ϕ)=μ(F¯i,ϕ)∀ϕ∈Vh,(5.4)(∂tEih,ζ)=(W¯ih,ζ)∀ζ∈Vh.
These are supplemented with discrete intial data, E0h≈Ĕ and W0h≈W̆, about which we will say more later (see later below (5.28)), and the discrete convolution, denoted by “⋆”, is defined by,(5.5)(ψ′′⋆X)i-1/2:=12∫ti-1tiψ′′(ti-s)dsXi+12∑n=1i-1∫tn-1tn(ψ′′(ti-s)+ψ′′(ti-1-s))dsXn
with empty sums set to zero and for any subset {Xn}n=1N⊂L2(Ω). Of course this is merely (ψ′′*X)¯i, the average of (ψ′′*X) at ti and ti-1, for a temporally piecewise constant X. Also, for mn:=ψ′(tn-1)-ψ′(tn+1)>0 (recall Assumptions 2.1), simple calculations give the alternative form(5.6)(ψ′′⋆X)i-1/2=-12(ψ′(0)-ψ′(k))Xi-12∑n=1i-1mi-nXn
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{mi}i=1N and {wi}i=1N of nonnegative real numbers we have for every j∈{2,3,…,N},
(5.7)(∑i=2j(∑n=1i-1mi-nwn)p)1/p⩽(∑n=1j-1mn)(∑i=1j-1wip)1/p,
for 1⩽p⩽∞ and the obvious interpretation if p=∞.Proof.
Obviously,∑n=1i-1mi-n=∑n=1i-1mn and so by two Hölder inequalities,
(5.8)∑i=2j(∑n=1i-1mi-nwn)p⩽(∑n=1j-1mn)p-1∑i=2j∑n=1i-1mi-nwnp.
We now interchange the order of summation with
(5.9)∑i=2j∑n=1i-1mi-nwnp=∑i=1j-1wip∑n=1j-imn,
apply Hölder's inequality again and take pth 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 Lemma4.4 but, before this, we need a bound on the discrete convolution.Lemma 5.2.
If Assumptions2.1 hold, then
(5.10)|∑i=1j((ψ′′⋆w)i-1/2,wi)|⩽(ψ′(0)-ψ¯j′)∑i=1j‖wi‖02
for every {wi}i=1N⊂L2(Ω) and every j∈{1,2,…,N}.Proof.
We have from the definition in (5.5) and (5.6) that
(5.11)∑i=1j((ψ′′⋆w)i-1/2,wi)=-12(ψ′(0)-ψ′(k))∑i=1j∥wi∥02-12∑i=2j∑n=1i-1mi-n(wn,wi),
and so applying Cauchy-Schwarz and then Lemma 5.1 gives
(5.12)|∑i=1j((ψ′′⋆w)i-1/2,wi)|⩽12(ψ′(0)-ψ′(k)+∑n=1j-1mn)∑i=1j∥wi∥02.
Now observe that ψ′(0)-ψ′(k)+∑n=1j-1mn=2ψ′(0)-2ψ¯j′.We can now state the discrete version of Lemma4.4.Lemma 5.3 (coercivity of discrete creep).
If Assumptions2.1 hold, then
(5.13)∑i=1j(ψ′(0)wi+(ψ′′⋆w)i-1/2,wi)⩾ψ¯j′∑i=1j∥wi∥02
for every {wi}i=1N⊂L2(Ω) and every j∈{1,2,…,N}.Proof.
Noting that(5.14)∑i=1j(ψ′(0)wi+(ψ′′⋆w)i-1/2,wi)⩾ψ′(0)∑i=1j∥wi∥02-|∑i=1j((ψ′′⋆w)i-1/2,wi)|,
the result follows from Lemma 5.2.Defining the norm(5.15)∥w∥∞,j:=max⁡{∥wi∥0:0⩽i⩽j}∀j⩽N,∀{wi}i=0N⊂L2(Ω),
we can now state the main result of this subsection.Theorem 5.4 (discrete well-posedness).
If Assumptions2.1 hold and we are given approximations (E0h,W0h)∈Vh×Vh to (Ĕ,W̆) then there exists a unique solution {(Eih,Wih)}i=1N⊂Vh×Vh to the discrete problem (5.3) and (5.4). Moreover,
(5.16)μɛ2∥Wh∥∞,j2+∥∇×Ejh∥02+2kμ(σ+ɛψ¯j′)∑i=1j∥∂tEih∥02⩽2μɛ∥W0h∥02+2∥∇×E0h∥02+8μɛ(k∑i=1j∥F¯i∥0)2
for each j∈{1,2,…,N}.Proof.
Since uniqueness implies existence for this finite dimensional problem, we need only prove the estimate and then takeE0h=W0h=F¯1=⋯=0. We select ϕ=2k∂tEih in (5.3) and ζ=2k∂tWih in (5.4), merge them using the common term 2k(∂tWih,∂tEih), and then sum the result over i=1,…,j. Noting the usual telescoping sum that arises from products like w¯i∂twi, we then use Lemma 5.3 to obtain
(5.17)μɛ∥Wjh∥02+∥∇×Ejh∥02+2kμ(σ+ɛψ¯j′)∑i=1j∥∂tEih∥02⩽μɛ∥W0h∥02+∥∇×E0h∥02+2kμ∑i=1j(F¯i,∂tEih).
It is clear from (5.4) that ∂tEih=W¯ih and so if i⩽j we easily conclude that ∥∂tEih∥0⩽∥Wh∥∞,j. Doubling the right-hand side to account for the ∥·∥∞,j norm a Young's inequality,
(5.18)|4kμ∑i=1j(F¯i,∂tEih)|⩽μɛ2∥W∥∞,j2+8μɛ(k∑i=1j∥F¯i∥0)2,
this completes the proof.The stability estimate given above is an exact analogue of the second estimate given for the continuous problem in Theorem4.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 theH(curl;Ω) projection as the linear continuous map, ℛ:V→Vh such that, for every u∈V, ℛu∈Vh is the unique solution to(5.19)(Ru,v)H(curl;Ω)=(u,v)H(curl;Ω)∀v∈Vh.
It is clear that the best approximation property,(5.20)∥u-Ru∥H(curl;Ω)⩽∥u-v∥H(curl;Ω)∀v∈Vh
is satisfied. Introducing also the L2(Ω) projection, x𝒫:L2(Ω)→Vh, through (𝒫u,v)=(u,v) for all v∈Vh, we then (as usual) split the errors according to(5.21)θi=Eih-RE(ti)∈Vh,λi=Wih-PW(ti)∈Vh,Θ(t)=E(t)-RE(t)∈V,Λ(t)=W(t)-PW(t)∈V,
so that Eih-E(ti)=θi-Θi and Wih-W(ti)=λi-Λi along with(5.22)(Λ(t),w)=0,(∇×Θ(t),∇×w)=-(Θ(t),w)∀w∈Vh.
The first result is a representation formula for the spatially discrete components of the error and for this we define, (5.23)Gi:=μɛΔiW+μ(σ+ɛψ′(0))(ΔiE+∂tΘi)-Θ¯i+μɛ(ψ′′⋆∂tΘ)i-1/2+μɛ(ψ′′*(Ė-∂tE))i¯,(5.24)Hi:=μɛΔiE+μɛ∂tΘi.Lemma 5.5.
For eachj∈{1,2,…,N}, we have
(5.25)2kμɛ∑i=1j(ψ′(0)∂tθi+(ψ′′⋆∂tθ)i-1/2,∂tθi)+μɛ∥λj∥02+∥∇×θj∥02+2kμσ∑i=1j∥∂tθi∥02=μɛ∥λ0∥02+∥∇×θ0∥02+2k∑i=1j((Gi,∂tθi)-(Hi,∂tλi)).Proof.
We average the continuous problem, (4.21), at the times ti and ti-1 and subtract the result from the discrete problem (5.3). Simultaneously adding and subtracting like terms in ℛE and 𝒫W 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 (∂tΛi,ϕ) is zero while (∇×Θ¯i,∇×ϕ) can be replaced by -(Θ¯i,ϕ). A simple manipulation provides
(5.26)(ψ′′*Ė)¯i-(ψ′′⋆∂tRE)i-1/2=(ψ′′⋆∂tΘ)i-1/2+(ψ′′*(Ė-∂tE))i¯
and so choosing ϕ=2k∂tθi generates the 𝒢i terms on the right and with ζ=2kμɛ∂tλi we can use (5.4), with (Λ¯i,ζ)=0, to get 2kμɛ(∂tθi,∂tλi)=kμɛ∂t∥λi∥02+2k(ℋi,∂tλi). Lastly we merge these and sum over i=1,2,…,j.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,(5.27)k∑i=1j(wi,∂tvi)=(wj,vj)-(w1,v0)-k∑i=2j(∂twi,vi-1)
and, secondly, that on taking ζ=μɛ∂tθi we can use (5.4) to get,(5.28)μɛ∥∂tθi∥0⩽μɛ∥λ∥∞,i+∥Hi∥0.
Thirdly, if the discrete intial data are defined by E0h:=ℛĔ and W0h:=𝒫W̆, then it follows that θ0=λ0=0.We can now give aformal  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 Assumptions2.1 hold and the discrete intial data are determined by E0h:=ℛĔ and W0h:=𝒫W̆, we have
(5.29)μɛ2‖Wj-Wjh‖02+‖∇×(Ej-Ejh)‖02+2kμ(σ+ɛψ¯j′)∑i=1j‖∂t(Ei-Eih)‖02⩽μɛ∥Λj∥02+2∥Θj∥H(curl;Ω)2+4kμ(σ+ɛψ¯j′)∑i=1j‖∂tΘi‖02+48μɛ‖Hj‖02+48μɛ(k∑i=2j‖∂tHi‖0)2+48μɛ(k∑i=1j‖Gi‖0)2+8kμϵ∑i=1j‖Gi‖0‖Hi‖0
for every j∈{1,2,…,N}, where the 𝒢i and ℋi are given in (5.23) and (5.24).Proof.
Lemmas5.5 and 5.3 give
(5.30)μɛ∥λ∥∞,j2+‖∇×θj‖02+2kμ(σ+ɛψ¯j′)∑i=1j∥∂tθi∥02⩽|4k∑i=1j((Gi,∂tθi)-(Hi,∂tλi))|
and using (5.27) and Young's inequality then gives,
(5.31)|4k∑i=1j(Hi,∂tλi)|⩽4μɛϵ∥λ∥∞,j2+2ϵμɛ∥Hj∥02+2ϵμɛ(k∑i=2j∥∂tHi∥0)2
for all ϵ>0. From (5.28), we find
(5.32)|4k∑i=1j(Gi,∂tθi)|⩽4kμɛ∑i=1j∥Gi∥0∥Hi∥0+2μɛϵ∥λ∥∞,j2+2ϵμɛ(k∑i=1j∥Gi∥0)2
for all ϵ>0. In fact, choosing ϵ=12 in the last two estimates gives
(5.33)μɛ2∥λ∥∞,j2+∥∇×θj∥02+2kμ(σ+ɛψ¯j′)∑i=1j∥∂tθi∥02⩽24μɛ∥Hj∥02+24μɛ(k∑i=2j∥∂tHi∥0)2+24μɛ(k∑i=1j∥Gi∥0)2+4kμɛ∑i=1j∥Gi∥0∥Hi∥0,
and so by applying the triangle inequality in the form ∥a±b∥2⩽2∥a∥2+2∥b∥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.(5.34)∥v̇(s)-∂tvi∥0⩽k∥v̈∥L1(ti-1,ti;L2(Ω)),fors∈[ti-1,ti],(5.35)∥Δiv∥0⩽k2∥vttt∥L1(ti-1,ti;L2(Ω)),(5.36)∥∂tvi∥0⩽k-1∥v̇∥L1(ti-1,ti;L2(Ω)),(5.37)∥∂tΔiv∥0⩽k∥vtttt∥L1(ti-2,ti;L2(Ω)).
We will also need the following estimates for the continuous and discrete convolutions.Lemma 5.7 (error in convolutions).
If Assumptions2.1 hold then we have
(5.38)k∑i=1j‖(ψ′′⋆∂tv)i-1/2‖0⩽(ψ′(0)-ψ¯j′)∥vt∥L1(0,tj;L2(Ω)),k∑i=1j‖(ψ′′*(vt-∂tv))¯i‖0⩽(ψ′(0)-ψ¯j′)k2∥vtt∥L1(0,tj;L2(Ω)),
for every i∈{1,2,…,j} and every j∈{1,2,…,N}.Proof.
First, from the definition in (5.6), we get after using (5.36) that
(5.39)k∥(ψ′′⋆∂tv)i-1/2∥0⩽12(ψ′(0)-ψ′(k))∥vt∥L1(ti-1,ti;L2(Ω))+12∑n=1i-1mi-n∥vt∥L1(tn-1,tn;L2(Ω)),
and then summing over i=1,…,j, we use Lemma 5.1 and arrive at
(5.40)k∑i=1j∥(ψ′′⋆∂tv)i-1/2∥0⩽12(ψ′(0)-ψ′(k)+∑n=1j-1mn)∥vt∥L1(0,tj;L2(Ω))
which is simplified to give the first result. The proof of the second is similar and begins with,
(5.41)(ψ′′*(vt-∂tv))¯i=12∫ti-1tiψ′′(ti-s)(v̇(s)-∂tvi)ds+12∑n=1i-1∫tn-1tn(ψ′′(ti-s)+ψ′′(ti-1-s))(v̇(s)-∂tvn)ds.
From (5.34) we define wn:=k∥vtt∥L1(tn-1,tn;L2(Ω))⩾∥v̇(s)-∂tvn∥0 for s∈[tn-1,tn], and then after a little work, we get
(5.42)∥(ψ′′*(vt-∂tv))¯i∥0⩽12(ψ′(0)-ψ′(k))wi+12∑n=1i-1mi-nwn.
Again we sum over i and apply Lemma 5.1 to arrive at,
(5.43)k∑i=1j∥(ψ′′*(vt-∂tv))¯i∥0⩽k22(ψ′(0)-ψ′(k)+∑n=1j-1mn)∥vtt∥L1(0,tj;L2(Ω)).
The same simplification completes the proof.Corollary5.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 “a≲b” to mean that “a⩽Cb” where C is a generic constant. In our working C will depend only on μ,ɛ,ψ′(0)>0 and σ⩾0 (but not on T or E).Lemma 5.8 (auxiliary bounds for𝒢i and ℋi).
If Assumptions2.1 hold then for every j∈{1,2,…,N},
(5.44)48μɛ∥Hj∥02+48μɛ(k∑i=2j∥∂tHi∥0)2+48μɛ(k∑i=1j∥Gi∥0)2+8kμɛ∑i=1j∥Gi∥0∥Hi∥0≲k4(∥Ettt∥L∞(0,tj;L2(Ω))2+∥Ett∥W12(0,tj;L2(Ω))2)+(∑i=1jk∥Θ¯i∥0)2+∥Θt∥L∞(0,tj;L2(Ω))2+∥Θt∥W11(0,tj;L2(Ω))2.Proof.
First, from (5.24), (5.35), and (5.36)
(5.45)∥Hi∥0⩽μɛk2∥Ettt∥L1(ti-1,ti;L2(Ω))+μɛk∥Θ̇∥L1(ti-1,ti;L2(Ω)),
and with a Young's inequality, it follows from this that
(5.46)48μɛ∥Hj∥02+8kμɛ∑i=1j∥Gi∥0∥Hi∥0≲μɛk4∥Ettt∥L∞(0,tj;L2(Ω))2+(k2∥Ettt∥L∞(0,tj;L2(Ω))+∥Θt∥L∞(0,tj;L2(Ω)))2+(∑i=1jk∥Gi∥0)2+μɛ∥Θt∥L∞(0,tj;L2(Ω))2,
where the hidden constant is merely a positive real number. Now, with the observation
(5.47)∂t2Θi=1k∂t∫ti-1tiΘ̇(s)ds=1k2∫ti-1tiΘ̇(s)-Θ̇(s-k)ds=1k2∫ti-1ti∫s-ksΘ̈(η)dηds
we see that
(5.48)∥∂t2Θi∥0⩽1k∫ti-2ti∥Θ̈(η)∥0dη.
Using this with (5.24) and (5.37) then gives
(5.49)∥∂tHi∥0⩽μɛk∥Etttt∥L1(ti-2,ti;L2(Ω))+μɛk∥Θtt∥L1(ti-2;ti;L2(Ω))
which, in turn, produces
(5.50)48μɛ(k∑i=2j∥∂tHi∥0)2⩽384μɛ(k4∥Etttt∥L1(0,tj;L2(Ω))2+∥Θtt∥L1(0,tj;L2(Ω))2).
The next task is to estimate the terms in ∥𝒢i∥0. To begin we get from (5.23), (5.35), and (5.36) that
(5.51)∥Gi∥0⩽μɛ∥ΔiW∥0+μ(σ+ɛψ′(0))(∥ΔiE∥0+∥∂tΘi∥0)+∥Θ¯i∥0+μɛ∥(ψ′′⋆∂tΘ)i-1/2∥0+μɛ∥(ψ′′*(Et-∂tE))¯i∥0,⩽μɛk2∥Wttt∥L1(ti-1,ti;L2(Ω))+μ(σ+ɛψ′(0))k2∥Ettt∥L1(ti-1,ti;L2(Ω))+μ(σ+ɛψ′(0))k∥Θ̇i∥L1(ti-1,ti;L2(Ω))+∥Θ¯i∥0+μɛ∥(ψ′′⋆∂tΘ)i-1/2∥0+μɛ∥(ψ′′*(Et-∂tE))¯i∥0.
Introducing the results from Lemma 5.7 into this provides,
(5.52)k∑i=1j∥Gi∥0≲k2∥Ett∥W12(0,tj;L2(Ω))+∥Θt∥L1(0,tj;L2(Ω))+∑i=1jk∥Θ¯i∥0,
where C depends only upon μ,ɛ,ψ′(0)>0 and σ⩾0 but not on T or E. Squaring this yields,
(5.53)(k∑i=1j∥Gi∥0)2≲k4∥Ett∥W12(0,tj;L2(Ω))2+∥Θt∥L1(0,tj;L2(Ω))2+(∑i=1jk∥Θ¯i∥0)2.
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,Vh, 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 L2(Ω) bounds) curl-conforming Nédélec spaces in [30], we find the following approximation result. Let the mesh be tetrahedral and regular; then if E∈Hℓ+1(Ω) there is an interpolator Π:V→Vh such that(5.54)∥E-ΠE∥0+∥∇×(E-ΠE)∥0⩽Chl∥E∥Hl+1(Ω).
Hence, in view of (5.20) we have(5.55)‖∂nΘ∂tn‖0+‖∇×∂nΘ∂tn‖0⩽Chl‖∂nE∂tn‖Hl+1(Ω)
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 E∈W14(0,T;L2(Ω))∩W12(0,T;Hℓ+1(Ω))∩W∞1(0,T;Hℓ+1(Ω)) with Ĕ,W̆∈Hℓ+1(Ω) then, under the same conditions as in Theorem 5.6, we have
(5.56)∥Wj-Wjh∥0+∥∇×(Ej-Ejh)∥0+((σ+ɛψ¯j′)∑i=1jk‖∂t(Ei-Eih)‖02)1/2≲hl(1+tj+((σ+ɛψ¯j′)tj)1/2)∥E∥W∞1(0,tj;Hl+1(Ω))+hl∥Et∥W11(0,tj;Hl+1(Ω))+k2(∥Ettt∥L∞(0,tj;L2(Ω))+∥Ett∥W12(0,tj;L2(Ω)))
for every j∈{1,2,…,N}.Proof.
First use Lemma5.8 to estimate the terms in 𝒢i and ℋi in Theorem 5.6, and remove the squares by using norm equivalence. This gives
(5.57)∥Wj-Wjh∥0+‖∇×(Ej-Ejh)‖0+((σ+ɛψ¯j′)∑i=1jk‖∂t(Ei-Eih)‖02)1/2≲∥Λj∥0+∥Θ(tj)∥H(curl;Ω)+∥Θt∥L∞(0,tj;L2(Ω))+∥Θt∥W11(0,tj;L2(Ω))+∑i=1jk∥Θ¯i∥0+((σ+ɛψ′¯j)∑i=1jk∥∂tΘi∥02)1/2+k2(∥Ettt∥L∞(0,tj;L2(Ω))+∥Ett∥W12(0,tj;L2(Ω))).
Now ∥Λj∥0⩽∥W-ΠW∥L∞(0,tj;L2(Ω))≲hℓ∥E∥W∞1(0,tj;Hℓ+1(Ω)) by the best approximation property of the L2(Ω) projection and with (5.20) and (5.36), we get
(5.58)∥Θ(tj)∥H(curl;Ω)+∥Θt∥L∞(0,tj;L2(Ω))+∥Θt∥W11(0,tj;L2(Ω))+∑i=1jk∥Θ¯i∥0+((σ+ɛψ¯j′)∑i=1jk∥∂tΘi∥02)1/2≲(1+tj)∥Θ∥L∞(0,tj;H(curl;Ω))+∥Θt∥W11(0,tj;L2(Ω))+(1+((σ+ɛψ¯j′)tj)1/2)∥Θ∥W∞1(0,tj;L2(Ω)).
Applying (5.55) and the triangle inequality then completes the proof.


## 5.1. The Numerical Scheme

Our numerical approximation to the electric field wave equation, (4.21), is the problem: for each i=1,2,…,N in turn, find a pair (Eih,Wih)∈Vh×Vh such that,(5.3)μɛ(∂tWih,ϕ)+(∇×E¯ih,∇×ϕ)+μσ(∂tEih,ϕ)+μɛ(ψ′(0)∂tEih+(ψ′′⋆∂tEh)i-1/2,ϕ)=μ(F¯i,ϕ)∀ϕ∈Vh,(5.4)(∂tEih,ζ)=(W¯ih,ζ)∀ζ∈Vh.
These are supplemented with discrete intial data, E0h≈Ĕ and W0h≈W̆, about which we will say more later (see later below (5.28)), and the discrete convolution, denoted by “⋆”, is defined by,(5.5)(ψ′′⋆X)i-1/2:=12∫ti-1tiψ′′(ti-s)dsXi+12∑n=1i-1∫tn-1tn(ψ′′(ti-s)+ψ′′(ti-1-s))dsXn
with empty sums set to zero and for any subset {Xn}n=1N⊂L2(Ω). Of course this is merely (ψ′′*X)¯i, the average of (ψ′′*X) at ti and ti-1, for a temporally piecewise constant X. Also, for mn:=ψ′(tn-1)-ψ′(tn+1)>0 (recall Assumptions 2.1), simple calculations give the alternative form(5.6)(ψ′′⋆X)i-1/2=-12(ψ′(0)-ψ′(k))Xi-12∑n=1i-1mi-nXn
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{mi}i=1N and {wi}i=1N of nonnegative real numbers we have for every j∈{2,3,…,N},
(5.7)(∑i=2j(∑n=1i-1mi-nwn)p)1/p⩽(∑n=1j-1mn)(∑i=1j-1wip)1/p,
for 1⩽p⩽∞ and the obvious interpretation if p=∞.Proof.
Obviously,∑n=1i-1mi-n=∑n=1i-1mn and so by two Hölder inequalities,
(5.8)∑i=2j(∑n=1i-1mi-nwn)p⩽(∑n=1j-1mn)p-1∑i=2j∑n=1i-1mi-nwnp.
We now interchange the order of summation with
(5.9)∑i=2j∑n=1i-1mi-nwnp=∑i=1j-1wip∑n=1j-imn,
apply Hölder's inequality again and take pth 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 Lemma4.4 but, before this, we need a bound on the discrete convolution.Lemma 5.2.
If Assumptions2.1 hold, then
(5.10)|∑i=1j((ψ′′⋆w)i-1/2,wi)|⩽(ψ′(0)-ψ¯j′)∑i=1j‖wi‖02
for every {wi}i=1N⊂L2(Ω) and every j∈{1,2,…,N}.Proof.
We have from the definition in (5.5) and (5.6) that
(5.11)∑i=1j((ψ′′⋆w)i-1/2,wi)=-12(ψ′(0)-ψ′(k))∑i=1j∥wi∥02-12∑i=2j∑n=1i-1mi-n(wn,wi),
and so applying Cauchy-Schwarz and then Lemma 5.1 gives
(5.12)|∑i=1j((ψ′′⋆w)i-1/2,wi)|⩽12(ψ′(0)-ψ′(k)+∑n=1j-1mn)∑i=1j∥wi∥02.
Now observe that ψ′(0)-ψ′(k)+∑n=1j-1mn=2ψ′(0)-2ψ¯j′.We can now state the discrete version of Lemma4.4.Lemma 5.3 (coercivity of discrete creep).
If Assumptions2.1 hold, then
(5.13)∑i=1j(ψ′(0)wi+(ψ′′⋆w)i-1/2,wi)⩾ψ¯j′∑i=1j∥wi∥02
for every {wi}i=1N⊂L2(Ω) and every j∈{1,2,…,N}.Proof.
Noting that(5.14)∑i=1j(ψ′(0)wi+(ψ′′⋆w)i-1/2,wi)⩾ψ′(0)∑i=1j∥wi∥02-|∑i=1j((ψ′′⋆w)i-1/2,wi)|,
the result follows from Lemma 5.2.Defining the norm(5.15)∥w∥∞,j:=max⁡{∥wi∥0:0⩽i⩽j}∀j⩽N,∀{wi}i=0N⊂L2(Ω),
we can now state the main result of this subsection.Theorem 5.4 (discrete well-posedness).
If Assumptions2.1 hold and we are given approximations (E0h,W0h)∈Vh×Vh to (Ĕ,W̆) then there exists a unique solution {(Eih,Wih)}i=1N⊂Vh×Vh to the discrete problem (5.3) and (5.4). Moreover,
(5.16)μɛ2∥Wh∥∞,j2+∥∇×Ejh∥02+2kμ(σ+ɛψ¯j′)∑i=1j∥∂tEih∥02⩽2μɛ∥W0h∥02+2∥∇×E0h∥02+8μɛ(k∑i=1j∥F¯i∥0)2
for each j∈{1,2,…,N}.Proof.
Since uniqueness implies existence for this finite dimensional problem, we need only prove the estimate and then takeE0h=W0h=F¯1=⋯=0. We select ϕ=2k∂tEih in (5.3) and ζ=2k∂tWih in (5.4), merge them using the common term 2k(∂tWih,∂tEih), and then sum the result over i=1,…,j. Noting the usual telescoping sum that arises from products like w¯i∂twi, we then use Lemma 5.3 to obtain
(5.17)μɛ∥Wjh∥02+∥∇×Ejh∥02+2kμ(σ+ɛψ¯j′)∑i=1j∥∂tEih∥02⩽μɛ∥W0h∥02+∥∇×E0h∥02+2kμ∑i=1j(F¯i,∂tEih).
It is clear from (5.4) that ∂tEih=W¯ih and so if i⩽j we easily conclude that ∥∂tEih∥0⩽∥Wh∥∞,j. Doubling the right-hand side to account for the ∥·∥∞,j norm a Young's inequality,
(5.18)|4kμ∑i=1j(F¯i,∂tEih)|⩽μɛ2∥W∥∞,j2+8μɛ(k∑i=1j∥F¯i∥0)2,
this completes the proof.The stability estimate given above is an exact analogue of the second estimate given for the continuous problem in Theorem4.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 theH(curl;Ω) projection as the linear continuous map, ℛ:V→Vh such that, for every u∈V, ℛu∈Vh is the unique solution to(5.19)(Ru,v)H(curl;Ω)=(u,v)H(curl;Ω)∀v∈Vh.
It is clear that the best approximation property,(5.20)∥u-Ru∥H(curl;Ω)⩽∥u-v∥H(curl;Ω)∀v∈Vh
is satisfied. Introducing also the L2(Ω) projection, x𝒫:L2(Ω)→Vh, through (𝒫u,v)=(u,v) for all v∈Vh, we then (as usual) split the errors according to(5.21)θi=Eih-RE(ti)∈Vh,λi=Wih-PW(ti)∈Vh,Θ(t)=E(t)-RE(t)∈V,Λ(t)=W(t)-PW(t)∈V,
so that Eih-E(ti)=θi-Θi and Wih-W(ti)=λi-Λi along with(5.22)(Λ(t),w)=0,(∇×Θ(t),∇×w)=-(Θ(t),w)∀w∈Vh.
The first result is a representation formula for the spatially discrete components of the error and for this we define, (5.23)Gi:=μɛΔiW+μ(σ+ɛψ′(0))(ΔiE+∂tΘi)-Θ¯i+μɛ(ψ′′⋆∂tΘ)i-1/2+μɛ(ψ′′*(Ė-∂tE))i¯,(5.24)Hi:=μɛΔiE+μɛ∂tΘi.Lemma 5.5.
For eachj∈{1,2,…,N}, we have
(5.25)2kμɛ∑i=1j(ψ′(0)∂tθi+(ψ′′⋆∂tθ)i-1/2,∂tθi)+μɛ∥λj∥02+∥∇×θj∥02+2kμσ∑i=1j∥∂tθi∥02=μɛ∥λ0∥02+∥∇×θ0∥02+2k∑i=1j((Gi,∂tθi)-(Hi,∂tλi)).Proof.
We average the continuous problem, (4.21), at the times ti and ti-1 and subtract the result from the discrete problem (5.3). Simultaneously adding and subtracting like terms in ℛE and 𝒫W 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 (∂tΛi,ϕ) is zero while (∇×Θ¯i,∇×ϕ) can be replaced by -(Θ¯i,ϕ). A simple manipulation provides
(5.26)(ψ′′*Ė)¯i-(ψ′′⋆∂tRE)i-1/2=(ψ′′⋆∂tΘ)i-1/2+(ψ′′*(Ė-∂tE))i¯
and so choosing ϕ=2k∂tθi generates the 𝒢i terms on the right and with ζ=2kμɛ∂tλi we can use (5.4), with (Λ¯i,ζ)=0, to get 2kμɛ(∂tθi,∂tλi)=kμɛ∂t∥λi∥02+2k(ℋi,∂tλi). Lastly we merge these and sum over i=1,2,…,j.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,(5.27)k∑i=1j(wi,∂tvi)=(wj,vj)-(w1,v0)-k∑i=2j(∂twi,vi-1)
and, secondly, that on taking ζ=μɛ∂tθi we can use (5.4) to get,(5.28)μɛ∥∂tθi∥0⩽μɛ∥λ∥∞,i+∥Hi∥0.
Thirdly, if the discrete intial data are defined by E0h:=ℛĔ and W0h:=𝒫W̆, then it follows that θ0=λ0=0.We can now give aformal  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 Assumptions2.1 hold and the discrete intial data are determined by E0h:=ℛĔ and W0h:=𝒫W̆, we have
(5.29)μɛ2‖Wj-Wjh‖02+‖∇×(Ej-Ejh)‖02+2kμ(σ+ɛψ¯j′)∑i=1j‖∂t(Ei-Eih)‖02⩽μɛ∥Λj∥02+2∥Θj∥H(curl;Ω)2+4kμ(σ+ɛψ¯j′)∑i=1j‖∂tΘi‖02+48μɛ‖Hj‖02+48μɛ(k∑i=2j‖∂tHi‖0)2+48μɛ(k∑i=1j‖Gi‖0)2+8kμϵ∑i=1j‖Gi‖0‖Hi‖0
for every j∈{1,2,…,N}, where the 𝒢i and ℋi are given in (5.23) and (5.24).Proof.
Lemmas5.5 and 5.3 give
(5.30)μɛ∥λ∥∞,j2+‖∇×θj‖02+2kμ(σ+ɛψ¯j′)∑i=1j∥∂tθi∥02⩽|4k∑i=1j((Gi,∂tθi)-(Hi,∂tλi))|
and using (5.27) and Young's inequality then gives,
(5.31)|4k∑i=1j(Hi,∂tλi)|⩽4μɛϵ∥λ∥∞,j2+2ϵμɛ∥Hj∥02+2ϵμɛ(k∑i=2j∥∂tHi∥0)2
for all ϵ>0. From (5.28), we find
(5.32)|4k∑i=1j(Gi,∂tθi)|⩽4kμɛ∑i=1j∥Gi∥0∥Hi∥0+2μɛϵ∥λ∥∞,j2+2ϵμɛ(k∑i=1j∥Gi∥0)2
for all ϵ>0. In fact, choosing ϵ=12 in the last two estimates gives
(5.33)μɛ2∥λ∥∞,j2+∥∇×θj∥02+2kμ(σ+ɛψ¯j′)∑i=1j∥∂tθi∥02⩽24μɛ∥Hj∥02+24μɛ(k∑i=2j∥∂tHi∥0)2+24μɛ(k∑i=1j∥Gi∥0)2+4kμɛ∑i=1j∥Gi∥0∥Hi∥0,
and so by applying the triangle inequality in the form ∥a±b∥2⩽2∥a∥2+2∥b∥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.(5.34)∥v̇(s)-∂tvi∥0⩽k∥v̈∥L1(ti-1,ti;L2(Ω)),fors∈[ti-1,ti],(5.35)∥Δiv∥0⩽k2∥vttt∥L1(ti-1,ti;L2(Ω)),(5.36)∥∂tvi∥0⩽k-1∥v̇∥L1(ti-1,ti;L2(Ω)),(5.37)∥∂tΔiv∥0⩽k∥vtttt∥L1(ti-2,ti;L2(Ω)).
We will also need the following estimates for the continuous and discrete convolutions.Lemma 5.7 (error in convolutions).
If Assumptions2.1 hold then we have
(5.38)k∑i=1j‖(ψ′′⋆∂tv)i-1/2‖0⩽(ψ′(0)-ψ¯j′)∥vt∥L1(0,tj;L2(Ω)),k∑i=1j‖(ψ′′*(vt-∂tv))¯i‖0⩽(ψ′(0)-ψ¯j′)k2∥vtt∥L1(0,tj;L2(Ω)),
for every i∈{1,2,…,j} and every j∈{1,2,…,N}.Proof.
First, from the definition in (5.6), we get after using (5.36) that
(5.39)k∥(ψ′′⋆∂tv)i-1/2∥0⩽12(ψ′(0)-ψ′(k))∥vt∥L1(ti-1,ti;L2(Ω))+12∑n=1i-1mi-n∥vt∥L1(tn-1,tn;L2(Ω)),
and then summing over i=1,…,j, we use Lemma 5.1 and arrive at
(5.40)k∑i=1j∥(ψ′′⋆∂tv)i-1/2∥0⩽12(ψ′(0)-ψ′(k)+∑n=1j-1mn)∥vt∥L1(0,tj;L2(Ω))
which is simplified to give the first result. The proof of the second is similar and begins with,
(5.41)(ψ′′*(vt-∂tv))¯i=12∫ti-1tiψ′′(ti-s)(v̇(s)-∂tvi)ds+12∑n=1i-1∫tn-1tn(ψ′′(ti-s)+ψ′′(ti-1-s))(v̇(s)-∂tvn)ds.
From (5.34) we define wn:=k∥vtt∥L1(tn-1,tn;L2(Ω))⩾∥v̇(s)-∂tvn∥0 for s∈[tn-1,tn], and then after a little work, we get
(5.42)∥(ψ′′*(vt-∂tv))¯i∥0⩽12(ψ′(0)-ψ′(k))wi+12∑n=1i-1mi-nwn.
Again we sum over i and apply Lemma 5.1 to arrive at,
(5.43)k∑i=1j∥(ψ′′*(vt-∂tv))¯i∥0⩽k22(ψ′(0)-ψ′(k)+∑n=1j-1mn)∥vtt∥L1(0,tj;L2(Ω)).
The same simplification completes the proof.Corollary5.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 “a≲b” to mean that “a⩽Cb” where C is a generic constant. In our working C will depend only on μ,ɛ,ψ′(0)>0 and σ⩾0 (but not on T or E).Lemma 5.8 (auxiliary bounds for𝒢i and ℋi).
If Assumptions2.1 hold then for every j∈{1,2,…,N},
(5.44)48μɛ∥Hj∥02+48μɛ(k∑i=2j∥∂tHi∥0)2+48μɛ(k∑i=1j∥Gi∥0)2+8kμɛ∑i=1j∥Gi∥0∥Hi∥0≲k4(∥Ettt∥L∞(0,tj;L2(Ω))2+∥Ett∥W12(0,tj;L2(Ω))2)+(∑i=1jk∥Θ¯i∥0)2+∥Θt∥L∞(0,tj;L2(Ω))2+∥Θt∥W11(0,tj;L2(Ω))2.Proof.
First, from (5.24), (5.35), and (5.36)
(5.45)∥Hi∥0⩽μɛk2∥Ettt∥L1(ti-1,ti;L2(Ω))+μɛk∥Θ̇∥L1(ti-1,ti;L2(Ω)),
and with a Young's inequality, it follows from this that
(5.46)48μɛ∥Hj∥02+8kμɛ∑i=1j∥Gi∥0∥Hi∥0≲μɛk4∥Ettt∥L∞(0,tj;L2(Ω))2+(k2∥Ettt∥L∞(0,tj;L2(Ω))+∥Θt∥L∞(0,tj;L2(Ω)))2+(∑i=1jk∥Gi∥0)2+μɛ∥Θt∥L∞(0,tj;L2(Ω))2,
where the hidden constant is merely a positive real number. Now, with the observation
(5.47)∂t2Θi=1k∂t∫ti-1tiΘ̇(s)ds=1k2∫ti-1tiΘ̇(s)-Θ̇(s-k)ds=1k2∫ti-1ti∫s-ksΘ̈(η)dηds
we see that
(5.48)∥∂t2Θi∥0⩽1k∫ti-2ti∥Θ̈(η)∥0dη.
Using this with (5.24) and (5.37) then gives
(5.49)∥∂tHi∥0⩽μɛk∥Etttt∥L1(ti-2,ti;L2(Ω))+μɛk∥Θtt∥L1(ti-2;ti;L2(Ω))
which, in turn, produces
(5.50)48μɛ(k∑i=2j∥∂tHi∥0)2⩽384μɛ(k4∥Etttt∥L1(0,tj;L2(Ω))2+∥Θtt∥L1(0,tj;L2(Ω))2).
The next task is to estimate the terms in ∥𝒢i∥0. To begin we get from (5.23), (5.35), and (5.36) that
(5.51)∥Gi∥0⩽μɛ∥ΔiW∥0+μ(σ+ɛψ′(0))(∥ΔiE∥0+∥∂tΘi∥0)+∥Θ¯i∥0+μɛ∥(ψ′′⋆∂tΘ)i-1/2∥0+μɛ∥(ψ′′*(Et-∂tE))¯i∥0,⩽μɛk2∥Wttt∥L1(ti-1,ti;L2(Ω))+μ(σ+ɛψ′(0))k2∥Ettt∥L1(ti-1,ti;L2(Ω))+μ(σ+ɛψ′(0))k∥Θ̇i∥L1(ti-1,ti;L2(Ω))+∥Θ¯i∥0+μɛ∥(ψ′′⋆∂tΘ)i-1/2∥0+μɛ∥(ψ′′*(Et-∂tE))¯i∥0.
Introducing the results from Lemma 5.7 into this provides,
(5.52)k∑i=1j∥Gi∥0≲k2∥Ett∥W12(0,tj;L2(Ω))+∥Θt∥L1(0,tj;L2(Ω))+∑i=1jk∥Θ¯i∥0,
where C depends only upon μ,ɛ,ψ′(0)>0 and σ⩾0 but not on T or E. Squaring this yields,
(5.53)(k∑i=1j∥Gi∥0)2≲k4∥Ett∥W12(0,tj;L2(Ω))2+∥Θt∥L1(0,tj;L2(Ω))2+(∑i=1jk∥Θ¯i∥0)2.
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,Vh, 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 L2(Ω) bounds) curl-conforming Nédélec spaces in [30], we find the following approximation result. Let the mesh be tetrahedral and regular; then if E∈Hℓ+1(Ω) there is an interpolator Π:V→Vh such that(5.54)∥E-ΠE∥0+∥∇×(E-ΠE)∥0⩽Chl∥E∥Hl+1(Ω).
Hence, in view of (5.20) we have(5.55)‖∂nΘ∂tn‖0+‖∇×∂nΘ∂tn‖0⩽Chl‖∂nE∂tn‖Hl+1(Ω)
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 E∈W14(0,T;L2(Ω))∩W12(0,T;Hℓ+1(Ω))∩W∞1(0,T;Hℓ+1(Ω)) with Ĕ,W̆∈Hℓ+1(Ω) then, under the same conditions as in Theorem 5.6, we have
(5.56)∥Wj-Wjh∥0+∥∇×(Ej-Ejh)∥0+((σ+ɛψ¯j′)∑i=1jk‖∂t(Ei-Eih)‖02)1/2≲hl(1+tj+((σ+ɛψ¯j′)tj)1/2)∥E∥W∞1(0,tj;Hl+1(Ω))+hl∥Et∥W11(0,tj;Hl+1(Ω))+k2(∥Ettt∥L∞(0,tj;L2(Ω))+∥Ett∥W12(0,tj;L2(Ω)))
for every j∈{1,2,…,N}.Proof.
First use Lemma5.8 to estimate the terms in 𝒢i and ℋi in Theorem 5.6, and remove the squares by using norm equivalence. This gives
(5.57)∥Wj-Wjh∥0+‖∇×(Ej-Ejh)‖0+((σ+ɛψ¯j′)∑i=1jk‖∂t(Ei-Eih)‖02)1/2≲∥Λj∥0+∥Θ(tj)∥H(curl;Ω)+∥Θt∥L∞(0,tj;L2(Ω))+∥Θt∥W11(0,tj;L2(Ω))+∑i=1jk∥Θ¯i∥0+((σ+ɛψ′¯j)∑i=1jk∥∂tΘi∥02)1/2+k2(∥Ettt∥L∞(0,tj;L2(Ω))+∥Ett∥W12(0,tj;L2(Ω))).
Now ∥Λj∥0⩽∥W-ΠW∥L∞(0,tj;L2(Ω))≲hℓ∥E∥W∞1(0,tj;Hℓ+1(Ω)) by the best approximation property of the L2(Ω) projection and with (5.20) and (5.36), we get
(5.58)∥Θ(tj)∥H(curl;Ω)+∥Θt∥L∞(0,tj;L2(Ω))+∥Θt∥W11(0,tj;L2(Ω))+∑i=1jk∥Θ¯i∥0+((σ+ɛψ¯j′)∑i=1jk∥∂tΘi∥02)1/2≲(1+tj)∥Θ∥L∞(0,tj;H(curl;Ω))+∥Θt∥W11(0,tj;L2(Ω))+(1+((σ+ɛψ¯j′)tj)1/2)∥Θ∥W∞1(0,tj;L2(Ω)).
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 Corollary5.9 does not depend upon T 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 ψ(t)~t as t→∞ (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 (6.1)ψ(t)=23.0459-21.3393e-t/17.67-0.7066e-t/0.9,φ(t)=0.0434+0.2433e-t/1.7969+0.7133e-t/0.3840;
and for methanol, (6.2)ψ(t)=13.5233-10.6487e-t/51.5-1.1183e-t/7.09-0.7563e-t/1.12,φ(t)=0.0739+0.1267e-t/14.2011+0.2502e-t/3.7374+0.5492e-t/0.5698.
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 (10-12s). The symbolic capability of matlab can be used to check the requirement that(6.3)(φ*ψ)(t)=∫0tφ(t-s)ψ(s)ds=t
and the values of t-φ*ψ 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.Graphical representation oft-φ*ψ for the numerically determined water and methanol creep-relaxation pairs discussed in Section 6.
(a)(b)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 (6.4)D(x,t)=ɛ0(x)E(x,t)+∫0tɛ1(x,t,s)E(x,s)ds,B(x,t)=μ0(x)H(x,t)+∫0tμ1(x,t,s)H(x,s)ds,J(x,t)=σ0(x)E(x,t)+∫0tσ1(x,t,s)E(x,s)ds,
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.



---

*Source: 923832-2010-12-27.xml*