Semi-discrete finite element approximation applied to Maxwell's equations in nonlinear media

TL;DR

This paper investigates a semi-discrete finite element method for Maxwell's equations in nonlinear Kerr media, providing error estimates for Nédélec element approximations.

Contribution

It introduces a priori error estimates for finite element approximations of Maxwell's equations in nonlinear media using Nédélec elements.

Findings

A priori error estimates are established for the finite element approximation.

The method is applicable to Maxwell's equations in Kerr-type nonlinear media.

The analysis confirms the effectiveness of Nédélec elements in this context.

Abstract

In this paper the semi-discrete finite element approximation of initial boundary value problems for Maxwell's equations in nonliear media of Kerr-type is investigated. For the case of N\'ed\'elec elements from the first family, a priori error estimates are established for the approximation.

Equations

$${\partial }_t\mathbf{D}-\nabla \times \mathbf{H}$$

$${\partial }_t\mathbf{B}+\nabla \times \mathbf{E}$$

$$\mathbf{B}$$

$$\mathbf{P}(\mathbf{E})$$

$$\mathbf{D}={\epsilon }_0{\epsilon }_s\mathbf{E}\text{with}{\epsilon }_s={\epsilon }_s(\mathbf{E}):=1+{\chi }^{(1)}+{\chi }^{(3)}|\mathbf{E}{|}^2,$$

$${\partial }_t\mathbf{D}=2{\epsilon }_0{\chi }^{(3)}(\mathbf{E}\cdot {\partial }_t\mathbf{E})\mathbf{E}+{\epsilon }_0{\epsilon }_s{\partial }_t\mathbf{E}={\epsilon }_0\left({\epsilon }_s\mathbf{I}+{\epsilon }_m\right){\partial }_t\mathbf{E}\text{with}{\epsilon }_m={\epsilon }_m(\mathbf{E}):=2{\chi }^{(3)}\mathbf{E}{\mathbf{E}}^{\top },$$

$$\epsilon (\mathbf{E}):={\epsilon }_0\left({\epsilon }_s(\mathbf{E})\mathbf{I}+{\epsilon }_m(\mathbf{E})\right),$$

$$\epsilon (\mathbf{E}){\partial }_t\mathbf{E}-\nabla \times \mathbf{H}=0\text{ in }{Q}_T,$$

$${\mu }_0{\partial }_t\mathbf{H}+\nabla \times \mathbf{E}=0\text{ in }{Q}_T.$$

$$\varepsilon_{0}^{-1}{\mathbf{\Phi}}^{\top}\varepsilon({\mathbf{\Psi}}){\mathbf{\Phi}}=\Big{(}1+\chi^{(1)}+\chi^{(3)}|{\mathbf{\Psi}}|^{2}\Big{)}|{\mathbf{\Phi}}|^{2}+2\chi^{(3)}|{\mathbf{\Psi}}\cdot{\mathbf{\Phi}}|^{2}\geq|{\mathbf{\Phi}}|^{2}.$$

$${\mathbf{C}}({\mathbf{E}}):=\big{(}\varepsilon({\mathbf{E}})\big{)}^{-1}.$$

$$\mathbf{C}(\mathbf{E})=\frac{1}{{\epsilon }_0}{\mathbf{C}}_m(\mathbf{E})$$

$${\mathbf{C}}_m(\mathbf{E})$$

$$=\frac{1}{{\epsilon }_s(\mathbf{E})}\left(\mathbf{I}-\frac{1}{1+{\chi }^{(1)}+3{\chi }^{(3)}|\mathbf{E}{|}^2}{\epsilon }_m(\mathbf{E})\right).$$

$${\epsilon }_0{\partial }_t\mathbf{E}-{\mathbf{C}}_m(\mathbf{E})\nabla \times \mathbf{H}=0\text{ in }{Q}_T,$$

$${\mu }_0{\partial }_t\mathbf{H}+\nabla \times \mathbf{E}=0\text{ in }{Q}_T.$$

$$\nu \times \mathbf{E}=0\text{on }(0,T)\times \partial \Omega ,$$

