Numerical validation of an explicit P1 finite-element scheme for Maxwell's equations in a polygon with variable permittivity away from its boundary

TL;DR

This paper numerically validates an explicit finite-element scheme for Maxwell's equations with variable permittivity in a polygon, confirming optimal convergence under the CFL condition through 2D experiments.

Contribution

It provides empirical validation of an explicit P1 finite-element scheme for Maxwell's equations with variable permittivity, extending prior theoretical analysis.

Findings

Convergence results are optimal when CFL condition is satisfied.

Numerical experiments confirm the theoretical reliability of the scheme.

The scheme effectively handles variable permittivity away from the boundary.

Abstract

This paper is devoted to the numerical validation of an explicit finite-difference scheme for the integration in time of Maxwell's equations in terms of the sole electric field, using standard linear finite elements for the space discretization. The rigorous reliability analysis of this numerical model was the object of another authors' arXiv paper. More specifically such a study applies to the particular case where the electric permittivity has a constant value outside a sub-domain, whose closure does not intersect the boundary of the domain where the problem is defined. Our numerical experiments in two-dimension space certify that the convergence results previously derived for this approach are optimal, as long as the underlying CFL condition is satisfied.

Tables

Table 1. Table 1: Maximum over the time steps of relative errors in the L 2 subscript 𝐿 2 L_{2} -norm, in the H 1 superscript 𝐻 1 H^{1} -seminorm and in the L 2 superscript 𝐿 2 L^{2} -norm of the time derivative for mesh sizes h l = 2 − l , l = 1 , … , 6 formulae-sequence subscript ℎ 𝑙 superscript 2 𝑙 𝑙 1 … 6 h_{l}=2^{-l},l=1,...,6 taking m = 2 𝑚 2 m=2 in ( 4.12 )

$l$	$nel$	$nno$	$e_l^1$	$e_{l-1}^1/e_l^1$	$e_l^2$	$e_{l-1}^2/e_l^2$	$e_l^3$	$e_{l-1}^3/e_l^3$
1	8	9	0.054247		0.2767		1.0789
2	32	25	0.013902	3.902100	0.1216	2.2755	0.4811	2.2426
3	128	81	0.003706	3.751214	0.0532	2.2857	0.2544	1.8911
4	512	289	0.000852	4.349765	0.0234	2.2735	0.1279	1.9891
5	2048	1089	0.000229	3.720524	0.0121	1.9339	0.0641	1.9953
6	8192	4225	0.000059	3.881356	0.0061	1.9836	0.0321	1.9969

Table 2. Table 2: Maximum over the time steps of relative errors in the L 2 subscript 𝐿 2 L_{2} -norm, in the H 1 superscript 𝐻 1 H^{1} -seminorm and in the L 2 superscript 𝐿 2 L^{2} -norm of the time derivative for mesh sizes h l = 2 − l , l = 1 , … , 6 formulae-sequence subscript ℎ 𝑙 superscript 2 𝑙 𝑙 1 … 6 h_{l}=2^{-l},l=1,...,6 taking m = 3 𝑚 3 m=3 in ( 4.12 )

$l$	$nel$	$nno$	$e_l^1$	$e_{l-1}^1/e_l^1$	$e_l^2$	$e_{l-1}^2/e_l^2$	$e_l^3$	$e_{l-1}^3/e_l^3$
1	8	9	0.043394		0.2784		1.0869
2	32	25	0.011451	3.789538	0.1098	2.5355	0.5305	2.0488
3	128	81	0.003343	3.425366	0.06	1.83	0.2586	2.0514
4	512	289	0.000781	4.385873	0.0248	2.4194	0.1306	1.9801
5	2048	1089	0.000202	3.866337	0.0119	2.0840	0.0654	1.9969
6	8192	4225	0.000052	3.884615	0.0059	2.0169	0.0327	2

Table 3. Table 3: Maximum over the time steps of relative errors in the L 2 subscript 𝐿 2 L_{2} -norm, in the H 1 superscript 𝐻 1 H^{1} -seminorm and in the L 2 superscript 𝐿 2 L^{2} -norm of the time derivative for mesh sizes h l = 2 − l , l = 1 , … , 6 formulae-sequence subscript ℎ 𝑙 superscript 2 𝑙 𝑙 1 … 6 h_{l}=2^{-l},l=1,...,6 taking m = 6 𝑚 6 m=6 in ( 4.12 )

