Convergence of Summation-by-Parts Finite Difference Methods for the Wave Equation

When using a finite difference method to solve a time dependent partial differential equation, the truncation error is often larger at a few grid points near a boundary or grid interface than in the interior. In computations, the observed convergence rate is often higher than the order of the large truncation error. In this paper, we develop techniques for analyzing this phenomenon, and particularly consider the second order wave equation. The equation is discretized by a finite difference operator satisfying a summation by parts property, and the boundary and grid interface conditions are imposed weakly by the simultaneous approximation term method. It is well-known that if the semi-discretized wave equation satisfies the determinant condition, that is the boundary system in Laplace space is nonsingular for all Re(s)≥0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(s)\ge 0$$\end{document}, two orders are gained from the large truncation error localized at a few grid points. By performing a normal mode analysis, we show that many common discretizations do not satisfy the determinant condition at s=0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s=0$$\end{document}. We then carefully analyze the error equation to determine the gain in the convergence rate. The result shows that stability does not automatically imply a gain of two orders in the convergence rate. The precise gain can be lower than, equal to or higher than two orders, depending on the boundary condition and numerical boundary treatment. The accuracy analysis is verified by numerical experiments, and very good agreement is obtained.


Introduction
In many physical problems arising in for example acoustics, seismology, and electromagnetism, the governing equations can be formulated as systems of second order time dependent hyperbolic partial differential equations (PDE). For wave propagation problems with smooth solutions, it has been shown that high order accurate time marching methods as well as high order spatial discretizations solve these problems more efficiently than low order methods [7,9,10]. Examples of high order accurate methods to discretize the wave equation include the discontinuous Galerkin method [5] and the spectral method [23]. In this paper, we in particular consider high order finite difference methods.
High order finite difference methods have been widely used for solving wave propagation problems. One major difficulty with high order spatial discretizations is the numerical treatment of boundary conditions and grid interface conditions. To achieve both stability and high accuracy, one candidate is the summation-by-parts simultaneous-approximation-term (SBP-SAT) finite difference method [4,22]. SBP finite difference operators were originally constructed in [11] for first derivatives. In this paper, we use the SBP finite difference operators in [15] to approximate second derivatives. Boundary conditions and grid interface conditions are imposed weakly by the SAT technique [2,3,13,14]. In the SBP-SAT framework, the energy method is often used to derive an energy estimate that guarantees stability.
A commonly used measure for accuracy is convergence rate, typically in L 2 norm. The convergence rate indicates how fast the error in the numerical solution approaches zero as the mesh size goes to zero. The error in the numerical solution is caused by truncation errors. The truncation error near a boundary is often larger than in the interior of the computational domain, but the large boundary truncation error is localized at a fixed number of grid points. As a consequence, its effect on the convergence rate may be weakened. This is often called gain in convergence. A rule of thumb says m orders are gained in convergence for a PDE with m th order spatial derivative, and such a gain is termed as optimal. The analysis of gain in convergence for different PDEs has been a long-standing research topic.
It is well-known that by directly applying the energy method to the error equation, 1/2 order is gained in the convergence rate compared with the largest truncation error. A noticeable exception is that in [1], it is proven that 3/2 orders are gained in the convergence rate for parabolic problems by a careful derivation of the energy estimate. However, numerical experiments in [1] show a gain of two orders in convergence. This indicates that although an energy estimate gives an upper bound of the error, it is often not sharp.
A more powerful tool for stability and accuracy analysis is normal mode analysis, which is used for example in [6] to prove that under reasonable conditions one order is gained in the convergence rate for first order hyperbolic systems. The idea is based on Laplace transforming the error equation in time, which leads to a system of linear equations referred to as the boundary system. The optimal gain straightforwardly follows if the boundary system is nonsingular for all Re(s) ≥ 0, where s is the dual variable of time in Laplace space. For such cases we use the same terminology as in [8, pp. 292] and say that the boundary system satisfies the determinant condition.
In [21] the concept of pointwise stability is put forward as a condition implying the determinant condition, and therefore leading to the optimal gain. In [22] a sufficient condition of pointwise stability for an initial-boundary-value problem with a first derivative in time is discussed. However, we find that for the wave equation, the determinant condition does not follow from pointwise stability as defined in [21], and such an example is presented in "Appendix 2".
In fact, there are many examples of discretized problems that violate the determinant condition at points along the imaginary axis, even though the discretization is stable by an energy estimate. In [17], the Schrödinger equation with a grid interface is considered and is shown to be of this type. Analysis in Laplace space is performed and yields sharper error estimates than the 1/2 order gain obtained by applying the energy method to the error equation in physical space.
In this paper, we consider the wave equation in the second order form. We find that also for this problem, the determinant condition is violated in many interesting settings even though there is an energy estimate. In particular, we consider problems with Dirichlet boundary conditions, Neumann boundary conditions and a problem with a grid interface, which are referred to as the Dirichlet problem, the Neumann problem and the interface problem, respectively. If the determinant condition is satisfied, we get the optimal convergence rate. If the determinant condition is not satisfied on the imaginary axis, the results in [8] are not valid. We then carefully derive an estimate for the solution of the boundary system to see how much is gained in the convergence rate. In addition, we have found by a careful analysis of the boundary system that in certain cases the gain in convergence is higher than the optimal gain, which we refer to as super-convergence.
The main contributions in this paper are: for the second order wave equation (1) an energy estimate does not imply an optimal convergence rate; (2) the determinant condition is not necessary for an optimal convergence rate; (3) if there is an energy estimate but the determinant condition is not satisfied, there can be an optimal gain of order 2, or a nonoptimal gain of order 1, or only order 1/2; (4) in certain cases it is possible to obtain a super-convergence, i.e. a gain of 2.5 orders.
The rest of the paper is organized as follows. We start in Sect. 2 with the SBP-SAT method applied to the one dimensional wave equation with Dirichlet boundary conditions. We apply the energy method and normal mode analysis, and derive results that correspond to the second, fourth and sixth order schemes. We then discuss the Neumann problem in Sect. 3. In Sect. 4, we consider the one dimensional wave equation on a grid with a grid interface. In Sect. 5, we perform numerical experiments. The computational convergence results support the accuracy analysis. In the end, we conclude the work in Sect. 6.

