Stability and convergence analysis of a domain decomposition FE/FD method for the Maxwell's equations in time domain

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



There are no comments yet.


page 22

page 28


Finite element analysis for a diffusion equation on a harmonically evolving domain

We study convergence of the evolving finite element semi-discretization ...

Sensitivity Analysis of Domain Decomposition method for 4D Variational Data Assimilation (DD-4DVAR DA)

We prove consistence, convergence and stability of the Domain Decomposit...

Analysis of algebraic flux correction for semi-discrete advection problems

We present stability and error analysis for algebraic flux correction sc...

Modelling of a Permanent Magnet Synchronous Machine Using Isogeometric Analysis

Isogeometric analysis (IGA) is used to simulate a permanent magnet synch...

Parallel domain decomposition solvers for the time harmonic Maxwell equations

The time harmonic Maxwell equations are of current interest in computati...

A FETI approach to domain decomposition for meshfree discretizations of nonlocal problems

We propose a domain decomposition method for the efficient simulation of...

Optimally Convergent Mixed Finite Element Methods for the Stochastic Stokes Equations

We propose some new mixed finite element methods for the time dependent ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

New computational techniques meet the needs of industry in developing efficient computational methods to simulate partial differential equations (PDEs). Especially, for simulations in higher dimensions and large computational domains. In this regard, certain

domain decomposition (DD) method leading to efficient schemes, in numerical investigations, has gained a lot of interest in numerical analysis community. This variant of the DD method is the subject of our current, and some related and ongoing, research. This type of schemes was previously studied, e.g. in [12, 30] for others than studied here problems.

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 electric field in non-conductive media, studied in [3, 4]. A stable, time-domain (TD), DDM scheme for Maxwell’s equations was proposed in [26, 27], and further verified in [13]. This method uses FDTD scheme on the, structured, FD part of the mesh, and edge elements on unstructured part. In applications, because of edge elements implementation, this method remains computationally expensive.

The domain decomposition FE/FD method for time-dependent Maxwell’s equations for electric field, assuming constant dielectric permittivity function in a finite difference domain, was considered in [3]

. This assumption simplifies the numerical schemes in both FE and FD domains and significantly reduces computational efforts for implementation of the whole DD method. Modified numerical scheme, energy estimate and numerical verifications of this method was presented in

[4]. However, the fully stability and convergence analysis with numerical implementations, in - and -norms of the developed FE and FD schemes, are not presented in the above studies. We fill this gap in the present work.

More specificaly, we present stability analysis for explicit schemes for both FEM and FDM in the DD hybrid FE/FD method. The 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. Similar to the DD approach of [4, 5], we decompose the computational domain such that FEM and FDM are used in different subdomains: FDM in simple geometry and FE in the subdomain where more detailed information is needed about the structure of this subdomain. This also allows application of adaptive FEM in such subdomain, see, e.g. [2, 3, 8, 9, 10, 28, 29].

Reliability and convergence of the domain decomposition method, studied in this work, are evident for solution of coefficient inverse problems (CIPs) in , see, e.g. [8, 9, 10, 28, 29]. For the case of CIPs, the computational domain is splitted into subdomains such that a simple discretization scheme can be used in a large region and more refined discretization scheme is applied in smaller, however more critical, part 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 is necessary on it’s boundary or in it’s neighborhood. In such cases it is convenient to condsider the numerical solution of time-dependent Maxwell’s equations in different subdomains with constant dielectric permittivity function in some subdomain and non-constant in the other ones. For the time-dependent Maxwell’s equations, the DD scheme of [4], which is analyzed in the present work, is used for solution of different CIPs to determine the dielectric permittivity function in non-conductive media using simulated and experimentally generated data, see [2, 8, 9, 10, 28, 29].

An outline of this paper is as follows. In Section 2 we introduce the mathematical model. In Section 3 we briefly present the domain decomposition FE/FD method and communication scheme between two methods. In Section 4 we describe the domain decomposition FE/FD method for solution of Maxwell’s equations and set up the finite element and finite difference schemes. Section 5 is devoted to the stability analysis. In Section 6, we derive optimal a priori error estimates in finite element method for the semi-discrete (spatial discretizations) problems. Finally, in Section 7 we present numerical implementations that justify the theoretical investigations of the paper. In what follows, will be a generic constant independent of all parameters, unless otherwise specifically specified, and not necessarliy the same at each occurance.

2 The mathematical model

The Cauchy problem for the electric field , , , of the Maxwell’s equations, under the assumptions that the dimensionless relative magnetic permeability of the medium is and the electric volume charges are zero, is given by


where, and are the dimensionless relative dielectric permittivity and electric conductivity functions, respectively. , and are the permittivity and permeability of the free space, respectively, and is the speed of light in free space. In this paper we consider the problem (1) in non-conductive media, i.e. , and hence study the initial value problem