$$\mathbf{E}(0,\mathbf{x})={\mathbf{E}}_0(\mathbf{x})\text{and }\mathbf{H}(0,\mathbf{x})={\mathbf{H}}_0(\mathbf{x})\text{for all }\mathbf{x}\in \Omega ,$$

$$\nabla \cdot ({\mu }_0{\mathbf{H}}_0)=0\text{in }\Omega ,{\mathbf{H}}_0\cdot \nu =0\text{on }\partial \Omega .$$

$$\nabla \cdot ({\mu }_0\mathbf{H})$$

$$(\mathbf{u},\mathbf{v}{)}_{\omega }:={\int }_{\Omega }\omega \mathbf{u}\cdot \mathbf{v}d\mathbf{x}\text{and}\Vert \mathbf{u}{\Vert }_{\omega }:=\Vert \mathbf{u}{\Vert }_{{\mathbf{L}}_{\omega }^2(\Omega )}:=\sqrt{(\mathbf{u},\mathbf{u}{)}_{\omega }}.$$

$$\Vert \mathbf{u}{\Vert }_{\omega }<\infty .$$

$$\Vert \mathbf{u}{\Vert }_{{\mathbf{L}}_{{\omega }_1}^2(\Omega )}\le \Vert \mathbf{u}{\Vert }_{{\mathbf{L}}_{{\omega }_2}^2(\Omega )}\text{for all }\mathbf{u}\in {\mathbf{L}}_{{\omega }_1}^2(\Omega )\cap {\mathbf{L}}_{{\omega }_2}^2(\Omega ).$$

$$\Vert \mathbf{u}{\Vert }_{{\mathbf{C}}^m(0,T;X)}:=\sum _{j=0}^m\underset{t\in (0,T)}{\sup }\Vert {\mathbf{u}}^{(j)}(t){\Vert }_X.$$

$${\int }_0^T\Vert \mathbf{u}(t){\Vert }_X^{p}dt<\infty$$

$$\Vert \mathbf{u}{\Vert }_{{\mathbf{L}}_p(0,T;X)}:={\left\{{\int }_0^T\Vert \mathbf{u}(t){\Vert }_X^{p}dt\right\}}^{1/p}.$$

$$\Vert \mathbf{u}{\Vert }_{{\mathbf{L}}_2(0,T;{\mathbf{L}}_{\omega }^2(\Omega ))}:={\left\{{\int }_0^T{\int }_{\Omega }|\mathbf{u}(t){|}^2\omega d\mathbf{x}dt\right\}}^{1/2}.$$

$$\mathbf{H}(\text{curl},\Omega )$$

$${\mathbf{H}}_0(\text{curl},\Omega )$$

$$\mathbf{H}(\text{div},\Omega )$$

$$\Vert \mathbf{u}{\Vert }_{\mathbf{H}(\text{curl},\Omega )}$$

$$\Vert \mathbf{u}{\Vert }_{\mathbf{H}(\text{div},\Omega )}$$

$$({\mathbf{E}},{\mathbf{H}})\in\big{(}{\mathbf{C}}^{1}(0,T;{\mathbf{L}}_{2,\varepsilon({\mathbf{E}})}(\Omega))\cap{\mathbf{C}}(0,T;{\mathbf{H}}_{0}(\textnormal{curl},\Omega))\big{)}\times\big{(}{\mathbf{C}}^{1}(0,T;{\mathbf{L}}_{2}(\Omega))\cap{\mathbf{C}}(0,T;{\mathbf{H}}(\textnormal{curl},\Omega))\big{)}$$

$$(\epsilon (\mathbf{E}){\partial }_t\mathbf{E},\mathbf{\Psi })-(\nabla \times \mathbf{H},\mathbf{\Psi })$$

Full text

Keywords: Semi-discrete finite element method, nonlinear Maxwell’s equations, error estimate

MSC 2010: 35Q61, 65M60, 65M15

1 Introduction