$l$	$nel$	$nno$	$e_l^1$	$e_{l-1}^1/e_l^1$	$e_l^2$	$e_{l-1}^2/e_l^2$	$e_l^3$	$e_{l-1}^3/e_l^3$
1	8	9	0.054228		0.2837		1.1120
2	32	25	0.012241	4.430030	0.0906	3.1313	0.4937	2.2524
3	128	81	0.002973	4.117389	0.0408	2.2206	0.2665	1.8525
4	512	289	0.000590	5.038983	0.0150	2.7200	0.1335	1.9963
5	2048	1089	0.000163	3.619631	0.0079	1.8987	0.0667	2.0015
6	8192	4225	0.000043	3.790698	0.0040	1.9750	0.0334	1.9970

Table 4. Table 4: Maximum over the time steps of relative errors in the L 2 subscript 𝐿 2 L_{2} -norm, in the H 1 superscript 𝐻 1 H^{1} -seminorm and in the L 2 superscript 𝐿 2 L^{2} -norm of the time derivative for mesh sizes h l = 2 − l , l = 1 , … , 6 formulae-sequence subscript ℎ 𝑙 superscript 2 𝑙 𝑙 1 … 6 h_{l}=2^{-l},l=1,...,6 taking m = 7 𝑚 7 m=7 in ( 4.12 )

$l$	$nel$	$nno$	$e_l^1$	$e_{l-1}^1/e_l^1$	$e_l^2$	$e_{l-1}^2/e_l^2$	$e_l^3$	$e_{l-1}^3/e_l^3$
1	8	9	0.054224		0.5710		1.1208
2	32	25	0.012483	4.343828	0.1505	3.7940	0.5024	2.2309
3	128	81	0.002751	4.537623	0.0686	2.1939	0.2688	1.8690
4	512	289	0.000627	4.387559	0.0240	2.8583	0.1339	2.0075
5	2048	1089	0.000158	3.968354	0.0114	2.1053	0.0669	2.0015
6	8192	4225	0.000040	3.949999	0.0057	2	0.0334	2.0030

Equations

$$\begin{array}[]{ll}\varepsilon\partial_{tt}{\bf e}+\nabla\times\nabla\times{\bf e}={\bf 0}&\mbox{ in }\Omega\times(0,T),\\ {\bf e}(\cdot,0)={\bf e}_{0}(\cdot),\mbox{ and }\partial_{t}{\bf e}(\cdot,0)={\bf e}_{1}(\cdot)&\mbox{ in }\Omega,\\ {\bf e}={\bf 0}&\mbox{ on }\partial\Omega\times(0,T),\\ \nabla\cdot(\varepsilon{\bf e})={\bf 0}&\mbox{ in }\Omega.\end{array}$$

$${\left({\partial }_{tt}\mathbf{e},\mathbf{v}\right)}_{\epsilon }+(\nabla \mathbf{e},\nabla \mathbf{v})+(\nabla \cdot \epsilon \mathbf{e},\nabla \cdot \mathbf{v})-(\nabla \cdot \mathbf{e},\nabla \cdot \mathbf{v})=0\forall t\in (0,T).$$

$$\begin{array}[]{l}\left(\partial_{tt}{\bf e}_{h},{\bf v}\right)_{\varepsilon}+(\nabla{\bf e}_{h},\nabla{\bf v})+(\nabla\cdot[\varepsilon{\bf e}_{h}],\nabla\cdot{\bf v})-(\nabla\cdot{\bf e}_{h},\nabla\cdot{\bf v})=0,\\ \\ {\bf e}_{h}(\cdot,0)={\bf e}_{0h}(\cdot)\mbox{ and }\partial_{t}{\bf e}_{h}(\cdot,0)={\bf e}_{1h}(\cdot)\mbox{ in }\Omega.\end{array}$$