The One Dimensional Wave Equation with Dirichlet Doundary Conditions
To describe the properties of the equation and numerical methods, we need the following definitions. Let w 1 (x) and w 2 (x) be real-valued functions in L 2 [0, ∞). An inner product is defined as (w 1 , w 2 ) = ∞ 0 w 1 w 2 dx with a corresponding norm w 1 2 = (w 1 , w 1 ), and w 1 is said to be in L 2 if w 1 < ∞.
The second order wave equation on a half-line in one space dimension takes the form We consider the equation up to some finite time t f . For boundary conditions, we use either the Dirichlet boundary condition: or the Neumann boundary condition The forcing function F, the initial data G,Ḡ and the boundary data M are compatible smooth functions with compact support. As a consequence, the true solution is smooth. For the problem (1) in a bounded domain, there is one boundary condition at each boundary. For the half line problem in consideration, the right boundary condition is substituted by requiring that the L 2 norm of the solution is bounded, i.e. U (·, t) < ∞. If the corresponding half line and Cauchy problems are well-posed, then the original initial-boundary-value problem is well-posed [8, pp. 256]. The problem (1) with (2) or (3) is well-posed if the solution can be bounded in terms of the data. Multiplying Eq. (1) by U t and integrating in space, with the homogeneous Dirichlet or Neumann boundary condition, it follows from the integration-by-parts principle where |||U ||| 2 = U t 2 + U x 2 is the continuous energy, a semi-norm of the solution U . The standard energy estimate follows from an integration in time of (4), Therefore, problem (1) is well-posed with either (2) or (3).

Remark 1
We are interested in the solution U , but the standard energy analysis above only shows that a semi-norm of U is bounded. To include U itself in the estimate, we use the relation An integration in time, together with (5), gives Therefore, U is bounded by the data in any bounded time interval t ∈ [0, t f ], and is in the space of L 2 .
To discretize the equation in space, we introduce an equidistant grid x i = ih, i = 0, 1, 2, · · · , and a grid function u i (t) ≈ U (x i , t). Furthermore, let u(t) = [u 0 (t), u 1 (t), · · · ] T where T denotes the transpose. We also define an inner product and norm for the grid functions a and b in R as (a, b) H = a T Hb and a 2 H = a T Ha, respectively, where H is a symmetric positive definite operator on the space of grid functions. The norm · H is equivalent to the standard discrete L 2 norm, and they are the same if H = h I with an identity operator I . An SBP operator has the standard central finite difference stencil in the interior and a special one-sided stencil near the boundary. The boundary closure is chosen such that the SBP operator mimics integration-by-parts via its associated inner product. The SBP operator approximating the second derivative D ≈ ∂ 2 /∂ x 2 has a structure where H is diagonal and positive definite, A is symmetric positive semi-definite and B = diag(−1, 0, 0, · · · ). S is a one sided approximation of the first derivative at the boundary.
In this paper, we use the diagonal norm SBP operators constructed in [15]. Although these operators are termed as 2 p th order accurate, the truncation error of D is O(h 2 p ) in the interior but only O(h p ) near the boundary. The truncation error of S is O(h p+1 ). Such operators are constructed for p = 1, 2, 3, 4 in [15], and p = 5 in [12].
An SBP operator itself does not impose any boundary condition. It is important that boundary conditions are imposed in a way such that an energy estimate can be derived to ensure stability. Such numerical boundary treatments for the wave equation include the projection method [16] and the ghost point approach [19], which impose boundary conditions strongly. In this paper, we instead use a weak enforcement technique, the SAT method [3], since it is easy to derive an energy estimate even in higher dimensions. The semi-discretization of (1) and (2) with the homogeneous boundary data reads: where T , E 0 = e 0 e T 0 and f (t) is the restriction of F(x, t) to the grid. On the right hand side, the first term Du is an SBP approximation of U x x , while the second term P D u imposes weakly the boundary condition U (0, t) = 0. It acts as a penalty term dragging the numerical solution at the boundary towards zero so that the boundary condition is simultaneously approximated, but the boundary condition is in general not satisfied exactly. The penalty parameter τ is to be determined so that an energy estimate of the discrete solution exists, which ensures stability.

