Error analysis of a fully discrete discontinuous Galerkin alternating direction implicit discretization of a class of linear wave-type problems

This paper is concerned with the rigorous error analysis of a fully discrete scheme obtained by using a central fluxes discontinuous Galerkin (dG) method in space and the Peaceman–Rachford splitting scheme in time. We apply the scheme to a general class of wave-type problems and show that the resulting approximations as well as discrete derivatives thereof satisfy error bounds of the order of the polynomial degree used in the dG discretization and order two in time. In particular, the class of problems considered includes, e.g., the advection equation, the acoustic wave equation, and the Maxwell equations for which a very efficient implementation is possible via an alternating direction implicit splitting.


Introduction
In this paper we consider the rigorous error analysis of a fully discrete alternating direction implicit (ADI) scheme for a general class of wave-type problems. In particular, we analyze the full discretization obtained by using a method-of-lines approach comprised of the Peaceman-Rachford scheme in time [37] and a central fluxes discontinuous Galerkin (dG) discretization in space [5,20] in a rather general framework. Given an initial state u 0 we seek the solution u of the wave-type problem M∂ t u = Lu + g, R + × , on an open, bounded and connected Lipschitz domain in ⊂ R d , d = 2, 3, . . ., supplemented with suitable boundary conditions. The latter will be incorporated into the domain of the spatial differential operator L. In (1), M ∈ L ∞ ( ) m×m is a symmetric and uniformly positive definite material tensor and g ∈ C(R + ; L 2 ( ) m ) is a source term. The operator L is a Friedrichs' operator [13] given by where the partial derivatives are to be understood in a distributional sense. The assumption that the coefficients L 0 , . . . , L d are constant is posed for the sake of presentation. However, we point out that all results hold true for more general coefficients fulfilling suitable regularity assumption, cf., [29] for details. To obtain a wellposed (and energy-dissipating) problem in this constant coefficient case, it is sufficient to assume that the symmetric part of L 0 is negative semi-definite, i.e., x T (L 0 + L 0 T )x ≤ 0 for all x ∈ R m , and that the remaining coefficients L 1 , . . . , L d are symmetric [5,23]. Note that it is possible to drop the condition on L 0 and still get a wellposed, albeit only shift-dissipative, problem. The class of problems described by (1) comprises, among others, the advection equation, the acoustic or elastic wave equations, and the Maxwell equations, cf., [2,3,5,29].
In this paper, we employ the framework of (abstract) evolution equations. This allows us to study the wellposedness of (1) within the well-understood and extensive theory of linear semigroups [36]. Moreover, one can then formally apply time integration schemes designed for ordinary differential equations, such as Runge-Kutta schemes or the Peaceman-Rachford scheme used here, in a straightforward way. We thus eliminate the material tensor M in front of the time derivative to obtain (2b) with f = M −1 g and Lu = M −1 Lu. By this reformulation, L inherits the structure of L and we thus refer to this operator as a Friedrichs' operator as well.
We analyze a fully discrete Peaceman-Rachford ADI scheme for the numerical solution of (2). ADI schemes are a class of splitting schemes first proposed for the heat equation discretized by finite differences in [37]. The main idea is to split the spatial differential operator L in terms of the spatial directions. Combining this splitting with a suitable time integration method, such as the Peaceman-Rachford scheme, results in an unconditionally stable integrator. The main feat of the resulting scheme is that each timestep can be performed at about the cost of an explicit integrator, which makes it extremely efficient. Roughly speaking, this efficiency is achieved by splitting in such a way that only partial derivatives in one spatial direction occur in each of the split operators, essentially decoupling the one-dimensional flows. Of course, this severely restricts the class of problems to which the method can be applied efficiently. One main restriction is that the problem has to be considered on cuboidal domains (or unions thereof). Further, up until recently, the method was basically restricted to finite difference discretizations of the heat equation or the acoustic wave equation, both in two dimensions, where the appropriate splitting can easily be deduced [37]. There are varieties, where more than two split operators are considered to tackle higher dimensions, see e.g., [15,27]. However, this results in a scheme that is only conditionally stable, hence losing a very favorable property of ADI schemes.
Rather surprisingly, in [33,45], a novel ADI splitting into two operators for the full 3D Maxwell equations was proposed for finite differences time domain (FDTD) discretizations on the Yee grid [43]. Consequently, this ADI-FDTD method received a lot of attention and gave rise to several variants. For instance, more efficient implementations can be found in [4,14,16,25] and energy preserving variants in [7,17,18], see also the references given in these papers.
Motivated by this new splitting, in [22], the authors of this paper identified a general class of problems of the form (2) and their space discretizations which admit a splitting that leads to an efficient ADI method with only two split operators, while preserving unconditional stability. Since this class only comprised diagonal material tensors, it was extended to more general tensors in [29].
Further, in [41] and the recent review [42], the Peaceman-Rachford and other, related splitting schemes were rewritten into a formulation referred to as fundamental implicit scheme. This reformulation avoids applications of the discrete operators on the right-hand side leading to an even more efficient scheme. One can directly apply these ideas to the scheme considered here and in [22], since it exhibits the exact same structures exploited in [41,42].
Unfortunately, finite differences in space prevent a rigorous error analysis under realistic regularity assumptions on the approximated solution. Hence, in [22], we combined the Peaceman-Rachford scheme with a dG discretization in space, cf., [5,20]. We showed that if a tensor-structured grid is used, the fully discrete method exhibits runtime behavior similar to the finite difference version. In particular, one step can be performed in linear complexity w.r.t. the number of elements used in the mesh. Here, we supplement these computational investigations by a rigorous error analysis of the scheme used in [22].
There are already some results on the stability and convergence of the Peaceman-Rachford scheme applied to finite difference space discretizations, cf., e.g., [12,18,25,26,31,38]. However, these results rely on strong regularity of the exact solution. So far, we are only aware of one publication that provides rigorous error bounds under realistic regularity assumptions, namely [19]. However, these results seem not to be applicable to the setting we consider here, because they rely on the assumption that norms of some concatenations of the discrete operators are bounded independently of the discretization parameters. To our knowledge, this is not the case here.
Besides full discretization results, there exist several rigorous semidiscrete analyses for the Peaceman-Rachford scheme applied to the Maxwell equations interpreted as an abstract Cauchy problem. In [21], the homogeneous and in [6,8,9], the inhomogeneous and damped equations were considered. Further, in [44], a detailed analysis of the abstract Maxwell problem with jumping material parameters was performed, providing error bounds of order 3/2 for suitable data. We point out that all these results pose assumptions only on the data, as opposed to the analysis carried out here, which assumes regularity of the exact solution. However, the regularity analysis carried out in [6,8,9,21] can be used to defer these assumptions to the data in the particular case of Maxwell equations.
The analysis we present here is based on the fact that the Peaceman-Rachford scheme can be interpreted as a perturbation of the Crank-Nicolson scheme. Hence, we were able to extend techniques established in [24,40]. These techniques were originally developed for a locally implicit time integration scheme but also apply to the full discretization obtained by the Crank-Nicolson and the leapfrog scheme, respectively. In fact, our extension of these techniques can also be applied to other perturbations of the Crank-Nicolson scheme in a straightforward manner.
Further, based on a stronger stability result for abstract evolution equations given in [28], we identify suitable approximations to the first-order derivatives in space and time of the solution of (2). We show that, if the solution is sufficiently regular, these approximations fulfill error bounds of the same order as the approximations to the solution given directly by the scheme. We point out that these results can be generalized in a straightforward way to other schemes that can be analyzed by the techniques from [24,40]. As far as we are aware of, such results were not given before in the literature. Nevertheless, we would like to mention that related error bounds in a discrete H 1 -norm were achieved for finite difference methods under much stronger regularity assumptions in [16,17].
In [23], as a preliminary work to the analysis of the full discretization of (2), we analyzed the error of the spatial semi-discretization by a central fluxes dG method [5,20]. However, in contrast to [23], in this paper, we directly define and formulate all our results in a weighted L 2 -inner product incorporating the material tensor M instead of the standard one. This approach appears to be more natural for our analysis. Further, by stating the required polynomial approximation results in this weighted inner product, we could eliminate Assumption 4.2 in [23]. As a consequence, all results now apply to more general material parameters, e.g., piecewise Lipschitz or piecewise smooth instead of piecewise constant ones. Since most proofs in [23] directly generalize to this weighted setting, we do not present them in detail here again.
Lastly, note that our analysis is not restricted to the class of problems which allow for an efficient ADI splitting cf., [22] for details. In fact, these restrictions are only necessary for the efficiency of the method, not for the numerical analysis of the discretization scheme. Hence, our results also hold for more general problems and geometries of the domain, cf., Sect. 4.1.
The paper is organized as follows. After providing some notation, in Sect. 2, we present the analytical framework used for our analysis and show wellposedness of (2). Subsequently, in Sect. 3, we briefly present the spatial discretization employed in our method-of-lines approach. Section 4 starts by providing suitable splittings, which are then used to fully discretize (2) by using the Peaceman-Rachford method in time. In this section we also show unconditional stability of the scheme and identify approximations to the first-order derivatives in space and time of the solution of (2). Then, in Sect. 5, we present our main results, which show rigorous error bounds for both the approximations directly provided by the scheme as well as the aforementioned approximations to the derivatives. We conclude with numerical experiments verifying our theoretical results in Sect. 6.

