Stability of Discontinuous Galerkin Spectral Element Schemes for Wave Propagation when the Coefficient Matrices have Jumps

We use the behavior of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_{2}$$\end{document}L2 norm of the solutions of linear hyperbolic equations with discontinuous coefficient matrices as a surrogate to infer stability of discontinuous Galerkin spectral element methods (DGSEM). Although the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_{2}$$\end{document}L2 norm is not bounded in terms of the initial data for homogeneous and dissipative boundary conditions for such systems, the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_{2}$$\end{document}L2 norm is easier to work with than a norm that discounts growth due to the discontinuities. We show that the DGSEM with an upwind numerical flux that satisfies the Rankine–Hugoniot (or conservation) condition has the same energy bound as the partial differential equation does in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_{2}$$\end{document}L2 norm, plus an added dissipation that depends on how much the approximate solution fails to satisfy the Rankine–Hugoniot jump.


Introduction
In wave propagation problems, it is natural to find interfaces where material properties like the wave propagation speeds or density abruptly change. Examples include interfaces between two dielectrics in electromagnetic wave propagation problems, or different rock layers in geophysics. At such interfaces the solutions can make discontinuous jumps, causing difficulties for many numerical methods.
One of the key features of discontinuous Galerkin (DG) methods is that the discontinuous approximation at element interfaces naturally allows jump discontinuities in the solution if element boundaries are placed along them. Consequently, DG spectral element methods have been used for over twenty years to solve problems with material discontinuities, both stationary [3,6,9,19] and moving [20]. Computations and theory in such works show that placing the discontinuities at element boundaries leads to exponentially convergent approximations.
In a paper on discontinuous interface problems, La Cognata and Nordström [13] noted that hyperbolic problems with discontinuous coefficients do not necessarily have their energy bounded in terms of the initial data when measured in the L 2 norm, even with constant coefficients and homogeneous and dissipative boundary conditions. Instead, the L 2 norm can increase or decrease, depending on the relative size of the wave speeds on either side of the discontinuity. The lack of a bound on the L 2 norm is not due to an instability in the usual sense, but is due to the fact that conservation at the interface, and the resulting jump in the solution, can increase the norm of the solution as a wave propagates across it. In an alternate norm, however, one that discounts the effect of the jump, the energy is bounded.
Here we propose a procedure where we use the L 2 norm as a surrogate to infer stability of discontinuous Galerkin spectral element methods (DGSEM) for the approximation of hyperbolic equations with discontinuous coefficient matrices. The L 2 norm is easier to work with since it does not require finding the discount factors, which are difficult to compute in general configurations of elements and interfaces. We show that the DGSEM with an upwind numerical flux that satisfies the Rankine-Hugoniot (or conservation) condition behaves as the partial differential equation (PDE) does in the L 2 norm, plus an added dissipation that depends on how much the approximate solution fails to satisfy the Rankine-Hugoniot jump.