Stability
In [2,13,14], it is shown that the operator A in (7) can be written as where α 2 p > 0 is a constant independent of h andÃ is symmetric positive semi-definite. This is often called the borrowing trick. The values of α 2 p are computed for p = 1, 2 and 3 in [13]. We use the same technique to compute α 2 p for p = 4 and 5, and list the results in Table 1.
To derive an energy estimate, we introduce the energy where τ ≥ 1/α 2 p . Multiplying Eq. (8) by u T t H from the left, it follows from (7) and (9) The discrete energy estimate is an analogue of the continuous energy estimate (5). We conclude that the semidiscretization (8) is energy stable if τ ≥ 1/α 2 p . The term u H can be included in the energy estimate in the same way as in Remark 1.

Accuracy Analysis by the Energy Method
The error equation for the pointwise error where ζ(t) h < ∞, T 2 p is the truncation error and 2 p is the accuracy order of the SBP operator D. The solution has compact support during the entire time interval in consideration, so does the truncation error. Therefore, T In the analysis, we only consider the leading term of the truncation error. We introduce the interior truncation error T 2 p,I , and the boundary truncation error T 2 p,B such that T 2 p = h 2 p T 2 p,I + h p T 2 p,B , where T 2 p,I and T 2 p,B are independent of h, but depend on the derivatives of U . Since the number of grid points with the large truncation error O(h p ) is finite and independent of h, we have ).
For 2 p = 2, 4, 6, 8, 10 and small h, the first term is much smaller than the second one. Thus, we have T 2 p h ≤ K h p+1/2 . Here and in the rest of the paper, we use the capital letter K in the estimate to denote some constant independent of the grid spacing. The constant K can be made precise by Taylor expansions, but we do not distinguish them from one estimate to another for a sake of simplified notations.
By applying the energy method to the error equation (11), we obtain for the discrete energy defined by (10) |||ζ This means that by the energy method we can only prove a gain in accuracy of 1/2 order. Therefore, the convergence rate is at least p + 1/2 if the numerical scheme is stable, that is if the penalty parameter τ ≥ 1/α 2 p .

Normal Mode Analysis for the Boundary Truncation Error
To derive a sharper estimate, we partition the pointwise error into two parts, the interior error I and the boundary error , such that ζ = I + . I is the error due to the interior truncation error, and is the error due to the boundary truncation error. I can be estimated by the energy method, yielding As mentioned above, I H can be bounded similarly. We perform a normal mode analysis to analyze for the second, fourth and sixth order SBP-SAT schemes. The boundary error equation is where h < ∞. In the analysis of , we take the Laplace transform in time of (14), for Re(s) > 0, whereˆdenotes Laplace transform. After multiplying by h 2 on both sides of (15), we obtain with the notations = sh Note that the coefficients in the first two terms in the right hand side of (16) are h-independent, and thatT 2 p,B is nonzero only at a few points near the boundary.
There are essentially two steps in the normal mode analysis. In the first step, by a detailed analysis of the error equation (16), we derive an estimate of the solution to (16) in terms of the forcing. We use Parseval's relation in the second step to derive an estimate for the error in the physical space, which involves an inverse Laplace transform. As will be seen later, the integration is performed along the vertical line Re(s) = ηh with η > 0 a fixed constant independent of h. It is therefore important to derive a sharp error estimate in Laplace space when Re(s) goes to zero. The final estimate in the physical space is in the form where K depends only on η, the final time t f and the derivatives of the true solution U . Convergence rate is an asymptotic property, and we need the solution to be smooth and well resolved. By the smoothness assumption |T 2 p (η + iξ)| decreases fast for large |ξ |, making contributions from |ξ | >K insignificant for some constantK . In the analysis we will only consider s = η + iξ , with η positive and |ξ | <K . Correspondingly we only consider |s| 1, Re(s) = O(h). In particular it suffices to considers = 0 when checking the determinant condition. We formalize this discussion by making the following assumption.
Assumption 1 There is a constantK < ∞ such that contributions from s = η +iξ , |ξ | >K to the estimate (17) do not influence the rate q.
We summarize the convergence result for the Dirichlet problem in Theorem 1.  Table 2.
Remark 2 The penalty parameter τ has an effect on the convergence rate for a stable scheme. More precisely, with the same SBP operator, the convergence rate is always higher with τ > 1/α 2 p than with τ = 1/α 2 p . We therefore should always choose τ > 1/α 2 p in computations, but bearing in mind that a moderate value of τ is appropriate because a very large τ has a negative impact on the CFL condition. The convergence rates are optimal for the second and fourth order schemes, while superconvergence is obtained with the sixth order scheme. We also comment that the overall convergence rate for the second order scheme is dominated by the interior truncation error.

