Stability and Convergence Analysis of a Domain Decomposition FE/FD Method for Maxwell’s Equations in the Time Domain

Abstract

Stability and convergence analyses for the domain decomposition finite element/finite difference (FE/FD) method are presented. The analyses are developed for a semi-discrete finite element scheme for time-dependent Maxwell’s equations. The explicit finite element schemes in different settings of the spatial domain are constructed and a domain decomposition algorithm is formulated. Several numerical examples validate convergence rates obtained in the theoretical studies.

1. Introduction

In industrial applications, developing efficient numerical methods to solve the underlying partial differential equations (PDEs) requires simulations in higher dimensions and large computational domains. Such domains often consist of subdomains with different geometrical shapes and material properties. In this regard, domain decomposition (DD) methods, leading to efficient schemes, have gained considerable interest in the numerical analysis community, see [1,2]. A certain variant of the DD method, described below, is the subject of our current research on Maxwell’s equations, which is closely related to the schemes that were previously studied in [3,4]. More specifically, the present work is a further development of the DD hybrid finite element/finite difference (FE/FD) method for time-dependent Maxwell’s equations for an electric field in non-conductive media.

A stable, time-domain (TD), DD scheme for Maxwell’s equations was proposed in [5,6] and further verified in [7]. This method uses the Yee scheme [8] on the structured FD part of the mesh, and finite elements, possessing edge elements, on the unstructured part. Because of the presence of edge elements, implementing this method remains computationally expensive. A domain decomposition FE/FD method for time-dependent Maxwell’s equations for an electric field, assuming constant dielectric permittivity function in finite difference domain, was considered in [3]. The constant permittivity assumption simplifies the numerical schemes in both FE and FD domains and significantly reduces computational efforts to implement the considered DD scheme. A modified numerical scheme, energy estimate, and numerical verification of this method were presented in [4]. However, stability and convergence analysis with numerical implementations, in $L_2$- and $H^1$-norms of the developed FE/FD schemes, are not presented in the above studies. This work is meant to fill the gap.

More specifically, we present stability analyses for explicit schemes for FE (stability of FD is as in [5,6]), in the DD hybrid FE/FD method. The domain decomposition method (DDM) is constructed such that FEM and FDM coincide on the common, structured, overlapping layer between the two subdomains. The resulting domain decomposition approach at the overlapping layers can be viewed as a FE scheme which avoids instabilities at the interfaces. We decompose the computational domain such that FEM and FDM are used in different subdomains: FDM in simple geometry and FEM in the subdomain where more detailed information is needed about the structure of this subdomain. This allows the application of adaptive FEM in such subdomain, see, e.g., [3,9], and the references therein.

Reliability and convergence of the domain decomposition method studied in this work are evident for the solution of coefficient inverse problems (CIPs) in ${\mathbb{R}}^3$, see, e.g., [9], and the references therein. For the case of CIPs, the computational domain is split into subdomains such that a simple discretization scheme can be used in a large region and more refined discretization scheme is applied in smaller, but more critical, parts of the domain. In most algorithms for solution of electromagnetic CIPs, to determine the dielectric permittivity function inside a computational domain, a qualitative collection of experimental measurements at the boundary or in a neighborhood of it is necessary. In such cases, it is convenient to consider the numerical solution for Maxwell’s equations in different subdomains with a constant dielectric permittivity function in some, adequately chosen ones, and non-constant in other ones. For time-dependent Maxwell’s equations, the DD scheme of [4], analyzed in the present work, is used for the solution of different CIPs in order to determine the dielectric permittivity function in non-conductive media using simulated experimental data.

There are many other studies relevant to this work, e.g., in [10,11,12,13,14,15] and references therein. However, our analysis of DDM significantly differs from these works. In addition to optimal convergence of the constructed scheme, the novelty and scientific impact of the paper is in tackling coercivity (or lack of coercivity). More specifically, a challenging issue here is the lack of coercivity for the bilinear form in the presence of the time-derivative term. To circumvent this, a common approach is to split the bilinear form, separating the time-derivative term, and deriving the coercivity for the remaining sub-bilinear form. Then one may derive the convergence of the entire scheme with a time-derivative term.

An outline of this paper is as follows. In Section 2 and Section 3 we introduce the mathematical model, and briefly present the structure of domain decomposition FE/FD schemes and communication between the two schemes. In Section 4, we describe the domain decomposition FE/FD method for the numerical solution of Maxwell’s equations and set up adequate finite element and finite difference schemes. Section 5 is devoted to the stability analysis for the original bilinear form, (with a norm involving the time derivative term). In Section 6, we derive optimal a priori error estimates in finite element approximation for the semi-discrete (spatial discretization) problem on the entire domain, which is easily reduced to the FE subdomain, having known boundary data (the restriction of FD solution to the boundary of the FE subdomain). Finally, in our conclusion, Section 7, we present numerical implementations that justify the theoretical investigations in the paper.

2. The Mathematical Model

The Cauchy problem for the electric field $E\left(x,t\right)=\left(E_1,E_2,E_3\right)\left(x,t\right)$, $x\in {\mathbb{R}}^3$, $t\in [0,T]$, of Maxwell’s equations, under the assumptions that the dimensionless relative magnetic permeability of the medium is ${\mu }_r\equiv 1$, and the electric volume charges are zero, is given by

$$\begin{array}{cc}\hfill \frac{1}{c^2}{\epsilon }_r\left(x\right)\frac{{\partial }^{2}E}{\partial t^2}+\nabla \times \nabla \times E& =-{\mu }_0\sigma \left(x\right)\frac{\partial E}{\partial t},\hfill \\ \hfill \nabla ·\left(\epsilon E\right)& =0,\hfill \end{array}$$

(1)

where ${\epsilon }_r\left(x\right)=\epsilon \left(x\right)/{\epsilon }_0$ and $\sigma \left(x\right)$ are the dimensionless relative dielectric permittivity and electric conductivity functions, respectively. ${\epsilon }_0$ and ${\mu }_0$ are the permittivity and permeability of the free space, respectively, and $c=1/\sqrt{{\epsilon }_0{\mu }_0}$ is the speed of light in free space. In this paper we consider problem (1) in non-conductive media, i.e., $\sigma \equiv 0$, and hence, to begin with, we shall consider the following initial value problem

