Convergence of summation-by-parts finite difference methods for the wave equation

In this paper, we consider finite difference approximations of the second order wave equation. We use finite difference operators satisfying the summation-by-parts property to discretize the equation in space. Boundary conditions and grid interface conditions are imposed by the simultaneous-approximation-term technique. Typically, the truncation error is larger at the grid points near a boundary or grid interface than that in the interior. Normal mode analysis can be used to analyze how the large truncation error affects the convergence rate of the underlying stable numerical scheme. If the semi-discretized equation satisfies a determinant condition, two orders are gained from the large truncation error. However, many interesting second order equations do not satisfy the determinant condition. We then carefully analyze the solution of the boundary system to derive a sharp estimate for the error in the solution and acquire the gain in convergence rate. The result shows that stability does not automatically yield a gain of two orders in convergence rate. The accuracy analysis is verified by numerical experiments.

1. Introduction. In many physical problems, such as acoustics, seismology, and electromagnetism, the underlying equations can be formulated as systems of second order time dependent hyperbolic partial differential equations. For wave propagation problems, it has been shown that high order accurate time marching methods as well as high order spatial discretization are more efficient to solve these problems on smooth domains [5,7]. The major difficulty with high order spatial discretization is the numerical treatment of boundary conditions and grid interface conditions. High order finite difference methods have been widely used for wave propagation problems. We use the SBP finite difference operators [8,11,15] to approximate the spatial derivatives. In order to guarantee strict stability and maintain high accuracy, the boundary conditions and grid interface conditions are imposed weakly by the SAT technique [1,2,9,10]. A summary of SBP-SAT methodology is made in [3]. Another candidate for the approximation of spatial derivatives is discontinuous Galerkin (dG) method, and it has been successfully used in [4] for the wave equation. dG can handle complex geometry of the computational domain, but often has more degrees of freedom than finite difference method for the same accuracy. An efficient implementation of dG might be challenging.
In the SBP-SAT framework, accuracy analysis has drawn less attention than stability. When using finite difference method to solve first order hyperbolic partial differential equations, it is in many cases possible to show that p + 1 convergence could be obtained if the boundary truncation error is O(h p ), where h is the mesh size [5]. In other words, we gain one order in convergence rate. For the scalar wave equation in the second order form, the hypothesis in [14] is that we gain two orders in convergence rate. This is often observed in numerical experiments. The hypothesis is proved for Schrödinger equation in [12,13]. In this paper, we use the same technique and show that for the second order wave equation, the penalty parameter in SAT plays an important role for the convergence rate. Two orders in convergence are gained if the penalty parameter is well chosen. Even if the numerical scheme is stable, a bad choice of the penalty parameter results in significant loss of accuracy.
The structure of the paper is as follows. We start with the SBP-SAT method applied to the one dimensional wave equation with Dirichlet boundary condition. We derive error estimates by the energy method. These estimates are weaker than the estimates indicated by the hypothesis, and are valid as long as the penalty parameter in the boundary term is chosen so that stability is ensured. We then use normal mode analysis for the second, fourth and sixth order methods to derive stronger error estimates. In §3, we consider the one dimensional wave equation on a grid with a grid interface. We extend the analysis to two dimension in §4. In §5, we perform numerical experiments. The computational convergence results support the accuracy analysis. The conclusions are drawn in §6.