Solution to the Error Equation
To begin with, we note that sufficiently far away from the boundary, the penalty terms and the boundary truncation error in (16) are not present. The coefficients in D correspond to standard central finite difference schemes: The corresponding characteristic equations are It is easy to verify by the von Neumann analysis that the interior numerical scheme is stable when applied to the corresponding periodic problem. From Lemma 12.1.3 in [8, pp. 379], it is straightforward to prove that there is no root of (19) with |κ| = 1 for Re(s) > 0. We call a root κ an admissible root if |κ| < 1. Since the problem in consideration is a half line problem, any κ with |κ| > 1 is not admissible because it results in an unbounded solution.
We will need the following specifics for the roots.
2 p = 4: Eq. (19b) has four roots: We find by Taylor expansion ats = 0 that Thus, the admissible roots are There is no closed form of the solution to a sixth order equation like (19c). However, (19c) withs = 0 can be factorized to a second order polynomial multiplied by a fourth order polynomial, thus an analytical solution can be obtained. We then perform a perturbation analysis to the six roots, and find analytical expressions for their dependence oñ s in a neighbourhood ofs = 0. The three admissible roots are given by (20c). We show the numerical values with four digits since the analytical expressions are very lengthy.
The pointwise error away from the boundary can be expressed as where the coefficients σ are determined by the boundary closure. The error equation (16) corresponding to a few grid points near the boundary also determines the pointwise errors that are not determined by 21. The general form for the L 2 norm ofˆ can be written as where d = 1, 1, 2 for 2p = 2, 4, 6, respectively. Note that by Lemma 1 we have that in the second term the factor 1 1−|κ 1 | 2 cannot be bounded independent of h whens = O(h).
to the leading order.
Proof Lets = x + iy where x, y are real numbers. Then |x| ≥ ηh.
The desired estimate follows.

By Lemma 2, (22) becomes
In the following sections, we analyze the error equation (16) corresponding to the grid points near the boundary, and derive bounds forˆ i and σ m in (23). The final estimate in the physical space of the form (17) follows by Parseval's relation. To keep a balance between clarity and paper length, we give a very detailed analysis for the second order scheme, and only show the main steps of the proof for the fourth and sixth order schemes.

Proof of Theorem 1 for the Second Order Scheme
The first three rows of (16) are affected by the boundary closure. They arẽ By Taylor expansion, it is straightforward to deriveT 2,B 0 = −Û x x x (0, s) to the leading order. We obtain the boundary system by rewriting (24) in a matrix-vector multiplication form If the determinant condition is satisfied, i.e. C 2D is nonsingular for Re(s) ≥ 0, then |Σ 2D | is bounded by a constant multiplying h 3 and a gain of two orders from the boundary truncation error follows. In other words, we obtain the optimal convergence rate if the determinant condition is satisfied. Below we demonstrate that by analyzing the components of Σ 2D , it is in certain cases possible to obtain an even higher order gain from the boundary truncation error, which is referred to as super-convergence.
In Lemma 2, we use Re(s) = ηh to estimate 1/(1 − |κ 1 | 2 ), with η a constant independent of h. When analyzing the boundary system, we need to use the sames = ηh because that is where the inverse Laplace transform is performed later. We write the solution to (25) as for Re(s) = ηh. Since η is a constant and h is arbitrarily small, it is important to check how Σ 2D (s, τ ) behaves ass approaches zero. We start by setting s = 0, and find that when τ = 2.5 the determinant condition is satisfied, i.e. C 2D (0, τ ) is nonsingular; but the determinant condition is not satisfied when τ = 2.5. The stability condition τ ≥ 1/α 2 = 2.5 motivates us to consider these two cases separately.
In the case τ > 2.5, the determinant condition is satisfied and we can expect at least an optimal gain of 2 orders in convergence. A perturbation analysis of Σ 2D ats = 0 shows that which together with (23) leads to Remember that |s| 1 corresponds to |s| <K for some constantK and h 1. By Assumption 1 and Parseval's relation, we have We have derived the estimate forˆ in the vicinity ofs = 0, but Assumption 1 allows us to use this estimate when integrating ξ from −∞ to ∞. By arguing that future cannot affect past [8, pp. 294], we can replace the upper limit of the integrals on both sides by a finite time t f . Since we obtain Thus, the boundary error is O(h 3.5 ), which is 2.5 orders higher than the boundary truncation error. In practical computations, this super-convergence is not seen because the dominating source of error is the interior error O(h 2 ) given by (13). We conclude that the overall convergence rate is 2, and this proves that q B = 3.5 and q = 2 in Theorem 1.

