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

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)\ge 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 non–optimal 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.

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,\infty )\). An inner product is defined as \((w_1,w_2)=\int _0^{\infty } w_1w_2 dx\) with a corresponding norm \(\Vert w_1\Vert ^2=(w_1,w_1)\), and \(w_1\) is said to be in L\(_2\) if \(\Vert w_1\Vert <\infty \).

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, \(\bar{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. \(\Vert U(\cdot ,t)\Vert < \infty \). 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 \({\left| \left| \left| U \right| \right| \right| }^2=\Vert U_t\Vert ^2+\Vert U_x\Vert ^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

and

to obtain

An integration in time, together with (5), gives

Therefore, \(\Vert U\Vert \) is bounded by the data in any bounded time interval \(t\in [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,\cdots ,\) and a grid function \(u_i(t)\approx U(x_i,t)\). Furthermore, let \(u(t)=[u_0(t),u_1(t),\cdots ]^T\) where \(^T\) denotes the transpose. We also define an inner product and norm for the grid functions a and b in \(\mathcal {R}\) as \((a,b)_H=a^THb\) and \(\Vert a\Vert _H^2=a^THa\), respectively, where H is a symmetric positive definite operator on the space of grid functions. The norm \(\Vert \cdot \Vert _H\) is equivalent to the standard discrete L\(_2\) norm, and they are the same if \(H=hI\) 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\approx \partial ^2/\partial x^2\) has a structure

where H is diagonal and positive definite, A is symmetric positive semi-definite and \(B=\text {diag}(-1,0,0,\cdots )\). 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 \(2p^{th}\) order accurate, the truncation error of D is \(\mathcal {O}(h^{2p})\) in the interior but only \(\mathcal {O}(h^{p})\) near the boundary. The truncation error of S is \(\mathcal {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 \(P_D=-H^{-1}S^TE_0-\tau h^{-1}H^{-1}E_0\), \(e_0=[1,0,0,\cdots ]^T\), \(E_0=e_0e_0^T\) 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_{xx}\), while the second term \(P_Du\) 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 \(\tau \) is to be determined so that an energy estimate of the discrete solution exists, which ensures stability.

2.1 Stability

In [2, 13, 14], it is shown that the operator A in (7) can be written as

where \(\alpha _{2p}>0\) is a constant independent of h and \(\tilde{A}\) is symmetric positive semi-definite. This is often called the borrowing trick. The values of \(\alpha _{2p}\) are computed for \(p=1,2\) and 3 in [13]. We use the same technique to compute \(\alpha _{2p}\) for \(p=4\) and 5, and list the results in Table 1.

To derive an energy estimate, we introduce the energy

where \(\tau \ge 1/\alpha _{2p}\). Multiplying Eq. (8) by \(u_t^TH\) 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 semi–discretization (8) is energy stable if \(\tau \ge 1/\alpha _{2p}\). The term \(\Vert u\Vert _H\) can be included in the energy estimate in the same way as in Remark 1.

2.2 Accuracy Analysis by the Energy Method

The error equation for the pointwise error \(\zeta _j=U(x_j,t)-u_j(t)\) is

where \(\Vert \zeta (t)\Vert _h<\infty \), \(T^{2p}\) is the truncation error and 2p 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^{2p}_j=0,\ j\ge J\) for some \(J\sim \mathcal {O}(h^{-1})\). The first m components of \(T^{2p}\) are of order \(\mathcal {O}(h^p)\), and all the other nonzero components are of order \(\mathcal {O}(h^{2p})\). For \(2p=2,4,6,8,10\), the corresponding values of m are 1, 4, 6, 8, 11. In the analysis, we only consider the leading term of the truncation error. We introduce the interior truncation error \(T^{2p,I}\), and the boundary truncation error \(T^{2p,B}\) such that \(T^{2p}=h^{2p}T^{2p,I}+h^pT^{2p,B}\), where \(T^{2p,I}\) and \(T^{2p,B}\) are independent of h, but depend on the derivatives of U. Since the number of grid points with the large truncation error \(\mathcal {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 \(\Vert T^{2p}\Vert _h\le 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 \(\tau \ge 1/\alpha _{2p}\).

2.3 Normal Mode Analysis for the Boundary Truncation Error

To derive a sharper estimate, we partition the pointwise error into two parts, the interior error \(\epsilon ^I\) and the boundary error \(\epsilon \), such that \(\zeta =\epsilon ^I+\epsilon \). \(\epsilon ^I\) is the error due to the interior truncation error, and \(\epsilon \) is the error due to the boundary truncation error. \(\epsilon ^I\) can be estimated by the energy method, yielding

As mentioned above, \(\Vert {\epsilon ^I}\Vert _H\) can be bounded similarly.

We perform a normal mode analysis to analyze \(\epsilon \) for the second, fourth and sixth order SBP-SAT schemes. The boundary error equation is

where \(\Vert \epsilon \Vert _h<\infty \). In the analysis of \(\epsilon \), we take the Laplace transform in time of (14),

for \(\text {Re}(s)>0\), where \(\hat{}\) denotes Laplace transform. After multiplying by \(h^2\) on both sides of (15), we obtain with the notation \(\tilde{s}=sh\)

Note that the coefficients in the first two terms in the right hand side of (16) are h-independent, and that \(\hat{T}^{2p,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\((\tilde{s})=\eta h\) with \(\eta >0\) a fixed constant independent of h. It is therefore important to derive a sharp error estimate in Laplace space when Re\((\tilde{s})\) goes to zero. The final estimate in the physical space is in the form

where K depends only on \(\eta \), 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 \(|\hat{T}^{2p} (\eta +i\xi )|\) decreases fast for large \(|\xi |\), making contributions from \(|\xi |>\bar{K}\) insignificant for some constant \(\bar{K}\). In the analysis we will only consider \(s=\eta +i\xi \), with \(\eta \) positive and \(|\xi |<\bar{K}\). Correspondingly we only consider \(|\tilde{s}|\ll 1\), Re\((\tilde{s})=\mathcal {O}(h)\). In particular it suffices to consider \(\tilde{s}=0\) when checking the determinant condition. We formalize this discussion by making the following assumption.

Assumption 1

There is a constant \(\bar{K}<\infty \) such that contributions from \(s=\eta +i\xi \), \(|\xi |>\bar{K}\) to the estimate (17) do not influence the rate q.

We summarize the convergence result for the Dirichlet problem in Theorem 1.

Theorem 1

Consider the second, fourth and sixth order stable SBP-SAT approximations (8) of the second order wave equation (1) with the Dirichlet boundary condition. With Assumption 1, the rates q in (17) depend on \(\tau \):

1. 1.

   If \(\tau =1/\alpha _{2p}\), then \(q=p+1/2\), i.e. the gain in convergence is only half an order.

2. 2.

   If \(\tau >1/\alpha _{2p}\), then the interior and boundary truncation errors lead to errors \(\mathcal {O}(h^{2p})\) and \(\mathcal {O}(h^{q_B})\) in the solution, respectively. The values \(q_B\) and the overall rates \(q=\min (2p,q_B)\) are listed in Table 2.

Remark 2

The penalty parameter \(\tau \) 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 \(\tau >1/\alpha _{2p}\) than with \(\tau =1/\alpha _{2p}\). We therefore should always choose \(\tau >1/\alpha _{2p}\) in computations, but bearing in mind that a moderate value of \(\tau \) is appropriate because a very large \(\tau \) has a negative impact on the CFL condition.

The convergence rates are optimal for the second and fourth order schemes, while super–convergence 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.

After some preliminaries we will prove Theorem 1 for \(2p=2,4\) and 6 in Sects. 2.3.2, 2.3.3 and 2.3.4, respectively.

2.3.1 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:

where

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 \(|\kappa |=1\) for \(\text {Re}(s)>0\). We call a root \(\kappa \) an admissible root if \(|\kappa |<1\). Since the problem in consideration is a half line problem, any \(\kappa \) with \(|\kappa |>1\) is not admissible because it results in an unbounded solution. We will need the following specifics for the roots.

Lemma 1

For \(2p=2,4,6\), the number of admissible roots of (19) satisfying \(|\kappa |<1\) for \(\text {Re}(\tilde{s})>0\) is 1,2,3, respectively. In the vicinity of \(\tilde{s}=0\), they are given by

Proof

\(2p=2\): Eq. (19a) has two roots: \(\kappa _{1,2}=1+\frac{1}{2}\tilde{s}^2\pm \frac{1}{2}\sqrt{\tilde{s}^4+4\tilde{s}^2}\). We find by Taylor expansion at \(\tilde{s}=0\) that \(\kappa _1=1-\tilde{s}+\mathcal {O}(\tilde{s}^2)\) and \(\kappa _2=1+\tilde{s}+\mathcal {O}(\tilde{s}^2)\). Thus, in the vicinity of \(\tilde{s}=0\) the admissible root is \(\kappa _1=1-\tilde{s}+\mathcal {O}(\tilde{s}^2)\).

\(2p=4\): Eq. (19b) has four roots:

We find by Taylor expansion at \(\tilde{s}=0\) that

Thus, the admissible roots are \(\kappa _1=1-\tilde{s}+\mathcal {O}(\tilde{s}^2)\) and \(\kappa _2=7-4\sqrt{3}+\mathcal {O}(\tilde{s}^2)\).

\(2p=6\): There is no closed form of the solution to a sixth order equation like (19c). However, (19c) with \(\tilde{s}=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 on \(\tilde{s}\) in a neighbourhood of \(\tilde{s}=0\). The three admissible roots are given by (20c). We show the numerical values with four digits since the analytical expressions are very lengthy. \(\square \)

The pointwise error away from the boundary can be expressed as

where the coefficients \(\sigma \) 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 \(\hat{\epsilon }\) 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 \(\frac{1}{1-|\kappa _1|^2}\) cannot be bounded independent of h when \(\tilde{s}=\mathcal {O}(h)\).

Lemma 2

Consider \(\text {Re}(\tilde{s})=\eta h>0\), where \(\eta \) is a constant independent of h. For \(2p=2,4,6\) the admissible root \(\kappa _1(\tilde{s})\) in 20 satisfies

to the leading order.

Proof

Let \(\tilde{s}=x+iy\) where x, y are real numbers. Then \(|x|\ge \eta h\).

The desired estimate follows. \(\square \)

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 \(\hat{\epsilon }_i\) and \(\sigma _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.

2.3.2 Proof of Theorem 1 for the Second Order Scheme

The first three rows of (16) are affected by the boundary closure. They are

By Taylor expansion, it is straightforward to derive \(\hat{T}^{2,B}_0=-\hat{U}_{xxx}(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\((\tilde{s})\ge 0\), then \(|\varSigma _{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 \(\varSigma _{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\((\tilde{s})=\eta h\) to estimate \(1/(1-|\kappa _1|^2)\), with \(\eta \) a constant independent of h. When analyzing the boundary system, we need to use the same \(\tilde{s}=\eta h\) because that is where the inverse Laplace transform is performed later. We write the solution to (25) as \(\varSigma _{2D}(\tilde{s},\tau )=h^3 C_{2D}^{-1} \hat{T}^{2,B}_v\) for Re\((\tilde{s})=\eta h\). Since \(\eta \) is a constant and h is arbitrarily small, it is important to check how \(\varSigma _{2D}(\tilde{s},\tau )\) behaves as \(\tilde{s}\) approaches zero. We start by setting \(\tilde{s}=0\), and find that when \(\tau \ne 2.5\) the determinant condition is satisfied, i.e. \(C_{2D}(0,\tau )\) is nonsingular; but the determinant condition is not satisfied when \(\tau =2.5\). The stability condition \(\tau \ge 1/\alpha _{2}=2.5\) motivates us to consider these two cases separately.

In the case \(\tau >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 \(\varSigma _{2D}\) at \(\tilde{s}=0\) shows that

which together with (23) leads to

Remember that \(|\tilde{s}|\ll 1\) corresponds to \(|s|<\bar{K}\) for some constant \(\bar{K}\) and \(h\ll 1\). By Assumption 1 and Parseval’s relation, we have

We have derived the estimate for \(\hat{\epsilon }\) in the vicinity of \(\tilde{s}=0\), but Assumption 1 allows us to use this estimate when integrating \(\xi \) from \(-\infty \) to \(\infty \). 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 \(\mathcal {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 \(\mathcal {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 \(\sigma _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 \(\tau =2.5\), the last two components of \(\varSigma _{2D}\) are infinite when \(\tilde{s}=0\) since \(C_{2D}(0,\tau )\) is singular and the determinant condition is not satisfied. We again perform a perturbation analysis to the solution \(\varSigma _{2D}\) and obtain

This, in the same way as the case \(\tau >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.

2.3.3 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\((\tilde{s})=\eta h\) and write the solution to the boundary system (27) as \(\varSigma _{4D}=h^4C_{4D}^{-1}(\tilde{s},\tau )\hat{T}_v^{4,B}\) where \(\varSigma _{4D}=[\sigma _1, \sigma _2, \hat{\epsilon }_0, \hat{\epsilon }_1]^T\), and

The matrix \(C_{4D}(0,\tau )\) is presented in “Appendix 1”.

In the case \(\tau >1/\alpha _4\), \(C_{4D}(0,\tau )\) is nonsingular and the determinant condition is satisfied. All the four components of \(\varSigma _{4D}\) are order \(\mathcal {O}(h^4)\) with \(\sigma _1\) independent of \(\tau \). 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 \(\mathcal {O}(h^4)\). In this case, the interior error is also \(\mathcal {O}(h^4)\) given by (13). Therefore, the convergence rate is 4. This proves that \(q_B=q=4\) in Theorem 1.

In the case \(\tau =1/\alpha _4\), \(C_{4D}(0,\tau _4)\) is singular and the determinant condition is not satisfied. By the energy estimate (12), the convergence rate is \(p+1/2=2.5\).

2.3.4 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 \(\varSigma _{6D}=[\sigma _1,\sigma _2,\sigma _3,\hat{\epsilon }_0,\hat{\epsilon }_1,\hat{\epsilon }_2]^T\):

where

The matrix \(C_{6D}(0,\tau )\) is presented in “Appendix 1”. It is singular if \(\tau =1/\alpha _{6}\), and nonsingular otherwise.

When \(\tau >1/\alpha _6\), we write the solution to (29) as \(\varSigma _{6D}(\tilde{s},\tau )=h^5C_{6D}^{-1}\hat{T}_v^{6,B}\). A calculation of \(\varSigma _{6D}(\tilde{s},\tau )\) at \(\tilde{s}=0\) shows that \(\sigma _1=0\) and the other five components \(\mathcal {O}(h^5)\), and a further perturbation analysis gives \(\sigma _1\sim \mathcal {O}(h^7)\) when \(\tilde{s}=\mathcal {O}(h)\). As a consequence, from (23) we obtain that the boundary error is \(\mathcal {O}(h^{5.5})\). In this case, the interior error is \(\mathcal {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.

When \(\tau =1/\alpha _{6}\), the convergence rate \(p+1/2=3.5\) is given by the energy estimate (12).

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.

3.1 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_{xx}\) 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_0S)\). As a consequence, the semi–discretization (31) can be written as

Multiplying Equation (32) by \(u_t^TH\) from the left, we obtain

The discrete energy \({\left| \left| \left| u \right| \right| \right| }_{N,h}^2= \Vert u_t\Vert _H^2+\Vert u\Vert _A^2 \) is bounded as

Here again, \(\Vert u\Vert _H\) can be included in the energy estimate as in Remark 1.

3.2 Accuracy

With the pointwise error defined as before, the error equation is

where \(\Vert \zeta \Vert _h<\infty \), \(T^{2p,N}\) is the truncation error and 2p 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 \(\zeta \) into two parts, the error \(\epsilon ^I\) due to the interior truncation error and the error \(\epsilon \) due to the boundary truncation error. \(\epsilon ^I\) can be estimated by the energy method, yielding \({\left| \left| \left| \epsilon ^I \right| \right| \right| }_{N,h}\le K h^{2p}\). The boundary error equation is

where \(h^pT^{2p,N,B}\) is the boundary truncation error. \(T^{2p,N,B}\) depends on the derivatives of the true solution U, but not h. Moreover, only the first m elements of \(T^{2p,N,B}\) are nonzero, where m depends on p but not h.

We Laplace transform Equation (35) in time. With the notation \(\tilde{s}=sh\), we obtain

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.

3.2.1 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 at \(\tilde{s}=0\), but we straightforwardly obtain

for \(\tilde{s}=\mathcal {O}(h)\). 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 \(\mathcal {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 \(\mathcal {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.

3.2.2 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

for the sixth order scheme. The truncation errors \(\hat{T}_u^{4,N,B}\), \(\hat{T}_u^{6,N,B}\) are the same as the Dirichlet case (28), (30) except the coefficients of the first component are 43/204 and -53845/163788, respectively. The matrices \(C_{4N}(0)\) and \(C_{6N}(0)\) presented in “Appendix 1” are singular, indicating that the determinant condition is violated. We consider \(\tilde{s}=\eta h\) with a constant \(\eta \) independent of h, and analyze the solution to the boundary system. For the fourth order case, \(\sigma _1\) in (23) is \(\sigma _1\sim \mathcal {O}(\tilde{s} h^4)\sim \mathcal {O}(\eta h^5)\), and the other components are bounded by \(\mathcal {O}(h^4)\). For the sixth order case \(\sigma _1,\hat{\epsilon }_0\sim \mathcal {O}(\tilde{s} h^5)\sim \mathcal {O}(\eta h^6)\), and the other components are bounded by \(\mathcal {O}(h^5)\). Both cases lead to a gain of 2.5 orders from the boundary truncation error. We summarize the results on the Neumann problem in the following theorem.

Theorem 2

Consider the second, fourth and sixth order stable SBP-SAT approximations (31) of the second order wave equation (1) with the Neumann boundary condition. With Assumption 1, the boundary truncation error lead to error \(\mathcal {O}(h^{q_B})\) in the solution, where \(q_B\) and the overall rates \(q=\min (2p, q_B)\) are listed in Table 3.

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 \(\bar{G}\) 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,\cdots \), and the grid on the right \(x_{j}=jh_R,\ j=0,1,2,\cdots \), with the grid size \(h_L\) and \(h_R\), respectively. The grid functions are \(u_j^L(t)\approx U(x_{-j},t)\) and \(u_j^R(t)\approx U(x_j,t)\). Denote

Both \(u_0^L=e_{0L}^Tu^L\) and \(u_0^R=e_{0R}^Tu^R\) approximate U(0, t). Both \((S_Lu^L)_0=e_{0L}^TS_Lu^L\) and \((S_Ru^R)_0=e_{0R}^TS_Ru^R\) approximate \(U_x(0,t)\). To simplify notation, we define \(\{u\}=u_0^L-u_0^R\) and \(\{Su\}=(S_Lu^L)_0-(S_Ru^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_{IL}\{u\}\) and \(P_{IR}\{u\}\) impose the interface conditions weakly and have the form

The penalty parameters \(\tau \) 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 \(\tilde{\tau }_{2p}=1/(2\alpha _{2p}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 \(\mathcal {O}(h^{2p})\) in the solution. In the remaining part of this section we will only consider the effect of the interface truncation error. Denote

The error equation in Laplace space reads

where \(\tilde{s}_L=sh_L\), \(\tilde{s}_R=sh_R\). Here the interface truncation errors \(\hat{T}^{2p,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 \(\tau =\tilde{\tau }_{2p}\) in (40), the rate \(q=p+1/2\) is given by the energy estimate. When \(\tau >\tilde{\tau }_{2p}\), the rate q is listed in Table 4.

Remark 5

In practical computations, we should always choose \(\tau >\tilde{\tau }_{2p}\) 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 (2p,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 \(\tau >\tilde{\tau }_{2p}\). 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 when \(\tilde{s}=\mathcal {O}(h)\).

4.1 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

where

The determinant condition is not satisfied because \(C_{2I}(0,\tau )\) shown in “Appendix 1” is singular. We find that when \(r=1\),

which leads to a boundary error \(\mathcal {O}(h^{3.5})\), i.e. a gain of 2.5 orders in convergence and \(q_B^{(h_L=h_R)}=3.5\), in the same way as in Sect. 2.3.2. When \(r\ne 1\), all components of \(\varSigma \sim \mathcal {O}(\tilde{s}^{-1}h^3)\sim \mathcal {O}(\tilde{h}^2)\). In this case, the gain in convergence is only 1 order, i.e. \(q_B^{(h_L\ne h_R)}=2\). In both cases, the overall convergence rate is \(q=2\).

4.2 Proof of Theorem 3 for the Fourth Order Scheme

In this case, the eight by eight boundary system takes the form

where

The boundary system is also singular at \(\tilde{s}=0\), and the determinant condition is not satisfied. The matrix \(C_{4I}(0,\tau )\) is presented in “Appendix 1”. When \(r=1\),

which leads to a boundary error \(\mathcal {O}(h^{4.5})\), i.e. a gain of 2.5 orders in convergence and \(q_B^{(h_L=h_R)}=4.5\). When \(r\ne 1\), all components of \(\varSigma \sim \mathcal {O}(h^4)\). In this case, we gain two orders in convergence, i.e. \(q_B^{(h_L\ne h_R)}=4\). In both cases, the overall convergence rate is \(q=4\).

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 \(\varDelta t\) in time is chosen so that the temporal error has a negligible impact on the spatial convergence result. The value of \(\varDelta 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

5.1 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 \(\varDelta 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.

A convergence test for the one dimensional wave equation with a Dirichlet and b Neumann boundary conditions

Full size image

With a \(2p^{th}\ (2p=2,4,6)\) order accurate method, the convergence rate is \(p+1/2\) in L\(_2\) norm if the penalty parameter \(\tau \) equals its limit. With a larger penalty parameter (increased by 20%), the accuracy of the numerical solution is improved significantly. The observed convergence rates by (42) for the second, fourth and sixth order accurate schemes are 2.00, 3.97 and 5.56, which are very close to the theoretical rates 2, 4 and 5.5 in 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 \(\tau >1/\alpha _6\), the components of the solution to the boundary system are of different magnitude: \(\sigma _1\sim \mathcal {O}(h^7)\) while the other components are of order \(\mathcal {O}(h^5)\). Therefore, the boundary truncation error only leads to \(\mathcal {O}(h^7)\) interior error. The large error \(\mathcal {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.

Error plot for the. a Dirichlet and. b Neumann problem with the sixth order accurate scheme

Full size image

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 \(\varOmega =[0,1]\), and the grid size ratio is 2:1. The step size in time is \(\varDelta 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.

Convergence tests for. a a grid interface problem. b a Neumann problem with schemes with added truncation errors

Full size image

According to the accuracy analysis in Sect. 4, for second and fourth order accurate methods with \(\tau =\tau _{2p}\), 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 \(\tau =\tau _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.

5.2 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 \(\mathcal {F}\) to be \((10\pi )^3\) and all the other components zero. We use the number \((10\pi )^3\) because the truncation error involves the third derivative of the true solution with the factor \((10\pi )^3\). The boundary truncation error of (44) has the same structure as in the original scheme, but increases to \(\mathcal {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 \(\mathcal {F}\) be \(h(10\pi )^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 \(\mathcal {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.

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.