To solve the problem (2) numerically, we consider it in a bounded domain (instead of whole ), with boundary , and employ a split scheme on : a hybrid, finite element/finite difference scheme, kind of domain decomposition, developed in [3, 4] and summarized in Algorithm 1. More specifically, we divide the computational domain into two subregions, and such that and is a subset of the convex hul of . The function is assumed to be constant in , and bounded and smooth in . The communication ß between and is arranged using an overlapping mesh structure through a two-element thick layer around as shown by blue and green common boundaries in Figure 1. The blue boundary is outer boundary of and inner boundary of . Similarly, the green boundary is the inner boundary of from which the solution is copied to the green boundary of .

The key idea with such a decomposition is to be able apply different numerical methods in different computational domains. For the numerical solution of (2) in we use the finite difference method on a structured mesh. In , we use finite elements on a sequence of unstructured meshes , with elements consisting of triangles in and tetrahedra in , both satisfying minimal angle condition. This approach combines the flexibility of the finite elements and the efficiency of the finite differences in terms of speed and memory usage and fits well for reconstruction algorithms presented below.

We assume that for some known constant , the function satisfies


Conditions (3) on and the relation


together with divergence free field , make equations in (2) independent of each others in and so that, in , we just need to solve the system of wave equations:



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 [21, 24]. However, the edge elements are less attractive for solution of time-dependent problems, since a linear system of equations should be solved at each time iteration step. In contrary, P1 elements can be efficiently used in a fully explicit finite element scheme with lumped mass matrix [11, 14, 18]. It is also well known that numerical solution of Maxwell’s equations using nodal finite elements is often unstable and results spurious oscillatory solutions [22, 25]. There are a number of techniques to overcome such instabilities, see, e.g. [15, 16, 17, 23, 25].

In [6, 7], a finite element analysis shows stability and consistency of the stabilized finite element method for the solution of (1) with . In the current study we show stability and convergence for the combined FEM/FDM scheme, under the condition (3) on , where the stabilized FEM is used for the numerical solution of (2) in and usual FDM discretization of (5) is applied in .


Here, we consider the case when for . Further, we assume that , and . Recall that we assumed non-conductive media: . In the presence of electric conductivity, additional -terms appear in the equations. They lead to more involved estimates and heavier implementations which we plan to perform in a forthcoming study.

Hence, in this note we study the following initial boundary value problem:


3 The structure of domain decomposition

a) b) c)  
Figure 1:   Domain decomposition and mesh discretization in . The mesh of is a combination of the quadrilateral finite difference mesh presented on c), and the finite element mesh presented on b). Domains and overlap by two layers of structured nodes such that they have common boundaries shown by green and blue colors.

We now describe the DD method between two domains, and , where the FEM is used for computation of the solution in , and FDM is used in . Communication between and is achieved letting overlapping of both meshes across a two-element thick layer around - see Figure 1. The common nodes of both and domains belong to either of the following boundaries (see Figure 1):

  • Nodes on the blue boundary - lie on the boundary of and are interior to ,

  • Nodes on the green boundary - lie on the inner boundary of and are interior to .

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

1:  On the mesh , where FDM is used, update the Finite Difference (FD) solution.
2:  On the mesh , where FEM is used, update the Finite Element (FE) solution.
3:  Copy FE solution obtained at nodes (nodes on the green boundary of Figure 1) as a boundary condition on the inner boundary for the FD solution in .
4:  Copy FD solution obtained at nodes (nodes on the blue boundary of Figure 1) as a boundary condition for the FE solution on of .
5:  Apply boundary condition at at the red boundary of .
Algorithm 1 Domain decomposition process for hybrid FE/FD scheme

By (3), at the overlapping nodes between and . 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, combined, finite element-finite difference schemes to solve the model problem (6). To do this, first we present the finite element scheme to solve (6) in entire . This induces the finite element scheme in . Then, we derive the finite difference scheme in , when the domain decomposition FE/FD structure is applied to solve (6) in .


The computational schemes derived in this section are explicit and therefore, for their converegence, the CFL condition below (see ,e.g. [6, 7]) should be satisfied


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

In the sequel, we denote the inner product of by , and the corresponding norm by . The scalar inner product in we denote by , and the associated norm by . Further, we let to be the boundary of , the inner boundary of and the outer boundary of .

4.1 Finite element discretization in

First we derive finite element scheme to solve the model problem (6) in whole . Next, we discretize in two steps: (i) the spatial discretization using a partition of into elements , where is a mesh function defined as , representing the local diameter of elements. We also denote by a partition of the boundary into boundaries of the elements such that, at least one of the vertices of these elements belong to . (ii) As for temporal discretization, we let be a uniform partition of the time interval into subintervals of length As usual, we also assume a minimal angle condition on elements in . To formulate the finite element method for (6) in , we introduce the finite element space for each component of the electric field defined by

where denote the set of piecewise-linear functions on . Setting we define (resp. ) to be the usual

-interpolant of

(resp. ) in (6). Also, because of the Dirichlet boundary data in (6) we need to choose the test function space as


Then, recalling (4), and the fact that the spatial semi-discrete problem in reads:

Find such that ,


We note that

implies that


Recalling (8), vanishing test functions at the boundary yields , and the final weak formulation for the semi-discrete problem in is: Find such that ,