In this paper we investigate the semi-discrete conforming finite element approximation to the solution of Maxwell’s equations for nonlinear media of Kerr-type. As a concrete example, we consider the (meanwhile classical) Nédélec elements from the so-called first family. To the best knowlege of the author, the nonlinear situation is not yet well investigated, most works dealing with nonlinear effects are computational or experimental (see, e.g., [AY18]). Here we derive energy (stability) estimates for the weakly formulated problem and error estimates for the semi-discretized problem.

Let $Q_{T}:=(0,T)\times\Omega$, where $\Omega\subset\mathbb{R}^{3}$ is a simply connected domain with a sufficiently smooth boundary $\partial\Omega$ and $T>0$ is the length of the time interval under consideration. Let ${\mathbf{D}},{\mathbf{B}},{\mathbf{E}},{\mathbf{H}}:\;Q_{T}\to\mathbb{R}^{3}$ represent the displacement field, the magnetic induction, the electric and magnetic field intensities, respectively. The time-dependent Maxwell’s equations in a nonlinear medium can be written in the form

[TABLE]

where the following constitutive relations hold:

[TABLE]

Here $\varepsilon_{0}>0$ and $\mu_{0}>0$ are the vacuum permittivity and the permeability, respectively. Often the constitutive relation for the polarization ${\mathbf{P}}={\mathbf{P}}({\mathbf{E}})$ is approximated by a truncated Taylor series [Boy03]. In the case of an isotropic material, it takes the form

[TABLE]

where $\chi^{(j)}:\;\Omega\to\mathbb{R}$ are the media susceptibility coefficients, $j=1,3$. Then from (3) we obtain the representation

[TABLE]

and it follows by a simple calculation that

[TABLE]

where ${\mathbf{I}}$ denotes the identy in $\mathbb{R}^{d}$. Setting

[TABLE]

the system (1)–(3) can be written as

[TABLE]

Next we state a simple result which in particular implies that the matrix $\varepsilon({\mathbf{E}})$ is regular for all electric field intensities ${\mathbf{E}}$ under consideration.

Lemma 1.1**.**

Let $\chi^{(1)},\chi^{(3)}\geq 0$ a.e. in $\Omega$. Then the symmetric matrix $\varepsilon({\mathbf{\Psi}})$ is uniformly positive definite a.e. for any ${\mathbf{\Psi}}\in\mathbb{R}^{3}$.

Proof.

For all ${\mathbf{\Phi}}\in\mathbb{R}^{3}$, it holds that

[TABLE]

∎

As in [PNTB09], we denote the inverse by

[TABLE]

By means of the Sherman-Morrison formula (see, e.g., [GvL96]), the matrix ${\mathbf{C}}({\mathbf{E}})$ can be given explicitely:

[TABLE]

with

[TABLE]

Therefore, if the formula (6) holds, the system (4)–(5) takes the form

[TABLE]

This is an appropriate formulation for the development of time-discrete numerical algorithms, see, e.g., [PNTB09].

A perfect conducting boundary condition on $\Omega$ is assumed so that

[TABLE]

where $\mathbf{\nu}$ denotes the outward unit normal on $\partial\Omega$. In addition, initial conditions have to be specified so that

[TABLE]

where ${\mathbf{E}}_{0}:\;\Omega\to\mathbb{R}^{3}$ and ${\mathbf{H}}_{0}:\;\Omega\to\mathbb{R}^{3}$ are given functions, and ${\mathbf{H}}_{0}$ satisfies

[TABLE]

The divergence-free condition in (11) together with (5) implies that

[TABLE]

2 Notation

For a real number $p\in[1,\infty]$, the symbol $L^{p}(\Omega)$ denotes the usual Lebesgue spaces equipped with the norm $\|\cdot\|_{L^{p}(\Omega)}$. The analogous spaces of vector fields $\mathbf{u}:\Omega\to\mathbb{R}^{3}$ are denoted by ${\mathbf{L}}_{p}(\Omega):=\left[L^{p}(\Omega)\right]^{3}$ with the norm $\|\cdot\|_{{\mathbf{L}}^{p}(\Omega)}$.