Linear Hyperbolic Systems with Discontinuous Coefficients
In this paper we establish the stability of a discontinuous Galerkin spectral element approximation to linear hyperbolic systems of equations of the form where u is the state vector, and ↔ f is the vector of fluxes, with coefficient matrices A j and unit coordinate vectorsx j . We assume throughout this paper that the coefficient matrices are piecewise constant, with discontinuities marking what we will refer to in this paper as material interfaces. To isolate the contribution of the discontinuities, we choose the coefficient matrices constant between material interfaces to avoid the complexity of possible exponential growth in the solution norm when the matrices vary over the domain.
We examine the problem defined in a domain Ω, as sketched in two space dimensions in Fig. 1. It is sufficient to consider two domains with a single material interface, so the domain is split into two subdomains Ω L and Ω R separated by an interface Γ . The external boundary is Γ b , along which we assume that proper, well-posed and dissipative boundary conditions are applied.
Since the system is hyperbolic, there exists a matrix of right eigenvectors, P, and a real diagonal matrix, We also assume that the matrices A j are simultaneously symmetrizable and that there exists a piecewise constant matrix S such that As a concrete example of the system (1), we pose the linear acoustic wave system where and where ρ is the density of the medium, c is the sound speed, and δ i j is the Kronecker delta. The state vector can be viewed as representing pressure, p, and three velocity components, u, v, w. The coefficient matrices are simultaneously symmetrizable by the matrix With jump discontinuities in the material parameters, ρ and c, the coefficient matrices and symmetrizer have jump discontinuities. We contrast the approximation of the system (1) with that of the approximation of systems that can be written in the form where E > 0 is diagonal and discontinuous at material interfaces while → B is continuous. The system (3), for example, can be re-written in the form (5) with symmetric matrices For equations of the form (5), there is a natural norm, in which the energy is bounded for homogeneous dissipative physical boundary conditions and nonconservative interface conditions, with that energy satisfying Stability of DG spectral approximations to equations in the form (5) has been shown specifically, for instance, for Maxwell's equations [6] and the elastic wave equations [19].
Remark 2. 1 The system (1) cannot in general be rewritten in the form (5). That would require that each A j can be written as A j = E −1 B j where E = E T > 0 and E contains all material properties. A counter example is the frozen coefficient compressible Euler equations [1].
As noted in [13], systems of the form (1) with discontinuous coefficient matrices do not necessarily have energy bounded by the initial data when measured in the L 2 norm, and we present an example here to motivate the situation. Figure 2 shows the p component of the analytic solution of acoustic wave reflection and transmission at a material boundary placed at x = 0 at three times: The initial incident wave, when the wave is interacting with the material discontinuity, and the reflected and transmitted waves after the interaction.
We plot the energy as a function of time, measured by the L 2 norm, in Fig. 3. We see that the L 2 energy is bounded, and even though the L 2 energy estimate does not show boundedness directly, energy is bounded in terms of the initial data in a norm that discounts the jump [13]. Note that there is a slight downturn in the energy in Fig. 3 as t → 3. The energy does decrease to zero after that time as the waves propagate out of the domain.
To establish the stability of the discontinuous Galerkin spectral element approximation of (1), we follow the roadmap presented in [16]. We first establish energy behavior of the PDE Fig. 3 Exact L 2 energy for the solution of the one dimensional acoustic wave equation for propagation of a wave across a material interface system, and then follow an equivalent discrete path to establish an equivalent behavior for the approximation. We begin with the study of the scalar one-dimensional advection problem, since it is easy to follow the steps, and then a symmetric system in one space dimension. Finally we use the symmetric system results to derive the energy bound for the general system in Sect. 2.3.

Energy Dynamics of the Scalar Problem in One Space Dimension
To motivate (and simplify) the general formulation, we start with the scalar advection equation with two domains as an introduction. Our discussion in this section restates that of [13], but introduces our notation used in succeeding sections.
We derive the energy dynamics of the solution to the scalar advection initial-boundaryvalue problem in the form (1) where a L , a R are constants, and a L = a R . The discussion that follows leads to the same types of conclusions if the wave speeds are both negative. We are interested here in problems where the domains couple and waves propagate from one side to the other. So we do not consider a L > 0, a R < 0, where the domains decouple as energy is dissipated at the interface, or a L < 0, a R > 0 where boundary conditions for both sides are required. We split the problem into two: Left, and Right where u * is the upwind specified interface condition chosen so that the Rankine-Hugoniot (or conservation) condition is satisfied. Thus, for the scalar equation, u * (t) = a L a R u(0 − , t). To find the energy equation, we multiply by the solution and integrate over the domains. Define the L 2 energy norms Adding together and re-arranging, where ||·|| 2 = ||·|| 2 L + ||·|| 2 R . Applying the homogeneous boundary condition on the left, 1 2 When we apply the interface condition, The quantity Q will be used later in this paper to define stability.
Finally, we substitute the interface value for u * , and rearrange so that 1 2 Equation (21) shows that the energy is dissipated by the interface only if a R > a L . Otherwise the interface generates energy, as illustrated in Fig. 3. In [13] it was shown that one can construct "discounted norms", in which the energy is bounded. If the second equation in (16) is multiplied by a constant α c > 0, then the weighted sum leads to 1 2 Then defining the new norm with the α c discount factor, we have