One dimensional wave equation.
To describe the properties of the equation and the numerical methods, we need the following definitions. Let w 1 (x) and w 2 (x) be real-valued functions in L 2 [a 1 , a 2 ]. The inner product is defined as (w 1 , w 2 ) = a2 a1 w 1 w 2 dx. The corresponding norm is w 1 2 = (w 1 , w 1 ). We consider the second order wave equation in one dimension with a forcing term F added on the right hand side: We consider the half line problem with Dirichlet boundary condition: The forcing function F , the initial data f 1 , f 2 and the boundary data g are compatible smooth functions with compact support. The L 2 norm of the solution is bounded, i.e. U (·, t) < ∞. With homogeneous Dirichlet boundary condition g(t) = 0, we multiply Equation (2.1) by U t and integrate by parts, and obtain is the continuous energy. The energy estimate follows from Gronwall's lemma Next, we introduce the 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 . We also define an inner product and norm for the vectors a and b in a Hilbert space over R as (a, b) H = a T Hb and a 2 H = a T Ha, respectively, where a T denotes the transpose of a and H is positive definite. An SBP operator mimics the property of the continuous integration-by-parts via the inner product and norm defined above.
In this paper, we use the diagonal norm SBP operator D approximating the second order derivative, i.e. D ≈ ∂ 2 ∂x 2 . The order of accuracy is denoted by 2p, p = 1, 2, 3, 4, 5. For a 2p th order diagonal norm SBP operator [8,11], the error for the approximation of the spatial derivative is O(h 2p ) in the interior, and O(h p ) near the boundary. The SBP operator can be written as D = H −1 (−A+BS), 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 with the order of accuracy p + 1.  [8,11]. We use the SAT technique to impose the Dirichlet boundary condition weakly. The semi-discretized equation corresponding to (2.1) and (2.2) with homogeneous boundary data reads: where e 0 = [1, 0, 0, · · · ] T , E 0 = e 0 e T 0 and F g is the grid function corresponding to F (x, t). The penalty parameter τ is to be determined so that an energy estimate of the discrete solution exists, which ensures stability. In a realistic setting, the problem would be posed on a finite interval with boundary conditions at both ends. In the half line case, the second boundary condition is replaced by u h < ∞.

Stability.
In [1,9,10], it is shown that the operator A can be written as where α 2p > 0 is as large as possible andÃ is symmetric positive semi-definite. This is often called the borrowing trick. We use the technique in [9] to compute the values of α 2p and list the results in Table 1 H + u 2Ã + Q is a discrete energy, and |||·||| h is the corresponding discrete energy norm. We then obtain the discrete energy estimate where E h,0 is the initial discrete energy. Clearly, stability is ensured if τ ≥ τ 2p := 1 α2p . Remark It is possible to include the term u 2 H in the discrete energy E h . A similar energy estimate holds [6]. This is useful when comparing the numerical results with the theoretical analysis.
2.2. Accuracy analysis by the energy estimate. Assuming that the continuous problem has a smooth solution U (x, t), we can derive the error equation for the pointwise error ξ j = U (x j , t) − u j (t). The error equation corresponding to (2.4) is Since the number of grid points with the large truncation error O(h p ) is finite and independent of h, we have For 2p = 2, 4, 6, 8, 10 and small h, the first term is much smaller than the second one. Thus, we have T 2p h ≤K B h p+1/2 . By applying energy method to the error equation (2.6), we obtain Therefore, the convergence rate is at least p + 1 2 if the numerical method is stable, that is if τ ≥ τ 2p .