Notation
In this section, we introduce the basic notation used throughout the paper.
The distributional partial derivative in ith coordinate direction of R d is denoted by We further denote the 1 -norm of α ∈ N d 0 by |α|. For vectors and matrices these derivatives act componentwise. Given two real Hilbert spaces X , · · X and Y , · · Y , we denote the set of all bounded operators from X to Y by B(X , Y ). The identity operator on a Hilbert space is denoted by I. We write the dual space of X as X and denote the canonical dual pairing between X and its dual by · · : X × X → R.
In Sect. 2 we will restrict a Friedrichs' operator A to its domain D(A), which incorporates the boundary conditions in such a way that the restriction becomes maximal dissipative. Then, we define the domain of the concatenation of two Friedrichs' operators A and B as D(AB) := v ∈ D(B) | Bv ∈ D(A) and extend this definition recursively to more than two operators.
For the remainder of this section, let K ⊂ R d be open and F ⊂ ∂ K . The set of polynomials of degree at most k ∈ N on K is denoted by Q k d (K ). Given two vector-valued functions u, v ∈ L 2 (K ) m , we denote the (standard) L 2 (K )-inner product by and for F ⊂ ∂ K and u | F , v | F ∈ L 2 (F) m , we denote the surface integral over F as For the analysis, we also need the following inner product weighted by the material tensor M given by By the properties of M, the weighted and L 2 -inner product are equivalent and thus, the space L 2 (K ) m , · · K is a Hilbert space. The norms induced by these inner products are denoted by · L 2 ,K , · K and · F , respectively.
The L 2 -Sobolev spaces on K are denoted by H q (K ), q ∈ N 0 , and we equip them with the norms where α ∈ N d 0 is a multi-index. Note that we use the weighted L 2 -norm, not the standard one, to define the Sobolev-norms. However, by equivalence, these spaces can be identified. Further, note that q = 0 yields the corresponding (weighted) L 2 -space.
We denote the spectral norm of a matrix A ∈ R m×m by A . Given a square matrixvalued field A ∈ L ∞ (K ) m×m , the essential supremum of the spectral norm of A is denoted by Lastly, by W q,∞ (K ), q ∈ N 0 , we denote the L ∞ -Sobolev spaces on K and for A ∈ W q,∞ (K ) m×m , we write