$$\left\{\begin{array}{cc}\hfill \epsilon \left(x\right)\frac{{\partial }^{2}E}{\partial t^2}+\nabla \times \nabla \times E& =0,x\in {\mathbb{R}}^3,t\in (0,T],\hfill \\ \hfill \nabla ·\left(\epsilon E\right)& =0,\mathrm{in}\Omega \times (0,T],\hfill \\ \hfill E(x,0)=f_0\left(x\right),\frac{\partial E}{\partial t}(x,0)& =f_1\left(x\right),x\in {\mathbb{R}}^3,t\in (0,T].\hfill \end{array}\right.$$

(2)

To solve problem (2) numerically, we need to consider it in a bounded domain $\Omega \subset {\mathbb{R}}^n,n=2,3$, with boundary $\partial \Omega$, and employ a split scheme on $\Omega$: a hybrid, finite element/finite difference scheme, a kind of domain decomposition described in Algorithm 1. In this approach, we divide the computational domain $\Omega$ into two subdomains, ${\Omega }_{FEM}$ and ${\Omega }_{FDM}$ such that $\Omega ={\Omega }_{FEM}\cup {\Omega }_{FDM}$   (note that ${\Omega }_{FEM}\cap {\Omega }_{FDM}\ne Ø$), and ${\Omega }_{FEM}$ is a subset of the convex hull of ${\Omega }_{FDM}$. The function $\epsilon \left(x\right)$ is assumed to be constant in ${\Omega }_{FDM}$, and bounded and smooth in ${\Omega }_{FEM}$. The communication between ${\Omega }_{FEM}$ and ${\Omega }_{FDM}$ is arranged using an overlapping mesh structure through a two-element thick layer around ${\Omega }_{FEM}$, as shown by the blue and green common boundaries in Figure 1. The blue boundary is the boundary of ${\Omega }_{FEM}$, the green boundary is the inner boundary of ${\Omega }_{FDM}$, and the red boundary is the outer boundary of $\Omega$ (which is also the boundary of the original domain $\Omega$).

Algorithm 1 Domain decomposition process for the hybrid FE/FD scheme

1: On the mesh ${\Omega }_{FDM}$, where FDM is used, update the Finite Difference (FD) solution. 2: On the mesh ${\Omega }_{FEM}$, where FEM is used, update the finite element (FE) solution. 3: Copy the FE solution obtained at nodes ${\omega }_{\diamond }$ (nodes on the green boundary of Figure 1) as a boundary condition on the inner boundary for the FD solution in ${\Omega }_{FDM}$. 4: Copy the FD solution obtained at nodes ${\omega }_{\mathrm{o}}$ (nodes on the blue boundary of Figure 1) as a boundary condition for the FE solution on $\partial {\Omega }_{FEM}$ of ${\Omega }_{FEM}$. 5: Apply boundary condition at $\partial \Omega$ at the red boundary of ${\Omega }_{FDM}$.

The key idea with such a decomposition is that it allows applying different numerical methods in different computational domains. For the numerical solution of (2) in ${\Omega }_{FDM}$, we use the finite difference method on a structured mesh. In ${\Omega }_{FEM}$, we use finite elements on a sequence of unstructured meshes $K_h=\left\{K\right\}$, with elements K consisting of triangles in ${\mathbb{R}}^2$ and tetrahedron’s in ${\mathbb{R}}^3$, both satisfying the minimal angle condition. This approach combines the flexibility of the finite element and the efficiency of the finite difference in terms of speed and memory usage and fits well for the reconstruction algorithm presented below.

Assumption 1.We assume that for some known constant $d>1$, the function $\epsilon \in C^2\left({\mathbb{R}}^3\right)$ satisfies   

$$\begin{array}{cc}\hfill \epsilon \left(x\right)\in \left[1,d\right],& \mathrm{in}\Omega \backslash {\Omega }_{FDM}\hfill \\ \hfill \epsilon \left(x\right)=1,& \mathrm{in}{\Omega }_{FDM}.\hfill \end{array}$$

(3)

As an immediate consequence of this assumption we have that

$$\epsilon \equiv 1\mathrm{in}{\Omega }_{FDM}\cap {\Omega }_{FEM}\mathrm{also}\epsilon \equiv 1\mathrm{in}\partial {\Omega }_{FDM}.$$

(4)

Assumption 1 has a physical meaning. The model problem (2) has a lot of applications in the field of coefficient inverse problems what we have pointed out in the introduction. One of the application of the model problem (2) is determination of dielectric permittivity function from measurements of the electric field at the boundary of the investigated domain, and this is a typical CIP. We refer for applications of such CIPs in [9,16,17,18,19] and references therein. This CIP is an ill-posed problem and needs several assumptions in order to be able to solve it, see details in [20,21,22]. One of the assumptions is the introduction of the set of admissible parameters for the dielectric permittivity function, which we define in the Condition (3). The lower bound in this set corresponds to the homogeneous media and the upper bound d is chosen as maximal values of the dielectric permittivity function and depends on the concrete physical problem. We refer to [16,18,23,24] where one can find values of dielectric permittivity function for different real-life applications with applications in medical imaging. Using these values can determine the upper and lower bounds in the set of admissible parameters for this function. See works [9,17,19] for examples of such sets for some real physical problems.

Note further that $\partial {\Omega }_{FEM}\subset {\Omega }_{FDM}$. Condition (3) on $\epsilon$ and the relation

$$\nabla \times \nabla \times E=\nabla (\nabla ·E)-\nabla ·(\nabla E),$$

(5)

together with divergence-free field E, make the equations in (2) independent of each other in ${\Omega }_{FEM}$ and ${\Omega }_{FDM}$ so that, in ${\Omega }_{FDM}$ (with $\epsilon =1$), we just need to solve the system of wave equations:

$$\frac{{\partial }^{2}E}{\partial t^2}-\Delta E=0,(x,t)\in {\Omega }_{FDM}\times (0,T].$$

(6)

Remark 1.

There are different discretizations of Maxwell’s equations in the time domain such as the FDTD Yee scheme [8,25], FE Lee–Madsen scheme [26], node-based least squares FE schemes [27,28], edge elements [29,30], and a number of discontinuous Galerkin FE schemes [31,32,33,34,35]. It is well known that, for stable implementation of the finite element solution of Maxwell’s equations, divergence-free edge elements are the most satisfactory ones from a theoretical point of view [29,30,36]. However, the edge elements are less attractive for the solution of time-dependent problems, since a linear system of equations should be solved at each time iteration step. On the contrary, P1 elements can be efficiently used in a fully explicit finite element scheme with lumped mass matrix, see, e.g., [37,38,39]. It is also well known that the numerical solution of Maxwell’s equations using nodal finite elements is often unstable and results in spurious oscillatory solutions [40,41]. There are a number of techniques to overcome such instabilities, see, e.g., [27,28,41,42,43].

In [44], a finite element analysis shows stability and consistency of the stabilized finite element method for the solution of (1) with $\sigma \left(x\right)=0$. In the current study we show stability and convergence for the combined FEM/FDM scheme, under the condition (3) on $\epsilon$ (the convergence would require an additional condition on $\epsilon$: namely (4)) where the stabilized FEM is used for the numerical solution of (2) in ${\Omega }_{FEM}$ and the usual FDM discretization of (6) is applied in ${\Omega }_{FDM}$.

We shall consider the following boundary and initial conditions:

$$\left\{\begin{array}{c}E(x,t)=0,x\in \partial \Omega ,t\in (0,T]\hfill \\ E(x,0)\in {\left[H^1(\Omega )\right]}^3,{\partial }_{t}E(x,0)\in \mathbf{H}(div,\Omega )\hfill \\ \nabla ·\left(\epsilon E(x,0)\right)=0,\nabla ·\left(\epsilon {\partial }_{t}E(x,0)\right)=0.\hfill \end{array}\right.$$

(7)

Recall that we assumed non-conductive media: $\sigma \equiv 0$. In the presence of electric conductivity, additional $\sigma$-terms appear in the equations. Then the $\sigma {\partial }_{t}E$-term adds an analytically advantageous, parabolic character to the equation.

Hence, in this note, we assume (3) and study the following stabilized version of the initial boundary value problem (2), with sufficiently smooth initial data $f_0$ and $f_1$:

$$\left\{\begin{array}{cc}\epsilon {\partial }_{tt}E-\Delta E+\nabla (\nabla ·(1-\epsilon )E)=0\hfill & \mathrm{in} \Omega \times (0,T),\hfill \\ E(·,0)=f_0(·),{\partial }_{t}E(·,0)=f_1(·)\hfill & \mathrm{in} \Omega ,\hfill \\ E=0\hfill & \mathrm{on} \partial \Omega \times (0,T).\hfill \end{array}\right.$$

(8)

3. The Structure of Domain Decomposition

We now describe the DD method between two subdomains, ${\Omega }_{FEM}$ and ${\Omega }_{FDM}$, where the FEM is used for computation of the solution in ${\Omega }_{FEM}$, and FDM is used in ${\Omega }_{FDM}$. Communication between ${\Omega }_{FEM}$ and ${\Omega }_{FDM}$ is achieved by allowing overlapping of both meshes across a two-element thick layer between $\partial {\Omega }_{FEM}$ and the inner boundary of ${\Omega }_{FDM}$, see Figure 1. The common nodes of both ${\Omega }_{FEM}$ and ${\Omega }_{FDM}$ domains belong to either of the following boundaries (see Figure 1):

• Nodes on the blue boundary ${\omega }_{\mathrm{o}}$—lie on the boundary $\partial {\Omega }_{FEM}$ of ${\Omega }_{FEM}$ and are interior to ${\Omega }_{FDM}$,
• Nodes on the green boundary ${\omega }_{\diamond }$—lie on the inner boundary $\partial {\Omega }_{FDM}$ of ${\Omega }_{FDM}$ and are interior to ${\Omega }_{FEM}$.

Then the main loop in time for the explicit hybrid FEM/FDM scheme, that solves (2), associated with appropriate boundary conditions, at each time step k, is described in Algorithm 1.

By (3), $\epsilon =1$ at the overlapping nodes between ${\Omega }_{FEM}$ and ${\Omega }_{FDM}$. Thus, FEM and FDM schemes coincide on the common, structured, overlapping layer. Hence, we avoid instabilities at interfaces.

4. Derivation of Computational Schemes

In this section, we construct a combined, finite element/finite difference scheme to solve the model problem (8). In this approach, first, we present the finite element scheme to solve (8) in entire $\Omega$. This induces the finite element scheme in ${\Omega }_{FEM}$. Then, we derive the finite difference scheme in ${\Omega }_{FDM}$, when the domain decomposition FE/FD structure is applied to solve (8) in $\Omega$.

In what follows, C will be a generic constant independent of all parameters, unless otherwise specifically specified, and not necessarily the same at each occurrence.

Remark 2.

The computational schemes derived in this section are explicit and therefore, for their convergence, the CFL condition below (see, e.g., [44]) should hold

$$\tau \le \frac{h}{\eta },\eta =C\sqrt{{1+3\parallel \epsilon -1\parallel }_{\infty }}.$$

(9)

Here C is a mesh independent constant, τ is the time step, and h is the mesh size.

In the sequel, we denote the inner product of ${\left[L^2(\Omega )\right]}^M,M\in \{1,2,3\},$ by $(·,·)$, and the corresponding $L_2$-norm by $\Vert ·\Vert$. The scalar inner product in ${\left[L^2\left({\Omega }_{FEM}\right)\right]}^M$ we denote by ${(·,·)}_{{\Omega }_{FEM}}$, and its associated $L_2$-norm by ${\parallel ·\parallel }_{{\Omega }_{FEM}}$. Further, we let $\partial {\Omega }_{FEM}$ be the boundary of ${\Omega }_{FEM}$, $\partial {\Omega }_{FDM}$ the inner boundary of ${\Omega }_{FDM}$, $\partial \Omega$ the outer boundary of ${\Omega }_{FDM}$, and ${\langle ·,·\rangle }_{\partial \Omega }$ denotes scalar product at the boundary $\partial \Omega$.

4.1. Finite Element Discretization in $\Omega$

First, we derive the finite element scheme to solve the model problem (8) in whole $\Omega$. The restriction of this FE solution to the inner $\partial {\Omega }_{FDM}$ is the approximation for ${E|}_{\partial }{\Omega }_{FDM}$. Next, we discretize ${\Omega }_{{FEM}_T}={\Omega }_{FEM}\times (0,T)$ in two steps: (i) the spatial discretization using a partition ${\mathcal{T}}_h$ of ${\Omega }_{FEM}$ into elements K, where $h=h\left(x\right)$ is a mesh function defined as ${h|}_K=h_K$, representing the local diameter of elements. We also denote by $\partial {\mathcal{T}}_h$ a partition of the boundary $\partial {\Omega }_{FEM}$ into boundaries $\partial K$ of the elements K such that in $2D$, two of the vertices, and in $3D$, three of the vertices of these elements belong to $\partial {\Omega }_{FEM}$. As usual, we also assume a minimal angle condition on elements K in $K_h$, see details in [45]. (ii) As for temporal discretization, we let $J_{\tau }$ be a uniform partition of the time interval $(0,T)$ into N subintervals $J=(t_{k-1},t_k]$ of length $\tau =T/N.$ Our implementations are for fully discrete, space-time scheme. However, temporal discretization (ii), being well-studied in the literature, is not considered in the analysis here. To formulate the finite element method for (8) in $\Omega$, we introduce the finite element space $W_h^E(\Omega )$ for each component of the electric field E defined by

$$W_h^E(\Omega ):=\{w\in H^1(\Omega ):w{|}_K\in P_1\left(K\right),\forall K\in {\mathcal{T}}_h\},$$

where $P_1\left(K\right)$ denote the set of piecewise-linear functions on K. Setting ${\mathbf{W}}_h^E(\Omega ):={\left[W_h^E(\Omega )\right]}^3$, and because of the homogeneous Dirichlet boundary data in (8) we need to choose the test function space as

$${\mathbf{W}}_{h,0}^E:=\{\mathbf{v}\in {\mathbf{W}}_h^E:\mathbf{v}=0,\mathrm{on} \partial \Omega \}.$$

(10)

Then, recalling (5) and the Dirichlet boundary condition in (8) the spatial semi-discrete problem in $\Omega$ reads:

Find$E_h\in {\mathbf{W}}_{h,0}^E(\Omega )$such that$\forall \mathbf{v}\in {\mathbf{W}}_{h,0}^E(\Omega )$,

$$\begin{array}{cc}\hfill \left(\epsilon {\partial }_{tt}{E}_h,\mathbf{v}\right)& +(\nabla E_h,\nabla \mathbf{v})+(\nabla ·\left(\epsilon E_h\right),\nabla ·\mathbf{v})-(\nabla ·E_h,\nabla ·\mathbf{v})\hfill \\ \hfill & -{\langle {\partial }_n{E}_h,\mathbf{v}\rangle }_{\partial \Omega }-{\langle n·(\nabla ·(\epsilon -1)E_h),\mathbf{v}\rangle }_{\partial \Omega }:=\sum _{j=1}^6{T}_j=0,\hfill \\ \hfill & E_h(·,0)=f_{0,h}(·)and{\partial }_t{E}_h(·,0)=f_{1,h}(·)\mathrm{in} \Omega .\hfill \end{array}$$

(11)

where $f_{0,h}$ and $f_{1,h}$ are suitable approximations for the initial data, specified below. Further, we note that

$$T_3=(\nabla ·\left(\epsilon E_h\right),\nabla ·\mathbf{v})=(\nabla \epsilon ·E_h,\nabla ·\mathbf{v})+(\epsilon \nabla ·E_h,\nabla ·\mathbf{v}),$$

implies

$$\begin{array}{cc}\hfill T_3+T_4& =(\nabla ·\left(\epsilon E_h\right),\nabla ·\mathbf{v})-(\nabla ·E_h,\nabla ·\mathbf{v})\hfill \\ \hfill & =(\nabla \epsilon ·E_h,\nabla ·\mathbf{v})+\left((\epsilon -1)\nabla ·E_h,\nabla ·\mathbf{v}\right).\hfill \end{array}$$

(12)

Recalling (10), vanishing test functions at the boundary yields $T_5=T_6\equiv 0$ ( $\epsilon =1$ on $\partial \Omega$ yields $T_6=0$), and the final weak formulation for the semi-discrete problem in $\Omega$ is:

Find$E_h\in {\mathbf{W}}_{h,0}^E(\Omega )$such that$\forall \mathbf{v}\in {\mathbf{W}}_{h,0}^E(\Omega )$,

$$B_{\Omega }(E_h,\mathbf{v}):=\left(\epsilon {\partial }_{tt}{E}_h,\mathbf{v}\right)+(\nabla E_h,\nabla \mathbf{v})+(\nabla \epsilon ·E_h,\nabla ·\mathbf{v})+((\epsilon -1)\nabla ·E_h,\nabla ·\mathbf{v})=0.$$

(13)

For the reflexive, inhomogeneous boundary condition, see the FE scheme for ${\Omega }_{FEM}$.

To get a fully discrete scheme for (8) we apply time discretization to (11). Then $E_h\left(k\tau \right)$, denoted by $E_h^k$, and the central difference scheme, for $k=1,\dots ,N-1$, yields

$$\begin{array}{c}\left(\epsilon \frac{E_h^{k+1}-2E_h^k+E_h^{k-1}}{{\tau }^2},\mathbf{v}\right)+(\nabla E_h^k,\nabla \mathbf{v})+(\nabla ·\left(\epsilon E_h^k\right),\nabla ·\mathbf{v})\hfill \\ -(\nabla ·E_h^k,\nabla ·\mathbf{v})=0\forall \mathbf{v}\in {\mathbf{W}}_{h,0}^E(\Omega ),\hfill \\ \hfill \\ E_h^0=f_{0,h}\mathrm{and}{E}_h^1=E_h^0+\tau f_{1,h}\mathrm{in} \Omega .\hfill \end{array}$$

(14)

Multiplying both sides of (14) by ${\tau }^2/\epsilon$ we get, $\forall \mathbf{v}\in {\mathbf{W}}_{h,0}^E(\Omega ),$ that

$$\begin{array}{c}\left(E_h^{k+1}-2E_h^k+E_h^{k-1},\mathbf{v}\right)+{\tau }^2(1/\epsilon \nabla E_h^k,\nabla \mathbf{v})+{\tau }^2(1/\epsilon \nabla ·\left(\epsilon E_h^k\right),\nabla ·\mathbf{v})\hfill \\ -{\tau }^2(1/\epsilon \nabla ·E_h^k,\nabla ·\mathbf{v})=0,\hfill \\ \hfill \\ E_h^0=f_{0,h}\mathrm{and}{E}_h^1=E_h^0+\tau f_{1,h},\mathrm{in} \Omega .\hfill \end{array}$$

(15)

Rearranging terms in (15) we get the following scheme: Given the initial data, $f_0$ and $f_1$, find $E_h\in {\mathbf{W}}_h^E(\Omega )$ such that $\forall \mathbf{v}\in {\mathbf{W}}_{h,0}^E(\Omega )$,

$$\begin{array}{cc}\hfill \left(E_h^{k+1},\mathbf{v}\right)& =\left(2E_h^k,\mathbf{v}\right)-\left(E_h^{k-1},\mathbf{v}\right)-{\tau }^2(1/\epsilon \nabla E_h^k,\nabla \mathbf{v})-{\tau }^2(1/\epsilon \nabla ·\left(\epsilon E_h^k\right),\nabla ·\mathbf{v})\hfill \\ \hfill & +{\tau }^2(1/\epsilon \nabla ·E_h^k,\nabla ·\mathbf{v})\hfill \\ \hfill \\ \hfill & E_h^0=f_{0,h}\mathrm{and}{E}_h^1=E_h^0+\tau f_{1,h}\mathrm{in} \Omega .\hfill \end{array}$$

(16)

4.2. Finite Element Discretization in ${\Omega }_{FEM}$

To solve the model problem (8) via the domain decomposition FE/FD method, we use the split $\Omega ={\Omega }_{FEM}\cup {\Omega }_{FDM}$, see Figure 1. Thus, in ${\Omega }_{FEM}$ we use FEM to solve the equation

$$\begin{array}{cc}\epsilon {\partial }_{tt}E-\Delta E+\nabla (\nabla ·(1-\epsilon )E)=0\hfill & \mathrm{in} {\Omega }_{FEM}\times (0,T),\hfill \\ E(·,0)=f_0(·),\mathrm{and} {\partial }_{t}E(·,0)=f_1(·)\hfill & \mathrm{in} {\Omega }_{FEM},\hfill \\ E=\left(E{|}_{\partial {\Omega }_{FEM}}\right)=g,\hfill & \mathrm{on} \partial {\Omega }_{FEM}\times (0,T),\hfill \\ \nabla ·\left(\epsilon E\right)=0\hfill & \mathrm{in} {\Omega }_{FEM}.\hfill \end{array}$$

(17)

Here, g corresponds to the restriction of the exact solution to $\partial {\Omega }_{FEM}$. Therefore test and trial functions are not vanishing at this boundary, hence the term corresponding to $T_5$ in (11) will appear in the weak formulation for ${\Omega }_{FEM}$.

To formulate the finite element method for (17) in ${\Omega }_{FEM}$, mimicking (11), we introduce the finite element space $W_h^E\left({\Omega }_{FEM}\right)$ for each component of the electric field E defined by

$$W_h^E\left({\Omega }_{FEM}\right):=\{w\in H^1\left({\Omega }_{FEM}\right):w{|}_K\in P_1\left(K\right),\forall K\in {\mathcal{T}}_h\}.$$

Setting ${\mathbf{W}}_h^E\left({\Omega }_{FEM}\right):={\left[W_h^E\left({\Omega }_{FEM}\right)\right]}^3$ we define ${\mathcal{P}}_{h}g$ to be the usual $L_2$-projection of g in $\partial {\Omega }_{FEM}$. Then, similar to the FE scheme for $\Omega$, we get the following finite element scheme for ${\Omega }_{FEM}$: Given approximate values $f_{0,h}$ and $f_{1,h}$ of $f_0$ and $f_1$, respectively, and ${\mathcal{P}}_{h}g$, find ${\tilde{E}}_h\in {\mathbf{W}}_h^E\left({\Omega }_{FEM}\right)$ such that

$$\begin{array}{cc}\hfill \left(\epsilon {\partial }_{tt}{\tilde{E}}_h,\mathbf{v}\right)& +(\nabla {\tilde{E}}_h,\nabla \mathbf{v})+(\nabla ·\left(\epsilon {\tilde{E}}_h\right),\nabla ·\mathbf{v})-(\nabla ·{\tilde{E}}_h,\nabla ·\mathbf{v})\hfill \\ \hfill & ={\langle n·(\nabla ·(\epsilon -1){\mathcal{P}}_{h}g),\mathbf{v}\rangle }_{\partial {\Omega }_{FEM}}+\langle {{\partial }_n(\mathcal{P}}_{h}g){,\mathbf{v}\rangle }_{\partial {\Omega }_{FEM}}\hfill \\ \hfill & =\{\epsilon \equiv 1\mathrm{on}\partial {\Omega }_{FEM}\}={\langle {\partial }_n\left({\mathcal{P}}_{h}g\right),\mathbf{v}\rangle }_{\partial {\Omega }_{FEM}},\forall \mathbf{v}\in {\mathbf{W}}_h^E\left({\Omega }_{FEM}\right)\hfill \\ \hfill & \hfill \\ \hfill & {\tilde{E}}_h(·,0)=f_{0,h}(·)\mathrm{and}{\partial }_t{\tilde{E}}_h(·,0)=f_{1,h}(·)\mathrm{in} {\Omega }_{FEM}\hfill \end{array}$$

(18)

where we set ${\tilde{E}}_h{{|}_{\partial {\Omega }_{FEM}}=E_h\Omega |}_{\partial {\Omega }_{FEM}}$, i.e., the restriction of the finite element solution on whole $\Omega$ to $\partial {\Omega }_{FEM}$, (which, here, contrary to the whole domain where $E_h{|}_{\partial \Omega }=0$, will be non-zero).

A corresponding fully discrete problem in ${\Omega }_{{FEM}_T}$ reads as follows:

Given $f_{0,h}$, $f_{1,h}$, ${\mathcal{P}}_{h}g$, $E_h^k$, and $E_h^{k-1}$; find $E_h^{k+1}$ such that $\forall \mathbf{v}\in {\mathbf{W}}_h^E\left({\Omega }_{FE{M}_T}\right)$,

$$\begin{array}{cc}\hfill \left({\tilde{E}}_h^{k+1},\mathbf{v}\right)=& \left(2{\tilde{E}}_h^k,\mathbf{v}\right)-\left({\tilde{E}}_h^{k-1},\mathbf{v}\right)-{\tau }^2(1/\epsilon \nabla {\tilde{E}}_h^k,\nabla \mathbf{v})+{\tau }^2(1/\epsilon \nabla ·{\tilde{E}}_h^k,\nabla ·\mathbf{v})\hfill \\ \hfill & -{\tau }^2(1/\epsilon \nabla ·\left(\epsilon {\tilde{E}}_h^k\right),\nabla ·\mathbf{v})+{\tau }^2{\langle {\partial }_n\left({\mathcal{P}}_{h}g\right)/\epsilon ,\mathbf{v}\rangle }_{\partial {\Omega }_{FEM}},\hfill \\ \hfill & \\ \hfill E_h^0=f_{0,h}& E_h^1=E_h^0+\tau f_{1,h}\mathrm{in} {\Omega }_{{FEM}_T}.\hfill \end{array}$$

(19)

4.3. Finite Difference Formulation

We recall now that from conditions (3) it follows that in ${\Omega }_{FDM}$ the function $\epsilon \left(x\right)=1$. This means that in ${\Omega }_{FDM}$ for the model problem (2) the forward problem will be

$$\begin{array}{cc}\hfill \frac{{\partial }^{2}E}{\partial t^2}-\Delta E=0& \mathrm{in} {\Omega }_{FDM}\times (0,T),\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill E(x,0)=f_0\left(x\right),E_t(x,0)=f_1\left(x\right)& \mathrm{in} {\Omega }_{FDM},\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill E=0& \mathrm{on} \partial \Omega \times (0,T),\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill {\partial }_{n}E={\partial }_n{E}_{FEM}& \mathrm{on} \partial {\Omega }_{FDM}\times (0,T).\hfill \end{array}$$

(23)

Using standard finite difference discretization of the Equation (20) in ${\Omega }_{FDM}$ we obtain the following explicit scheme for the solution of the forward problem:

$$E_{l,j,m}^{k+1}={\tau }^2\Delta E_{l,j,m}^k+2E_{l,j,m}^k-E_{l,j,m}^{k-1}.$$

(24)

In the system of equations above, $E_{l,j,m}^k$ is the solution at the time iteration step k at the discrete point $(l,j,m)$, $\tau$ is the time step, and $\Delta E_{l,j,m}^k$ is the discrete Laplacian.

Note that (24) is associated with the Neumann boundary condition above, where the Dirichlet boundary conditions $E=E_{FEM}$ can be considered as well.

5. Stability

5.1. Stability Estimate for Continuous Problem in $\Omega$

In this section, we derive stability estimates for the semi-discrete approximations. For stability in $\Omega$, these estimates are extensions of the stability approach derived for the continuous problem in [4]. As for the stability in ${\Omega }_{FEM}$, we get slightly different norms involving contributions corresponding to the reflexive boundary, $\partial {\Omega }_{FEM}$. We define the discrete version of a triple norm induced by the weak variational formulation of (8), where we also use the relation (12) (in ${\Omega }_{FDM},\epsilon =1$, therefore, the second term on the right-hand side in (12) is not present in here). Note that in the continuous case $\nabla ·\left(\epsilon E\right)=0$, (whereas, in general $\nabla ·\left(\epsilon E_h\right)\ne 0$), and the continuous variational formulation in whole $\Omega$, reads as:

Find$E\in {\mathbf{W}}_0^E(\Omega )$such that for all$\mathbf{v}\in {\mathbf{W}}_0^E(\Omega )$

$$\left\{\begin{array}{c}\left(\epsilon {\partial }_{tt}E,\mathbf{v}\right)+(\nabla E,\nabla \mathbf{v})+((\nabla \epsilon )·E,\nabla ·\mathbf{v})+((\epsilon -1)\nabla ·E,\nabla ·\mathbf{v})=0,\hfill \\ \\ E(·,0)=f_0(·)\mathrm{and}{\partial }_{t}E(·,0)=f_1(·)\mathrm{in} \Omega .\hfill \end{array}\right.$$

(25)

Remark 3.

We emphasize that the bilinear form induced by (25), containing the time derivative term, is not coercive (whereas, without the $\epsilon {\partial }_{tt}E$ term, both (25) and its spatially discrete version are coercive). Further $H^1(\Omega )$-conforming finite element may result in spurious oscillatory solutions. A remedy is through modifying the equation by adding a gauge constraint of Coulomb-type, see, e.g., [36,43]. This is supplied by the “zero”-term: $\nabla ·\left(\epsilon E\right)=0$, in (8), which we add in the continuous variational formulation in (11). This, however, is not necessarily true in the discrete forms, e.g., in (18), where most likely $\nabla ·\left(\epsilon E_h\right)\ne 0.$

Taking $\mathbf{v}={\partial }_{t}E$ in (25), (with the boundary condition $E=0$ on $\partial \Omega$), yields   

$$\begin{array}{cc}\hfill \left(\epsilon {\partial }_{tt}E,{\partial }_{t}E\right)& +(\nabla E,\nabla {\partial }_{t}E)+((\nabla \epsilon )·E,\nabla ·{\partial }_{t}E)\hfill \\ \hfill & +((\epsilon -1)\nabla ·E,\nabla ·{\partial }_{t}E)\equiv 0,\hfill \end{array}$$

(26)

which, due to the fact that $\epsilon$ is independent of t, can be rewritten as

$$\begin{array}{cc}\hfill \frac{1}{2}\frac{d}{dt}\left(\epsilon {\partial }_{t}E,{\partial }_{t}E\right)& +\frac{1}{2}\frac{d}{dt}\left(\nabla E,\nabla E\right)+\left((\nabla \epsilon )·E,\nabla ·{\partial }_{t}E\right)\hfill \\ \hfill & +\frac{1}{2}\frac{d}{dt}\left((\epsilon -1)\nabla ·E,\nabla ·E\right)\equiv I_1+I_2+I_3+I_4=0.\hfill \end{array}$$

(27)

Proposition 1.

Let $\Omega \subset {\mathbb{R}}^n,n=2,3$ be a bounded domain with piecewise linear boundary $\partial \Omega$. Let Assumption 1 for function ε hold. Assume that there exists a solution $E\in H^2\left({\Omega }_T\right)$ of problem (8). Then, Equation (8) has a unique solution $E\in H^2\left({\Omega }_T\right)$. Further let $f_1\in L_{\epsilon }^2(\Omega )$ and $f_0\in L_2(\Omega )\cap H_{\left|\nabla \epsilon \right|+\epsilon -1}^1(\Omega )$, then there is a constant C, independent of ε and t, such that, $\forall t\in (0,T]$, the following stability estimate holds true

$$\begin{array}{cc}\hfill {\left|\right|\left|E\right|\left|\right|}_{\epsilon }^2\left(t\right):=& {∥{\partial }_{t}E∥}_{\epsilon }^2\left(t\right)+{∥\nabla E∥}^2\left(t\right)+{∥\nabla ·E∥}_{\left|\nabla \epsilon \right|+\epsilon -1}^2\left(t\right)+{∥E∥}_{\left|\nabla \epsilon \right|}^2\left(t\right)\hfill \\ \hfill & +{\int }_0^t{∥E_t∥}_{\left|\nabla \epsilon \right|}^{2}ds+{\int }_0^t{∥\nabla ·E∥}_{\left|\nabla \epsilon \right|}^{2}ds\hfill \\ \hfill \le & C\left({∥f_1∥}_{\epsilon }^2+{∥\nabla f_0∥}^2+{∥\nabla ·f_0∥}_{\left|\nabla \epsilon \right|+\epsilon -1}^2+{∥f_0∥}_{\left|\nabla \epsilon \right|}^2\right).\hfill \end{array}$$

(28)

Proof. 

A version of the estimate (28) could be proved, as in [4], Theorem 4.1, by setting $s=1$ and $j\equiv 0$, where the nonsymmetric term $I_3$ is the source of challenge.

A rather involved proof of a similar result is derived in [46], using Nitsche’s penalty method. A direct approach using Fubini’s theorem (through a change of integration order in x and t, followed by partial integration in t), yields the estimate of the non-symmetric term $I_3$. This, however, is with the price of the appearance of the integral terms on the left-hand side, viz.

$${\int }_0^t{I}_{3}ds=-{\int }_{\Omega }\left({\int }_0^t((\nabla \epsilon )·E_t)(\nabla ·E)ds-{\left[(\nabla \epsilon ·E)(\nabla ·E)\right]}_0^t\right):=I_{3,1}+I_{3,2}.$$

(29)

Then, we can derive the following inequalities:

$$I_{3,1}\le \frac{1}{2}{\int }_0^t{∥E_t∥}_{\left|\nabla \epsilon \right|}^{2}ds+\frac{1}{2}{\int }_0^t{∥\nabla ·E∥}_{\left|\nabla \epsilon \right|}^{2}ds,$$

(30)

and

$$I_{3,2}\le \frac{1}{2}{∥E∥}_{\left|\nabla \epsilon \right|}^2\left(t\right)+\frac{1}{2}{∥\nabla ·E∥}_{\left|\nabla \epsilon \right|}^2\left(t\right)+\frac{1}{2}{∥E∥}_{\left|\nabla \epsilon \right|}^2\left(0\right)+\frac{1}{2}{∥\nabla ·E∥}_{\left|\nabla \epsilon \right|}^2\left(0\right).$$

(31)

Combining the contributions from the symmetric terms: $I_1,I_2$, and $I_4$ with (29)–(31), an adequate choice (large enough) of the constant C yields the desired result.  ☐

Below we translate this stability to the semi-discrete problem.

5.2. Stability Estimate for the Semi-Discrete Problem in $\Omega$

The stability for the semi-discrete problem in $\Omega$ is basically as in the continuous case above where all E:s are replaced by $E_h\in {\mathbf{W}}_{h,0}^E$, i.e., $E_h{|}_{\partial \Omega }\equiv 0$ and some relevant assumptions in the discrete data, viz.

Lemma 1.

Assume that the approximations for the initial-data $f_0$ and $f_1$: $f_{0,h}$ and $f_{1,h}$ satisfy the regularity conditions $f_{1,h}\in L_{\epsilon }^2(\Omega )$ and $f_{0,h}\in H^1(\Omega )\cap H_{\epsilon -1}^1(\Omega )$, then for each $t\in (0,T]$,

$$\left|\right||E_h{\left|\right||}_{\epsilon }^2\left(t\right)\le C\left({∥f_{1,h}∥}_{\epsilon }^2+{∥\nabla f_{0,h}∥}^2+{∥\nabla ·f_{0,h}∥}_{\left|\nabla \epsilon \right|+\epsilon -1}^2+{∥f_{0,h}∥}_{\left|\nabla \epsilon \right|}^2\right).$$

(32)

where

$$\begin{array}{cc}\hfill \left|\right||E_h{\left|\right||}_{\epsilon }^2\left(t\right):=& {∥{\partial }_t{E}_h∥}_{\epsilon }^2\left(t\right)+{∥\nabla E_h∥}^2\left(t\right)+{∥\nabla ·E_h∥}_{\left|\nabla \epsilon \right|+\epsilon -1}^2\left(t\right)+{∥E_h∥}_{\left|\nabla \epsilon \right|}^2\left(t\right)\hfill \\ \hfill & +{\int }_0^t{∥{\partial }_t{E}_h∥}_{\left|\nabla \epsilon \right|}^{2}ds+{\int }_0^t{∥\nabla ·E_h∥}_{\left|\nabla \epsilon \right|}^{2}ds.\hfill \end{array}$$

(33)

5.3. Stability of the Semi-Discrete Problem in ${\Omega }_{FEM}$

The stability of the semi-discrete problem in ${\Omega }_{FEM}$, relying on the variational formulation (18), and due to the appearance of the data function g, is slightly different from (32). We reformulate (18), in view of (12), and with $\mathbf{v}={\mathbf{v}}_h={\partial }_t{E}_h$ as: given $E_h(·,0)={f_0}_h$, ${\partial }_t{E}_h(·,0)={f_1}_h$, and ${\mathcal{P}}_{h}g$ (the $L_2$-projection of the boundary data g in (17)), find $E_h\in {\mathbf{W}}_h^E\left({\Omega }_{FEM}\right)$ such that

$$\begin{array}{cc}\hfill \left(\epsilon {\partial }_{tt}{E}_h,{\partial }_t{E}_h\right)& +(\nabla E_h,\nabla {\partial }_t{E}_h)+((\nabla \epsilon )·E_h,\nabla ·{\partial }_t{E}_h)+((\epsilon -1)\nabla ·E_h,\nabla ·{\partial }_t{E}_h)\hfill \\ \hfill & ={\langle {\partial }_n\left({\mathcal{P}}_{h}g\right),{\partial }_t{E}_h\rangle }_{\partial {\Omega }_{FEM}},\hfill \end{array}$$

(34)

Note that $\epsilon =1$ at the boundary of ${\Omega }_{FEM}$, hence the boundary term with the coefficient $\epsilon -1$ does not appear in hear.

To deal with the $(\nabla \epsilon )·E_h$-term, we use the discrete version of the $I_3$-estimate above and get

$$\begin{array}{cc}\hfill {\int }_0^t{I_h}_{3}ds& =-{\int }_{{\Omega }_{FEM}}\left({\int }_0^t((\nabla \epsilon )·{\partial }_t{E}_h)(\nabla ·E_h)ds-{\left[(\nabla \epsilon ·E_h)(\nabla ·E_h)\right]}_0^t\right)\hfill \\ \hfill & :={I_h}_{3,1}+{I_h}_{3,2}.\hfill \end{array}$$

(35)

Then

$${I_h}_{3,1}\le \frac{1}{2}{\int }_0^t{∥{\partial }_t{E}_h∥}_{\left|\nabla \epsilon \right|,{\Omega }_{FEM}}^{2}ds+\frac{1}{2}{\int }_0^t{∥\nabla ·E_h∥}_{\left|\nabla \epsilon \right|,{\Omega }_{FEM}}^{2}ds,$$

(36)

and

$$\begin{array}{cc}\hfill {I_h}_{3,2}\le & \frac{1}{2}\left({∥E_h∥}_{\left|\nabla \epsilon \right|,{\Omega }_{FEM}}^2+{∥\nabla ·E_h∥}_{\left|\nabla \epsilon \right|,{\Omega }_{FEM}}^2\right)\left(t\right)\hfill \\ \hfill & +\frac{1}{2}\left({∥E_h∥}_{\left|\nabla \epsilon \right|,{\Omega }_{FEM}}^2+{∥\nabla ·E_h∥}_{\left|\nabla \epsilon \right|,{\Omega }_{FEM}}^2\right)\left(0\right).\hfill \end{array}$$

(37)

As for the boundary term on the right-hand side of (34), with ${\mathcal{P}}_{h}g=E_{h,FDM}\sim E_h^{\Omega }{|}_{\partial {\Omega }_{FEM}}$, we have that

$$\langle {{\partial }_n(\mathcal{P}}_{h}g),{\partial }_t{\mathcal{P}}_h{g\rangle }_{\partial {\Omega }_{FEM}}\le \frac{1}{2}{∥{\partial }_n\left({\mathcal{P}}_{h}g\right)∥}_{\partial {\Omega }_{FEM}}^2+\frac{1}{2}{∥{\partial }_t\left({\mathcal{P}}_{h}g\right)∥}_{\partial {\Omega }_{FEM}}^2.$$

(38)

Now, in view of (34)–(38), and taking the time integral, as in the case of stability in $\Omega$, we end up with the following $L_2\left({\Omega }_{FEM}\right)$-type stability estimate:   

$$\begin{array}{cc}\hfill \left|\right||E_h{\left|\right||}_{\epsilon }^{2,FEM}\left(t\right):=& \left({∥{\partial }_t{E}_h∥}_{\epsilon }^2+{∥\nabla E_h∥}^2+{∥E_h∥}_{\left|\nabla \epsilon \right|}^2+{∥\nabla ·E_h∥}_{\left|\nabla \epsilon \right|+\epsilon -1}^2\right)\left(t\right)\hfill \\ \hfill & +{\int }_0^t{∥{\partial }_t{E}_h∥}_{\left|\nabla \epsilon \right|,{\Omega }_{FEM}}^{2}ds+{\int }_0^t{∥\nabla ·E_h∥}_{\left|\nabla \epsilon \right|,{\Omega }_{FEM}}^{2}ds\hfill \\ \hfill \le & C_1\left({∥{\partial }_t{E}_h∥}_{\epsilon }^2+{∥\nabla E_h∥}^2+{∥E_h∥}_{\left|\nabla \epsilon \right|}^2+{∥\nabla ·E_h∥}_{\left|\nabla \epsilon \right|+\epsilon -1}^2\right)\left(0\right)\hfill \\ \hfill & +C_2{\int }_0^t{∥{\partial }_n\left({\mathcal{P}}_{h}g\right)∥}_{\partial {\Omega }_{FEM}}^{2}ds+C_3{\int }_0^t{∥{\partial }_t\left({\mathcal{P}}_{h}g\right)∥}_{\partial {\Omega }_{FEM}}^2.\hfill \end{array}$$

(39)

Note the appearance of the time integrals of spatial norms of the two boundary terms on the right-hand side of the estimate for $\left|\right||E_h{\left|\right||}_{\epsilon }^{2,FEM}\left(t\right)$ in (39).

Summing up, we have the following stability estimate for the semi-discrete problem in ${\Omega }_{FEM}$:

Lemma 2.

Under the following regularity assumptions on the interpolants for initial conditions: $f_{1,h}\in L_{\epsilon }^2\left({\Omega }_{FEM}\right)$, $f_{0,h}\in H^1\left({\Omega }_{FEM}\right)\cap L_{\left|\nabla \epsilon \right|}^2\left({\Omega }_{FEM}\right)$, $\nabla ·f_{0,h}\in L_{\left|\nabla \epsilon \right|+\epsilon -1}^2\left({\Omega }_{FEM}\right)$, and with both boundary data: $g_h$, and ${\partial }_t{E}_h\in L_1\left((0,T);L_2(\partial {\Omega }_{FEM})\right)$, we have that there are suitably chosen positive constants $C_i,i=1,2,3$, such that for all $t\in (0,T]$, the following stability estimate for the semi-discrete ${\Omega }_{FEM}$ problem holds true

$$\begin{array}{cc}\hfill \left|\right||E_h{\left|\right||}_{\epsilon }^{2,FEM}\left(t\right)\le & C_1\left({∥f_{1,h}∥}_{\epsilon }^2+{∥\nabla f_{0,h}∥}^2+{∥f_{0,h}∥}_{\left|\nabla \epsilon \right|}^2+{∥\nabla ·f_{0,h}∥}_{\left|\nabla \epsilon \right|+\epsilon -1}^2\right)\hfill \\ \hfill & +C_2{\int }_0^t{∥{\partial }_n\left({\mathcal{P}}_{h}g\right)∥}_{\partial {\Omega }_{FEM}}^{2}ds+C_3{\int }_0^t{∥{\partial }_t\left({\mathcal{P}}_{h}g\right)∥}_{\partial {\Omega }_{FEM}}^2.\hfill \end{array}$$

(40)

6. Error Estimates: Semi-Discrete (SD) Problems

In what follows, assuming certain regularity of the data set, $E\in {\mathbf{W}}^E\left(\tilde{\Omega }\right)\cap H^s\left(\tilde{\Omega }\right)$, and with the spectral polynomial order p chosen as the assumed regularity degree $s\ge 1$ we can prove semi-discrete spatial error estimates of the form

$$\left|\right||E-E_h{\left|\right||}_{\tilde{\Omega }}\le C{h}^p\sim C{h}^s,\mathrm{for}\tilde{\Omega }=\Omega \mathrm{or}\tilde{\Omega }={\Omega }_{FEM}.$$

(41)

Below we use a very similar argument to derive (41) for the spatial domains $\Omega$ and ${\Omega }_{FEM}$. To do so, we rely on an approach (see, e.g., [46]) based on a split of the bilinear form $B(·,·)$ to a temporal term, plus a time-harmonic sub-bilinear form below denoted $a(·,·)$.

Splitting of B. Here we extract from B the bilinear form $a(·,·)$, by suppressing the time-derivative term, and establish well-posedness for a. Then we deal with the convergence of the original, time-dependent scheme of B. Thus

$$a(E,\mathbf{v}):=(\nabla E,\nabla \mathbf{v})+\left(\right(\nabla \epsilon )·E,\nabla ·\mathbf{v})+\left(\right(\epsilon -1)\nabla ·E,\nabla ·\mathbf{v}),$$

(42)

Due to the homogeneous Dirichlet boundary condition, the inherited Galerkin orthogonality yields:

$$a(E-E_h,\chi )=0,\forall \chi \in {\mathbf{W}}_{0,h}^E.$$

(43)

Now we can write the equation (25) as Find $E\in {\mathbf{W}}^E(\Omega )$ such that for all $\mathbf{v}\in {\mathbf{W}}_{0,h}^E$,

$$(\epsilon {\partial }_{tt}E,\mathbf{v})+a(E,\mathbf{v})=0.$$

(44)

The consistency of $a(·,·)$ is a consequence of (44).

We also define the norm associated to $a(·,·)$:

$${\left[∥E∥\right]}^2:={∥E∥}^2+{∥\nabla E∥}^2+{∥\nabla ·E∥}_{\epsilon -1+\left|\nabla \epsilon \right|}^2,$$

(45)

which, as an immediate concern, requires the condition below   

$$\left|\nabla \epsilon \right|<min(1/2,\epsilon -1).$$

(46)

The coercivity of $a(·,·)$ is due to vanishing boundary terms and assuming (for simplicity) that $|\Omega |=1$, we use the Poincare inequality $∥E∥\le ∥\nabla E∥$, and the second condition in (46) to write   

$$\begin{array}{cc}\hfill a(E,E)=& (\nabla E,\nabla E)+\left(\right(\nabla \epsilon )·E,\nabla ·E)+\left(\right(\epsilon -1)\nabla ·E,\nabla ·E)\hfill \\ \hfill \ge & \frac{1}{2}{∥E∥}^2+\frac{1}{2}{∥\nabla E∥}^2-\frac{1}{2}{∥E∥}_{\left|\nabla \epsilon \right|}^2-\frac{1}{2}{∥\nabla ·E∥}_{\left|\nabla \epsilon \right|}^2\hfill \\ \hfill & +\frac{1}{4}{∥\nabla ·E∥}_{\epsilon -1}^2+{∥\nabla ·E∥}_{\left|\nabla \epsilon \right|}^2\hfill \\ \hfill \ge & \frac{1}{2}{∥E∥}^2+\frac{1}{4}{∥\nabla E∥}^2+\frac{1}{4}{∥\nabla ·E∥}_{(\epsilon -1+2|\nabla \epsilon \left|\right)}\hfill \\ \hfill \ge & \frac{1}{4}\left({∥E∥}^2+{∥\nabla E∥}^2+{∥\nabla ·E∥}_{\epsilon -1+\left|\nabla \epsilon \right|}^2\right)=\frac{1}{4}{\left[∥E∥\right]}^2,\hfill \end{array}$$

(47)

where the second inequality is a consequence of reusing Poincare inequality for the term $-\frac{1}{2}{∥E∥}_{\left|\nabla \epsilon \right|}^2$, and we have a desired coercivity for $a(·,·)$.

The continuity of $a(·,·)$ is evident from the coercivity (47) and (46), which yields

$$max\left(∥E∥,∥\nabla E∥,{∥\nabla ·E∥}_{\epsilon -1+\left|\nabla \epsilon \right|}\right)\le \frac{1}{2}\left[∥E∥\right]$$

(48)

$$\begin{array}{cc}\hfill a(\eta ,\xi )& =(\nabla \eta ,\nabla \xi )+\left(\right(\nabla \epsilon )·\eta ,\nabla ·\xi )+\left(\right(\epsilon -1)\nabla ·\eta ,\nabla ·\xi )\hfill \\ \hfill & \le ∥\nabla \eta ∥∥\nabla \xi ∥+{∥\eta ∥}_{\left|\nabla \epsilon \right|}{∥\nabla ·\xi ∥}_{\left|\nabla \epsilon \right|}+{∥\nabla ·\eta ∥}_{\epsilon -1}{∥\nabla ·\xi ∥}_{\epsilon -1}\hfill \\ \hfill & \le \frac{1}{4}\left[∥\eta ∥\right]\left[∥\xi ∥\right]+\frac{1}{4}\left[∥\eta ∥\right]\left[∥\xi ∥\right]+\frac{1}{4}\left[∥\eta ∥\right]\left[∥\xi ∥\right]\hfill \\ \hfill & =\frac{3}{4}\left[∥\eta ∥\right]\left[∥\xi ∥\right]\le \left[∥\eta ∥\right]\left[∥\xi ∥\right].\hfill \end{array}$$

(49)

6.1. Error Estimates: SD Problem in $\Omega$

We define a Ritz–Galerkin type projection ${\mathcal{Q}}_h:H(div,\Omega )\to {\mathbf{W}}_{0,h}^E(\Omega )$ by

$$a({\mathcal{Q}}_{h}E,\mathbf{v})=a(E,\mathbf{v})\forall \mathbf{v}\in {\mathbf{W}}_{0,h}^E(\Omega ).$$

Then using the estimates (see [46] for the details):

$${\left[∥\xi ∥\right]}^2\le -Ca(\xi ,\xi )\le C\left[∥\xi ∥\right]\left[∥\eta ∥\right],$$

(50)

and

$$\left[∥\eta ∥\right]\le C{h}^{\zeta -1}{∥E∥}_{H^{\zeta }(\Omega )},1\le \zeta \le 2.$$

(51)

we can deduce that (see Lemma 5.4 in [46], $H^1$-estimate and interpolation),

$${∥E-{\mathcal{Q}}_{h}E∥}_{\Omega }+h^{\zeta -1}\left[∥E-{\mathcal{Q}}_{h}E∥\right]\le C{h}^{\zeta }{∥E∥}_{H^{\zeta }(\Omega )},1\le \zeta \le 2.$$

(52)

Finally, we recall the Galerkin orthogonality:

$$a(E-E_h,\chi )=0,\forall \chi \in {\mathbf{W}}_{0,h}^E(\Omega ).$$

(53)

Our main result is the following theorem.

Theorem 1.

Let E and $E_h$ be the solutions for the continuous and semi-discrete, weakly formulated, problems (25) and (13), respectively. Assume that $E\left(t\right),{\partial }_{t}E\left(t\right),{\partial }_{tt}^{2}E\in H^{\zeta }(\Omega ),1\le \zeta \le 2$. Then, there exists a constant C, independent of $x,t,h$, but dependant on ε, such that for $t\in [0,T]$, the error $e\left(t\right)=E(·,t)-E_h(·,t)$ is estimated as

$$\begin{array}{cc}\hfill {∥{\partial }_t^{r}e\left(t\right)∥}_{\Omega }& \le C{h}^{\zeta }{∥{\partial }_t^{r}E\left(t\right)∥}_{H^{\zeta }(\Omega )}+C{h}^{\zeta }{t}^{1-r}{\int }_0^t{∥{\partial }_{ss}^{2}E∥}_{H^{\zeta }(\Omega )}ds,r=0,1,\hfill \\ \hfill ∥e\left(t\right)∥{]}_{\Omega }& \le C{h}^{\zeta -1}{∥E\left(t\right)∥}_{H^{\zeta }(\Omega )}+C{h}^{\zeta -1}{\int }_0^t{∥{\partial }_{ss}^{2}E∥}_{H^{\zeta }(\Omega )}ds.\hfill \end{array}$$

(54)

Proof. 

We shall prove the theorem for $\zeta =2$. The case $\zeta =1$ is a simpler version, and $\zeta \in [1,2]$ follows using classical interpolation of the Sobolev spaces, see [47]. For $\zeta =2$, we use a split of error via the Ritz projection:

$$e:=E-E_h=(E-{\mathcal{Q}}_{h}E)+({\mathcal{Q}}_{h}E-E_h):=\rho +\theta .$$

To bound $\theta$ we note that using the Galerkin orthogonality (53),

$$\begin{array}{cc}\hfill (\epsilon {\theta }_{tt},\chi )+a(\theta ,\chi )& ={(\epsilon {\mathcal{Q}}_h{E}_{tt},\chi )}_{\Omega }+a{({\mathcal{Q}}_{h}E,\chi )}_{\Omega }-{(\epsilon {\partial }_{tt}^2{E}_h,\chi )}_{\Omega }-a{(E_h,\chi )}_{\Omega }\hfill \\ \hfill & ={(\epsilon {\mathcal{Q}}_h{E}_{tt},\chi )}_{\Omega }+a{({\mathcal{Q}}_{h}E,\chi )}_{\Omega }-{(\epsilon {\partial }_{tt}^{2}E,\chi )}_{\Omega }-a{(E,\chi )}_{\Omega }\hfill \\ \hfill & =-{(\epsilon {\rho }_{tt},\chi )}_{\Omega },\forall \chi \in {\mathbf{W}}_{0,h}^E.\hfill \end{array}$$

Then, for $\chi ={\theta }_t$,

$${(\epsilon {\theta }_{tt},{\theta }_t)}_{\Omega }+a(\theta ,{\theta }_t)=-{(\epsilon {\rho }_{tt},{\theta }_t)}_{\Omega }.$$

This yields (noting that $\epsilon$ is a function only in space)

$$\frac{d}{dt}\left(\epsilon {∥{\theta }_t∥}_{\Omega }^2+a(\theta ,\theta )\right)\le 2{∥\epsilon {\rho }_{tt}∥}_{\Omega }{∥{\theta }_t∥}_{\Omega }.$$

Note that, by the definition $\theta \left(0\right)={\theta }_t\left(0\right)=0$. Hence, integrating over $[0,t]$ we get

$$\begin{array}{cc}\hfill \epsilon {∥{\theta }_t\left(t\right)∥}_{\Omega }^2+a(\theta \left(t\right),\theta \left(t\right))& \le 2{\int }_0^t{∥\epsilon {\rho }_{ss}∥}_{\Omega }{∥{\theta }_s∥}_{\Omega }ds\hfill \\ \hfill & \le 2\epsilon \underset{s\in [0,t]}{max}{∥{\theta }_s∥}_{\Omega }{\int }_0^t{∥{\rho }_{ss}∥}_{\Omega }ds\hfill \\ \hfill & \le 2{\left({\int }_0^t{∥{\rho }_{ss}∥}_{\Omega }\right)}^2+\frac{1}{2}{\left(\epsilon \underset{s\in [0,t]}{max}{∥{\theta }_s∥}_{\Omega }\right)}^2.\hfill \end{array}$$

(55)

Now, using the fact that $a\left(\theta \right(t),\theta (t\left)\right)\ge 0$, and since (55) holds for all $t\in [0,T]$,

$$\frac{1}{2}\epsilon {\left(\underset{s\in [0,t]}{max}{∥{\theta }_s∥}_{\Omega }\right)}^2\le 2{\left({\int }_0^T{∥{\rho }_{ss}∥}_{\Omega }ds\right)}^2,$$

which inserting into (55) and using (52) noting that $\epsilon \ge 1$ and $h\sim \tau$ we get

$$\begin{array}{cc}\hfill \epsilon {∥{\theta }_t\left(t\right)∥}_{\Omega }^2+a(\theta \left(t\right),\theta \left(t\right))& \le 4{\left({\int }_0^T{∥{\rho }_{ss}∥}_{\Omega }\right)}^2\hfill \\ \hfill & \le 4{\left(C{h}^2{\int }_0^T{∥E_{ss}∥}_{H^2(\Omega )}ds\right)}^2\le 4\epsilon {\left(C{h}^2{\int }_0^T{∥E_{ss}∥}_{H^2(\Omega )}ds\right)}^2.\hfill \end{array}$$

(56)

Consequently,

$${∥{\theta }_t\left(t\right)∥}_{\Omega }\le C{h}^2{\int }_0^T{∥E_{ss}∥}_{H^2(\Omega )}ds,$$

and

$${\left[∥\theta \left(t\right)∥\right]}_{\Omega }\le C{h}^2{\int }_0^T{∥E_{ss}∥}_{H^2(\Omega )}ds\le C{h}^2{\int }_0^T{∥E_{ss}∥}_{H^2(\Omega )}ds.$$

Note further that

$$\begin{array}{cc}\hfill 2{∥\theta ∥}_{\Omega }\frac{d}{dt}{∥\theta ∥}_{\Omega }& =\frac{d}{dt}{∥\theta ∥}_{\Omega }^2=\frac{d}{dt}{\int }_{\Omega }{\left|\theta \right|}^{2}dx\hfill \\ \hfill & =2{\int }_{\Omega }\theta ·{\theta }_{t}dx\le 2{∥\theta ∥}_{\Omega }{∥{\theta }_t∥}_{\Omega }.\hfill \end{array}$$

(57)

Which integrating yields   

$${∥\theta \left(t\right)∥}_{\Omega }\le {\int }_0^t{∥{\theta }_s\left(s\right)∥}_{\Omega }ds\le C{h}^{2}t{\int }_0^T{∥E_{ss}∥}_{H^2(\Omega )}ds.$$

(58)

Summing up, we get the desired result.    ☐

Remark 4.

Note that the errors in initial data, $f_0-f_{0,h}$ and $f_1-f_{1,h}$, do not appear in the above estimates. In general, and for sufficiently regular $f_0$ and $f_1$, the initial data errors can always be dominated by other error terms in appropriate norms.

6.2. Error Estimates: SD Problem in ${\Omega }_{FEM}$

Here the non-vanishing boundary data causes additional asymmetry (the first one: $\left(\right(\nabla \epsilon )·E,\nabla ·\mathbf{v})$ was tackled above), which may be symmetrized using a Nitsche-type approach of [46]. We circumvent Nitsche invoking additional conditions on the Ritz projection at the boundary $\partial {\Omega }_{FEM}$. The variational form and the finite element formulation in ${\Omega }_{FEM}$ are

$$\begin{array}{cc}\hfill {(\epsilon {\partial }_{tt}E,\mathbf{v})}_{{\Omega }_{FEM}}& +{(\nabla E,\nabla \mathbf{v})}_{{\Omega }_{FEM}}+{((\nabla \epsilon )·E,\nabla ·\mathbf{v})}_{{\Omega }_{FEM}}+{((\epsilon -1)\nabla ·E,\nabla ·\mathbf{v})}_{{\Omega }_{FEM}}\hfill \\ \hfill & -{\langle {\partial }_{n}E,\mathbf{v}\rangle }_{\partial {\Omega }_{FEM}}=0,\forall \mathbf{v}\in {\mathbf{W}}^E\left({\Omega }_{FEM}\right),\hfill \end{array}$$

(59)

and

$$\begin{array}{cc}\hfill {(\epsilon {\partial }_{tt}{E}_h,\mathbf{v})}_{{\Omega }_{FEM}}& +{(\nabla E_h,\nabla \mathbf{v})}_{{\Omega }_{FEM}}+{((\nabla \epsilon )·E_h,\nabla ·\mathbf{v})}_{{\Omega }_{FEM}}+{((\epsilon -1)\nabla ·E_h,\nabla ·\mathbf{v})}_{{\Omega }_{FEM}}\hfill \\ \hfill & -{\langle {\partial }_n{E}_h,\mathbf{v}\rangle }_{\partial {\Omega }_{FEM}}=0,\forall \mathbf{v}\in {\mathbf{W}}_h^E\left({\Omega }_{FEM}\right),\hfill \end{array}$$

(60)

respectively. Hence, subtracting (60) from (59), where in (59) we restrict $\mathbf{v}\in {\mathbf{W}}_h^E\left({\Omega }_{FEM}\right)$, we end up with the error equation for $\tilde{e}:=E-E_h$, in ${\Omega }_{FEM}$:

$$\begin{array}{cc}\hfill {(\epsilon {\partial }_{tt}\tilde{e},\mathbf{v})}_{{\Omega }_{FEM}}& +{(\nabla \tilde{e},\nabla \mathbf{v})}_{{\Omega }_{FEM}}+{((\nabla \epsilon )·\tilde{e},\nabla ·\mathbf{v})}_{{\Omega }_{FEM}}+{((\epsilon -1)\nabla ·\tilde{e},\nabla ·\mathbf{v})}_{{\Omega }_{FEM}}\hfill \\ \hfill & -{\langle {\partial }_n\tilde{e},\mathbf{v}\rangle }_{\partial {\Omega }_{FEM}}=0,\forall \mathbf{v}\in {\mathbf{W}}_h^E\left({\Omega }_{FEM}\right).\hfill \end{array}$$

(61)

To be concise, for ${\Omega }_{FEM}$, we suppress the subscript ${\Omega }_{FEM}$ from the inner products, and to distinguish from $\Omega$, use $\tilde{e}$, $\tilde{B}$, and $\tilde{a}$ for the error and the bilinear forms in ${\Omega }_{FEM}$. Now we define a somewhat modified Ritz–Galerkin projection: ${\tilde{\mathcal{Q}}}_h:H(div,{\Omega }_{FEM})\to {\mathbf{W}}_h^E\left({\Omega }_{FEM}\right)$ by

$$\tilde{a}({\tilde{\mathcal{Q}}}_h\mathbf{u},\mathbf{v})-{\langle {\partial }_n\left({\tilde{\mathcal{Q}}}_h\mathbf{u}\right),\mathbf{v}\rangle }_{\partial {\Omega }_{FEM}}=\tilde{a}(\mathbf{u},\mathbf{v})-{\langle {\partial }_n\mathbf{u},\mathbf{v}\rangle }_{\partial {\Omega }_{FEM}}\forall \mathbf{v}\in {\mathbf{W}}_h^E\left({\Omega }_{FEM}\right),$$

(62)

with $\tilde{a}$ defined as a, but now over ${\Omega }_{FEM}$, and with the additional property of satisfying weak Neumann boundary condition that behaves like Galerkin approximation:

$${\langle {\partial }_n\left({\tilde{\mathcal{Q}}}_h\mathbf{u}\right),\mathbf{v}\rangle }_{\partial {\Omega }_{FEM}}={\langle {\partial }_n{\mathbf{u}}_h,\mathbf{v}\rangle }_{\partial {\Omega }_{FEM}}\forall \mathbf{v}\in {\mathbf{W}}_h^E\left({\Omega }_{FEM}\right),$$

(63)

Note that denoting $\tilde{B}:=B_{{\Omega }_{FEM}}$, (61) yields a modified Galerkin orthogonality, viz.

$$\tilde{B}(E,\mathbf{v})=\tilde{B}(E_h,\mathbf{v}),\forall \mathbf{v}\in {\mathbf{W}}_h^E\left({\Omega }_{FEM}\right),$$

(64)

which can be rewritten as:

$$(\epsilon E_{h,tt},\mathbf{v})+\tilde{a}(E_h,\mathbf{v})-{\langle {\partial }_n{E}_h,\mathbf{v}\rangle }_{\partial {\Omega }_{FEM}}=(\epsilon E_{tt},\mathbf{v})+\tilde{a}(E,\mathbf{v})-{\langle {\partial }_{n}E,\mathbf{v}\rangle }_{\partial {\Omega }_{FEM}}$$

(65)

Now, we follow the procedure above (using the very similar estimates of $\rho$ and $\theta$ for $\tilde{\rho }$ and $\tilde{\theta }$), and let

$$\tilde{e}:=E-E_h:=(E-{\tilde{\mathcal{Q}}}_{h}E)+({\tilde{\mathcal{Q}}}_{h}E-E_h):=\tilde{\rho }+\tilde{\theta }.$$

(66)

Then by this definition   

$$\begin{array}{cc}\hfill {(\epsilon {\partial }_{tt}\tilde{\theta },\mathbf{v})}_{{\Omega }_{FEM}}& -{\langle {\partial }_n\tilde{\theta },\mathbf{v}\rangle }_{\partial {\Omega }_{FEM}}+\tilde{a}{(\tilde{\theta },\mathbf{v})}_{{\Omega }_{FEM}}\hfill \\ \hfill & ={(\epsilon {\tilde{\mathcal{Q}}}_h{E}_{tt},\mathbf{v})}_{{\Omega }_{FEM}}-{(\epsilon E_{h,tt},\mathbf{v})}_{{\Omega }_{FEM}}-{\langle {\partial }_n{\tilde{\mathcal{Q}}}_{h}E,\mathbf{v}\rangle }_{\partial {\Omega }_{FEM}}\hfill \\ \hfill & +{\langle {\partial }_n{E}_h,\mathbf{v}\rangle }_{\partial {\Omega }_{FEM}}+\tilde{a}{({\tilde{\mathcal{Q}}}_{h}E,\mathbf{v})}_{{\Omega }_{FEM}}-\tilde{a}{(E_h,\mathbf{v})}_{{\Omega }_{FEM}}\hfill \end{array}$$

(67)

Thus recalling (65) and (62), we end up with

$$\begin{array}{cc}\hfill {(\epsilon {\partial }_{tt}\tilde{\theta },\mathbf{v})}_{{\Omega }_{FEM}}& -{\langle {\partial }_n\tilde{\theta },\mathbf{v}\rangle }_{\partial {\Omega }_{FEM}}+\tilde{a}{(\tilde{\theta },\mathbf{v})}_{{\Omega }_{FEM}}\hfill \\ \hfill & ={(\epsilon {\tilde{\mathcal{Q}}}_h{E}_{tt},\mathbf{v})}_{{\Omega }_{FEM}}-{\langle {\partial }_n{\tilde{\mathcal{Q}}}_{h}E,\mathbf{v}\rangle }_{\partial {\Omega }_{FEM}}+\tilde{a}{({\tilde{\mathcal{Q}}}_{h}E,\mathbf{v})}_{{\Omega }_{FEM}}\hfill \\ \hfill & -{(\epsilon E_{tt},\mathbf{v})}_{{\Omega }_{FEM}}+{\langle {\partial }_{n}E,\mathbf{v}\rangle }_{\partial {\Omega }_{FEM}}-\tilde{a}{(E,\mathbf{v})}_{{\Omega }_{FEM}}\hfill \\ \hfill & =-{(\epsilon {\tilde{\rho }}_{tt},\mathbf{v})}_{{\Omega }_{FEM}}.\hfill \end{array}$$

(68)

The remaining part is as in the proof of the previous theorem for the error estimates for $\Omega$. Thus summing up we have proved the following error estimate in ${\Omega }_{FEM}$:

Theorem 2.

Let E and $E_h$ be the solutions for the continuous and semi-discrete, weakly formulated, problems (25) and (13), respectively. Further, assume that $E\left(t\right),{\partial }_{t}E\left(t\right),{\partial }_{tt}^{2}E\in H^{\zeta }(\Omega )$, $1\le \zeta \le 2$. Then, there exists a constant C depending on ε but independent of $x,t,h$ such that for $t\in [0,T]$ and $\tilde{e}$ defined as in (66);

$$\begin{array}{cc}\hfill {∥{\partial }_t^r\tilde{e}\left(t\right)∥}_{{\Omega }_{FEM}}& \le C{h}^{\zeta }{∥{\partial }_t^{r}E\left(t\right)∥}_{H^{\zeta }\left({\Omega }_{FEM}\right)}+C{h}^{\zeta }{t}^{1-r}{\int }_0^t{∥{\partial }_{ss}^{2}E∥}_{H^{\zeta }\left({\Omega }_{FEM}\right)}ds,r=0,1\hfill \\ \hfill ∥\tilde{e}\left(t\right)∥{]}_{{\Omega }_{FEM}}& \le C{h}^{\zeta -1}{∥E\left(t\right)∥}_{H^{\zeta }\left({\Omega }_{FEM}\right)}+C{h}^{\zeta -1}{\int }_0^t{∥{\partial }_{ss}^{2}E∥}_{H^{\zeta }\left({\Omega }_{FEM}\right)}ds.\hfill \end{array}$$

(69)

Corollary 1.

Under the conditions in the theorem (2) and recalling (63) we have

$${\left\{∥\tilde{e}∥\right\}}_{{\Omega }_{FEM}}:={∥\tilde{e}∥}_{{\Omega }_{FEM}}+{∥\tilde{e}∥}_{\partial {\Omega }_{FEM}}=\mathcal{O}{h}^{\zeta -1/2},1\le \zeta \le 2.$$

(70)

Proof. 

Using the split (66) and (63) we have

$${(\tilde{e},\tilde{\theta })}_{\partial {\Omega }_{FEM}}={(\rho ,\theta )}_{\partial {\Omega }_{FEM}}.$$

(71)

As an immediate consequence of (71) we get

$${∥\tilde{e}∥}_{\partial {\Omega }_{FEM}}\le {∥\tilde{\rho }∥}_{\partial {\Omega }_{FEM}}.$$

(72)

Further, by the trace theorem

$${∥\tilde{\rho }∥}_{\partial {\Omega }_{FEM}}\le {∥h^{-1/2}\tilde{\rho }∥}_{{\Omega }_{FEM}}=\mathcal{O}\left(h^{\zeta -1/2}\right),1\le \zeta \le 2,$$

(73)

which yields the desired result.    ☐

7. Numerical Examples

In this section, we present numerical examples justifying the theoretical results of the previous two sections. For convergence tests the domain decomposition algorithm (see Algorithm 2), implemented in the software package WavES [48], was used. We note that because of using explicit FE and FD schemes in ${\Omega }_{FEM}$ and ${\Omega }_{FDM}$, correspondingly, we need to choose a time step $\tau$ according to the CFL stability condition (9) derived in [44] so that the whole hybrid scheme remains stable.

Algorithm 2 Domain decomposition algorithm for the solution of Maxwell’s equation (Equation (2)). At every time step k we performed the following operations:

1: Compute $E_h^{k+1}$ in ${\Omega }_{FDM}$ using the explicit finite difference scheme (24) with known $E_h^k$, and $E_h^{k-1}$-values. 2: Compute $E_h^{k+1}$ in ${\Omega }_{FEM}$ by using the finite element scheme (19) with known $E_h^k,E_h^{k-1}$. 3: For the finite difference method in ${\Omega }_{FDM}$, use the values of the function $E_h^{k+1}$ at nodes ${\omega }_{\diamond }$ (green boundary of Figure 1), which are computed using the finite element scheme (19), as a boundary condition at the inner boundary of ${\Omega }_{FDM}$. 4: Apply appropriate boundary condition at the outer boundary of ${\Omega }_{FDM}$. 5: For the finite element method in ${\Omega }_{FEM}$, use the values of the functions $E_h^{k+1}$ at nodes ${\omega }_{\mathrm{o}}$ (blue boundary of the Figure 1), which are computed using the finite difference scheme (24) as a boundary condition. 6: Apply swap of the solutions for the computed function $E_h^{k+1}$ to be able to perform the algorithm on a new time level k.

Numerical tests are performed in time interval $(0,T)=(0,0.25)$ and in the spatial dimensionless computational domain

$$\Omega =\left\{(x,y):x\in [0,1],y\in [0,1]\right\},$$

(74)

which is split into the finite element domain

$${\Omega }_{FEM}=\left\{(x,y):x\in [0.25,0.75],y\in [0.25,0.75]\right\}$$

(75)

and the finite difference domain ${\Omega }_{FDM}$, thus $\Omega ={\Omega }_{FEM}\cup {\Omega }_{FDM}$, see Figure 1.

The model problem in all our tests that is stated for the electric field $E=(E_1,E_2)$ is as follows:

$$\begin{array}{cc}\epsilon {\partial }_{tt}E+\nabla (\nabla ·E)-▵E-\nabla \nabla ·\left(\epsilon E\right)=F\hfill & \mathrm{in} \Omega \times (0,T),\hfill \\ E(·,0)=0 \mathrm{and} {\partial }_{t}E(·,0)=0\hfill & \mathrm{in} \Omega ,\hfill \\ E=0\hfill & \mathrm{on} \partial \Omega \times (0,T).\hfill \end{array}$$

(76)

We have the functions

$$\begin{array}{cc}\hfill E_1& =\frac{1}{\epsilon }2\pi {sin}^2\pi xcos\pi ysin\pi y\frac{t^2}{2},\hfill \\ \hfill E_2& =-\frac{1}{\epsilon }2\pi {sin}^2\pi ycos\pi xsin\pi x\frac{t^2}{2},\hfill \end{array}$$

(77)

as the exact solution $E=(E_1,E_2)$ of the model problem (76) with the source data $F=(F_1,F_2)$ which corresponds to this exact solution.

The function $\epsilon$ in (76) is defined as

$$\epsilon (x,y)=\left\{\begin{array}{cc}1+{sin}^m\pi (2x-0.5)·{sin}^m\pi (2y-0.5)\hfill & \mathrm{in} {\Omega }_{FEM},\hfill \\ 1\hfill & \mathrm{in} {\Omega }_{FDM}.\hfill \end{array}\right.$$

(78)

We choose $m=2,4,6,8$ in our numerical examples, see Figure 2, for these functions in the domain ${\Omega }_{FEM}$. We note that the exact solution (77) satisfies the divergence-free condition $\nabla ·\left(\epsilon E\right)=0$ for $\epsilon$ defined by (78), the homogeneous initial conditions, as well as the homogeneous Dirichlet conditions for all times.

The computational domain ${\Omega }_{FEM}\times (0,T)$ was discretized into triangular elements with mesh sizes $h_l=2^{-l},l=3,4,5,6$, and the mesh in ${\Omega }_{FDM}\times (0,T)$ was decomposed into squares of the same mesh sizes as described in Section 3, see Figure 1. The time step was chosen corresponding to the stability criterion (9) as ${\tau }_l=0.025·2^{-l}$ for $l=3,4,5,6$. Convergence results of the proposed finite element scheme computed in $L_2$ and $H^1$ norms are presented in

Table 1. Relative errors $e_l^1$ and $e_l^2$ in the $L_2$ and $H^1$ norms, respectively, for mesh sizes $h_l=2^{-l}$, $l=3,\dots ,6$, for $m=2$ in (78).

l	$\mathit{nel}$	$\mathit{nno}$	${\mathit{e}}_{\mathit{l}}^1$	$\frac{{\mathit{el}}_{\mathit{h}}^1}{{\mathit{el}}_{2\mathit{h}}^1}$	${\mathit{r}}_1$	${\mathit{e}}_{\mathit{l}}^2$	$\frac{{\mathit{el}}_{\mathit{h}}^2}{{\mathit{el}}_{2\mathit{h}}^2}$	${\mathit{r}}_2$
3	128	81	$4.878·{10}^{-2}$	−	−	$3.902·{10}^{-1}$	−	−
4	512	289	$1.222·{10}^{-2}$	$3.992$	$1.997$	$1.955·{10}^{-1}$	$1.996$	$0.997$
5	2048	1089	$2.654·{10}^{-3}$	$4.604$	$2.203$	$8.492·{10}^{-2}$	$2.302$	$1.203$
6	8192	4225	$5.15·{10}^{-4}$	$5.151$	$2.365$	$3.297·{10}^{-2}$	$2.575$	$1.365$

,

Table 2. Relative errors $e_l^1$ and $e_l^2$ in the $L_2$ norm and in the $H^1$ norm, respectively, for mesh sizes $h_l=2^{-l},l=3,\dots ,6$, for $m=4$ in (78).

l	$\mathit{nel}$	$\mathit{nno}$	${\mathit{e}}_{\mathit{l}}^1$	$\frac{{\mathit{el}}_{\mathit{h}}^1}{{\mathit{el}}_{2\mathit{h}}^1}$	${\mathit{r}}_1$	${\mathit{e}}_{\mathit{l}}^2$	$\frac{{\mathit{el}}_{\mathit{h}}^2}{{\mathit{el}}_{2\mathit{h}}^2}$	${\mathit{r}}_2$
3	128	81	$5.54·{10}^{-2}$	−	−	$4.432·{10}^{-1}$	−	−
4	512	289	$8.438·{10}^{-3}$	$6.565$	$2.7148$	$1.35·{10}^{-1}$	$3.283$	$1.715$
5	2048	1089	$1.581·{10}^{-3}$	$5.337$	$2.4160$	$5.06·{10}^{-2}$	$2.668$	$1.416$
6	8192	4225	$3.35·{10}^{-4}$	$4.722$	$2.2394$	$2.143·{10}^{-2}$	$2.361$	$1.239$

,

Table 3. Relative errors $e_l^1$ and $e_l^2$ in the $L_2$ norm and in the $H^1$ norm, respectively, for mesh sizes $h_l=2^{-l},l=3,\dots ,6$, for $m=6$ in (78).

l	$\mathit{nel}$	$\mathit{nno}$	${\mathit{e}}_{\mathit{l}}^1$	$\frac{{\mathit{el}}_{\mathit{h}}^1}{{\mathit{el}}_{2\mathit{h}}^1}$	${\mathit{r}}_1$	${\mathit{e}}_{\mathit{l}}^2$	$\frac{{\mathit{el}}_{\mathit{h}}^2}{{\mathit{el}}_{2\mathit{h}}^2}$	${\mathit{r}}_2$
3	128	81	$1.856·{10}^{-2}$	−	−	$1.485·{10}^{-1}$	−	−
4	512	289	$5.168·{10}^{-3}$	$3.592$	$1.845$	$8.268·{10}^{-2}$	$1.796$	$0.778$
5	2048	1089	$1.594·{10}^{-3}$	$3.243$	$1.697$	$5.099·{10}^{-2}$	$1.621$	$0.697$
6	8192	4225	$3.6·{10}^{-4}$	$4.432$	$2.148$	$2.301·{10}^{-2}$	$2.216$	$1.148$

and

Table 4. Relative errors $e_l^1$ and $e_l^2$ in the $L_2$ norm and in the $H^1$ norm, respectively, for mesh sizes $h_l=2^{-l},l=3,\dots ,6$, for $m=8$ in (78).

l	$\mathit{nel}$	$\mathit{nno}$	${\mathit{e}}_{\mathit{l}}^1$	$\frac{{\mathit{el}}_{\mathit{h}}^1}{{\mathit{el}}_{2\mathit{h}}^1}$	${\mathit{r}}_1$	${\mathit{e}}_{\mathit{l}}^2$	$\frac{{\mathit{el}}_{\mathit{h}}^2}{{\mathit{el}}_{2\mathit{h}}^2}$	${\mathit{r}}_2$
3	128	81	$1.131·{10}^{-2}$	−	−	$1.045·{10}^{-1}$	−	−
4	512	289	$5.669·{10}^{-3}$	$2.304$	$1.204$	$9.071·{10}^{-2}$	$1.152$	$0.2$
5	2048	1089	$1.711·{10}^{-3}$	$3.314$	$1.728$	$5.475·{10}^{-2}$	$1.657$	$0.728$
6	8192	4225	$3.83·{10}^{-4}$	$4.468$	$2.16$	$2.451·{10}^{-2}$	$2.234$	$1.16$

for $m=2,4,6,8$ in (78). Relative norms in these tables were computed as

$$\begin{array}{cc}\hfill e_l^1& {\displaystyle =\frac{{\displaystyle \underset{1\le k\le N}{max}\parallel E^k-E_h^k\parallel }}{{\displaystyle \underset{1\le k\le N}{max}\parallel E^k\parallel }},}\hfill \\ \hfill e_l^2& {\displaystyle =\frac{{\displaystyle \underset{1\le k\le N}{max}\parallel \nabla (E^k-E_h^k)\parallel }}{{\displaystyle \underset{1\le k\le N}{max}\parallel \nabla E^k\parallel }}.}\hfill \end{array}$$

(79)

E and $E_h$ are the exact and computed FE solutions in ${\Omega }_{FEM}\times (0,T)$, respectively, and $N=T/{\tau }_l$. Logarithmic convergence rates $r_1,r_2$ in these tables are computed, viz.

$$\begin{array}{cc}\hfill r_1& =\frac{\left|log\left(\frac{e{l}_h^1}{e{l}_{2h}^1}\right)\right|}{|log(0.5\left)\right|},\mathrm{and}{r}_2=\frac{\left|log\left(\frac{e{l}_h^2}{e{l}_{2h}^2}\right)\right|}{|log(0.5\left)\right|},\hfill \end{array}$$

(80)

where $e{l}_h^{1,2},e{l}_{2h}^{1,2}$ are relative norms computed via (79) on the mesh ${\mathcal{T}}_h$ with the mesh size h and $2h$, respectively. Figure 3 shows convergence of the relative $L_2$- and $H^1$-norms computed via (79) and compared with exact behavior of h and $h^2$.

Computed hybrid FE/FD versus exact solutions with $m=8$ in (78) at the time $t=0.25$ are presented in Figure 4, where $|E_h|$ is computed using the domain decomposition algorithm (Algorithm 2) for different meshes with sizes $h_l=2^{-l},l=3,4,5,6$. The top figures of Figure 4 present hybrid FE/FD meshes which were used for computations; common hybrid FE/FD solution in $\Omega$ is presented in the middle figures; the bottom figures show only FD solution as part of the common hybrid solution in ${\Omega }_{FDM}$. Interpreting these figures, we observe smooth behavior of the hybrid solution across finite element/finite difference boundary, as was predicted in theory.

Furthermore, through these figures, as well as tables and Figure 3, we observe that with increasing l in $h_l=2^{-l},l=3,4,5,6$, the computational errors approach the second order convergence in $L_2$- and first order in $H^1$-norm for $m=2,4,6,8$. Therefore, we can conclude that the finite element scheme in the hybrid FE/FD method, considered in ${\Omega }_{FEM}$, behaves like a first-order method in the $H^1$-norm and a second-order method in the $L_2$-norm. These results are all in good agreement with the analytic estimates derived in Section 5 and Section 6, as well as with results presented for finite element method in [44] for the whole $\Omega$.

8. Conclusions

In this paper, we present stability and convergence analysis for the domain decomposition FE/FD method for time-dependent Maxwell’s equations developed in [3,4]. The convergence is optimal due to the assumed maximal available regularity of the exact solution in a Sobolev space in correlation with the optimal interpolation rates in the power s of the mesh function h, see e.g., (41) and (69).

The analysis is performed for the semi-discrete (spatial discretization) problem for the constructed finite element schemes in two different settings: in $\Omega$ and ${\Omega }_{FEM}$. The temporal discretization algorithms are constructed using the CFL condition (9) derived in [44]. We have implemented several numerical examples that validate the robustness of the theoretical studies.