Normal mode analysis for the boundary truncation error.
To derive a sharp 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 We perform a normal mode analysis to analyze the boundary truncation error ǫ for second, fourth and sixth order SBP-SAT scheme. The boundary error equation is where ǫ h < ∞. In the analysis of ǫ, we take the Laplace transform in time of (2.9), for Re(s) > 0. Lets = sh. After multipling h 2 on both sides of (2.10), we obtain Note that the first three terms in the right hand side of (2.11) are h-independent. We assume that the true solution U (x, t) is a smooth function. Then |T 2p,B (s)| decreases fast for large |s|. In order to observe the asymptotic convergence behaviour, h must be sufficiently small. Thus, we only need to consider a vicinity of |s| = 0. In the analysis, we assume |s| ≤ K and consider Re(s) ≥ η > 0. If the semi-discretized equation is 2p q(τ > τ 2p ) q(τ = τ 2p ) 2 2 1.5 4 4 2.5 6 5 3.5 where C depends only on η, the final time t f and the derivatives of the true solution U . The results are summarized in Theorem 2.1. Theorem 2.1. For the second, fourth and sixth order stable SBP-SAT approximations of the second order wave equation (2.1-2.2) on a half line with a Dirichlet boundary condition, the rates q in (2.12) depend on τ , and are listed in Table 2.
In §5, we present results from numerical experiments, which agree well with the theoretical results in Table 2. After some preliminaries we will prove Theorem 2.1 for 2p = 2, 4 and 6 in §2.3.1, 2.3.2 and 2.3.3, respectively.
To begin with, we note that sufficiently far away from the boundary, the two penalty terms and the boundary truncation error in (2.11) are not present. The coefficients in D correspond to the standard central finite difference scheme. We have 2p = 2 : s 2ǫ j = D + D −ǫj , j = 3, 4, 5, · · · , 2p = 4 : s 2ǫ The corresponding characteristic equations are 2p = 2 : κ 2 − (2 +s 2 )κ + 1 = 0, If the numerical scheme is stable, then there is no root with |κ| = 1 for Re(s) > 0. The proof of this can be found in [6]. We will need the following specifics for the roots.
In the following analysis, we use Lemma 2.3 with |s| ≤ K and Re(s) ≥ η, where K and η are positive numbers. That is, 2ηh to the leading order.
2.3.1. Proof of Theorem 2.1 for the second order scheme. Away from the boundary, the solution of (2.11) is The first three rows of (2.11) are affected by the penalty terms. They arẽ (2.16) By Taylor expansion, it is straightforward to deriveT 2,B 0 = −Û xxx (0, s) to the leading order. We write (2.16) in the matrix vector multiplication form with the help of (2.15), and obtain We use |·| to denote the standard Euclidean norm of vectors and matrices. We want to derive an estimate for |Σ|, which requires to investigate the invertibility of C 2D (s, τ ) in a vicinity ofs = 0. We note that The determinant of C 2D (0, τ ) is given by det(C 2D (0, τ )) = 5−2τ . Clearly, C 2D (0, τ ) is singular if τ = 2.5 and nonsingular otherwise. The stability condition is τ ≥ τ 2 = 2.5. We consider these two relevant cases separately.
In the case τ > τ 2 , C 2D (0, τ ) is nonsingular. By the continuity of C 2D (s, τ ), for any γ > 1 and sufficiently small h we have Thus, to the leading order, Note that in the above, we use Parseval's relation twice. By arguing that the future can not affect the past, 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 ). In this case, the interior error is O(h 2 ) by (2.8), which is the dominating source of error. Therefore, the convergence rate is 2.