Remark 3
By only checking the determinant condition, we can prove the optimal convergence rate, i.e. a gain of two orders from the boundary truncation error. The above analysis demonstrates that it is important to analyze the components of the solution to the boundary system. Super-convergence is obtained when σ 1 is much smaller than the other components of the solution to the boundary system when h is sufficiently small. In this case, the error in the solution is larger at a few grid points near the boundary than in the interior.
In the case τ = 2.5, the last two components of Σ 2D are infinite whens = 0 since C 2D (0, τ ) is singular and the determinant condition is not satisfied. We again perform a perturbation analysis to the solution Σ 2D and obtain This, in the same way as the case τ > 2.5, leads to a gain of 0.5 order from the boundary truncation error and the overall convergence rate p + 1/2 = 1.5, which is the same as predicted by the energy estimate.

Proof of Theorem 1 for the Fourth Order Scheme
The first four rows of (16) are affected by the penalty terms, which can be written as the boundary system As in the second order case in Sect. 2.3.2, we again consider Re(s) = ηh and write the solution to the boundary system (27) The matrix C 4D (0, τ ) is presented in "Appendix 1". In the case τ > 1/α 4 , C 4D (0, τ ) is nonsingular and the determinant condition is satisfied. All the four components of Σ 4D are order O(h 4 ) with σ 1 independent of τ . This gives the estimate to the leading order By Parseval's relation and by arguing that future cannot affect past, we obtain Thus, the boundary error is O(h 4 ). In this case, the interior error is also O(h 4 ) given by (13). Therefore, the convergence rate is 4. This proves that q B = q = 4 in Theorem 1.

Proof of Theorem 1 for the Sixth Order Scheme
Similar to the previous two cases, we consider the six-by-six boundary system and analyze its solution Σ 6D = [σ 1 , σ 2 , σ 3 ,ˆ 0 ,ˆ 1 ,ˆ 2 ] T : wherê The matrix C 6D (0, τ ) is presented in "Appendix 1". It is singular if τ = 1/α 6 , and nonsingular otherwise. When τ > 1/α 6 , we write the solution to (29) as Σ 6D (s, τ ) = h 5 C −1  (23) we obtain that the boundary error is O(h 5.5 ). In this case, the interior error is O(h 6 ) given by (13). Thus, the convergence rate is 5.5, which is a super-convergence and is half order higher than the optimal convergence. This proves that q B = q = 5.5 in Theorem 1.

The One Dimensional Wave Equation with Neumann Boundary Conditions
We consider Eq. (1) with the Neumann boundary condition (3), and use the same assumption of the data as for the Dirichlet problem. As will be seen later, the determinant condition is not satisfied. Our analysis gives sharp error estimates for such problems as well.

Stability
To discretize the equation in space, we use the grid and grid functions introduced for the Dirichlet problem. The semi-discretized equation of the Neumann problem is On the right hand side of (31), the first term approximates U x x and the second term imposes weakly the boundary condition U x (0, t) = 0. Since we consider a half line problem, the SBP operator D is in the form D = H −1 (−A − E 0 S). As a consequence, the semidiscretization (31) can be written as

Multiplying Equation (32) by u T t H from the left, we obtain
Here again, u H can be included in the energy estimate as in Remark 1.

Accuracy
With the pointwise error defined as before, the error equation is where ζ h < ∞, T 2 p,N is the truncation error and 2 p is the accuracy order of the SBP operator. In the same way as for the Dirichlet problem, by applying the energy method to the error equation (34), we obtain an estimate This means that the convergence rate of (32) is at least p + 1/2. In the following, we use normal mode analysis to derive a sharper bound of the error, which agrees well with the results from numerical experiments. As in the Dirichlet problem, we partition the error ζ into two parts, the error I due to the interior truncation error and the error due to the boundary truncation error. I can be estimated by the energy method, yielding I N ,h ≤ K h 2 p . The boundary error equation is Since the operator A in (36) is only symmetric positive semi-definite, we expect that the determinant condition is not satisfied for this problem. The characteristic equations of the interior error equations are the same as for the Dirichlet problem (19), and Lemma 1 and 2 are also applicable to the Neumann problem. Below we analyze the second, fourth and sixth order accurate cases.