Analytical setting
In this section, we briefly present some properties of Friedrichs' operators and systems of the form (2). We refer to [2,5,11,23] for most of the details.
The graph space of a Friedrichs' operator L is defined as As a straightforward consequence of [5,Lem. 7.2], the graph space H (L) together with the (weighted) graph norm form a Hilbert space and by definition we have L ∈ B(H (L), L 2 ( ) m ). The notation H (L) for the graph space is inspired by H (div) and H (curl), which are important, since they arise for the wave and Maxwell equations, respectively. In this weak setting, the meaning of boundary values is not clear. In order to still get hold of them, we need to introduce a more general concept of boundary values. For this, we first need to define the formal adjoint of L.
The formal adjoint L allows us to define the following boundary operator, which is the aforementioned generalization of boundary values.

Definition 2.2
We call L ∂ : H (L) → H (L) defined by the boundary operator associated with L.
Basically, (3) is a generalized integration-by-parts formula. Further, note that the boundary operator L ∂ is symmetric and bounded, i.e., L ∂ ∈ B(H (L), H (L) ), cf., [11,Lem. 2.2]. Boundary conditions in this framework are implemented by incorporating them into the domain of the operator L. To this end, we follow the procedure in [11,Sec. 2.1] and assume the existence of the following abstract boundary operator.

Assumption 2.3 There exists a bounded operator
Since the kernel of a bounded operator on a Hilbert space is closed, both ker(L ∂ −L ) and ker(L ∂ + L ) are Hilbert spaces if endowed with the graph norm of L.
To show wellposedness of (2), we use semigroup theory. Hence, we define the domain of L as D(L) := ker(L ∂ − L ). In fact, if restricted to this space, L is maximal dissipative [23,Thm. 3.5]. Note that D(L) can be seen as the subspace of H (L) in which the boundary conditions defined by L are incorporated.
As a direct consequence of the maximal dissipativity, the famous Lumer-Phillips Theorem [10, Thm. II.3.15, Cor. II. 3.20] implies that the restriction L | D(L) is the generator of a contraction semigroup, which we denote by e tL t≥0 . This immediately yields the following wellposedness and stability result, see [36] and [28,Lem. 2.4] for the second stability bound.
Then, for a given initial value u 0 ∈ D(L), there exists a unique solution u ∈ (2) given by the variation-of-constants formula Further, we have the stability bounds

Spatial discretization
We use the dG method with central fluxes to discretize the Friedrichs' operator L. The content of this section is a brief presentation of previously established results (mostly from [11,23,29]), hence we omit the rather technical and lengthy proofs.
To avoid taking domain approximation into account, in the following, we assume to be polyhedral. Given a mesh T h of , we denote the diameter of an (open) mesh element or cell K ∈ T h by h K . The index h = max K ∈T h h K denotes the maximal diameter of all elements in T h or meshsize of T h . To keep the notation of mesh-dependent norms concise and intuitive, we further define the piecewise constant Since we analyze the approximation error w.r.t. h, we consider a sequence of meshes T H = T h h∈H , where H is a countable collection of positive numbers h < 1 with 0 as only accumulation point. To investigate the aforementioned error, we have to ensure a certain quality of the meshes as the meshsize approaches zero and we therefore assume that T H is admissible in the sense of [5,Def. 1.57]. This means that the mesh sequence is shape and contact regular and has optimal polynomial approximation properties. We denote the regularity parameter of T H by ρ.
The faces of a mesh T h are collected in the set F h = F int h ∪ F bnd h , which is decomposed into the set of interior faces or interfaces F int h and the set of boundary faces F bnd h . Further, for each K ∈ T h we introduce the subset of faces that compose the boundary of said element as F K h and denote the maximum number of faces per element by Note that this number is bounded independently of h ∈ H due to the admissibility of T H , cf., [5,Lemma 1.41].
We denote the outward unit normal vector to an element K ∈ T h by n K . Additionally, we define a face normal vector to each face F ∈ F h denoted by n F in the following way. Given a boundary face F ∈ F bnd h , we simply define n F as the outward unit normal vector to . For each interface F ∈ F int h , we denote the two neighboring elements w.r.t. F arbitrarily by K F 1 and K F 2 and fix this choice. The face normal vector n F is then defined as the outward unit normal vector to K F 1 . Next, we denote the discrete approximation space consisting of broken polynomials of degree at most k in each variable by Throughout the rest of the paper, we denote objects contained in or mapping into the discrete space V h in bold face. In particular, we slightly deviate from the notation introduced in the beginning by denoting the identity on V h by I I I.