For the, reflexive, inhomogeneous boundary condition, see the FE scheme for .

To get fully discrete scheme for (6) we apply time discretization to (9) approximating , denoted by , where we use the central difference scheme, for :


Multiplying both sides of (12) by we get, that


Rearranging terms in (13) we get the following scheme: Given the approximate initial data, and , find such that ,


4.2 Finite element discretization in

To solve the model problem (6) via the domain decomposition FE/FD method, we use the split , see Figure 1. Thus in we use FEM to solve the equation


Here, is the restriction of the solution obtained by the FDM in to and therefore the test functions are not vanishing at the boundary and hence the term corresponding to in (9) will be appearing in the weak formulation.

To formulate the finite element method for (15) in , mimiking (9), we introduce the finite element space for each component of the electric field defined by

Setting we define , , and to be the usual -interpolants of , , and , respectively, in . Then, similar to the FE scheme for , we get the following finite element scheme for : Given , , and , find such that


A corresponding fully discrete problem in reads as follows:

Given , , , , and ; find such that



Note that, in (15), Dirichlet boundary condition can be considered as well.

4.3 Fully discrete FE scheme for the electric field in

We expand the functions in terms of the standard continuous piecewise linear functions in space as


where denote unknown coefficients at and the spatial mesh point . Then, substituting of (18) in (14), and setting , we obtain the linear system of equations:


Note that, unlike (14), now the contributions at boundary of the element appear in . Here, are the block mass matrices in space, are the block stiffness matrices in space, denote the nodal values of , and is a uniform time step. Now we define the mapping from the reference element onto such that and let be the piecewise linear local basis function on such that . Then, the explicit formulas for the entries in system of equations (19), for each element , can be written as:


where , and , denote the , and , scalar products on and , respectively. Note that here is only the part of the boundary of element that lies at .

To obtain fully explicit scheme we approximate with the lumped mass matrix , (see [14, 18, 6] for the details corresponding to the Maxwell’s system (2)). Next, we multiply (20) by , and get the following explicit, fully discrete method in :


4.4 Fully discrete scheme for the electric field in

As in the fully discrete FE scheme (19) in , we obtain fully discrete FE scheme in in the domain decomposition setting: Expanding the functions of via the continuous piecewise linear functions in space as in (18), and then substituting them in (17), (with , and ), we get the linear system of equations:


Here, is the block mass matrices in space, restricted to , otherwise the same as in (19), are the block stiffness matrices in space as in (19),

is the assembled load vector ,

denote the nodal values of , is the time step. All quantities are for . Defining the mapping for the reference element in the mesh generated in as in the previous section, the formulas for entries of all matrices in the system (22) are the same as those in (20), and the entries of load vector are computed as


Again, approximating with the lumped mass matrix , we obtain the following fully explicit scheme:


4.5 Finite difference formulation

We recall now that from conditions (3) it follows that in the function . This means that in for the model problem (2) the forward problem will be


Using standard finite difference discretization of the equation (25) in we obtain the following explicit scheme for the solution of the forward problem:


In the, system of, equations above, is the solution at the time iteration at the discrete point , is the time step, and is the discrete Laplacian.

Note that, in (29), the Dirichlet boundary consitions can be considered as well.

5 Stability

In this section we derive stability estimates for the semi-discrete approximations. For stability in these estimates are extensions of the stability approach derived for the continuous problem in [4]. As for the stability in we get slightly different norms involving contributions corresponding to the reflexive boundary: . We use discrete version of a triple norm induced by the weak variational formulation of (6), where we use the relation (10) (which is not necessary in the continuous case where , however, in general ):

Find such that



In general, in non-divergent free case, the bilinear form induced by (30) is not coercive. Further -conforming finite element may result in spurious solutions. A remedy is through modifying the equation by adding a gauge constrain of Coulomb-type, see, e.g. [21] and [23]. This is supplied by the ”zero”-term: , in (6), which we add in the continuous variational formulation in (9). This, however, is not necessarily true in the discrete forms, e.g. in (16), where most likely Taking in (30), (we used the boundary condition on ), yields


which, due to the fact that is independent of , can be rewritten as


Proposition [4]

Let be a bounded domain with piecewise linear boundary . Then, the equation (6) has a unique solution . Further Let and , then there is a constant such that, , the following stability estimate holds true


The estimate (33) is proved in [4], Theorem 4.1 by setting and . Integrating (33) over the time interval we get the desired result. We omit the details. ∎

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

5.1 Stability estimate for the semi-discrete problem in

The stability for the semi-discrete problem in is basically as in the continuous case above where all :s are replaced by with some relevant assumptions in the discrete data, viz.


Assume that the interpolants of the data and : and satisfy the regularity conditions and , then for each




5.2 Stability of the semi-discrete problem in

The stability of the semi-discrete problem in , relying on the variational formulation (16), and due to the appearance of the data function , is slightly different from (34). We rewrite (16), in view of (10), and with as: Given , , and , find such that


To deal with the -term we rewrite (36) in its original form as (9) for and with :