In what follows we have to deal with weighted function spaces. Given a weight ${\mathbf{\omega}}:\;\Omega\to\mathbb{R}$, where the values of ${\mathbf{\omega}}$ are positive a.e. on $\Omega,$ we define a weighted inner product and a weighted norm by

[TABLE]

The space ${\mathbf{L}}_{{\mathbf{\omega}}}^{2}(\Omega)$ consists of vector fields $\mathbf{u}:\Omega\to\mathbb{R}^{3}$ with Lebesgue-measurable components and such that

[TABLE]

In the case ${\mathbf{\omega}}=1,$ the subscript is omitted. An elementary property of weighted spaces, which we will apply at different places without special emphasis, is the monotonicity w.r.t. the weight: If ${\mathbf{\omega}}_{1}$, ${\mathbf{\omega}}_{2}$ are two weights such that ${\mathbf{\omega}}_{1}\leq{\mathbf{\omega}}_{2}$ a.e. on $\Omega,$ then

[TABLE]

As transient problems are addressed, we will work with functions that depend on time and have values in certain Banach spaces. If $\mathbf{u}=\mathbf{u}(t,\mathbf{x})$ is a vector field of the space variable $\mathbf{x}$ and the time variable $t$, it is suitable to separate these variables in such a way that $\mathbf{u}(t)=\mathbf{u}(t,\cdot)$ is considered as a function of $t$ with values in a Banach space, say $X$, with the norm $\|\cdot\|_{X}$. That is, for any $t\in(0,T)$, the mapping $\mathbf{x}\mapsto\mathbf{u}(t,\mathbf{x})$ is interpreted as a parameter-dependent element $\mathbf{u}(t)$ of $X$. In this sense we will write ${\mathbf{E}}(t)={\mathbf{E}}(t,\cdot)$, ${\mathbf{H}}(t)={\mathbf{H}}(t,\cdot)$ and so on.

The space ${\mathbf{C}}^{m}(0,T;X)$, $m\in\mathbb{N}\cup\{0\}$, consists of all continuous functions $\mathbf{u}:\;(0,T)\to X$ that have continuous derivatives up to order $m$ on $(0,T).$ It is equipped with the norm

[TABLE]

For the sake of consistency in the notation we will write ${\mathbf{C}}(0,T;X):={\mathbf{C}}^{0}(0,T;X)$.

The space ${\mathbf{L}}_{p}(0,T;X)$ with $p\in[1,\infty)$ contains (equivalent classes of) strongly measurable functions $\mathbf{u}:\;(0,T)\to X$ such that

[TABLE]

(for the definition of strongly measurable functions we refer to [KJF77]). The norm on ${\mathbf{L}}_{p}(0,T;X)$ is defined by

[TABLE]

These spaces can be equipped with a weight, too. In particular, we will write

[TABLE]

Finally, all the above definitions can be extended to the standard Sobolov spaces of functions with weak spatial derivatives of maximal order $r\in\mathbb{N}$ in $L^{p}(\Omega)$: $W^{r,p}(\Omega)$ with norm $\|\cdot\|_{W^{r,p}(\Omega)}$. If $p=2$, we write $H^{r}(\Omega):=W^{r,2}(\Omega)$ and $\|\cdot\|_{H^{r}(\Omega)}:=\|\cdot\|_{W^{r,2}(\Omega)}$.

The space $H_{0}^{1}(\Omega)$ is defined as the closure of $C_{0}^{\infty}(\Omega)$ with respect to the norm $\|\cdot\|_{H^{1}(\Omega)}$, where $C_{0}^{\infty}(\Omega)$ denotes the space of all arbitrarily often differentiable functions with compact support on $\Omega$. It is well knwon that $H_{0}^{1}(\Omega)$ is a closed subspace of $H^{1}(\Omega)$ and consists of elements $u$ such that $u=0$ on $\partial\Omega$ in the sense of traces [AF03]. As in the case of the $L^{p}$-spaces, we shall write ${\mathbf{W}}^{r,p}(\Omega):=[W^{r,p}(\Omega)]^{3}$ and so on.