$$\begin{array}[]{l}\displaystyle\left(\frac{{\bf e}_{h}^{k+1}-2{\bf e}_{h}^{k}+{\bf e}_{h}^{k-1}}{\tau^{2}},{\bf v}\right)_{\varepsilon}+(\nabla{\bf e}_{h}^{k},\nabla{\bf v})+(\nabla\cdot\varepsilon{\bf e}_{h}^{k},\nabla\cdot{\bf v})-(\nabla\cdot{\bf e}_{h}^{k},\nabla\cdot{\bf v})=0\;\forall{\bf v}\in{\bf V}_{h},\\ \\ {\bf e}_{h}^{0}={\bf e}_{0h}\mbox{ and }{\bf e}_{h}^{1}={\bf e}_{h}^{0}+\tau{\bf e}_{1h}\mbox{ in }\Omega.\end{array}$$

$${\int }_K\epsilon {\mathbf{u}}_{|K}\cdot {\mathbf{v}}_{|K}d\mathbf{x}\approx \epsilon (G_K)area(K)\sum _{i=1}^3\frac{\mathbf{u}(S_{K,i})\cdot \mathbf{v}(S_{K,i})}{3},$$

$$\begin{array}[]{l}\displaystyle\left(\frac{{\bf e}_{h}^{k+1}-2{\bf e}_{h}^{k}+{\bf e}_{h}^{k-1}}{\tau^{2}},{\bf v}\right)_{\varepsilon_{h},h}+(\nabla{\bf e}_{h}^{k},\nabla{\bf v})+(\nabla\cdot\varepsilon{\bf e}_{h}^{k},\nabla\cdot{\bf v})-(\nabla\cdot{\bf e}_{h}^{k},\nabla\cdot{\bf v})=0\;\forall{\bf v}\in{\bf V}_{h},\\ \\ {\bf e}_{h}^{0}={\bf e}_{0h}\mbox{ and }{\bf e}_{h}^{1}={\bf e}_{h}^{0}+\tau{{\bf e}_{1}}_{h}\mbox{ in }\Omega.\end{array}$$

$$\eta :=2+|\epsilon {|}_{1,\infty }+2|\epsilon {|}_{2,\infty };$$

$$\Vert \nabla v\Vert \le C{h}^{-1}\Vert v{\Vert }_{{\epsilon }_h,h}\text{mbox}forallv\in V_h,$$

$$\tau \le h/\nu \text{mbox}with\nu =C(1+3\Vert \epsilon -1{\Vert }_{\infty }{)}^{1/2}.$$

$$\boxed{\begin{array}[]{l}\displaystyle\displaystyle\max_{1\leq m\leq N-1}\left\|[\partial_{t}({\bf e}_{h}-{\bf e})]^{m+1/2}\right\|+\displaystyle\displaystyle\max_{2\leq m\leq N}\|\nabla({\bf e}_{h}^{m}-{\bf e}^{m})\|\\ \leq{\mathcal{C}}(\tau+h+h^{2}/\tau)\displaystyle\left\{\|{\bf e}\|_{H^{4}[\Omega\times(0,T)]}+|{\bf e}_{0}|_{2}+|{\bf e}_{1}|_{2}\right\}.\mbox{ \rule{5.69054pt}{5.69054pt}}\end{array}}$$

$$\begin{array}[]{ll}\varepsilon\partial_{tt}{\bf e}-\nabla^{2}{\bf e}-\nabla\nabla\cdot(\varepsilon-1){\bf e}={\bf f}&\mbox{ in }\Omega\times(0,T),\\ {\bf e}(\cdot,0)={\bf 0}\mbox{ and }\partial_{t}{\bf e}(\cdot,0)={\bf 0}&\mbox{ in }\Omega,\\ {\bf e}={\bf 0}&\mbox{ on }\partial\Omega\times(0,T).\end{array}$$

$$\begin{array}{rl}{e}_1& =\frac{1}{\epsilon }2\pi {\sin }^2\pi x\cos \pi y\sin \pi y\frac{t^2}{2},\\ e_2& =-\frac{1}{\epsilon }2\pi {\sin }^2\pi y\cos \pi x\sin \pi x\frac{t^2}{2}\end{array}$$