Second Order Accurate Scheme
Only the first row of the error equation is affected by the boundary closure. In this case, the boundary system is the scalar equation with and the error is in the form (21a) with j starting at 0. Clearly the determinant condition is violated ats = 0, but we straightforwardly obtain In the same way as for the Dirichlet problem, an estimate in physical space of the type (17) follows with q = 2. Since the interior error is also O(h 2 ), the convergence rate is 2.

Remark 4
In this case, the gain from the boundary truncation error is only one order, and does not follow the p + 2 optimal gain. The interior error O(h 2 ) restricts the overall convergence rate to second order. Therefore, a gain of more than one order can normally not be observed.
In the numerical experiments, we design a stable numerical scheme to verify that the gain in convergence is indeed only one order.

Fourth and Sixth Order Accurate Schemes
We follow the above approach to derive the four by four boundary system for the fourth order scheme, and the sixth by sixth boundary system   Table 3.

The One Dimensional Wave Equation with a Grid Interface
In this section, we consider another example that does not satisfy the determinant condition. It is the Cauchy problem for the second order wave equation in one space dimension where the forcing function F , the initial data G andḠ are compatible smooth functions with compact support. It is straightforward to derive the energy estimate of the form (5) for Eq. (38). We solve the equation on a grid with a grid interface at x = 0. With the assumption that the true solution is smooth, it is natural to impose two interface conditions at x = 0: continuity of the solution and continuity of the first derivative of the solution.
We introduce the grid on the left x − j = − jh L , j = 0, 1, 2, · · · , and the grid on the right x j = jh R , j = 0, 1, 2, · · · , with the grid size h L and h R , respectively. The grid functions Both u L 0 = e T 0L u L and u R 0 = e T 0R u R approximate U (0, t). Both (S L u L ) 0 = e T 0L S L u L and (S R u R ) 0 = e T 0R S R u R approximate U x (0, t). To simplify notation, we define {u} = u L 0 − u R 0 and {Su} = (S L u L ) 0 − (S R u R ) 0 . The semi-discretized equation reads where f L and f R are the restriction of F(x, t) on the grid. The penalty terms P I L {u} and P I R {u} impose the interface conditions weakly and have the form The penalty parameters τ is chosen so that the semi-discretization is stable. By the energy method, the numerical scheme is stable if The stability proof is found in [13,Lemma 2.4,pp. 215]. If h L = h R := h, thenτ 2 p = 1/(2α 2 p h). By the energy estimate, the convergence rate is at least p+1/2 if the semi-discretization is stable. In order to derive a sharper estimate, we follow the same approach as in the previous sections. The interior truncation error results in an error O(h 2 p ) in the solution. In the remaining part of this section we will only consider the effect of the interface truncation error. Denote , j = 0, 1, 2, · · · . The error equation in Laplace space reads wheres L = sh L ,s R = sh R . Here the interface truncation errorsT 2 p,L/R,B are nonzero at only a few grid points near the interface. The errors are estimated for the second and fourth order schemes by the normal mode analysis, with the results summarized in the following theorem.

Theorem 3 For the second and fourth order stable SBP-SAT approximation (39) of the wave equation with a grid interface, the numerical solution converges to the true solution at rate q.
When τ =τ 2 p in (40), the rate q = p + 1/2 is given by the energy estimate. When τ >τ 2 p , the rate q is listed in Table 4.

Remark 5
In practical computations, we should always choose τ >τ 2 p so that the optimal convergence rate is obtained. Note that in this case, the gain in convergence depends on whether the grid spacing h L is equal to h R as shown in Table 4. The overall rate q = min(2 p, q B ), however, is not affected by the relation of h L and h R .
We only need to analyze the values q B when τ >τ 2 p . We use the notation h = h L = rh R where r is a fixed mesh size ratio. The analysis follows in the same way as before. We construct the characteristic equation on each side and solve them to obtain the general solution, which are given by (19) and (20), respectively. The boundary system is formulated by substituting the general solution to the error equation. The accuracy order is determined by how the solution of this system behaves with respect to h whens = O(h).

Proof of Theorem 3 for the Second Order Scheme
The last three rows of (41a) and the first three rows of (41b) are affected by the penalty terms, which lead to the six by six boundary system The determinant condition is not satisfied because C 2I (0, τ ) shown in "Appendix 1" is singular. We find that when r = 1, = 2. In both cases, the overall convergence rate is q = 2.

Proof of Theorem 3 for the Fourth Order Scheme
In this case, the eight by eight boundary system takes the form  (0, s).
The boundary system is also singular ats = 0, and the determinant condition is not satisfied. The matrix C 4I (0, τ ) is presented in "Appendix 1". When r = 1,