Remark 3.1
To keep the presentation concise, we consider polynomials of the same degree on each element K ∈ T h . Note, however, that the dG method can easily handle different polynomial degrees on different elements and all our results can be generalized to this case in a straightforward manner. Further, there is more freedom in choosing the discrete spaces, e.g., polynomials of total degree at most k. More details on this can be found in As a consequence of the admissibility of T H , some important properties of the discrete spaces can be inferred. In particular, since admissibility implies shape and contact regularity, the inverse inequality and the discrete trace inequality respectively, hold for all We point out that the result on the elements is originally stated in the standard L 2 -inner product, not the weighted version used here. However, by equivalence, it also holds in this setting. The inverse inequality (5) further implies a similar inequality for the Friedrichs' operator L instead of the gradient, namely where To identify a best approximation in V h to some L 2 -function (in some sense), we define the (weighted) L 2 -projection π h : The error caused by this projection is then denoted as Since the mesh sequence T H is admissible and hence possesses optimal polynomial approximation properties in the sense of [5, Def. 1.55], we are able to infer bounds that measure the quality of this approximation. However, since we use a weighted L 2projection and the necessary results [5, Lem. 1.58 & 1.59] are given for the standard one, we can not immediately apply them. Still, the next lemma shows that these results also hold in our setting.

Lemma 3.2 For all
where C π and C π,∂ are independent of K and h.

Proof
To prove this, one shows that, by the properties of M, the projection error π h w.r.t. the weighted inner product yields the same quality of approximation as the standard L 2 -projection w.r.t. · L 2 ,K . Using the equivalence of the weighted and standard norm then yields the claim. We omit the details.
We aim to approximate functions in space by elements of the discrete approximation space V h . These discrete functions can be discontinuous along element borders but are polynomials on the elements. Thus, they are apt to approximate functions that are sufficiently smooth on the elements. The broken Sobolev spaces defined by contain such functions. They are Hilbert spaces if endowed with the norm Since neither functions contained in the discrete approximation space V h , nor those contained in the broken Sobolev spaces H q (T h ) need to be continuous across the boundaries of mesh elements, their evaluation at these boundaries is not well-defined. Thus, for an interface F ∈ F int h , we write v | K F 1 and v | K F 2 for the limit of a function v approaching F from K F 1 or K F 2 , respectively. With this, we define the average and the jump of a function v across an interior face F ∈ F int h as respectively. These values serve as a measure for the discontinuities of the discrete functions and are used to establish coupling between the otherwise decoupled elements in the discrete problem. Lastly, analogously to the broken Sobolev spaces H q (T h ), we define the spaces

Friedrichs' operators in the discrete setting
To define the discrete operator, it is convenient to have a more concrete representation of boundary values than the abstract boundary operator L ∂ from Definition 2.2. To achieve this, we restrict ourselves to spaces whose elements admit square-integrable traces -in contrast to those of H (L). In particular, we consider the intersection of H (L) and the broken Sobolev space H 1 (T h ). Then, the abstract boundary operator L ∂ can be expressed by means of the usual integration-by-parts formula (on each element K ∈ T h ). Namely, by [23,Lem. 4.3] for Following the same line of reasoning, it is convenient to also have a representation of the abstract boundary operator L in terms of boundary fields. This is achieved by making the following assumption, which is in fact fulfilled in many practical situations. We refer to [11,Sec. 5] for some examples.

Assumption 3.3
The boundary operator L is associated with a matrix-valued boundary field L ∈ L ∞ ( ) m×m such that for v, w sufficiently smooth we have Lastly, we give an auxiliary result needed to bound defects occurring in the error analysis of the fully discrete scheme. It can be proven by multiple applications of the product rule applied to L = M −1 L and exploiting that L is composed of constant coefficients and first-order partial derivatives.

Discrete Friedrichs' operators
We now define the central fluxes dG discretization of L. Instead of defining it only on the discrete approximation space V h , which would be sufficient for formulating the (semi-)discrete problem, we extend its domain of definition to the space D(L) ∩ H 1 (T h ) m . This will prove useful later for the error analysis. Hence, we define the discrete operator domain associated with L as The discrete operator is then defined as follows.
Note that L L L is well-defined as an easy consequence of the Riesz representation theorem, cf., [23,Sec. 4.3].
In the following, we state some important properties of the discrete operator needed for the analysis of the fully discrete method. Proofs of the first results can be found in [23,Sec. 4.3]. Note again that, despite the fact that the results therein are given w.r.t. the standard L 2 -inner products, the proofs can be carried out completely analogously in the weighted setting. The first proposition shows that the discrete operator is consistent in the sense that its application to a sufficiently smooth function yields the projection of the continuous operator applied to this function. Further, it states that the discrete operator inherits the dissipativity on the approximation space.
Further, its restriction to V h is dissipative, i.e., we have Note that this immediately implies that L L L is also maximal dissipative, since V h is finite dimensional.
We formulate the following two results as a slightly more general version than the respective ones given and proved in [23]. We need these more general results for the last property of the discrete operator stated in this section. However, as the proofs are completely analogous to the ones given in [23,Sec. 4.3 & Appendix] by using the broken weighted L 2 -norms h p · instead of the non-broken version, we omit them here. The first result states that the discrete operator L L L satisfies an inverse inequality, reflecting the corresponding property (6) of the original operator.
Then, for all p ∈ Z, the discrete Friedrichs' operator L L L fulfills the inverse inequality and, in particular, The constants are given by C inv,L, The next result gives a bound on the broken weighted L 2 -norm of L L L applied to the projection error e v π of a sufficiently smooth function v defined in (7). In some sense, this quantifies how well the discrete operator L L L approximates the continuous one, as The constants are given by C π,L, p = 1 2 N ∂ C tr C π,∂ C ,L + C L (1 + ρ p+1/2 ) and C π,L = C π,L,0 .
The last result is another approximation property and generalizes Proposition 3.8 to more than one discrete operator. We point out that this result and its proof are generalizations of [35,Thm. 6.3]. There, the assumption of a quasi-uniform mesh sequence is needed and the result is restricted to powers of one particular operator (namely the Maxwell operator and its central and upwind fluxes dG discretization). Note that for each operator that is applied, we lose one order of h in the approximation. The proof is given in the appendix.
and, in particular, where the constants are independent of h.