$$\varepsilon(x,y)=\left\{\begin{array}[]{ll}1+\sin^{m}\pi(2x-0.5)\cdot\sin^{m}\pi(2y-0.5)&\textrm{in $[0.25,0.75]\times[0.25,0.75]$},\\ 1&\textrm{otherwise},\end{array}\right.$$

$$\begin{array}{rl}{e}_l^1& =\frac{\underset{1\le k\le N}{\max }\Vert {\mathbf{e}}^k-{\mathbf{e}}_h^k\Vert }{\underset{1\le k\le N}{\max }\Vert {\mathbf{e}}^k\Vert },\\ e_l^2& =\frac{\underset{1\le k\le N}{\max }\Vert \nabla ({\mathbf{e}}^k-{\mathbf{e}}_h^k)\Vert }{\underset{1\le k\le N}{\max }\Vert \nabla {\mathbf{e}}^k\Vert },\\ e_l^3& =\frac{\underset{1\le k\le N-1}{\max }\Vert \{{\partial }_t(\mathbf{e}-{\mathbf{e}}_h){\}}^{k+1/2}\Vert }{\underset{1\le k\le N-1}{\max }\Vert \{{\partial }_t\mathbf{e}{\}}^{k+1/2}\Vert }.\end{array}$$

Full text

1 Motivation

The purpose of this article is to provide a numerical validation of an explicit $P_{1}$ finite element solution scheme of hyperbolic Maxwell’s equations for the electric field with constant dielectric permittivity in a neighborhood of the boundary of the computational domain. This numerical model was thoroughly studied in [2] from the theoretical point of view. The focus of that paper was on the case of three-dimensional domains on whose boundary absorbing boundary conditions are enforced. However as pointed out therein, all the underlying analytical results trivially extend to the two-dimensional case and to other types of boundary conditions such as Dirichlet and Neumann conditions. Actually in our experiments we consider only two-dimensional test-problems, in which Dirichlet boundary conditions are prescribed.

The standard continuous $P_{1}$ FEM is a tempting possibility to solve Maxwell’s equations, owing to its simplicity. It is well known however that, for different reasons, this method is not always well suited for this purpose. The first reason is that in general the natural function space for the electric field is not the Sobolev space ${\bf H}^{1}$, but rather in the space ${\bf H}(curl)$. Another issue difficult to overcome with continuous Lagrange finite elements is the prescription of the zero tangential-component boundary conditions for the electric field, which hold in many important applications. All this motivated the proposal by Nédélec about four decades ago of a family of ${\bf H}(curl)$-conforming methods to solve these equations (cf. [24]). These methods are still widely in use, as much as other approaches well adapted to such specific conditions (see e.g. [1], [13] and [27]). A comprehensive description of finite element methods for Maxwell’s equations can be found in [21].

There are situations however in which the $P_{1}$ finite element method does provide an inexpensive and reliable way to solve the Maxwell’s equations. In this work we address one of such cases, characterized by the fact that the electric permittivity is constant in a neighborhood of the whole boundary of the domain of interest. This is because, at least in theory, whenever the electric permittivity is constant, the Maxwell’s equations simplify into as many wave equations as the space dimension under consideration. More precisely here we show by means of numerical examples that, in such a particular case, a space discretization with conforming linear elements, combined with a straightforward explicit finite-difference scheme for the time integration, gives rise to optimal approximations of the electric field, as long as a classical CFL condition is satisfied.

Actually this work is strongly connected with studies presented in [3, 4] for a combination of the finite difference method in a sub-domain with constant permittivity with the finite element method in the complementary sub-domain. As pointed out above, the Maxwell’s equations reduces to the wave equation in the former case. Since the analysis of finite-difference methods for this type of equation is well established, only an explicit $P_{1}$ finite element scheme for Maxwell’s equations is analyzed in this paper.

In [3, 4] a stabilized domain-decomposition finite-element/finite-difference approach for the solution of the time-dependent Maxwell’s system for the electric field was proposed and numerically verified. In these works [3, 4] different manners to handle a divergence-free condition in the finite-element scheme were considered. The main idea behind the domain decomposition methods in [3, 4] is that a rectangular computational domain is decomposed into two sub-domains, in which two different types of discretizations are employed, namely, the finite-element domain in which a classical $P_{1}$ finite element discretization is used, and the finite-difference domain, in which the standard five- or seven-point finite difference scheme is applied, according to the space dimension. The finite element domain lies strictly inside the finite difference domain, in such a way that both domains overlap in two layers of structured nodes. First order absorbing boundary conditions [16] are enforced on the boundary of the computational domain, i.e. on the outer boundary of the finite-difference domain. In [3, 4] it was assumed that the dielectric permittivity function is strictly positive and has a constant value in the overlapping nodes as well as in a neighborhood of the boundary of the domain. An explicit scheme was used both in the finite-element and finite-difference domains.

We recall that for a stable finite-element solution of Maxwell’s equation divergence-free edge elements are the most satisfactory from a theoretical point of view  [24, 21]. However, the edge elements are less attractive for solving time-dependent problems, since a linear system of equations should be solved at every time iteration. In contrast, $P_{1}$ elements can be efficiently used in a fully explicit finite element scheme with lumped mass matrix [15, 20]. On the other hand it is also well known that the numerical solution of Maxwell’s equations with nodal finite elements can result in unstable spurious solutions [22, 25]. Nevertheless a number of techniques are available to remove them, and in this respect we refer for example to [17, 18, 19, 23, 25]. In the current work, similarly to [3, 4], the spurious solutions are removed from the finite element scheme by adding the divergence-free condition to the model equation for the electric field. Numerical tests given in [4] demonstrate that spurious solutions are removable indeed, in case an explicit $P_{1}$ finite-element solution scheme is employed.

Efficient usage of an explicit $P_{1}$ finite-element scheme for the solution of coefficient inverse problems (CIPs), in the particular context described above was made evident in [5]. In many algorithms aimed at solving electromagnetic CIPs, a qualitative collection of experimental measurements is necessary on the boundary of a computational domain, in order to determine the dielectric permittivity function therein. In this case, in principle the numerical solution of the time-dependent Maxwell’s equations is required in the entire space $\mathbb{R}^{3}$ (see e.g. [5, 6, 7, 8, 9, 10], but instead it can be more efficient to consider Maxwell’s equations with a constant dielectric permittivity in a neighborhood of the boundary of a computational domain. The explicit $P_{1}$ finite-element scheme considered in this work was numerically tested in the solution of the time-dependent Maxwell’s system in both two- and three-dimensional geometry (cf. [4]). It was also combined with a few algorithms to solve different CIPs for determining the dielectric permittivity function in connection with the time-dependent Maxwell’s equations, using both simulated and experimentally generated data (see [6, 7, 8, 9, 10]). In short, the formal reliability analysis of such a method conducted in this work, corroborates the previously observed adequacy of this numerical approach.

An outline of this paper is as follows: In Section 2 we describe in detail the model problem being solved, and give its equivalent variational form. In Section 3 we set up the discretizations of the model problem in both space and time, and recall the main results of the reliability analysis conducted in [2] for the underlying numerical model. Section 4 is devoted to the numerical experiments that validate such results. We conclude in Section 5 with a few comments.

2 The model problem

The particular form of Maxwell’s equations for the electric field ${\bf e}=(e_{1},e_{2})$ in a bounded domain $\Omega$ of $\Re^{2}$ with boundary $\partial\Omega$ that we deal with in this work is as follows. First we consider that $\Omega=\bar{\Omega}_{in}\cup\Omega_{out}$, where $\Omega_{in}$ is an interior open set whose boundary does not intersect $\partial\Omega$ and $\Omega_{out}$ is the complementary set of $\bar{\Omega}_{in}$ with respect to $\Omega$. Now in case ${\bf e}$ satisfies (homogeneous) Dirichlet boundary conditions, we are given ${\bf e}_{0}\in[H^{1}(\Omega)]^{2}$ and ${\bf e}_{1}\in{\bf H}(div,\Omega)$ satisfying $\nabla\cdot(\varepsilon{\bf e}_{0})=\nabla\cdot(\varepsilon{\bf e}_{1})=0$ where $\varepsilon$ is the electric permittivity. $\varepsilon$ is assumed to belong to $W^{2,\infty}(\Omega)$ and to fulfill $\varepsilon\equiv 1$ in $\Omega_{out}$ and $\varepsilon\geq 1$. Incidentally, throughout this article we denote the standard semi-norm of $C^{m}(\bar{\Omega})$ by $|\cdot|_{m,\infty}$ for $m>0$ and the standard norm of $C^{0}(\bar{\Omega})$ by $\|\cdot\|_{0,\infty}$.

In doing so, the problem to solve is:

[TABLE]

Remark 2.1**.**

As pointed out above the analysis carried out in [2] applies to the case where absorbing boundary conditions $\partial_{n}{\bf e}=-\partial_{t}{\bf e}$ are prescribed, where $\partial_{n}{\bf e}$ represents the outer normal derivative of ${\bf e}$ on $\partial\Omega$. This choice was motivated by the fact that they correspond to practical situations addressed in [6, 7, 8, 9, 10].

We next set (4.11) in variational form. With this aim we denote the standard inner product of $[L^{2}(\Omega)]^{2}$ by $(\cdot,\cdot)$ and the corresponding norm by $\parallel\{\cdot\}\parallel$. Further, for a given non-negative function $\omega\in L^{\infty}(\Omega)$ we introduce the weighted $L^{2}(\Omega)$-semi-norm $\|\{\cdot\}\|_{\omega}:=\sqrt{\int_{\Omega}|\omega||\{\cdot\}|^{2}d{\bf x}}$, which is actually a norm if $\omega\neq 0$ everywhere in $\bar{\Omega}$. We also introduce, the notation $({\bf a},{\bf b})_{\omega}:=\int_{\Omega}\omega{\bf a}\cdot{\bf b}d{\bf x}$ for two fields ${\bf a},{\bf b}$ which are square integrable in $\Omega$. Notice that if $\omega$ is strictly positive this expression defines an inner product associated with the norm $\|\{\cdot\}\|_{\omega}$.

Then requiring that ${\bf e}_{|t=0}={\bf e}_{0}$ and $\{\partial_{t}{\bf e}\}_{|t=0}={\bf e}_{1}$ and ${\bf e}=0$ on $\partial\Omega\times[0,T]$, we write for all ${\bf v}\in[H_{0}^{1}(\Omega)]^{2}$,

[TABLE]

We recall that the equivalence of problem (2.2) with Maxwell’s equations (4.11) was established in [2].

3 The numerical model

Henceforth we restrict our studies to the case where $\Omega$ is a polygon.

3.1 Space semi-discretization

Let $V_{h}$ be the usual $P_{1}$ FE-space of continuous functions related to a mesh ${\mathcal{T}}_{h}$ fitting $\Omega$, consisting of triangles with maximum edge length $h$, belonging to a quasi-uniform family of meshes (cf. [12]).

Setting ${\bf V}_{h}:=[V_{h}\cap H^{1}_{0}(\Omega)]^{2}$ we define ${\bf e}_{0h}$ (resp. ${\bf e}_{1h}$) to be the usual ${\bf V}_{h}$-interpolate of ${\bf e}_{0}$ (resp. ${\bf e}_{1}$). Then the semi-discretized problem is space we wish to solve writes,

*Find *${\bf e}_{h}\in{\bf V}_{h}$ * such that $\forall{\bf v}\in{\bf V}_{h}$ *

[TABLE]

3.2 Full discretization

To begin with we consider a natural centered time-discretization scheme to solve (3.3), namely: Given a number $N$ of time steps we define the time increment $\tau:=T/N$. Then we approximate ${\bf e}_{h}(k\tau)$ by ${\bf e}_{h}^{k}\in{\bf V}_{h}$ for $k=1,2,...,N$ according to the following scheme for $k=1,2,\ldots,N-1$:

[TABLE]

Owing to its coupling with ${\bf e}_{h}^{k}$ and ${\bf e}_{h}^{k-1}$ on the left hand side of (3.4), ${\bf e}_{h}^{k+1}$ cannot be determined explicitly by (3.4) at every time step. In order to enable an explicit solution we resort to the classical mass-lumping technique. We recall that for a constant $\varepsilon$ this consists of replacing on the left hand side the inner product $({\bf u},{\bf v})_{\varepsilon}$ by an inner product $({\bf u},{\bf v})_{\varepsilon.h}$, using the trapezoidal rule to compute the integral of $\int_{K}\varepsilon{\bf u}_{|K}\cdot{\bf v}_{|K}d{\bf x}$ (resp. $\int_{K\cap\partial\Omega}{\bf u}_{|K}\cdot{\bf v}_{|K}dS$), for every element $K$ in ${\mathcal{T}}_{h}$, where ${\bf u}$ stands for ${\bf e}_{h}^{k+1}-2{\bf e}_{h}^{k}+{\bf e}_{h}^{k-1}$. It is well-known that in this case the matrix associated with $(\varepsilon{\bf e}_{h}^{k+1},{\bf v})_{h}$ for ${\bf v}\in{\bf V}_{h}$, is a diagonal matrix. In our case $\varepsilon$ is not constant, but the same property will hold if we replace in each element $K$ the integral of $\varepsilon{\bf u}_{|K}\cdot{\bf v}_{|K}$ in a triangle $K\in{\mathcal{T}}_{h}$ as follows:

[TABLE]

where $S_{K,i}$ are the vertexes of $K$, $i=1,2,3$, $G_{K}$ is the centroid of $K$.

Before pursuing we define the auxiliary function $\varepsilon_{h}$ whose value in each $K\in{\mathcal{T}}_{h}$ is constant equal to $\varepsilon(G_{K})$. Then still denoting the approximation of ${\bf e}_{h}(k\tau)$ by ${\bf e}_{h}^{k}$, for $k=1,2,...,N$ we determine ${\bf e}_{h}^{k+1}$ by,

[TABLE]

3.3 Convergence results

Recalling the assumption that $\varepsilon\in W^{2,\infty}(\Omega)$ we first set

[TABLE]

Next we recall the classical inverse inequality (cf. [12]) together with a result in [11] according to which,

[TABLE]

where $C$ is a mesh-independent constant.

Now we assume that $\tau$ satisfies the following CFL-condition:

[TABLE]

. We further assume that the solution ${\bf e}$ to equation (4.11) belongs to $[H^{4}\{\Omega\times(0,T)\}]^{2}$.

Let us define a function ${\bf e}_{h}$ in $\bar{\Omega}\times[0,T]$ whose value at $t=k\tau$ equals ${\bf e}_{h}^{k}$ for $k=1,2,\ldots,N$ and that varies linearly with $t$ in each time interval $([k-1]\tau,k\tau)$, in such a way that $\partial_{t}{\bf e}_{h}({\bf x},t)=\displaystyle\frac{{\bf e}^{k}_{h}({\bf x})-{\bf e}^{k-1}_{h}({\bf x})}{\tau}$ for every ${\bf x}\in\bar{\Omega}$ and $t\in([k-1]\tau,k\tau)$. We also define ${\bf a}^{m+1/2}(\cdot)$ for any field ${\bf a}(\cdot,t)$ to be ${\bf a}(\cdot,[m+1/2]\tau)$.

Then denoting by $|\cdot|_{m}$ the standard semi-norm of Sobolev space $H^{m}(\Omega)$ for $m\in{\mathcal{N}}$, according to [2] we have:

Provided the CFL condition (3.8) is fulfilled and $\tau$ also satisfies $\tau\leq 1/[2\eta]$, under the above regularity assumption on ${\bf e}$, there exists a constant ${\mathcal{C}}$ depending only on $\Omega$, $\varepsilon$ and $T$ such that,

[TABLE]

(3.9) means that, as long as $\tau$ varies linearly with $h$, first order convergence of scheme (3.5) in terms of either $\tau$ or $h$ holds in the sense of the norms on the left hand side of (3.9).

4 Numerical validation

We perform numerical tests in time $(0,T)=(0,0.5)$ in the computational domain $\Omega=[0,1]\times[0,1]$ for the model problem in two-dimension space, namely

[TABLE]

for the electric field ${\bf e}=(e_{1},e_{2})$.

The source data ${\bf f}$ (the right hand side) is chosen such that the function

[TABLE]

is the exact solution of the model problem (4.10) In (4.11) the function $\varepsilon$ is defined to be,

[TABLE]

where $m$ is an integer greater than one. In Figure 1 the function $\varepsilon$ are illustrated for different values of $m$.

The solution given by (4.11) satisfies homogeneous initial conditions together with homogeneous Dirichlet conditions on the boundary $\partial\Omega$ of the square $\Omega$ for every time $t$. In our computations we used the software package Waves [31] only for the finite element method applied to the solution of the model problem (4.10). We note that this package was also used in [4] to solve the the same model problem (4.10) by a domain decomposition FEM/FDM method.

We discretized the computational domain $\Omega\times(0,T)$ denoting by ${K_{h}}_{l}=\{K\}$ a partition of the spatial domain $\Omega$ into triangles $K$ of sizes $h_{l}=2^{-l},l=1,...,6$. We let $J_{\tau_{l}}$ be a partition of the time domain $(0,T)$ into time intervals $J=(t_{k-1},t_{k}]$ of uniform length $\tau_{l}$ for a given number of time intervals $N,l=1,...,6$. We choose the time step $\tau_{l}=0.025\times 2^{-l},\;l=1,...,6$, which provides numerical stability for all meshes.

We performed numerical tests taking $m=2,...,9$ in (4.12) and computed the maximum value over the time steps of the relative errors measured in the function $L_{2}$-norm and the $H^{1}$-semi-norm and in the $L_{2}$ norm for the time-derivative, respectively represented by

[TABLE]

Here ${\bf e}$ is the exact solution given by (4.11) and ${\bf e}_{h}$ is the computed solution, while $N=T/\tau_{l}$.

In Tables 1-4 method’s convergence in these three senses is observed taking $m=2,3,6,7$ in (4.12).

Figure 2 shows convergence of our numerical scheme based on a $P_{1}$ space discretization, taking the function $\varepsilon$ defined by (4.12) with $m=2$ (on the left) and $m=3$ (on the right) for $\varepsilon(x)$. Similar convergence results are presented in Figures 3 and 4 taking $m=6,7,8,9$ in (4.12).

Observation of these tables and figures clearly indicates that our scheme behaves like a first order method in the (semi-)norm of $L^{\infty}[(0,T);H^{1}(\Omega)]$ for ${\bf e}$ and in the norm of $L^{\infty}[(0,T);L^{2}(\Omega)]$ for $\partial_{t}{\bf e}$ for all the chosen values of $m$. As far as the values $m=6$ and $m=8$ are concerned this perfectly conforms to the a priori error estimates established in [2]. However those tables and figures also show that such theoretical predictions extend to cases not considered in our analysis such as $m=2$ and $m=3$, in which the regularity of the exact solution is lower than assumed, or yet in the cases $m=7$ and $m=9$, in which the minimum of $\varepsilon$ is not attained on the boundary. Otherwise stated some of our assuptions seem to be of academic interest only and a lower regularity of the solution such as $\{H^{2}[\Omega\times(0,T)]\}^{2}$ should be sufficient to attain optimal first order convergence in both senses. On the other hand second-order convergence can be expected from our scheme in the norm of $L^{\infty}[(0,T);L^{2}(\Omega)]$ for ${\bf e}$, according to Tables 1-4 and Figures 2-4.

5 Summary

In this work we validated the reliability analysis conducted in [2] for a numerical scheme to solve Maxwell’s equations of electromagnetism, combining an explicit finite difference time discretization with a lumped-mass $P1$ finite element space discretization. The scheme is effective in the particular case where the dielectric permittivity is constant in a neighborhood of the boundary of the spacial domain. After presenting the problem under consideration for the electric field we supplied the detailed description of such a scheme and recalled the a priori error estimates that hold for the latter under suitable regularity assumptions specified in [2]. Then we showed by means of numerical experiments performed for a test-problem in two-dimension space with known exact solution, that the convergence results given in [2] are confirmed in practice. Furthermore we presented convincing evidence that such theoretical predictions extend to solutions with much lower regularity than the one assumed in our analysis. Similarly optimal second-order convergence is observed in a norm other than those in which convergence was formally established. In short we undoubtedly indicated that Maxwell’s equations can be efficiently solved with classical conforming linear finite elements in some relevant particular cases, among which lie the model problem (4.10).

Acknowledgment: The research of the first author is supported by the Swedish Research Council grant VR 2018-03661. The second author gratefully acknowledges the financial support provided by CNPq/Brazil through grant 307996/2008-5.