Numerical Experiments
In this section, we perform numerical experiments to verify the accuracy analysis. We investigate how the accuracy of the numerical solution is affected by the large truncation error localized near boundaries and grid interfaces. In the analysis, we use the half line problem. However, in the numerical experiments, we use a bounded domain and impose boundary conditions on all boundaries weakly. For the integration in time, we write the semi-discretized equation as a first order system in time and use the classical fourth order Runge-Kutta method. The step size Δt in time is chosen so that the temporal error has a negligible impact on the spatial convergence result. The value of Δt is given in each numerical experiment.
The L 2 error at a given time point are computed as the norm of the difference between the exact solution u ex and the numerical solution u h according to where d is the spatial dimension. The convergence rate is computed by

The One Dimensional Wave Equation
We choose the analytical solution and use it to impose the Dirichlet boundary conditions at x = 0 and x = 1. The step size in time Δt = 0.1h is much smaller than the step size restricted by the stability condition. We use it because then the temporal error has a negligible impact on the spatial convergence rate. The analytical solution and its derivatives do not vanish at the boundaries at the final time t = 2. The errors in discrete L 2 norm are plotted versus the number of grid points in Fig. 1(a), and the convergence rates in L 2 norm are shown in the same figure.
With a 2 p th (2 p = 2, 4, 6) order accurate method, the convergence rate is p + 1/2 in L 2 norm if the penalty parameter τ equals its limit. With a larger penalty parameter (increased by 20%), the accuracy of the numerical solution is improved significantly. The observed  Table 2, Theorem 1.
The super-convergence of the sixth order accurate scheme is also observed when solving the Schrödinger equation [17,18].
To test the Neumann problem, we again use (43) as the analytical solution and impose the Neumann boundary condition at two boundaries x = 0 and x = 1. In Fig. 1b we show the L 2 errors and the convergence rates in L 2 norm, which corroborate the accuracy analysis very well: for the second, fourth and sixth order accurate schemes, the observed convergence rates 2.01, 4.04 and 5.54 are quite close to the theoretical values 2, 4 and 5.5.
For the sixth order accurate scheme with τ > 1/α 6 , the components of the solution to the boundary system are of different magnitude: σ 1 ∼ O(h 7 ) while the other components are of order O(h 5 ). Therefore, the boundary truncation error only leads to O(h 7 ) interior error. The large error O(h 5 ) in the solution is located only at a few points near the boundary, which can be seen in Fig. 2a. The same is true for the Neumann problem with the sixth order accurate scheme, and the error plot in Fig. 2b also shows that the error in the solution is localized.
Next, we test how the accuracy and convergence are affected by the large truncation error localized near a grid interface. We choose the same analytical solution (43). The grid interface is located at x = 0.5 in the computational domain Ω = [0, 1], and the grid size ratio is 2:1. The step size in time is Δt = 0.1h, where h is the smaller grid size. Note that the analytical solution does not vanish at the grid interface. We use the SAT technique to impose the outer boundary condition weakly and choose the boundary penalty parameters strictly larger than their limits. The interface conditions are also imposed by the SAT technique. We use different interface penalty parameters to see how they affect the accuracy and convergence. The results are shown in Fig. 3a.
According to the accuracy analysis in Sect. 4, for second and fourth order accurate methods with τ = τ 2 p , the expected convergence rates are 1.5 and 2.5 in L 2 norm. This is clearly observed in Fig. 3a. With a larger penalty parameter, a much better convergence result is obtained. For the second and fourth order methods, we get the optimal (second and fourth, respectively) order of convergence in L 2 norm. For the sixth order accurate scheme, the convergence rate is about 3.4 with τ = τ 6 , which is in line with the p + 1/2 convergence rate. With a larger penalty parameter, the convergence rate is about 5.4. Here we again observe a super-convergence phenomenon.

Non-Optimal Convergence
We have shown several cases where the gain in convergence does not follow the optimal p + 2 rule. For example, the sixth order accurate SBP-SAT scheme has a super-convergence property, which is both proved theoretically and observed in the numerical experiments. There are also cases where the accuracy analysis predicts a lower or higher order gain. Examples include the second and fourth order schemes for the Neumann problem, where analysis indicates a gain of 1 and 2.5 orders, respectively. In these cases, the interior error also plays an important role in the overall convergence rate. It is therefore unclear from the computation what the precise gain is in computations. We design the following two experiments to confirm that our accuracy analysis indeed gives sharp error estimates. We use (43) as the analytical solution and impose the Neumann boundary condition at two boundaries x = 0 and x = 1. With this analytical solution, the equation does not have a forcing term, and f = 0 in (32). We modify the numerical scheme (32) to In the first experiment, we use the second order accurate scheme and choose the first component of F to be (10π) 3 and all the other components zero. We use the number (10π) 3 because the truncation error involves the third derivative of the true solution with the factor (10π) 3 . The boundary truncation error of (44) has the same structure as in the original scheme, but increases to O(1), i.e. order 0. Note that energy stability still holds for (44). The L 2 error versus the number of grid points is plotted in Fig. 3b. Clearly, the convergence rate is 1, and is a gain of 1 order from the boundary truncation error. This corroborates the accuracy analysis very well.
In the second experiment, we use the fourth order accurate scheme and let the first four components of F be h(10π) 4 [43/204, −1/12, 5/516, 11/588] T , and the other components zero. The boundary truncation error of (44) has the same structure as in the original scheme, but increases to O(h). We plot the L 2 error versus the number of grid points in Fig. 3b, and observe a 3.5 orders convergence rate, i.e. a gain of 2.5 orders from the boundary truncation error. This again demonstrates that the error estimate is sharp.