Full discretization
This section is concerned with deriving the full discretization of the wave-type problem (2). To this end, we use a method-of-lines approach and first spatially discretize before we discretize in time. We achieve the first step by replacing the continuous operator L by its discrete counterpart L L L and projecting the data onto the discrete space V h . This yields the spatially semi-discrete wave-type problem (9b) with f f f π := π h f and initial value u u u 0 π := π h u 0 . The spatially discrete problem (9) is an evolution equation posed on the finitedimensional approximation space V h . Hence, and since the discrete operator is dissipative (and maximal due to the finite-dimensionality), its wellposedness and stability can be analyzed completely analogously to the original problem.
Full discretization is now achieved by discretizing (9) in time via the Peaceman-Rachford scheme. Before we can do so, we have to introduce a splitting of the (discrete) spatial operator, since the Peaceman-Rachford scheme is a splitting scheme. There are two ways to achieve this goal, which, in this case, lead to the same result. Namely, we either split the continuous operator, leading to two Friedrichs' operators, which we then discretize via the dG method. Or we discretize the full operator and then split the discrete one. As we deem it to be the more systematic way, we present the first approach in the next section.

Splitting
We split the Friedrichs' operator L = A + B into two operators by splitting its coefficients. More precisely, we split L i = A i + B i for all i = 0, . . . , d such that A 0 , B 0 ∈ R m×m are negative semi-definite and A 1 , B 1 , . . . , A d and analogously to L we define A = M −1 A and B = M −1 B, respectively. We extend all concepts introduced in Sect. 2 to the split operators, using analogous notation.
Additionally, we assume that we can split the boundary operator L accordingly, meaning that we have L = A + B with A and B fulfilling Assumption 2. As we are interested in a splitting of the discrete operator, we also discretize these split operators. Hence, we define the discrete operator domains associated with A and B, respectively, and the central fluxes dG discretizations

the Friedrichs' operators
A and B analogously to Definition 3.5. Using bilinearity of the inner product and the properties of the coefficients and boundary fields discussed above yields i.e., A A A and B B B indeed constitute a splitting of L L L.

Peaceman-Rachford scheme
We now use the Peaceman-Rachford scheme to obtain the fully discrete scheme that we analyze in the remainder of the paper. Let τ > 0 be the timestep size, t n := nτ , f n := f (t n ) and f f f n π := f f f π (t n ). The Peaceman-Rachford scheme applied to the semidiscrete wave-type problem (9) is given by where u u u n τ ≈ u(t n ) is the fully discrete approximation at time t n , n ∈ N, and u u u n+1/2 τ is an intermediate value.
We point out that (11) is in fact equivalent to the scheme used in [9,34], which relies on a different intermediate value u u u n+1/2 τ (for f = 0) but produces the same u u u n τ .
In particular, one can show that this change results in u u u n+1/2 τ being an approximation to u(t n+1/2 ) of the same order as u u u n τ is to u(t n ). This can be shown by using exactly the same arguments we use later.
To analyze the scheme, it is convenient to recast (11) into the equivalent form and Note that the resolvents of A A A and B B B exist because they are maximal dissipative. This readily implies the following result, which states that the fully discrete scheme is wellposed and provides a closed solution formula for the approximation.
Next, we show that the fully discrete scheme is unconditionally stable, which is a well-known fact for the Peaceman-Rachford scheme in this maximal dissipative setting. We follow the proof in [21] to show stability of the full discretization. However, some adaptions are necessary for the fully discrete version to get a bound that is independent of the discretization parameters. We proceed in two steps and first give a bound on arbitrary powers of the system operator S S S pr .  M-norm), see e.g., [39] for a proof of the latter. This proves the claim.

Theorem 4.2 Let v v v ∈ V h and q ∈ N. Then, for all h ∈ H and all τ > 0, we have
Note that the occurrence of the discrete operator B B B on the right-hand side means that the bound provided in Theorem 4.2 is not uniform in h. However, the next result shows that the scheme is unconditionally stable nonetheless.
Proof The discrete variation-of-constants formula (15)  Using the triangle inequality shows the first bound.
To show the second one, we add and subtract τ 2 π h Bu 0 and use Proposition 3.6 (consistency) and the contractivity of π h to obtain The bound on this term follows by Proposition 3.8 (approximation property) and proceeding analogously for the term involving A A A proves the claim.
Next, we show the discrete analogon to the second bound in Corollary 2.4. Since this stability result involves discrete derivatives of the approximations, we define the discrete time derivative and the discrete space derivative at time t n+1/2 as respectively. While the definition of the discrete time derivative is rather intuitive, the discrete space derivative is more involved. Its structure is due to the splitting nature of the Peaceman-Rachford scheme. In fact, if one considers the semi-discrete analysis in time, one sees that the full-step approximations at time t n are only contained in D(B), whereas the half-step approximations at time t n+1/2 are only contained in D(A), cf., [29] for details. Hence, one has to use the split operators to obtain a reasonable approximation on the whole family of meshes.