Furthermore, we need the following Hilbert spaces that are related to the (weak) rotation and divergence operators:

[TABLE]

These spaces are equipped with the norms (resp. induced norms)

[TABLE]

We refer to [DL76], [RT77], [GR86] and [Cia02] for details about these spaces.

3 Weak formulations

We assume that a solution

[TABLE]

of the nonlinear Maxwell’s equations (4)–(5) exists and is unique.

We multiply equation (4) by a test function ${\mathbf{\Psi}}\in{\mathbf{L}}_{2}(\Omega)$ and integrate over $\Omega$. Similarly we multiply (5) by a test function ${\mathbf{\Phi}}\in{\mathbf{H}}(\textnormal{curl},\Omega)$, integrate the result over $\Omega$ and integrate by parts the second term. This shows that it is natural to look for a weak solution ({\mathbf{E}},{\mathbf{H}})\in\big{(}{\mathbf{C}}^{1}(0,T;{\mathbf{L}}_{2,\varepsilon({\mathbf{E}})}(\Omega))\cap{\mathbf{C}}(0,T;{\mathbf{L}}_{2}(\Omega))\big{)}\times\big{(}{\mathbf{C}}^{1}(0,T;{\mathbf{L}}_{2}(\Omega))\cap{\mathbf{C}}(0,T;{\mathbf{H}}(\textnormal{curl},\Omega))\big{)} of (4)–(5) such that

[TABLE]

Alternatively, the use of test functions ${\mathbf{\Psi}}\in{\mathbf{H}}_{0}(\textnormal{curl},\Omega)$ and ${\mathbf{\Phi}}\in{\mathbf{L}}_{2}(\Omega)$ and the integration by parts in the equation (4) leads to the notion of a weak solution ({\mathbf{E}},{\mathbf{H}})\in\big{(}{\mathbf{C}}^{1}(0,T;{\mathbf{L}}_{2,\varepsilon({\mathbf{E}})}(\Omega))\cap{\mathbf{C}}(0,T;{\mathbf{H}}_{0}(\textnormal{curl},\Omega))\big{)}\times\big{(}{\mathbf{C}}^{1}(0,T;{\mathbf{L}}_{2}(\Omega))\big{)} of (4)–(5) such that

[TABLE]

In both cases, the initial conditions (10) have to be satisfied at least in the sense of

${\mathbf{C}}(0,T;{\mathbf{L}}_{2}(\Omega))$.

Remark 1*.*

As a consequence of the embedding (as sets)

[TABLE]

and of the fact that $C_{0}^{\infty}(\Omega)$ is dense in $L^{2}(\Omega)$ we see that ${\mathbf{H}}(\textnormal{div},\Omega)$ is a dense subset of ${\mathbf{L}}_{2}(\Omega)$ [AF03]. Therefore the test space in (17) can be reduced to ${\mathbf{H}}(\textnormal{div},\Omega)$. Also note that ${\mathbf{H}}\in{\mathbf{C}}(0,T;{\mathbf{H}}(\textnormal{div},\Omega))$ due to (12).

Remark 2*.*

In the case where $\mu_{0}$ is not a constant but a highly variable function $\mu=\mu(\mathbf{x})$ it is more convenient to use the magnetic flux density ${\mathbf{B}}=\mu{\mathbf{H}}$ instead of ${\mathbf{H}}$ as a dependent variable [MM95]. In such a case, the formulation (16)–(17) is replaced by

[TABLE]

where

[TABLE]

Next we will formulate a stability result for the problem (14)–(15). For this, we extend this problem to the case of a nontrivial right-hand side for a moment.

Theorem 3.1**.**

Let ${\mathbf{J}}_{e}\in{\mathbf{L}}_{2}(0,T;{\mathbf{L}}_{2,\varepsilon_{0}^{-1}(1+\chi^{(1)})^{-1}}(\Omega))$, ${\mathbf{J}}_{m}\in{\mathbf{L}}_{2}(0,T;{\mathbf{L}}_{2,\mu_{0}^{-1}}(\Omega))$ – the electric and magnetic current densities, respectively – be given and assume that the system