Remark 2.2
The weighted norm discounts the effect of the jump, with the result that viewed in the discounted norm, the energy no longer appears to increase.

Remark 2.3
The use of the discounted norm scales to multiple material interfaces and multiple space dimensions by choosing α c to be the minimum over all the ratios of downwind to upwind wave speed ratios.

Remark 2.4
Alternatively, unlike the general case noted in Remark 2.1, the scalar equation (10) can be recast to the form (5) by dividing by the wave speed. Let ε = 1/a > 0. Then If the nonconservative boundary condition at the interface, u * = u(0 , is used, then following the same procedure as (16)-(21), where the weighted energy norm is given by Using the norm (27), but with the conservative interface condition (14), the solution still has a bound like (21), namely So when the conservative interface condition is used, the weighted energy norm is also bounded only when a L /a R ≤ 1 .
The discounted norm is equivalent to the L 2 norm. In the discounted situation, where α c < 1, and Therefore, Equivalence of the norms means that the L 2 norm is actually bounded in terms of the initial data, even though the energy method does not show it directly through (21). From (23) and (31), Thus, for the constant coefficient problem.

Energy Dynamics for Hyperbolic Systems in One Space Dimension
We now increase the complexity and extend the scalar one-dimensional analysis to the general system (1) in one space dimension. We derive the energy equation for the one-dimensional hyperbolic system where the coefficient matrices are for now assumed to be symmetric. Under this assumption, there is a matrix P such that A = PΛP −1 satisfying P −1 = P T . For the moment, let us assume that A has no zero eigenvalues. We also assume that the number of positive and negative eigenvalues does not change across the interface. In other words, there is no eigenvalue that changes sign at the jump. Depending on the sign change, boundary/interface conditions are either lost or gained. More general conditions where the sign of the eigenvalues changes in multi-physics applications are considered in [5]. Finally, we assume that appropriate boundary and initial data are applied.
To find the interface condition at x = 0 for the system (34), we split the system into right and left going waves. The characteristic variables for the system (34) are where w + is associated with the positive eigenvalues of A and w − is associated with the negative ones. They are chosen upwind at the interface according to where here and in the following, the subscripts R and L correspond to the values at x = 0 + and x = 0 − , respectively. The w ± * are computed so that the Rankine-Hugoniot condition is satisfied at the stationary interface. Let us write Then (37) can be written as Let us put the unknowns on the left, and the knowns on the right, giving Equation (40) provides a system of equations for the unknowns. The matrices on the left of (40) have a special structure since P is the matrix of right eigenvectors and Λ is a diagonal matrix. Let n be the number of positive eigenvalues out of a total of m. Then define and where → p j is the eigenvector associated with the eigenvalue λ j and the eigenvalues are ordered in decreasing order, largest to smallest with λ j > 0 for j ≤ n. Then we can write (40) as Given the structure of the M ± matrices, the equations can be combined to produce a single system for the unknowns where Existence and uniqueness of the inflow characteristic vectors w ± * therefore depends on the existence of the inverse of the matrix M L R . That matrix is comprised of eigenvectors of the coefficient matrix evaluated on the left and eigenvectors evaluated on the right. On the one hand, if the eigenvectors of the coefficient matrix do not change across the material discontinuity, then, since the eigenvectors are independent, M −1 L R M L R is diagonal. As an example, the eigenvectors of the acoustic wave system (3) are constant, being independent of the material properties on either side. On the other hand, if the eigenvectors change across the interface and the matrix M −1 L R M RL is not diagonal, then the problem is ill-posed [5]. We therefore require that the eigenvectors be preserved across the jumps so that (46)

Remark 2.5
The Rankine-Hugoniot (conservation) condition (37) limits the form of the interface condition significantly. If only boundedness is desired, more general coupling conditions are allowed [5].

Remark 2.6
One can see that the assumption that there are no zero eigenvalues is not a restriction. If there are zero eigenvalues, then the associated characteristic variables w 0 are multiplied by the zero matrix and have no contribution to the system. Therefore those quantities can be eliminated, leaving (43) and what followed. The w 0 vector is determined by the initial data.
Going back to the original equations, (34), we compute the energy equation by multiplying by the state and integrating over the domain, giving where, now, ||u|| 2 = u, u and PBT represents the terms coming from the physical boundary conditions on the left and right. Since we are only interested here in the interface conditions, we will assume that the physical boundary conditions are well posed so that PBT ≥ 0. In that case, 1 2 where is the interface contribution to the energy. Following the steps in the scalar analysis, we now apply the interface boundary conditions on Q. We decompose the system into characteristic variables. Then we use the fact that A is symmetric, making P −1 = P T . With this decomposition, taking into account the upwinding of the characteristic variables, (36). We now gather the right-going and left-going wave contributions (c.f. (19)), and then use the fact thatΛ − < 0, to get the final form of the interface contribution, which we write in terms of its characteristic components, Equation (52) is the system version of the scalar interface condition seen in (19). As in the scalar problem, one can construct a discounted norm ||·|| B for which the associated interface term Q B is non-positive and the discounted norm is bounded when the coupling matrix, M, exists and is diagonal. For instance, one simple choice is to let where the entries with μ ± > 0 are counted according to the number of positive and negative eigenvalues of A. Then multiplying the system on x > 0 by B from the left, defining the norm ||u|| B = u, Bu 1 2 , and following the same steps leading to (52), the interface contribution to the energy is One then only needs to find μ + small enough and μ − large enough to ensure that Q B ≤ 0, in which case the new energy ||u|| 2 L + ||u|| 2 B,R is bounded in terms of the initial data. Since the coupling matrix M is diagonal, let us split it as

Extension to Non-symmetric Equations and an Arbitrary Domain
The results of the previous two sections extend to general geometries and non-symmetric coefficient matrices. In preparation for the generalization of (52), we note that within each subdomain the coefficient matrices are constant, and therefore we can re-write (1) in a split form as We also define the inner product and norm over a subdomain D = Ω L or Ω R as To form the energy, we take the inner product of (57) with the vector S −1 T S −1 u, giving Let us define u s = S −1 u to be the symmetric system state. Then since S is constant within the subdomains and A s = S −1 AS, We then apply multidimensional integration by parts and symmetry to the divergence term where → n is the outward normal at the boundary of D, and note that the volume term cancels the third term in (61), leaving only the boundary integral, Then over the domain Ω, where L/R represent the states on either side of the interface with respect to the normal → n. We can now get a bound for the multidimensional system similar to (48). The integrand in the interface integral is identical to that in (47), with A ← → A s L ·n and u ← u s . Therefore, if the boundary conditions along Γ b are properly posed and dissipative, Therefore, Q is still given by (52), but now formulated in the new symmetrized variables. Note that the norm defined by ||u s || 2 = S −1 T S −1 u, u is equivalent to the norm ||u|| since S −1 T S −1 > 0.

Stability
In summary, for hyperbolic systems of the form (1), with discontinuities in the coefficient matrices and homogeneous, dissipative boundary conditions, the L 2 norm of the solution obeys (65). The integrand of the interface contribution, Q, is of the form (52), where the characteristic variables are evaluated from the upwind side and satisfy the Rankine-Hugoniot condition. It is not necessarily non-negative, depending on the relative wave speeds from either side of the interface, so the L 2 norm of the solution is not bounded in general by the initial data. An example of such behavior was shown in Fig. 3. Although the L 2 norm (or, for that matter, weighted norms, see Remark 2.4) is not always bounded in terms of the initial data, there exists an energy in a discounted norm that is bounded in the usual way provided that the coupling matrix between the upwind and downwind states is diagonal.
Thus, we have two views of stability at our disposal, which we will call direct and inferred: -Direct Stability When the L 2 norm is bounded, we directly have L 2 stability. This is seen in scalar problems if a L /a R ≤ 1 in (21). For the system, the equivalent is when -Inferred Stability Nonetheless, even if the L 2 norm is not directly bounded, we have seen that one can construct a discounted norm in which it is, e.g. (23) for scalar problems and for systems when (56) is satisfied. Stability in some discounted norm is therefore inferred, or implicit, if Q is given by (52).
In general geometries it may not be easy to find the discounted norm in which the solution is bounded. Finding the precise coefficients requires satisfying conditions like (56). When multiple subdomains exist in multiple space dimensions, T -type intersections between materials are possible. The discount factors must then take into account all subdomain boundaries and be adjusted globally so that at each interface (56) is still satisfied. For these reasons, it is easier to monitor the behavior (65) of the simpler L 2 norm as a surrogate to infer wellposedness of the system. Further insights into how choosing the norm affects how the energy is bounded or not can be found in [14].
Stability of a numerical approximation of a system follows that of the PDE, and so we state the stability condition for the approximation as: Definition 2.1 A scheme approximating the discontinuous coefficient problem (1) is said to have inferred stability if the discrete approximation of the standard L 2 norm is bounded as in (65) and the approximation to the integrand, Q N ≈ Q, satisfies where W is the approximation of w.

The Discontinuous Galerkin Spectral Element Discretization
In this section, we briefly summarize the important discretization steps. For a detailed description and derivation of the scheme, we refer to e.g. [4,10,21].
The first step is to divide the computational domain into a mesh of non-overlapping, possibly curved, hexahedral (quadrilateral in 2D) elements, {e l } K l=1 . Each hexahedron is mapped from physical space to a reference space cube To retain spectral accuracy and exponential convergence in the presence of jump discontinuities, we require element faces to be aligned with the material interface, Γ , so that polynomial approximation is not made across the discontinuity.
From the mapping, we can compute the metric terms Note that we need to carefully evaluate the metric terms to get a discretely divergence-free contravariant basis J → a i , which is necessary to guarantee free-stream preservation of the discretization [7] and stability of the volume terms [4,12].
The second step of the discetization process is to transform the problem (57) from physical to reference space. In reference space, (57) becomes where we collect the metric terms in the block matrix with the identity matrix, I, having the size as the state vector u. The third step is the variational Galerkin formulation. We first approximate the solution with an interpolatory polynomial of degree N , and denote polynomial approximations with capital letters u ≈ U = I N (u), where I N denotes the interpolation operator. In the spectral collocation framework, one typically uses a nodal basis for the interpolation. Furthermore, for hexahedral/quadrilateral elements, we use a tensor-product of one-dimensional nodal Lagrange basis functions spanned on the Legendre-Gauss-Lobatto nodes. The same polynomial approximation is used for all quantities, e.g. for the contravariant flux function To get the variational formulation, we multiply the transformed PDE (67) by polynomial test functions ϕ, which are linear combinations of the nodal basis functions. Then we integrate over the reference element E and use integration-by-parts to arrive at wheren is the reference space outward pointing normal vector to the face ∂ E. Finally, we replace the integration in (69) by quadrature and cubature rules, collocated with the Legendre-Gauss-Lobatto interpolation. Note, that the Legendre-Gauss-Lobatto nodes include the boundary nodes and hence surface and volume integration nodes partially coincide with the interpolation ansatz and I N (·) can be dropped. Furthermore, we introduce the yet to be defined numerical flux function n, which depends on the two states U L,R at the interface and approximates the normal flux through the interface. Note that we assume the coefficients A are mostly constant, but when they jump, the mesh is aligned so that an element interface is at the jump. Hence, the numerical flux function at the coefficient jump interface depends not only on the solutions left and right, but also on the coefficients left and right: F * n = F * n (U L,R ; A L,R ). Applying quadrature, we get the formal statement of the DGSEM, where ·, · N and ∂ E,N represent the volume and surface quadratures, see [11]. The right hand side of (70) is written in terms of the normal covariant fluxes and is equivalent to that written in terms of the contravariant ones [21]. The resulting high-order semi-discretization is integrated with a proper high-order accurate explicit Runge-Kutta time integrator, which is stable under the typical CFL-type time step restriction.

Stability of the Discontinuous Galerkin Approximation
We establish the stability bound from the weak form of the equation, (70). We then follow the path taken in Sect. 2 for the continuous problem to examine the discontinuous interface term: We examine the scalar problem for insights, then the symmetric one-dimensional system, and finally the general problem for the DGSEM approximation.

Discrete Stability Estimate
For a detailed derivation of the discrete stability estimate, which parallels the continuous analysis, we refer to [4,21]. Here, we will only sketch some important intermediate steps.
To get the stability estimate, we replace the test function ϕ with the approximate solution polynomial and the symmetrizer matrices, writing where we define the symmetrized discrete flux F s n that uses the symmetric coefficient matrices Using the fact that the symmetrizer matrix S commutes with the metric block matrix M (see e.g. [4,21]) we see that the volume terms cancel, leaving only surface terms, When we sum over all elements, inner surface terms appear twice (with different normal vectors), whereas element surfaces that are at the physical domain boundary appear only once and are denoted as physical boundary terms ( PBT). The interior element surface contributions split into two parts: Surfaces that fall on the material interface Γ , and those across which the coefficient matrices are the same, which we call smooth interface boundary terms, SIBT. The sum over all elements can then be written as written in terms of the jump operator, U = U R − U L . Assuming that the discrete physical boundary terms are dissipative, the discrete L 2 norm satisfies which mimics the continuous stability (65) if SIBT ≤ 0. We thus need a proper numerical flux function F s, * n to control discrete stability, i.e. to guarantee that the integrand satisfies pointwise at each node on element faces along the discretization of Γ , and SIBT ≤ 0. The dissipativity of the SIBT for the upwind numerical flux has been shown elsewhere, e.g. [8], [21]. Therefore, in the following we will assume SIBT ≤ 0 and concern ourselves only with the discontinuous interface terms.

Stability for the Scalar Problem
In the DG approximation, the Rankine-Hugoniot condition and the inflow boundary condition are enforced weakly with the upwind numerical flux If summation by parts is applied again to the second term in (70), one gets the strong form of the approximation, in which the integrand of the boundary term is [21] As the solution converges, this difference goes to zero, and the Rankine-Hugoniot condition is satisfied. Furthermore, so that when the approximation converges, the analytical inflow boundary condition, U R = a L a R U L is approached, as required, c.f. (14). With (76), the interface contribution for the scalar problem is Factoring the quadratic, The quadraticQ(η; a L a R ) is concave up and has a minimum when η * = a L /a R , sincẽ When η * = a L /a R , the Rankine-Hugoniot condition is satisfied by the states on either side. The value of that minimum isQ(η * ; a L a R ) = 1 − a L a R . It then follows that the contribution to the energy in the numerical approximation matches that of the PDE, (21), plus a dissipation term dependent on how much the Rankine-Hugoniot condition is not satisfied by the approximate solution. If we define β = a L /a R , and note that the minimum value ofQ is 1 − β, we can separate out that term giving Re-writing the interface contribution in the final form of (82) will be a key step in showing inferred stability of the approximation for the more complex case of a system of equations.
When we substitute for η and β, Let us compare: In the continuous case, we have (21), with whereas discretely, Thus, according to the definition of stability, Definition 2.1, the DGSEM approximation of the scalar problem with the upwind numerical flux has inferred stability. (84) and (85) shows explicitly what is interpreted as stability. The first term in (85) can be positive or negative depending on a L /a R , but matches that of the PDE, (84). The approximation is therefore directly stable if a L /a R ≤ 1, just like the PDE. The second term is always non-positive and represents dissipation of the energy by the approximation.

Remark 4.2
For the scalar problem it is straightforward, as for the continuous problem, to show energy boundedness in a discounted norm by scaling the downwind domain contributions before summing over the elements. When the global sum (in this case, over two elements) is formed, the parenthetical term in the second line of (80) becomes Like the original quantity,Q N in (80),Q N ,α c is concave up, with minimum at the same point, η * , with minimum valueQ Since one can always show bounded energy in the new discounted norm by choosing α c to match the analytical value for any (positive) wavespeeds, the condition (85) infers stability. The amount of numerical dissipation in that norm depends on the particular choice of α c , however.

Stability for the One-Dimensional Symmetric System
We now parallel Sect. 2.2 and extend the analysis to a symmetric PDE system in one space dimension. For the system, the DG approximation has the interface contribution The upwind numerical flux is now with the equalities between the left and right representations arising by virtue of the Rankine-Hugoniot condition. The key observation is that But U T = PW T = W T P T and for the symmetric system P T P = I, so Now, and Therefore, Looking at the second jump term in (89), so Therefore, forming Q N and gathering right and left going wave components, Terms cancel, leaving Following (82), we now add and subtract terms to match the PDE form, which is to write Now, let Then To show that the approximation is stable according to Definition 2.1, then, we just need to show that R ± ≥ 0, since the other terms match those of the PDE. To do so, letW ± = Λ± W ± . Then Similarly, Thus, the interface contribution matches that of the PDE plus an additional dissipation and has inferred stability, satisfying Definition 2.1 with

Stability of the General Problem
As in the continuous problem, we use the analysis of the one-dimensional problem to imply stability of the multidimensional one. As before, replace U ← U s and A ← → A s · → n. Then Q N is given by (106), with the eigenvalues (and eigenvectors to construct the characteristic variables) coming from → A s · → n. Therefore the approximation to the general multidimensional problem is stable according to Definition 2.1.

Remark 4.3
The key features of the stability analysis are the use of summation by parts, and a stable implementation of the boundary terms. As such, the analysis extends to other methods that have the summation by parts property and allow discontinuities at subdomain interfaces, such as summation by parts finite difference techniques, e.g. as used in [15,17,18].

Example
As an example, we consider the scattering of a plane wave off a plane material interface, approximating the system of equations (1) with the state vector and coefficient matrices (3) reduced to two space dimensions. The problem has exact incident, transmitted and reflected plane wave solutions of the form where ψ is a given wavefunction, a is the amplitude, → k is the wavevector, ω is the frequency. For the incident wavevector the reflected and transmitted wavevectors are with amplitudes a r a i = where For the wavefunction, we choose the Gaussian with σ 2 = −(MT ) 2 /(4 ln(10 −4 )), M = 4 and period T = 2π/ω. We compute the problem on the square domain [−5, 5] 2 with 400 square elements and the material interface at x = 0. The solution parameters are provided in Table 1.
The results are shown in Figs. 4 and 5. Figure 4 shows the contours of the p component of the solution at time t = 5.0, which is near the time of the maximum L 2 energy, computed with sixth order polynomials. Clearly seen is the jump discontinuity at the interface. The L 2 energy is plotted as a function of time in Fig. 5, for polynomial degrees N = 2, 3 and 6. Although the L 2 energy initially grows, it reaches a maximum around time t = 4.5. Figure 5 shows that the computed energy converges from below to the exact as the polynomial order is increased. In fact, it converges exponentially with polynomial degree, as expected [2] for a spectral element method. Also, as expected due to the additional dissipation at physical, smooth and discontinuous interfaces, the computed energies fall below the exact curve and are worst for low order approximations.

Conclusions
We have shown that the interface treatment of the discontinuous Galerkin spectral element method with the upwind numerical flux is stable for hyperbolic systems with discontinuous coefficient matrices when the eigenvectors are preserved across the interface. Examples include systems like Maxwell's equations, or acoustic and elastic wave equations. The new feature of our approach was to show that the discrete L 2 norm of the approximate solution grows no faster than the same norm of the continuous solution. By matching the L 2 norm, we avoid having to find the precise conditions for a discounted norm in which the energy is bounded in terms of the initial data (for homogenous and dissipative boundary conditions). The numerical flux only weakly enforces the inflow boundary condition and the Rankine-Hugoniot condition. Viewing stability in terms of the L 2 norm shows that the dissipation introduced by the upwind numerical flux depends on the amount by which the approximate solution fails to satisfy the Rankine-Hugoniot condition.