Corollary 4.4
For all h ∈ H and all τ > 0, the approximations u u u n τ n≥0 and u u u n+1/2 τ n≥0 given by the fully discrete scheme (11) satisfy where C is independent of h, τ and n if t n ≤ T for some Proof First, note that by (14) we have Further, the definition (13) of S S S pr and the discrete splitting property (10) yields We point out that taking the limit τ → 0, we obtain the derivative of S S S pr w.r.t. τ evaluated at 0 on the left-hand and the discrete operator L L L on the right-hand side. This corresponds to the fact that S S S pr is a (time-)discrete version of the semigroup generated by L L L and this relation reflects the fact that the derivative of a semigroup evaluated at 0 yields its generator.
Having this in mind, we take the discrete derivative of the sequence {u u u n τ } n≥0 and apply the discrete variation-of-constants formula (15) and (17) where we further performed an index shift to combine the two sums. Note that this can be seen as applying a discrete version of the Leibnitz rule of integration. Taking norms and applying Theorem 4.2 yields an appropriate bound on ∂ τ u u u n+1/2 τ . Thus, it remains to show the bound on L L Lu u u n+1/2 τ to prove (16). To do so, we add the (half-step) Peaceman-Rachford iterations (11a) and (11b) to obtain Taking norms and using the already established bound on ∂ τ u u u n+1/2 τ yields (16). For the uniform bound, we use the fundamental theorem of calculus to obtain and analogously for the same term including A A A. The remaining terms can be bounded as in the last part of the proof of Corollary 4.3. For completeness, we state the concrete bound in the appendix.

Error analysis
In this section, we derive the main results of this paper, namely the error bounds for the fully discrete scheme (11) resulting from the Peaceman-Rachford scheme in time and a central fluxes dG discretization in space. More precisely, we show that the approximation u u u n τ gained from performing n steps of this scheme converges to the exact solution u(t n ) of the original problem (2) at time t n with order two in time (the classical order of the Peaceman-Rachford scheme) and order k in space, given that the solution is sufficiently smooth. Further, we show that discrete derivatives of the approximations w.r.t. space and time converge to their continuous counterparts exhibiting the same rates if a bit more regularity is assumed.
In order to perform the analysis, we split the full discretization error e n into e n = u(t n ) − u u u n τ = e n π + e e e n , where e n π = u(t n ) − π h u(t n ) is the projection error at time t n and e e e n = π h u(t n ) − u u u n τ is the discretization error after n steps. Note that by Lemma 3.2 we already have a bound on e n π of order k + 1 for sufficiently smooth solutions. Thus, it suffices to bound the discretization error e e e n . Besides this, we also show bounds for the time derivative error e e e n+1/2

Error recursion
We start by showing that the discretization error e e e n satisfies the fully discrete scheme (12) up to a defect. To keep the notation concise we define R R R pr :

abbreviating the concatenation of the resolvents of A A A and B B B.
The error analysis is based on the fact that the Peaceman-Rachford scheme can be interpreted as a perturbation of the well-known Crank-Nicolson scheme. Thus, we decompose the full defect into a defect already present in the Crank-Nicolson scheme (which corresponds to the quadrature error of the trapezoidal rule) and an additional defect caused by the perturbation. where the dG-Crank-Nicolson defect d d d n cn is given by

and the (dG-Peaceman-Rachford) perturbation defect d d d n pr is given by
Proof We begin by inserting the projected exact solution into the recursion (12a) of the fully discrete scheme. This causes an error, which we define as R R R pr (d d d n cn + d d d n pr ) and we thus obtain π h u(t n+1 ) = S S S pr π h u(t n ) Expanding the products and using the splitting property (10) proves the claim.
As a consequence, Theorem 4.1 yields the representation formula e e e n+1 = n j=0 S S S n− j for the discretization error. Further, by Theorem 4.2, this representation formula immediately implies an error bound for our scheme if the defects fulfill appropriate bounds. Note that the Crank-Nicolson defect already closely resembles the (projected) quadrature error δ n cn of the trapezoidal rule applied to ∂ t u, i.e., where we used the fundamental theorem of calculus and the differential equation (2a) for the second equality. The only difference is that the second term involves the operator L L Lπ h instead of π h L. This discrepancy results in the spatial error term. Moreover, after again applying the fundamental theorem, the perturbation defect is already of order three in τ , which is the order required to obtain order two globally. The only task still to be done is thus to show that the remaining factors can be bounded independently of the discretization parameters. In fact, these observations will be key to bounding the defect as we will see in the next section.

Bounds on the defects
Next, we derive appropriate bounds on the defects occurring in the error recursion (20). We start with the dG-Crank-Nicolson defect.
Proof We exploit the observation made at the end of Sect. 5.1 and rewrite d d d n cn as where we have added and subtracted the term involving π h L for the first and used Proposition 3.6 (consistency) for the second equality. The first term is just the projected quadrature error (22) of the trapezoidal rule. Using the associated Peano kernel representation, its M-norm can be bounded by the first term in (23). The M-norm of the second term can be bounded by using Proposition 3.8 (approximation property). Combining these bounds yields (23). For the second bound (24), we take the discrete derivative of d d d n cn , which yields We use the definition of δ n cn in (22) together with multiple applications of the fundamental theorem of calculus and the substitution s → s + τ to obtain The integrand w.r.t. s is again the quadrature error of the trapezoidal rule and the bound thus follows from its Peano kernel representation (and transformation of the integrals). For the second term in (25) we use the fundamental theorem of calculus to obtain Using Proposition 3.8 (approximation property) concludes the proof.