[TABLE]

together with the initial conditions (10) has a weak solution

[TABLE]

If

[TABLE]

denotes the nonlinear electromagnetic energy at the time $t$, the following energy law in differential form holds:

[TABLE]

Proof.

Taking ${\mathbf{\Psi}}={\mathbf{E}}$ and ${\mathbf{\Phi}}={\mathbf{H}}$ in (20)–(21) and adding the result gives

[TABLE]

An elementary calculation shows that

[TABLE]

hence

[TABLE]

Then it follows from equation (23) that

[TABLE]

∎

Corollary 3.2**.**

Under the assumptions of Thm. 3.1, for all $t\in(0,T)$, the following estimate is valid:

[TABLE]

Proof.

We estimate the right-hand side of (22) by means of the Cauchy-Bunyakovski-Schwarz’ inequality (twice – in the integral version and in the finite sum version):

[TABLE]

Since $\frac{d}{dt}\sqrt{W}=\frac{1}{2\sqrt{W}}\frac{dW}{dt}$ for $W(t)>0$ (the case $W(t)=0$ is trivial), it follows from (22) that

[TABLE]

Integrating this inequality w.r.t. $t$, we get

[TABLE]

Then

[TABLE]

∎

Remark 3*.*

Analogous results as in Thm. 3.1 and Cor. 3.2 can be obtained for the corresponding “non-homogeneous” version of (16)–(17) and for the subsequent semi-discretizations.

4 Spatial discretization

4.1 Semi-discretization of the weak formulations

Let ${\mathbf{W}}_{h}\subset{\mathbf{L}}_{2}(\Omega)$, ${\mathbf{U}}_{h}\subset{\mathbf{H}}(\textnormal{curl},\Omega)$, ${\mathbf{U}}_{0h}\subset{\mathbf{H}}_{0}(\textnormal{curl},\Omega)$, and ${\mathbf{V}}_{h}\subset{\mathbf{H}}(\textnormal{div},\Omega)$ be finite-dimensional subspaces.

The semi-discrete (in space) problem for the system (14)–(15) consists in determining elements $({\mathbf{E}}_{h},{\mathbf{H}}_{h})\in{\mathbf{C}}^{1}(0,T;{\mathbf{W}}_{h})\times{\mathbf{C}}^{1}(0,T;{\mathbf{U}}_{h})$ such that

[TABLE]

For the the equations (16)–(17), the semi-discrete problem involves the determination of elements $({\mathbf{E}}_{h},{\mathbf{H}}_{h})\in{\mathbf{C}}^{1}(0,T;{\mathbf{U}}_{0h})\times{\mathbf{C}}^{1}(0,T;{\mathbf{V}}_{h})$ such that

[TABLE]

The initial conditions for both problems read formally as

[TABLE]

where the concrete requirements to the particular choice of the discrete initial data $({\mathbf{E}}_{0h},{\mathbf{H}}_{0h})$ will be seen later (Thm. 5.1).

4.2 The choice of the finite element spaces

In the rest of the paper we will choose the so-called first family of Nédélec edge elements, usually denoted by ${\mathbf{D}}_{k}$ and ${\mathbf{R}}_{k}$ for $k\in\mathbb{N}$, for the construction of the concrete finite element spaces (for details see [Néd80] or [Mon03, Ch. 5]). That is, given an arbitrary member ${\mathcal{T}_{h}}$ of a family of triangulations of $\Omega$ consisting of open tetrahedra $K$, we set

[TABLE]

where ${\mathbf{P}}_{\ell}:=[\mathcal{P}_{\ell}]^{3}$ and $\mathcal{P}_{\ell}$ is the space of scalar real-valued polynomials in three variables of maximal degree $\ell\in\mathbb{N}\cup\{0\}$. To deal with the case of the Nédélec formulation (16)–(17), we still have to introduce the space ${\mathbf{U}}_{0h}:={\mathbf{U}}_{h}\cap{\mathbf{H}}_{0}(\textnormal{curl},\Omega)$. In the subsequent error analysis, we will make use of some projection operators. For the Lee-Madsen formulation (24)–(25), we need projections ${\mathbf{P}}_{h}:\;{\mathbf{L}}_{2}(\Omega)\to{\mathbf{W}}_{h}$, $\Pi_{h}:\;{\mathbf{H}}(\textnormal{curl},\Omega)\to{\mathbf{U}}_{h}$.