2.3.2.
Proof of Theorem 2.1 for the fourth order scheme. Away from the boundary, the solution of (2.11) is The first four rows of (2.11) are affected by the penalty terms. They arẽ (2.18) By Taylor expansion, we havê to the leading order. We investigate the solution of (2.18) in a vicinity ofs = 0 for Re(s) > 0. Fors = 0, we write (2.18) in the matrix vector multiplication form. By using (2.17), we obtain ] T and C 4D (0, τ ) is given in Appendix 1. We solve for det(C 4D (0, τ )) = 0, and obtain τ = 2(4834 Since the stability condition is τ ≥ τ 4 ≈ 3.986350342, we need to consider the two cases, τ > τ 4 and τ = τ 4 , separately. In the case τ > τ 4 , C 4D (0, τ ) is nonsingular. By the continuity of C 4D (s, τ ), for any γ > 1 and sufficiently small h, we have In L 2 norm, we have By Lemma 2.3, the second term can be bounded as For small h, the first and third terms can be bounded as We then obtain, to the leading order, By Parseval's relation and by arguing that the future can not affect the past, we obtain Thus, the boundary error is O(h 4 ). In this case, the interior error is also O(h 4 ) given by (2.8). Therefore, the convergence rate is 4.

2.3.3.
Proof of Theorem 2.1 for the sixth order scheme. Away from the boundary, the solution of (2.11) is Near the boundary, the first six rows of (2.11) are affected by the penalty terms. In the same way as for the second and fourth order method, we obtain a linear system of equations in Σ = [ǫ 0 ,ǫ 1 ,ǫ 2 , σ 1 , σ 2 , σ 3 ] T : where C 6D (s, τ ) is given in Appendix 2. We find that C 6D (0, τ ) is singular if τ = τ 6 , and nonsingular otherwise. When τ > τ 6 , we use Taylor expansion to expressT 6,B v in terms of the derivatives of the true solution U , and Lemma 2.3 to derive an estimate for (2.19). The boundary error is O(h 5 ). In this case, the interior error is O(h 6 ) given by (2.8). Thus, the convergence rate is 5.
When τ = τ 6 , the convergence rate p + 1 2 = 3.5 is given by the energy estimate. 3. One dimensional wave equation with a grid interface. In this section, we consider the Cauchy problem for the second order wave equation in one dimension where the forcing function F , the initial data f 1 and f 2 are compatible smooth functions with compact support. The L 2 norm of the solution is bounded, i.e. U (·, t) < ∞. It is straightforward to derive the energy estimate of the form (2.3).
We solve the equation on a grid with a grid interface. 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 interface is located at x = 0. The grid functions are u L j (t) ≈ U (x −j , t) and u R j (t) ≈ U (x j , t). Denote e 0L = [· · · , 0, 0, 1] T , e 0R = [1, 0, 0, · · · ] T , where F gL and F gR are the grid functions corresponding to F (x, t). The penalty parameters τ L and τ R are chosen so that the semi-discretization is stable. By the energy method, the numerical scheme is stable if τ L = −τ R ≥ hL+hR 4α2phLhR :=τ 2p . The stability proof is found in [9, Lemma 2.4, p. 215].
By the energy estimate, the convergence rate is at least p + 1 2 if the semidiscretization is stable. In order to derive a sharper estimate, we follow the same approach as for the half line problem in §2. Denote 2p q(τ >τ 2p ) q(τ =τ 2p ) 2 2 1.5 4 4 2.5 Table 3: Theoretical convergence result for one dimensional wave equation with a grid interface.
and τ = τ L = −τ R . Assuming that the continuous problem has a smooth solution U (x, t), we obtain the error equation with the truncation error T 2p,L and T 2p,R as the forcing terms. We partition the truncation error into two parts such that L T 2p,L,I and h 2p R T 2p,R,I are the interior truncation error, h p L T 2p,L,B and h p R T 2p,R,B are the interface truncation error. Next, we partition the error into two parts such that ξ L = ǫ L + ǫ L,I and ξ R = ǫ R + ǫ R,I , where ǫ L,I and ǫ R,I are due to the interior truncation error, ǫ L and ǫ R are due to the interface truncation error. Denote h = h L = rh R , where r is the grid size ratio. ǫ L,I and ǫ R,I can be estimated by the energy method, yielding (3.3) ǫ L,I 2 h + ǫ R,I 2 h ≤ K I h 2p . We take the Laplace transform in time of the interface error equation, and obtaiñ wheres L = sh L ands R = sh R . With the same assumption of the true solution U as in the previous case, in the analysis we consider |s| ≤ K and Re(s) ≥ η > 0. We use normal mode analysis to derive the error estimates for the second and fourth order method. The result is summarized in the following theorem. Theorem 3.1. For second and fourth order stable SBP-SAT approximation (3.2) of the Cauchy problem with a grid interface, the numerical solution converges to the true solution at rate q. The values of q depend on τ , and are listed in Table 3.
After some preliminaries we will prove Theorem 3.1 for 2p = 2 and 4 in §3.1 and 3.2, respectively. To begin with, we note that sufficiently far away from the interface, the three penalty terms and the interface truncation errors in (3.4) and (3.5) are 0. Thus, the coefficients in D L and D R correspond to the standard central finite difference scheme. On each side, the corresponding characteristic equations are given by (2.13), and the solutions of the characteristic equations are given by (2.14).

Proof of Theorem 3.1 for the second order scheme. Away from the interface, the solutions of (3.4) and (3.5) arê
Near the interface, the last three rows of (3.4) and the first three rows of (3.5) are affected by the penalty terms. We write them in the matrix vector multiplication form, and obtain a linear system of equations: The Taylor expansion of C 2I (s, τ ) ats = 0 is given in Appendix 3. In order to derive an estimate for |Σ|, we investigate the invertibility of C 2I (s, τ ) in a vicinity of |s| = 0 with Re(s) > 0. We find that C 2I (0, τ ) is singular with rank 5 if τ >τ 2 , and singular with rank 4 if τ =τ 2 . We consider these two cases separately.
In the case τ >τ 2 , we use the Taylor expansion of C 2I (s, τ ) ats = 0, that is, , and Lemma 3.4 in [12, pp. 539] to derive an estimate for |C −1 2I (s, τ )|. We need to perform singular value decomposition (SVD) on C 2I (0, τ ), which is very difficult due to the variable τ in the matrix. Instead, we compute numerically the SVD of C 2I (0, τ ) for τ = 1.1τ 2 , 1.2τ 2 , · · · , 5τ 2 . In all cases, we find that (U * C ′ 2I (0, τ )V ) 66 = 0, andT 2,B v is not in the column space of C(0, τ ). Therefore, in a vicinity of |s| = 0, we have |C −1 2I (s, τ )| ≤ Kc|s| −1 . Thus, where Re(s) ≥ hη > 0 and K c is independent of h and s. Denote ǫ 2 h = ǫ L 2 hL + ǫ R 2 hR . In L 2 norm, we have By Parseval's relation and the argument that the future can not affect the past, we obtain Thus, the interface error is O(h 2 ). In this case, the interior error is also O(h 2 ) given by (3.3). Therefore, the convergence rate is 2.
In the case τ =τ 2 , the numerical scheme is still stable. By the energy estimate, the convergence rate is p + 1 2 = 1.5. In §5, the numerical experiment agrees with this result.
3.2. Proof of Theorem 3.1 for the fourth order scheme. Away from the interface, the solutions of (3.4) and (3.5) arê Near the interface, the last four rows of (3.4) and the first four rows of (3.5) are affected by the penalty terms. We write them in the matrix vector multiplication form, and obtain a linear system of equations: The Taylor expansion of C 4I (s, τ ) ats = 0 is given in Appendix 4. In order to derive an estimate for |Σ|, we investigate the invertibility of C 4I (s, τ ) in a vicinity ofs = 0. We find that C 4I (0, τ ) is singular with rank 7 if τ >τ 4 , and singular with rank 6 if τ =τ 4 . We consider these two cases separately.
In the case τ >τ 4 , Lemma 3.4 in [12, pp. 539] is applicable. Also, there exists a nontrivial vector y such that C 4I (0, τ )y =T 4,B v , i.e.T 4,B v is in the column space of C 4I (0, τ ). Therefore, we have where K Σ is independent of h and s. As for the second order case, it then follows that the interface error is O(h 4 ). In this case, the interior error is O(h 4 ) given by (3.3). Therefore, the convergence rate is 4.
In the case τ =τ 4 , the numerical scheme is still stable, and the energy estimate gives a sharp estimate. The convergence rate is p + 1 2 = 2.5. In §5, the numerical experiment agrees with this result.

Two dimensional wave equation.
In this section, we consider the second order wave equation in two dimension: We consider the half plane problem in x with Dirichlet boundary condition. We also assume that the solution is 2π-periodic in y: The functions F , f 1 , f 2 and g are compatible smooth functions with compact support in x and 2π-periodicity in y. The L 2 norm of the solution is bounded, i.e. U < ∞.
Similar to the wave equation in one dimension, with homogeneous Dirichlet boundary condition the continuous energy E = U t 2 + U x 2 + U y 2 is bounded by the data. Next, we introduce the grid and grid functions. In order to simplify notation, we assume that the mesh size is h in both x and y direction. The grid is x i = ih, i = 0, 1, 2, · · · , y j = jh, j = 0, 1, 2, · · · , N y , Ny , and the grid functions are u ij (t) ≈ U (x i , y j , t). We arrange the grid functions u ij (t) column wise in a vector u. The general notation used to expand an operator from one dimension to two dimension is A x = A x ⊗ I y , where ⊗ is Kronecker product and I y is the identity operator. With homogeneous boundary condition, the semi-discretized equation reads D xx is the standard SBP operator approximating ∂ 2 ∂x 2 and D yyp is the standard central finite difference operator approximating ∂ 2 ∂y 2 . In the same way as for the wave equation in one dimension, a discrete energy estimate is obtained if the penalty parameter τ is chosen so that τ ≥ 1 α2p , where the values of α 2p are listed in Table 1. 4.1. Accuracy analysis. We analyze the accuracy of the numerical scheme by normal mode analysis. Denote ǫ the pointwise error due to the boundary truncation error. Then the boundary error equation is where ǫ h < ∞. We take the Laplace transform in t and Fourier transform in y, and obtain wheres + = s 2 − h 2Dω yyp ,s = sh andD ω yyp is the Fourier transform of D yyp . For different order of accuracy, we havẽ Note thatD ω yyp is nonpositive. Lemma 4.1. For Re(s) > 0 andD ω yyp defined above, it follows that Re(s + ) > 0. Proof. By the definition of principal square root, Re(s + ) ≥ 0. Assume Re(s + ) = 0, thens 2 − h 2Dω yyp is real and nonpositive, i.e. Re(s 2 − h 2Dω yyp ) ≤ 0. We then obtain Re(s 2 ) ≤ Re(h 2Dω yyp ) ≤ 0. This contradicts Re(s) > 0. Therefore, Re(s + ) > 0. The truncation error T 2p can be written as i0 ; T 2p i1 ; · · · ; T 2p iNy ] is an N y + 1-by-1 vector. For 2p th order accurate method, we have to the leading order, where m = 1, 4, 6, 8, 11 for 2p = 2, 4, 6, 8, 10, c k is a constant independent of h, U (p+2) := ∂ p+2 U ∂x p+2 and 0 is a zero vector of size N y + 1-by-1 . We take the Fourier transform in y and Laplace transform in t, and obtaiñ Similar to the one dimensional case, we only need to consider a vicinity of |s| = 0 with Re(s) > 0. By Lemma 4.1, we have Re(s + ) > 0. The boundary part of Equation (4.4) can be written in the matrix vector multiplication form C(s + , τ )Σ ω = h p+2T 2p,B ω . We then need to investigate the invertibility of C(s + , τ ) and the dependence of |C −1 (s + , τ )| on h. We note that Equation (4.4) is analogous to the error equation (2.11) in one dimensional case, and the matrix-valued function C(·, τ ) is the same with only the argument changed. Therefore, based on the normal mode analysis for the one dimensional problem in §2.3.1-2.3.3, for 2p = 2, 4, 6 and τ > τ 2p the boundary error in the Fourier and Laplace space is bounded as where K ω depends on η but not h. Since the above estimate holds for every ω, we sum them up and obtain By using the Parseval's relation in the Fourier space, we obtain We then use Parseval's relation in the Laplace space, and obtain For 2p = 2, 4, 6, the boundary error is bounded in O(h p+2 ). We know from the energy estimate that the interior error is bounded in O(h 2p ). Therefore, we conclude that for second, fourth and sixth order methods, the convergence rates are 2, 4 and 5 only if τ > τ 2p . If τ = τ 2p , the convergence rate is p + 1 2 given by the energy estimate. The numerical experiments in §5 agree with this conclusion.

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 the boundary and grid interface. In the analysis, we use the half line and half plane problem. However, in the numerical experiments, we use a bounded domain and impose boundary conditions on all boundaries weakly. We discretize in time using the classical fourth order Runge-Kutta method. The step size ∆t in time is chosen so that the temporal error has a negligible impact of the spatial convergence result. The value of ∆t is given in each numerical experiment.
The L 2 error and maximum error 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 dimension of the equation. The convergence rates are computed by 5.1. 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 analytical solution does not vanish at the boundaries. The L 2 error, the convergence rates in L 2 norm and maximum norm are shown in Table 4. The convergence rates behave as expected. With a 2p th (2p = 2, 4, 6) order accuracy method, the convergence rate is p + 1 2 in L 2 norm if the penalty parameter τ equals its limit. In this case, the pointwise error in the Laplace space is bounded in O(h p ). We can hope for p th order convergence rate in the maximum norm. With an increased penalty parameter, the accuracy of the numerical solution is improved significantly, and the convergence rate p + 2 is obtained. For second and fourth order accuracy cases, the accuracy difference between τ = 1.2τ 2p and τ = 3τ 2p is almost invisible. For sixth order accuracy case, one can see the small improvement in accuracy between τ = 1.2τ 2p and τ = 3τ 2p . The step size in time is ∆t = 0.1h in all cases, except for the last two grid refinements of the sixth order scheme with τ = 3τ 2p where ∆t = 0.05h is used in order to observe the expected convergence behaviour. Note that these step sizes are much smaller than the step size restricted by the stability condition. With τ = 1.2τ 2p , the numerical method is stable if ∆t/h ≤ c, where the Courant number c = 1.3, 0.9 and 0.6 for second, fourth and sixth order method. For an increased penalty parameter τ = 3τ 2p , the corresponding Courant numbers are 0.8, 0.5 and 0.4.
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 (5.1). The grid interface is located at x = 0.5 in the computational domain Ω = [0, 1], and the grid size ratio is 2 : 1. 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 its limit. 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 Table 5, where N denotes the number of grid points in the coarse domain.
According to the accuracy analysis in §3, for second and fourth order accuracy methods with τ = τ 2p , the expected convergence rates are 1.5 and 2.5 in L 2 norm. This is clearly observed in Table 5. We also observe in Table 5 that the convergence rate in the maximum norm behaves as p. With an increased penalty parameter, a much better convergence result is obtained. For the second and fourth order methods, we get second and fourth order convergence in both L 2 norm and maximum norm. For the sixth order method, we expect fifth order convergence in L 2 norm. The numerical experiments give a higher convergence result. This has also been observed for Schrödinger equation in [12,13]. The step size in time is ∆t = 0.1h.    Remark The numerical experiments with 8 th and 10 th order SBP operators are also carried out. When the penalty parameter equals its lower limit, 4.5 and 5.5 convergence rate in L 2 norm, and 4 and 5 convergence in maximum norm are clearly observed. With a larger penalty parameter, we expect the convergence to be 6 and 7, respectively. However, due to the high accuracy, the error decreases fast to the machine precision. Therefore, the expected convergence rates are not clearly observed.