Conclusion
For the second order wave equation a stable numerical scheme does not automatically satisfy the determinant condition, nor does it imply an optimal gain in convergence rate. We have considered stable SBP-SAT finite difference schemes for the Dirichlet, Neumann and interface problems. In all these cases the boundary/interface truncation error is larger than the interior truncation error. We find that the schemes satisfy the determinant condition only for the Dirichlet problem when the penalty parameter is greater than its limit value, and a gain of two orders in convergence follows by the standard analysis. By a careful analysis of the solution to the boundary system, we prove that the gain in convergence for some schemes is in fact 2.5 orders, i.e. half an order higher than the optimal gain.
For all the other cases the determinant condition is not satisfied, but there is nonetheless a gain in convergence of 0.5, 1, 2, or even 2.5 orders. In those cases, only a detailed analysis of the boundary system reveals the precise gain.
We have performed numerical experiments by using both the standard SBP-SAT scheme, and a modified version to confirm that our accuracy analysis gives sharp error estimates.
In this paper, we have performed accuracy analysis for problems in one space dimension. The technique and results can be straightforwardly generalized to two space dimensions if the boundary condition is periodic in one of the spatial dimensions. In a coming paper, we will show how to perform accuracy analysis for problems in two space dimensions when the boundary conditions in both dimensions are non-periodic. We will in particular consider a corner problem when the large truncation error is located on a few grid points at a conner in two space dimensions.  The boundary systems for the interface problem are Together with the estimate (6) and the Sobolev inequality [8,Lemma 8.3.1], we obtain that U is bounded by the data in maximum norm where U (x, t) ∞ := sup{|U (x, t)| : x ≥ 0} for any fixed t. Next, let V = U t . Then V satisfies with initial condition F(x, 0), and homogeneous Neumann boundary condition V x (0, t) = 0.
The equation in V is in the same form as the equation in U , but with a different forcing and initial data. In the same way, V is bounded by the data in maximum norm Before analyzing the semi-discrete case, we recall the definition of pointwise stability, Definition 2.2, in [21].
Definition 2.2 in [21]: The approximation, v, is strongly pointwise stable if, for all h ≤ h 0 , the estimate v(t) 2 ∞ ≤ K (t)( f 2 + max 0≤τ ≤t holds. Here K (t) is a bounded function in any finite time interval and does not depend on the data. ( · denotes some norm.) The approximation is pointwise stable if (48) holds with g(t) = 0.
In the above definition, v(t) is the semi-discrete solution, f is the initial data, F is the forcing in the equation and g is the boundary data. We also recall Theorem 2.13 in [21].
Theorem 2.13 in [21]:If v and v t are pointwise stable discrete solutions to (24'), then with r = 2 p − 2 the global order of accuracy is 2 p.
In the above theorem, (24') is the SBP-SAT discretization of the wave equation, v is the semi-discrete solution, v t is the time derivative of v, and r ? is the order of boundary truncation error. This theorem says that for the wave equation if the semi-discrete solution and its time derivative are pointwise bounded, then two orders are gained in convergence.
In the following, we derive the pointwise estimates for the semi-discrete solution and its time derivative. The discrete energy estimate (33) gives where g andḡ are the restrictions of G(x) andḠ(x) to the grid. and The term ḡ 2 A can be bounded by using Lemma 3. For the term − H −1 Ag H , we note Dg H also appears in the term g A in Lemma 3 and we use the same bound here. Only the first component of H −1 E 0 Sg is nonzero and is where K depends on the particular form of S but not h. As a consequence, − H −1 Ag H is bounded by G(x) and its derivatives up to order 2 p + 2 in maximum norm. We can then again use Lemma 3.2 in [21] to derive a pointwise bound on u t .
In conclusion, with an SBP-SAT finite difference scheme for the wave equation with Neumann boundary condition, the numerical solutions u and u t are pointwise bounded according to Definition 2.2 in [21]. However, the determinant condition is not satisfied as shown in Sect. 3, nor is the p + 2 optimal gain obtained. In particular, for p = 1 the gain is only 1 order.