Fully discrete error bounds
We are now able to state the main results of this paper. Since the assumptions for all results in this section are similar, we encapsulate them in the following assumption.

Assumption 5.4
For ∈ N, the exact solution of the wave-type problem (2) satis- Further, the material tensor fulfills M ∈ W 1,∞ (T h ) m×m .
Our first main result gives bounds on the full discretization error of the scheme. Note that the regularity assumptions posed on the exact solution are in fact equivalent to the combination of the corresponding ones posed in Lemmas 5.2 and 5.3. Further, we point out that in the general setting of this paper, these convergence rates are expected, since they coincide with the classical order of the Peaceman-Rachford scheme and, in this general setting, the optimal convergence rate for the central fluxes dG method [32], respectively.

Theorem 5.5
Let h ∈ H, τ > 0, let k ≥ 1 be the polynomial degree of the dG method, and let Assumption 5.4 be fulfilled with = 1. Then, for all n ∈ N 0 , the dG-Peaceman-Rachford error satisfies where C only depends on t n+1 , C π , C π,L , C pr,u , C pr, f , Proof We decompose the error as in (19). The projection error e n+1 π is bounded by Lemma 3.2. Thus, it remains to bound the discretization error e e e n+1 . We use (21)  ).
The claim now follows by Lemmas 5.2 and 5.3.
Owing to the stability bound on the discrete derivatives in Corollary 4.3, we are also able to show that the errors e e e n+1/2 τ and e e e n+1/2 L converge with the same orders if we assume some additional regularity. For the sake of readability, we do not give the details of the full bounds as in Theorem 5.5 here. However, they can be found in the appendix.
We start with the time derivative error. Note that we again need the combined regularity already required for the corresponding bounds in Lemmas 5.2 and 5.3.

Theorem 5.6
Let h ∈ H, τ > 0, let k ≥ 1 be the polynomial degree of the dG method, and let Assumption 5.4 be fulfilled with = 2. Then, for all n ∈ N 0 , the full dG-Peaceman-Rachford time derivative error satisfies where C is independent of h, τ and n if t n ≤ T for some T > 0.
Proof As for the dG-Peaceman-Rachford error, we first split the full error into a projection error and the discretization error e e e Taking the norm and using the stability result Theorem 4.2 and the bounds on the defects in Lemmas 5.2 and 5.3 (together with transforming some integrals) yields the desired bound.
Lastly, we show the error bound for the approximate space derivative. Note that besides the combination of the corresponding assumptions in Lemmas 5.2 and 5.3 we need some additional regularity of the inhomogeneity.

Theorem 5.7
Let h ∈ H, τ > 0, let k ≥ 1 be the polynomial degree of the dG method, let Assumption 5.4 be fulfilled with = 2, and f ∈ C(R + ; H k+1 (T h ) m ) ∩ C 2 (R + ; L 2 ( ) m ). Then, for all n ∈ N 0 , the full dG-Peaceman-Rachford space derivative error satisfies where C is independent of h, τ and n if t n ≤ T for some T > 0.
Proof Again, we first split the full error into a projection and a discretization part. Note that by ∂ t u, f ∈ C(R + ; H k+1 (T h ) m ) we immediately obtain Lu ∈ C(R + ; H k+1 (T h ) m ) via the original differential equation (2a), allowing us to bound the projection error. Hence, it remains to bound the discretization error e e e n+1/2 L .
To do so, we compare (18) to the projected original differential equation (2a) at time t n+1/2 to obtain e e e n+1/2 L = e e e n+1/2 We have already bounded e e e n+1/2 τ in the proof of Theorem 5.6 and the second term can be bounded by Taylor expansion. The last term can be bounded by using the fundamental theorem of calculus and subsequently proceeding analogously to the corresponding part of the proof of Lemma 5.2.

Numerical experiments
In this section, we present numerical experiments to verify the theoretical results obtained in this paper. Three different cases are considered, all of them consisting of the full 3D Maxwell equations with perfectly conducting boundary conditions, each time supplied with different data. Each data set is chosen such that there is an analytical solution with which the approximations can be compared to obtain the exact approximation error.