Two dimensional wave equation.
For two dimensional simulation, we choose the analytical solution u = cos(10x + 1) cos(10y + 2) cos(10 √ 2t + 3), 0 ≤ x ≤ 1, 0 ≤ y ≤ 1, 0 ≤ t ≤ 2, and use it to impose the Dirichlet boundary conditions. The convergence results are shown in Table 6. We clearly observe that with a 2p th order method and τ = τ 2p the convergence rates in L 2 norm are p + 1 2 and in maximum norm are p. With an increased penalty parameter, the errors in L 2 norm are much smaller. The convergence rates are p + 2 in both L 2 norm and maximum norm. The results from the numerical experiments agree well with the analysis in §4. The step size in time is ∆t = 0.1h in all cases, except for the sixth order scheme with τ = 1.2τ 2p where ∆t = 0.05h is used in order to observe the expected convergence behaviour.
6. Conclusion. In this paper, we consider the accuracy of the finite difference solution of the second order wave equation by the SBP-SAT methodology. The large truncation error is localized near the boundary and grid interface. The Dirichlet boundary condition and interface condition are imposed weakly by the SAT technique. The penalty parameter plays an important role in the accuracy and convergence rate of the numerical method. For stability, the penalty parameter has a lower limit. If the Second Order Fourth Order Sixth Order     Table 6: Convergence for two dimensional wave equation.
lower limit is used, for 2p th order SBP approximation the convergence rate is p + 1 2 in L 2 norm, and p in maximum norm. With an increased penalty parameter, a much better convergence result is obtained. By using normal mode analysis, we point out that for second, fourth and sixth order method, the convergence rates in L 2 norm are 2, 4 and 5. Therefore, p + 2 hypothesis is true only if the penalty parameter is strictly larger than its limit required by stability. The numerical experiments agree well with the analysis. In practice, we should choose a penalty parameter strictly larger than its limit value so that the optimal convergence rate is obtained.  where k(τ ) = 3.165067038τ − 9.382128117, and the last two columns of the matrix