Let ${\mathbf{P}}_{h}$ be the standard ${\mathbf{L}}_{2}(\Omega)$-projection operator onto ${\mathbf{W}}_{h}$, i.e. for given $\mathbf{w}\in{\mathbf{L}}_{2}(\Omega)$ the image ${\mathbf{P}}_{h}\mathbf{w}\in{\mathbf{W}}_{h}$ is defined by

[TABLE]

For this operator the following standard error estimate holds: If $\mathbf{w}\in{\mathbf{H}}^{k}(\Omega)$, then

[TABLE]

Moreover, since $\nabla\times{\mathbf{U}}_{h}\subset{\mathbf{W}}_{h}$ (see the beginnings of the proofs of [Mon91, Thm. 3.3] or [Mon03, Lemma 5.40]), it holds that

[TABLE]

Next, for $\mathbf{v}\in{\mathbf{H}}(\textnormal{curl},\Omega)$ we define ${\mathbf{\Pi}}_{h}\mathbf{v}\in{\mathbf{U}}_{h}$ by

[TABLE]

where $S_{h}^{k}$ is defined as

[TABLE]

see [Mon91, Subsect. 4.2].

If $\mathbf{v}\in{\mathbf{H}}^{k+1}(\Omega)$ such that $\nabla\cdot\mathbf{v}=0$ in $\Omega$ and $\mathbf{\nu}\cdot\mathbf{v}=0$ on $\partial\Omega$, then there exists a constant $C>0$ such that

[TABLE]

(see [Mon91, Thm. 4.6]).

5 An error estimate for the semi-discrete problem

In this section we formulate and prove the main result.

Theorem 5.1** (Semi-discrete error estimate for the Lee-Madsen formulation).**

Let $k\in\mathbb{N}$, $\chi^{(1)},\chi^{(3)}\in L_{\infty}(\Omega)$, ${\mathbf{E}}_{0}\in{\mathbf{L}}_{\infty}(\Omega),$ ${\mathbf{H}}_{0}\in{\mathbf{L}}_{2,\mu_{0}}(\Omega)$ satisfying (11),

[TABLE]

be the weak solution of the system (14)–(15), and

[TABLE]

be the finite element solution of the system (24)–(25) respectively, where the inclusion is to be understood uniformly w.r.t. the mesh parameter $h$ in the sense that $\|{\mathbf{E}}_{h}\|_{{\mathbf{C}}^{1}(0,T;{\mathbf{L}}_{\infty}(\Omega))}$ is bounded by a constant independent of $h$. Then there exists a constant $C>0$ independent of $h$ such that the following error estimate holds:

[TABLE]

(the detailed structure of the bound is given at the end of the proof).

Proof.

We set ${\mathbf{\Psi}}:={\mathbf{\Psi}}_{h}\in{\mathbf{W}}_{h}$ in (14) and ${\mathbf{\Phi}}:={\mathbf{\Phi}}_{h}\in{\mathbf{U}}_{h}$ in (15):

[TABLE]

By means of the projection operators ${\mathbf{P}}_{h}$ and ${\mathbf{\Pi}}_{h}$ defined in (28) and (31)–(32), resp., from this we get

[TABLE]

The last term on the right-hand side of (34) vanishes thanks to the properties of ${\mathbf{\Pi}}_{h}$, see [Mon91, eq. (2.4)] and (31).

The second term on the right-hand side of (35) can be omitted because of the commutation property $\partial_{t}{\mathbf{\Pi}}_{h}{\mathbf{H}}={\mathbf{\Pi}}_{h}\partial_{t}{\mathbf{H}}$, which results from the continuity properties of the operator ${\mathbf{\Pi}}_{h}$. The last term on the right-hand side vanishes thanks to the property (30) of ${\mathbf{P}}_{h}$.