Linear Maxwell equations
For all cases, we consider the linear, isotropic and undamped Maxwell equations on = [0, 2] × [0, 1] 2 in the following form. Given a three-dimensional, open, bounded and connected Lipschitz domain ⊂ R 3 , we seek the electric field E : R + × → R 3 and the magnetic field H : R + × → R 3 solving the Maxwell system Here, the initial data E 0 , H 0 : → R 3 , external current J : R + × → R 3 and material parameters ε, μ : → R \ {0} (permittivity and permeability, respectively) are given. These equations are supplied with perfectly conducting boundary conditions on = ∂ , i.e., we require where n is the outward normal vector on . Defining we see that (27) can be written in the form of (1). Further, by [11,Lem. 3.5], the perfectly conducting boundary conditions (27d) define an operator fulfilling Assumptions 2.3 and 3.3. Hence, if we have ε, μ ∈ L ∞ ( ) and J ∈ C(R + ; L 2 ( ) 3 ) this problem fits the framework of this paper with d = 3 and m = 6. In order to employ the Peaceman-Rachford method, we use the splitting proposed in [33,45]. In particular, we first split the curl operator into and with this define the split operators It can be shown that these split operators fulfill all necessary assumptions for our theory and we refer to [29, Sec. 6.5.3] for the details.
The scheme was implemented with the help of the C++ finite element library deal.ii [1], which was used to discretize the operators in space. The code can be found at https://www.waves.kit.edu/dg-ADI.php. All experiments were performed on uniform tensorial meshes of various mesh widths and we used a polynomial degree of k = 2 on all elements for the dG method. Lastly, the simulation interval for all experiments was chosen to be [0, 2] and we plot the maximal error over all timesteps in all graphs.
Cavity solution. As a first example, we consider the well-known cavity solution of the homogeneous Maxwell system with constant material parameters. In particular, we set J ≡ 0, assume ε, μ ∈ R to be constant and consider the family of solutions to (27) given by for (t, x) ∈ R + × . Here, we denote by c = (εμ) −1/2 the speed of light, κ = (κ 1 , κ 2 , κ 3 ) ∈ R 3 + is the wave vector and = c κ is the angular frequency. Further, E 1 , E 2 , E 3 are preset amplitudes of the waves. For our numerical experiment we chose the parameter set which, in particular, leads to a solution satisfying the perfectly conducting boundary condition (27d). Varying material parameters. For the second example, we adapt the cavity solution such that it allows for varying material parameters. In fact, one can show that (28) is a solution to the homogeneous (i.e., J ≡ 0) Maxwell equations (27) for sufficiently smooth ε, μ : → R with μ ≡ 1 ε and (∇ε) × E = 0 on for all t ∈ R + . We exploit this and choose the solution given by (28) with the parameters for x ∈ . Besides satisfying the above mentioned criteria, note that, again, the resulting solution also fulfills the perfectly conducting boundary conditions (27d). Inhomogeneous Problem. As a last example we consider the following solution to the inhomogeneous Maxwell equations. Let ∈ C 2 ( ) with n × ∇ = 0 on and ∈ C 1 (R). Then, for x ∈ , t ∈ R + , the solution of (27) with J (t, x) = − (t) ∇ (x) and ε ≡ μ ≡ 1 is given by In particular, we choose for this example.

Results
First, note that all three considered examples are sufficiently regular to fulfill the requirements of our main results Theorems 5.5, 5.6 and 5.7. Hence, we expect at least the convergence orders presented in these theorems. Fig. 1 Errors produced by the fully discrete scheme (11) applied to the problem solved by the cavity solution plotted against the timestep τ . The employed mesh sizes are given in the legend Fig. 2 Errors produced by the fully discrete scheme (11) applied to the problem with varying material parameters plotted against the timestep τ . The employed mesh sizes are given in the legend Further, since the first two examples are homogeneous, both the discrete as well as the continuous objects considered in Theorems 5.6 and 5.7 coincide by the numerical scheme (12) and the continuous problem (2), respectively. Therefore, we expect the same errors for both, which was indeed the case in all numerical results.
We performed several simulations with varying meshsizes h and timesteps τ to verify our theoretical results. The resulting errors are displayed in Figs. 1, 2, 3, 4 and 5.
In all experiments, one clearly sees the second-order convergence in time (whenever the spatial error is small enough so that the temporal error is dominant). The spatial error contained in e e e n , however, behaves about an order better than predicted in the homogeneous examples, i.e., we see convergence of order about k + 1 in h. This is probably due to the high regularity of both the considered solutions and the employed meshes. Such an improved order of convergence on a regular rectangular grid was shown in [30] for the neutron transport equation in dimension two. For general meshes and less regular solutions this might not be the case, see e.g., [40, Fig. 3.3] for a counterexample on unstructured meshes. In addition, this is not the case for the errors in the derivatives, e e e n τ and e e e n L (which coincide here as explained before), where we see second-order convergence. In the inhomogeneous example we see no influence of the spatial error at all. However, this is to be expected, since we used k = 2 in the experiments and thus both the solution and the inhomogeneity (as well as their spatial derivatives) are contained in the approximation space V h (for all t ∈ R + ). Hence, it is Fig. 3 Errors produced by the fully discrete scheme (11) applied to the inhomogeneous problem plotted against the timestep τ . The employed mesh sizes are given in the legend Fig. 4 Errors produced by the fully discrete scheme (11) applied to the problem solved by the cavity solution plotted against the mesh width h. The approximations were produced by using 5120 timesteps for time integration easy to see that the defects derived in Sect. 5.1 (and consequently the derived bounds on the errors) are in fact independent of h, which explains this behavior.
Lastly, note that the errors in the space derivative of the solution corresponding to Theorem 5.7 are of approximately one to two orders of magnitude higher than those of the error in the solution. This can probably be explained by the fact that additional terms involving the inhomogeneity enter the error constant.  5 Errors produced by the fully discrete scheme (11) applied to the problem with varying material parameters plotted against the mesh width h. The approximations were produced by using 5120 timesteps for time integration To tackle the first term, we use Proposition 3.7 (inverse inequality), yielding h p L L L r +1 (L L L r . . . L L L 1 π h − π h L r . . . L 1 )v ≤ C inv,L r +1 , p h p−1 (L L L r . . . L L L 1 π h − π h L r . . . L 1 )v , which can be bounded by the induction hypothesis with p − 1 instead of p.
Applying r times Lemma 3.4 yields the asserted bound for this term.
Full bound for Corollary 4.4. Next, we explicitly state the constant C used in Corollary 4.4. It is defined by 0 ∂ t f (s) + τ 2 A∂ t f (s) + τ 2 C π,A |∂ t f (s)| 1,T h ds .