Therefore (34)–(35) simplify to

[TABLE]

Now, subtracting (36)–(37) from the system (24)–(25) and taking into consideration that $\mu_{0}$ is constant, we obtain:

[TABLE]

Now we will deal with the first term of (38), where we have in mind the choice ${\mathbf{\Psi}}_{h}={\mathbf{E}}_{h}-{\mathbf{P}}_{h}{\mathbf{E}}$ in what follows:

[TABLE]

The treatment of $\delta_{1}$ is quite obvious. With ${\mathbf{E}}_{h}-{\mathbf{E}}={\mathbf{\Psi}}_{h}+{\mathbf{P}}_{h}{\mathbf{E}}-{\mathbf{E}}$ we get

[TABLE]

The term $\delta_{2}$ is decomposed as follows:

[TABLE]

For $\delta_{3}$, we use the following decomposition:

[TABLE]

With these decompositions, equation (38) takes the form

[TABLE]

or, after some rearrangement,

[TABLE]

Then:

[TABLE]

Since

[TABLE]

it follows that

[TABLE]

From the estimates

[TABLE]

and, analogously,

[TABLE]

we conclude that

[TABLE]

For the right-hand side, we have:

[TABLE]

where the positive constants $C_{1}$, $C_{2}$, $C_{3}$ depend on certain norms of $\chi^{(1)}$, $\chi^{(3)}$, ${\mathbf{E}}$, and ${\mathbf{E}}_{h}$. Combining the estimates (41) and (42) with (40), we get

[TABLE]

This finally leads to

[TABLE]

where

[TABLE]

Now we consider (39) with ${\mathbf{\Phi}}_{h}={\mathbf{H}}_{h}-{\mathbf{\Pi}}_{h}{\mathbf{H}}$ and get

[TABLE]

Adding both inequalities and making use of the commutation property of ${\mathbf{P}}_{h}$, we arrive at

[TABLE]

The projection errors can be estimated by means of (29) and (33), that is, for ${\mathbf{E}},\partial_{t}{\mathbf{E}}\in{\mathbf{H}}^{k}(\Omega)$ and $\partial_{t}{\mathbf{H}}\in{\mathbf{H}}^{k+1}(\Omega)$, we have that

[TABLE]

In this way the above estimate can be written as

[TABLE]

Setting

[TABLE]

we get

[TABLE]

Integrating this inequality, we obtain

[TABLE]

By the monotonicity of the weighted norms w.r.t. the weight and the nonnegativity of the integral term on the left-hand side, we see that

[TABLE]

On the other hand, we have the estimates

[TABLE]

and

[TABLE]

Combining (44), (45), (46) with (43), we get

[TABLE]

or, equivalently,

[TABLE]

where $C_{6}^{2}:=\|1+\chi^{(1)}\|_{L_{\infty}(\Omega)}+3\|\chi^{(3)}\|_{L_{\infty}(\Omega)}\|{\mathbf{E}}(0)\|_{L_{\infty}(\Omega)}^{2}$.

In the paper [Daf79], a Gronwall-type lemma (Lemma 4.1) is specified which gives a bound on the value $w(T)$ provided an inequality like (47) is satisfied:

[TABLE]

From this and the triangle inequality in conjunction with (29) and (33) the statement follows. Indeed, since

[TABLE]

we get

[TABLE]

Then

[TABLE]

which implies the stated estimate. ∎

Remark 4*.*

Note that the constant $C$ in this estimate behaves as $Te^{C_{4}T}$ for large $T$.

6 Conclusion

In this paper we have investigated a semi-discrete conforming finite element approximation to the solution of Maxwell’s equations for nonlinear media of Kerr-type using Nédélec elements from the first family. We have demonstrated energy (stability) estimates for the weakly formulated problem and error estimates for the semi-discretized problem. The results can be extended to other conforming finite element methods provided the corresponding projection operators ${\mathbf{P}}_{h}$ and ${\mathbf{\Pi}}_{h}$ ((28), (31)–(32)) admit analogous properties.