Convergence of finite difference methods for the wave equation in two space dimensions

TL;DR

This paper extends normal mode analysis to two-dimensional wave equations solved by finite difference methods, providing a framework to understand boundary truncation errors and their impact on convergence, supported by numerical experiments.

Contribution

It introduces a general framework for analyzing convergence of finite difference methods for 2D wave equations, including boundary and corner truncation errors, extending previous 1D analyses.

Findings

Boundary truncation errors affect overall convergence.

The analysis applies to second order wave equations.

Numerical results confirm theoretical predictions.

Abstract

When using a finite difference method to solve an initial--boundary--value problem, the truncation error is often of lower order at a few grid points near boundaries than in the interior. Normal mode analysis is a powerful tool to analyze the effect of the large truncation error near boundaries on the overall convergence rate, and has been used in many previous literatures for different equations. However, existing work only concerns problems in one space dimension. In this paper, we extend the analysis to problems in two space dimensions. The two dimensional analysis is based on a diagonalization procedure that decomposes a two dimensional problem to many one dimensional problems of the same type. We present a general framework of analyzing convergence for such one dimensional problems, and explain how to obtain the result for the corresponding two dimensional problem. In particular, we consider two kinds of truncation errors in two space dimensions: the truncation error along an entire boundary, and the truncation error localized at a few grid points close to a corner of the computational domain. The accuracy analysis is in a general framework, here applied to the second order wave equation. Numerical experiments corroborate our accuracy analysis.

Tables

Table 1. Table 1: Result of accuracy analysis for the one dimensional problem. 2 ​ p : : 2 𝑝 absent 2p: order of the interior truncation error. p : : 𝑝 absent p: order of the boundary truncation error. q = min ⁡ ( 2 ​ p , p + Gain ) : : 𝑞 2 𝑝 𝑝 Gain absent q=\min(2p,p+\text{Gain}): overall convergence rate. ι , ι 0 : : 𝜄 subscript 𝜄 0 absent \iota,\iota_{0}: energy stability is obtained if the penalty parameter ι 𝜄 \iota satisfies ι ≥ ι 0 𝜄 subscript 𝜄 0 \iota\geq\iota_{0} .

	$2p=2$		$2p=4$		$2p=6$
	$p+$Gain	$q$	$p+$Gain	$q$	$p+$Gain	$q$
Dirichlet $\left(\iota ={\iota }_0\right)$	1+0.5	1.5	2+0.5	2.5	3+0.5	3.5
Dirichlet $\left(\iota >{\iota }_0\right)$	1+2	2	2+2	4	3+2.5	5.5
Neumann	1+1	2	2+2	4	3+2.5	5.5

Equations

$$U_{tt}=U_{xx}+U_{yy}+F,$$

$$0\le x<\infty ,0\le y\le 1,0\le t\le t_f,$$

$$\begin{array}{rl}& x_i=(i-1)h,i=1,\cdots ,\\ & y_j=(j-1)h,j=1,\cdots ,N_y,\end{array}$$

$$\Vert w{\Vert }_{2D}^2=h^2\sum _{i=1}^{N_y}\sum _{j=1}^{\infty }|w_{i,j}{|}^2.$$

$$\Vert w_{i,:}{\Vert }_{1D,x}^2=h\sum _{j=1}^{\infty }|w_{i,j}{|}^2,\Vert w_{:,j}{\Vert }_{1D,y}^2=h\sum _{i=1}^{N_y}|w_{i,j}{|}^2.$$

$$u_{tt}=\frac{1}{h^2}(Q_x\otimes I_y)u+\frac{1}{h^2}(I_x\otimes Q_y)u+F_h.$$

$$\frac{1}{h^2}{Q}_y\phi =-\lambda \phi ,$$

$${\zeta }_{tt}=\frac{1}{h^2}(Q_x\otimes I_y)\zeta +\frac{1}{h^2}(I_x\otimes Q_y)\zeta +h^{p}T,$$

$$T=[T_0;T_1;\cdots ,T_k;0;0;0;\cdots ],$$

$$s^2\hat{\zeta }=\frac{1}{h^2}(Q_x\otimes I_y)\hat{\zeta }+\frac{1}{h^2}(I_x\otimes Q_y)\hat{\zeta }+h^p\hat{T},$$

$$s^2\hat{\zeta }=\frac{1}{h^2}(Q_x\otimes I_y)\hat{\zeta }-(I_x\otimes \Phi \Lambda {\Phi }^{-1})\hat{\zeta }+h^p\hat{T},$$

$$s^2(I_x\otimes {\Phi }^{-1})\hat{\zeta }=\frac{1}{h^2}(Q_x\otimes {\Phi }^{-1})\hat{\zeta }-(I_x\otimes \Lambda {\Phi }^{-1})\hat{\zeta }+h^p(I_x\otimes {\Phi }^{-1})\hat{T}.$$

$$s^2\widehat{ϵ}=\frac{1}{h^2}(Q_x\otimes I_y)\widehat{ϵ}-(I_x\otimes \Lambda )\widehat{ϵ}+h^p(I_x\otimes {\Phi }^{-1})\hat{T},$$

$$\underset{({\tilde{s}}_{+}^{(r)}{)}^2}{\underbrace{({\tilde{s}}^2+h^2{\lambda }^{(r)})}}{\widehat{ϵ}}^{(r)}=Q_x{\widehat{ϵ}}^{(r)}+h^{p+2}{\hat{\tau }}^{(r)},$$

$${\hat{\tau }}^{(r)}=[{\hat{\tau }}_0^{(r)},{\hat{\tau }}_1^{(r)},\cdots ,{\hat{\tau }}_k^{(r)},0,0,0,\cdots {]}^T.$$

$$\Vert \widehat{ϵ}{\Vert }_{2D}^2=h\sum _{r=1}^{N_y}\Vert {\widehat{ϵ}}^{(r)}{\Vert }_{1D,x}^2.$$

$$\Vert {\widehat{ϵ}}^{(r)}{\Vert }_{1D,x}^2\le \frac{K{h}^{2g}}{{\eta }^m}\sum _{i=0}^k\left(|{\hat{\tau }}_i^{(r)}{|}^2+{\left|\widehat{\frac{{\partial }^b{\tau }_i^{(r)}}{\partial t^b}}\right|}^2\right),$$

$${\int }_0^{t_f}\Vert \zeta {\Vert }_{2D}^{2}dt\le \frac{K{e}^{2\eta t_f}{h}^{2g}}{{\eta }^m}{\int }_0^{t_f}{\chi }^2(\Phi )\sum _{i=0}^k\left(\Vert T_i{\Vert }_{1D,y}^2+{\Vert \frac{{\partial }^b{T}_i}{\partial t^b}\Vert }_{1D,y}^2\right)dt.$$

$$\Vert \widehat{ϵ}{\Vert }_{2D}^2\le \frac{K{h}^{2g}}{{\eta }^m}\sum _{i=0}^k\left(\Vert {\tau }_i{\Vert }_{1D,y}^2+{\Vert \widehat{\frac{{\partial }^b{\tau }_i}{\partial t^b}}\Vert }_{1D,y}^2\right),$$

$${\int }_0^{t_f}{e}^{-2\eta t}\Vert ϵ{\Vert }_{2D}^{2}dt\le \frac{K{h}^{2g}}{{\eta }^m}{\int }_0^{t_f}{e}^{-2\eta t}\sum _{i=0}^k\left(\Vert {\tau }_i{\Vert }_{1D,y}^2+{\Vert \frac{{\partial }^b{\tau }_i}{\partial t^b}\Vert }_{1D,y}^2\right)dt,$$

$${\int }_0^{t_f}{e}^{-2\eta t}\left(\Vert {\tau }_i{\Vert }_{1D,y}^2+{\Vert \frac{{\partial }^b{\tau }_i}{\partial t^b}\Vert }_{1D,y}^2\right)dt\le {\int }_0^{t_f}{e}^{-2\eta t}\Vert {\Phi }^{-1}{\Vert }^2\left(\Vert T_i{\Vert }_{1D,y}^2+{\Vert \frac{{\partial }^b{T}_i}{\partial t^b}\Vert }_{1D,y}^2\right)dt.$$

$$\begin{array}{rl}\Vert \Phi {\Vert }^2{\int }_0^{t_f}{e}^{-2\eta t}\Vert ϵ{\Vert }_{2D}^{2}dt& ={\int }_0^{t_f}{e}^{-2\eta t}\Vert \Phi {\Vert }^2\Vert (I_x\otimes {\Phi }^{-1})\zeta {\Vert }_{2D}^{2}dt\\ & \ge {\int }_0^{t_f}{e}^{-2\eta t}\Vert \zeta {\Vert }_{2D}^{2}dt.\end{array}$$

$${\int }_0^{t_f}{e}^{-2\eta t}\Vert \zeta {\Vert }_{2D}^{2}dt\le \frac{K{h}^{2g}}{{\eta }^m}{\int }_0^{t_f}{e}^{-2\eta t}{\chi }^2(\Phi )\sum _{i=0}^k\left(\Vert T_i{\Vert }_{1D,y}^2+{\Vert \frac{{\partial }^b{T}_i}{\partial t^b}\Vert }_{1D,y}^2\right)dt,$$

$$U_{tt}=U_{xx},0\le x<\infty ,t\ge 0,$$

$$u_{tt}=\frac{Q}{h^2}u,$$

$${\tilde{s}}^2\hat{\zeta }=Q\hat{\zeta }+h^{p+2}\hat{T},$$

$${\hat{T}}_i=a_i\frac{{\partial }^{p+2}}{\partial x^{p+2}}\hat{U}(0,s),i=0,\cdots ,k,$$

$$C(\tilde{s})\Sigma =h^{p+2}{\hat{T}}_C,$$

$$\Sigma =[{\hat{\zeta }}_1,{\hat{\zeta }}_2,\cdots ,{\hat{\zeta }}_d,{\sigma }_1,{\sigma }_2,\cdots ,{\sigma }_l{]}^T,$$

$$\hat{\zeta }=[{\hat{\zeta }}_1,{\hat{\zeta }}_2,\cdots ,{\hat{\zeta }}_d,\sum _{j=1}^l{\sigma }_j,\sum _{j=1}^l{\sigma }_j{\kappa }_j,\sum _{j=1}^l{\sigma }_j{\kappa }_j^2,\cdots {]}^T.$$

Full text

Keywords: Convergence rate, Accuracy, Two space dimensions, Normal mode analysis, Finite difference method, Second order wave equation

1 Introduction

Wave propagation problems can often be efficiently discretized with high order finite difference methods. Due to stability consideration, the formal accuracy order of the discretization scheme is typically considerably lower close to computational boundaries than that of the interior scheme. However, the numerical solution often converges at a rate higher than indicated by the boundary truncation error, a phenomenon termed as gain in convergence. This phenomenon can partly be understood from the fact that the number of grid points with the lower order stencil is independent of grid spacing. Analysis is needed to determine the precise order of gain in convergence.

There are two different methods for analyzing how much is gained in the convergence rate, the energy method and the normal mode analysis [10, 11]. Applying the energy method in a straightforward way indicates a half order gain in convergence. An exception is found in [1], where a careful energy analysis performed to the heat equation gives a gain of one and a half orders in convergence. However, the computations show a gain of two orders.

Sharp error estimates can be obtained by the normal mode analysis. In [28], it is noted by a normal mode approach that the gain in convergence can equal to the highest order of spatial derivatives in the equation. We refer to this gain as optimal. When a so–called determinant condition is satisfied, normal mode analysis can straightforwardly be used to prove that the gain in convergence is at least optimal. For first order hyperbolic equations, technical assumptions of finite difference schemes are given in [9] under which the determinant condition is satisfied, thus one order is gained in convergence. A detailed analysis in [24, 25] for a class of discretizations for the Schrödinger equation proves the gain is two orders, which is equal to the optimal gain. Results for the wave equation are presented in [30], where it is shown that the gain in convergence is not unified, but depends on boundary conditions and numerical boundary treatments. For both the Schrödinger equation and the wave equation, the determinant condition is not satisfied in many of the cases considered, even though the schemes are stable. We remark that the theoretical convergence rate obtained from the normal mode analysis is in a generalized sense, analogue to the concept of stability in the generalized sense [11, Chapter 12.3].

The accuracy analyses in the above mentioned references are limited to problems in one space dimension. However, the one dimensional analysis cannot always explain two dimensional numerical results. One example is two dimensional problems discretized in a multi–block setting, where numerical solutions in regions with different grid spacings are coupled using interpolation [13, 20, 26]. For cases like this, truncation errors located at grid points along the interface between two mesh blocks are of lower accuracy order compared to interior truncation errors. This is because of one–sided difference stencils on each side of the interface. At a few grid points on the edge of the interface, the accuracy of the truncation error is often one or two additional orders lower, caused by one–sided interpolations. We refer to such points as corner points, shown in Figure 1(A). The number of corner points depends on the particular discretization but is independent of grid spacing. Because the dominating truncation errors are localized in both spatial dimensions, it is a situation that does not occur in one dimensional problems. Numerical results in the literature indicate that such two dimensional cases may lead to higher convergence rates than what is predicted by the corresponding one dimensional analysis, see [14] for the advection equation, [24] for the Schrödinger equation and [31] for the wave equation. The numerical results also indicate that depending on the partial differential equation and the numerical interface treatment, the gain in convergence rates may be different. Accuracy analysis for two dimensional problems is needed in order to fully understand these results.

In this paper, we present a general accuracy analysis framework for semi–discrete partial differential equations in two space dimensions, to better understand the effect localized truncation errors have on the overall spatial discretization error. As a model problem, we consider the semi–discretization of the wave equation by finite difference methods. High order methods solve wave propagation problems more efficiently than low order methods on smooth domains [12, 15]. However, it is a challenging task to construct stable and high order accurate methods for wave equations in the presence of boundaries and interfaces. One way to do this is to combine the summation–by–parts (SBP) finite difference method with the simultaneous approximation term (SAT) to impose boundary conditions. The SBP–SAT finite difference method has been successfully used to solve many types of differential equations numerically, and is our choice for the spatial discretization of the wave equation. However, the technique we develop to analyze accuracy is not limited to this class of methods.

Since the equation is linear, we analyze the interior truncation error and boundary truncation error separately. The interior truncation error can be analyzed straightforwardly by the energy method. For the lower order boundary truncation error, we use again the superposition principle and analyze separately the effect of truncation errors along different parts of boundaries. Therefore, in the analysis we will focus on the truncation error along one boundary $x=0$ and consider two cases shown in Figure 1(B) and 1(C). In Figure 1(B), the dominating truncation error is located at only a few grid points on a boundary, i.e. localized in two space dimensions. This is a simplified model of the case shown in Figure 1(A). In addition, we also consider the case shown in Figure 1(C), which corresponds to when the dominating truncation error is located along an entire boundary. To further simplify the analysis, we will consider semi–infinite domains, where the boundaries at $x=1$ in Figure 1(B) and Figure 1(C) are moved to infinity. This separation of boundaries is justified by the same arguments which justify the separation in stability analysis [11].

The accuracy analysis relies on a transformation of the two dimensional problem to many one dimensional problems of the same type. Each one dimensional problem is then analyzed by the normal mode analysis, a technique involving Laplace transformation in time. We show how for a one dimensional problem the error is dictated by the behaviour of a boundary system in the origin of the $s$–plane, where $s$ is the Laplace dual variable of time. This was discussed previously in [24, 25, 30], and in this paper we make the arguments more precise. We also show how to use the one dimensional results to obtain results for both two dimensional cases.

The outline of the paper is as follows. In Section 2 we present the two dimensional accuracy analysis, as well as the accuracy analysis for the relevant one dimensional problems. Properties for the semi–discrete system discretized with SBP operators and weak treatment of boundary conditions using the SAT method are described in Section 3. Sections 4 and 5 contain numerical experiments and conclusions, respectively.

2 Accuracy analysis for the two dimensional wave equation

We consider the two dimensional wave equation

[TABLE]

on a domain

[TABLE]

with suitable initial and boundary conditions so that (1) is wellposed. We in particular consider the Dirichlet and Neumann boundary conditions. A comprehensive discussion on the wellposedness of second order hyperbolic equations is found in [16]. Because the focus in this paper is the accuracy analysis of numerical methods solving equation (1), we assume that the solution is sufficiently smooth and in L2 at any time. This gives rise to compatibility conditions between the initial and boundary data at the space–time corner. One way to guarantee that the compatibility conditions are satisfied is to assume that both spatial derivatives of the initial data on the boundary, and temporal derivatives of the boundary data at $t=0$, vanish up to sufficiently high order. We take this approach in the present work.

The domain is discretized by an equidistant grid with a grid spacing $h$ in both spatial directions

[TABLE]

with $h=1/(N_{y}-1)$. As we will measure errors in both one dimensional and two dimensional spaces, we distinguish between three different discrete norms. Let $w$ be a two dimensional grid function on the grid (2). We define the discrete norm as

[TABLE]

We will also need norms of restrictions of grid functions to lines with constant $x$ or $y$ coordinates. They are denoted by

[TABLE]

The spatial derivatives are approximated by finite difference operators yielding the semi–discretization

[TABLE]

In equation (3), $u$ is a grid function approximating the true solution $U$ on the grid and is arranged column–wise, i.e. the first $N_{y}$ components of $u$ are the numerical solutions at the grid points on the boundary $x=0$. The finite difference operators $Q_{x}/h^{2}$ and $Q_{y}/h^{2}$ approximate second derivatives in space, including an implementation of the boundary conditions, and $I_{x}$, $I_{y}$ are identity operators. The Kronecker product $\otimes$ is used to extend the operators from one space dimension to two space dimensions. The grid function $F_{h}$ is the projection of the forcing function $F(x,y,t)$ on the grid. Inhomogeneous boundary data would also appear in the right–hand side of (3), but is excluded here for a sake of simplified notation as it has no effect on the accuracy analysis.

To consider accuracy, we need the semi–discretization (3) to be stable. One way to prove stability is to show by the energy method that the numerical solution in some appropriate norm is bounded by the data, and such a scheme is called energy stable. Another way is to prove stability in a generalized sense by the Laplace transform technique in the normal mode analysis framework, see details in [11, Chapter 12]. In particular, we consider operators $Q_{y}$ that satisfies the following property.

Assumption 1**.**

There exists a symmetric positive definite operator $P$ such that with $Q_{y}$ in (3) the product $PQ_{y}$ is symmetric negative semi–definite, and the spectral norms of $P^{1/2}$ and $P^{-1/2}$ are uniformly bounded.

The spectral norm, denoted by $\|\cdot\|$, is induced from the standard Euclidean vector norm. For any value $N_{y}$, operators in the $y$–direction can be represented by matrices. Assumption 1 on the operator $Q_{y}$ ensures that an energy estimate can be obtained when $Q_{y}$ is used as the spatial discretization operator in the corresponding one dimensional problem. This assumption can be satisfied by many kinds of discretizations for which a standard energy estimate can be obtained. One example is a finite difference operator satisfying a summation–by–parts (SBP) property with a weak imposition of boundary conditions and properly chosen penalty parameters [2, 22]. In this case, $H=hP$ is the operator associated with the SBP norm with the grid spacing $h$, and the condition number of $P^{1/2}$ is independent of $h$. The SBP finite difference method is discussed in more detail in Section 3.

The following lemma describes an important property of the discretization operator $Q_{y}$ that will be needed for the two dimensional accuracy analysis.

Lemma 1**.**

Consider the eigenvalue problem

[TABLE]

with $Q_{y}$ in (3). Under Assumption 1, $Q_{y}$ is diagonalizable $Q_{y}/h^{2}=-\Phi\Lambda\Phi^{-1}$ by $\Phi=[\varphi_{1},\varphi_{2},\cdots,\varphi_{N_{y}}]$, where $\Lambda$ is a diagonal matrix with real and non–negative diagonal entries. In addition, $\|\Phi\|=\|P^{-1/2}\|$ and $\|\Phi^{-1}\|=\|P^{1/2}\|$ are uniformly bounded with respect to $N_{y}$.

The proof of Lemma 1 can be found in Appendix A. As will be seen later, with such an operator $Q_{y}$ the error equation for the two dimensional problem can be transformed by a diagonalization technique to a number of one dimensional problems of the same type.

Let $U_{h}$ be the true solution $U$ projected on the grid, and the point–wise error be $\zeta(t)=U_{h}(t)-u(t)$. As discussed in the introduction, we will only analyze the effect of the truncation error, denoted $\mathcal{O}(h^{p})$, caused by the one–sided stencil at the boundary $x=0$. The error equation is

[TABLE]

where $h^{p}T$ is the boundary truncation error and $T$ is independent of $h$. Only the components of $T$ corresponding to the grid points located at $x=ih,\ i=0,\cdots k$ are nonzeros, where $k$ is a small constant independent of $h$. We therefore write

[TABLE]

where $h^{p}T_{i},\ i=0,\cdots k$, is the truncation error at the grid points located at $x=ih$. Here, $T$ is a two dimensional grid function, $T_{i}$ is a one dimensional grid function in the $y$–direction, and $\mathbb{0}$ is a zero vector of length $N_{y}$.

2.1 Diagonalization of the error equation

To begin with, we consider the case with the dominating truncation error along the entire boundary shown in Figure 1(C), for which $\|T_{i}\|_{1D,y}^{2}=\mathcal{O}(1)$. In Section 2.3, we consider the localized case in Figure 1(B) with $\|T_{i}\|_{1D,y}^{2}=\mathcal{O}(h)$.

The next step of the normal mode analysis is to perform a Laplace transform in time of (5)

[TABLE]

where $s$ is the time dual in the Laplace space. We diagonalize $Q_{y}$ as in Assumption 1 and Lemma 1, and rewrite (7) as

[TABLE]

where $\Lambda$ is diagonal with diagonal entries $\lambda^{(r)}\geq 0$ in ascending order for $r=1,2,\cdots,N_{y}$. Multiplying the above equation by $(I_{x}\otimes\Phi^{-1})$ from the left, we obtain

[TABLE]

With the notation $\hat{\epsilon}=(I_{x}\otimes\Phi^{-1})\hat{\zeta}$, (8) becomes

[TABLE]

with the operator in the $y$–direction diagonalized. This is the spectrally decomposed form, which consists of $N_{y}$ scalar difference equations

[TABLE]

where $r=1,2,\cdots,N_{y}$ and $\tilde{s}=sh$. For every $r$, we have

[TABLE]

Note the close relation between $\hat{\tau}^{(r)}$ and $\hat{\tau}_{i}=\Phi^{-1}\hat{T}_{i}$: the $i^{th}$ entry $\hat{\tau}^{(r)}_{i}$ is the same as the $r^{th}$ entry of $\hat{\tau}_{i}$. In addition, $\hat{\epsilon}$ in (9) is related to $\hat{\epsilon}^{(r)}$ in (10) by

[TABLE]

With the notation $(\tilde{s}^{(r)}_{+})^{2}=\tilde{s}^{2}+h^{2}\lambda^{(r)}$, we have transformed the two dimensional error equation to $N_{y}$ one dimensional error equations in the Laplace space. This transformation can be understood as a variant of Fourier transform. A general result on the estimate of the error $\zeta$ for the two dimensional problem is stated in the following theorem.

Theorem 1**.**

If for all $s$ with Re$(s)=\eta>0$, $\hat{\epsilon}^{(r)}$ in (10) is bounded as

[TABLE]

with $\eta,g,m$ and $b$ independent of $h$, then

[TABLE]

Here, $\chi(\Phi)$ is the condition number of $\Phi$.

We note that in practice, to validate (12) we do not need to check for every $r$. It suffices to check the corresponding one dimensional problem, that is (12) without any shift when $h^{2}\lambda^{(r)}=0$. To see this, we need the following lemma describing an important property of $\tilde{s}^{(r)}_{+}$.

Lemma 2**.**

Let $\tilde{s}_{+}=\sqrt{\tilde{s}^{2}+\gamma}$ for some real $\gamma\geq 0$. If Re$(\tilde{s})\geq\delta\geq 0$, then Re$(\tilde{s}_{+})\geq\delta$.

The proof of Lemma 2 is given in Appendix B. This lemma implies that each error equation (10) can be seen as the corresponding one dimensional error equation perturbed in the direction to the right of the complex plane. Therefore, the estimate of (10) for $r=1,2,\cdots,N_{y}$ is no worse than the estimate for the corresponding one dimensional problem.

We will show in Section 2.2 how an estimate of the type (12) can be derived by normal mode analysis for a one dimensional problem of the type (10). In Theorem 1 and the rest of the paper, we use the capital letter $K$ in the estimates to denote some constant independent of the grid spacing $h$, where the precise value of $K$ could be different from one estimate to another. Below we give the proof of Theorem 1.

Proof of Theorem 1.

A sum of (12) in $r$ leads to

[TABLE]

and Parseval’s relation gives

[TABLE]

where $\eta$ is a constant independent of $h$, $\epsilon=(I_{x}\otimes\Phi^{-1})\zeta$ and $\tau_{i}=\Phi^{-1}T_{i}$. It is obvious by the property of the induced norm that

[TABLE]

In addition, we also have

[TABLE]

Finally, we obtain the following estimate by (14), (15) and (16),

[TABLE]

where $\chi(\Phi)=\|\Phi\|\|\Phi^{-1}\|$ is the condition number of $\Phi$, and by Lemma 1 it is uniformly bounded. The estimate (13) is obtained by multiplying (17) with $e^{2\eta t_{f}}$ on both sides, by using $e^{-2\eta t}<1$ and $e^{2\eta(t_{f}-t)}>1$ for $0<t<t_{f}$. ∎

Note that by Lemma 1, $\chi(\Phi)$ is uniformly bounded. Theorem 1 indicates that the convergence rate of a two dimensional problem is at least equal to that of the corresponding one dimensional problem. A special case for a two dimensional problem is when the dominating truncation error is localized also in the $y$–direction, i.e. the number of grid points with the dominating truncation error is independent of $h$, yielding $\|T_{i}\|_{1D,y}^{2}=\mathcal{O}(h),i=0,\cdots,k$, in (13). Straightforwardly, this leads to an additional half order gain. A more detailed analysis in Section 2.3 shows that the gain for a two dimensional problem with a localized truncation error could be a full order higher compared with the gain for a corresponding one dimensional problem.

2.2 Background: accuracy analysis for one dimensional problems

We would like to explore the relation between the one dimensional error equation (10) and the corresponding estimate (12). Therefore, accuracy analysis for the corresponding one dimensional problem is needed, and is presented in this section. Part of the analysis can also be found in [30], where it is argued that the error estimate in Laplace space is determined by its behavior at the origin. In the paper, we give a proof of that.

We consider the one dimensional wave equation

[TABLE]

with appropriate initial and boundary conditions so that the problem is well–posed and the true solution is in L2 and smooth. Similar to the two dimensional case, (18) is discretized by a finite difference operator and the boundary conditions are imposed so that the semi–discretization is energy stable. The semi–discretization is

[TABLE]

and the error equation in the Laplace space is

[TABLE]

where $s$ is the time dual in Laplace space and $\tilde{s}=sh$. Here $u$ and $\hat{\zeta}$ are vectors containing the numerical solution and point–wise error in the Laplace space, respectively. The operator $Q/h^{2}$ approximates the second derivative in space including an implementation of the boundary condition. The vector $\hat{T}$ has only a few non–zero components corresponding to the large truncation error close to the boundary $x=0$. Since the true solution is smooth, the non–zero components of $\hat{T}$ are to the leading order on the form

[TABLE]

where $a_{i}$ can be obtained by the Taylor expansion. The constants $a_{i}$ and $k$ are determined by the precise form of $Q$, and are independent of $h$.

Remark 1**.**

Equation (20) is in the same form as (10), with $\tilde{s}^{2}$ in the left–hand side of (20) and a perturbed $(\tilde{s}_{+}^{(r)})^{2}$ in the left–hand side of (10). We keep this in mind in the following analysis.

The difference equation (20) can be solved by first considering the components corresponding to the grid points away from the boundary where the forcing $h^{p+2}\hat{T}$ is zero and $Q/h^{2}$ is the standard central finite difference stencil. This step gives the characteristic equation with roots being functions of $\tilde{s}$. Since the problem under consideration is a half–line problem, in the estimate we only need to include the admissible roots $\kappa(\tilde{s})$, for which $|\kappa(\tilde{s})|<1$ for all Re$(\tilde{s})>0$. With a $2l^{th}$ order central finite difference stencil, there are $l$ admissible roots $\kappa_{1},\cdots,\kappa_{l}$.

The next step in solving (20) is to consider its first few rows, corresponding to the discretization close to the boundary, where $h^{p+2}\hat{T}$ does not vanish. The corresponding equations can be written in a matrix–vector multiplication form

[TABLE]

which is referred to as the boundary system. The matrix $C(\tilde{s})$ is called the determinant matrix and depends on the boundary stencil, and is often of small size. For example, for the second order discretization in the SBP–SAT framework, $C(\tilde{s})$ is a 3–by–3 matrix with the Dirichlet boundary condition, and a scalar with the Neumann boundary condition. For the precise form of these boundary systems, see equation (25) and (37) in [30], respectively. For higher order discretizations, the dimension of $C(\tilde{s})$ is in general $k+1$. The first $k+1$ components of $\hat{T}_{C}$ is taken from $T_{i},i=0,\cdots,k$ in (21), with zeros appending to the end if needed. The unknown vector $\Sigma$ is

[TABLE]

where $d$ depends on the numerical boundary treatment.

We can express the general solution of the difference equation (20) by using (23) as

[TABLE]

Then the L2 norm of $\hat{\zeta}$ is

[TABLE]

There are two steps to derive an estimate of (25). One step is to estimate $1/(1-|\kappa_{j}|^{2})$, and is derived in Section 2.2.2. In addition, we also need to estimate the other terms in (25), i.e. the components of $\Sigma$:

[TABLE]

As will be seen later, when $\tilde{s}=0$ the characteristic equation always has an admissible root equal to 1, denoted by $\kappa_{1}=1$. Therefore, we consider Re$(\tilde{s})=\eta h=\mathcal{O}(h)$ in the limit as $h$ approaches zero. The solution to (22) can be written as

[TABLE]

because of the following lemma.

Lemma 3**.**

If the discretization (19) is stable, then $C(\tilde{s})$ in the boundary system (22) is non–singular for all Re$(\tilde{s})>0$.

The proof of Lemma 3 is very similar to the proof of Lemma 12.1.1 in [11, pp. 378], but for completeness we include it in Appendix C.

In the following section, we derive estimates for each component of $\Sigma$, denoted by $\Sigma_{i}$, where $i=1,\cdots,l+d$. For simplified notation, we use the maximum norm of vectors

[TABLE]

where $\|\Sigma\|_{\max}=\max_{i}|\Sigma_{i}|$. We emphasize again that the dimension of the vector $\Sigma$ in (27) is finite and is independent of the grid spacing. In certain cases, we do not use the last inequality in (28) because some components of $\Sigma$ may vanish and computing $C^{-1}(\tilde{s})\hat{T}_{C}$ directly gives a sharper estimate.

2.2.1 Estimates of $|\Sigma_{i}|$

When $\tilde{s}$ is away from the imaginary axis, $C(\tilde{s})$ is non–singular by the energy stability. In many common stable semi–discretizations of the second order wave equation, the determinant matrix $C(\tilde{s})$ is singular for some $\tilde{s}$ on the imaginary axis. We will show below that singularities away from the origin typically have no influence on the order of accuracy, and can be handled separately. A singularity at $\tilde{s}=0$ on the other hand can introduce a loss of accuracy compared with the optimal gain. The following lemma makes the arguments precise.

Lemma 4**.**

Consider the boundary system (22) at $\tilde{s}=i\tilde{\xi}+\eta h$ where $\eta>0$ is a constant independent of $h$. Let $\delta$ be a small but $h$–independent constant.

• •

When $\tilde{s}$ is in the vicinity of the origin, i.e. $|\tilde{s}|\leq\delta$, there exists a non–negative integer $w$ such that

[TABLE]

If $C(0)$ is nonsingular, then $w=0$. Otherwise, $w$ can be determined by the precise form of the boundary system (22).

• •

When $\tilde{s}$ is away from the origin, i.e. $|\tilde{s}|>\delta$, there exists a non–negative integer $\alpha$ such that

[TABLE]

The value of $\alpha$ can be determined by the precise form of the boundary system (22).

The proof of Lemma 4 is found in Appendix D.1.

If the boundary system (22) is singular at the origin, then we need to use (29) to estimate $|\Sigma_{i}|$ leading to a $w^{th}$ order loss. If the singularity occurs on the imaginary axis away from the origin, (30) introduces no accuracy order loss, but temporal derivatives of the truncation error appear in the estimate of $|\Sigma_{i}|$ . As an example, with Neumann boundary conditions we have $w=1$ for a second order method and $w=0$ for a fourth order method, see [30].

2.2.2 Estimates of $1/(1-|\kappa_{j}|^{2})$

The characteristic equation and its solution only depend on the finite difference stencil in the interior of the computational domain. To discretize the wave equation (18), it is common to use the standard central finite difference stencil as the interior stencil. In this case, we have the following lemma for the roots.

Lemma 5**.**

Consider the characteristic equation corresponding to the difference equation $\tilde{s}^{2}\zeta=\tilde{Q}\zeta$, where $\tilde{Q}$ is a standard central difference operator approximating the second derivative. Let $\tilde{s}=i\tilde{\xi}+\eta h$ with a constant $\eta>0$ independent of $h$. The roots of the characteristic equation have the following properties.

• •

When $\tilde{s}$ is in the vicinity of the origin, there is one admissible root $\kappa_{1}(\tilde{s})$ to the characteristic equation satisfying

[TABLE]

The constant $K$ depends on the finite difference stencil, but not on $h$ or $\eta$. For the other admissible roots, $\frac{1}{1-|\kappa(\tilde{s})|^{2}}$ is bounded independently of $h$.

• •

When $\tilde{s}$ is away from the origin, there can be admissible roots satisfying

[TABLE]

for some $\beta>0$. Such roots do not introduce any accuracy order loss when estimating (25).

The proof of Lemma 5 can be found in Appendix E.1. When $\tilde{s}$ is in a vicinity of the origin, the bound on $1/(1-|\kappa_{1}(\tilde{s})|^{2})$ in (31) leads to a $h^{-1}$ factor in (25). When $\tilde{s}$ is away from the origin, $1/(1-|\kappa_{j}(\tilde{s})|^{2})$ does not introduce any $h$–dependent factor in (25). Together with Lemma 4, we therefore conclude that the order of accuracy is determined by the error in the vicinity of the origin in Laplace space.

2.2.3 Final estimates in the Laplace space

The final error estimate is determined by the combined effect of the estimates (29) and (31) from Lemma 4 and 5, respectively. In the worst scenario, we have

[TABLE]

which leads to

[TABLE]

where the constant $K$ is independent of $h$. A more detailed analysis might reveal that $\sigma_{1}/(1-|\kappa_{1}|^{2})$ satisfies a sharper bound than given in (33). The value $b=2\alpha+\beta$ in (34) is given by $\alpha$ in (30) and $\beta$ in (32). Comparing with the estimate (12) in Theorem 1 we have $2g=2p+4-2w$ and $m=1+2w+b$.

2.3 Accuracy analysis for two dimensional problems with corner truncation errors

We now consider a two dimensional problem where the large truncation error along a boundary is only localized at a few grid points at a corner, see Figure 1(B).

Without loss of generality we consider only non-zero elements in $T_{0}$, i.e. the truncation error has the form (6) with $k=0$ and $T_{0}=[\mathcal{O}(1),0,\cdots]^{T}$. An immediate consequence is that $\|T_{0}\|_{1D,y}\sim\mathcal{O}(h^{1/2})$, which together with the uniform boundedness of $\|\Phi^{-1}\|$ in Lemma 1, leads to

[TABLE]

A direct application of Theorem 1 to this problem thus leads to an additional gain of a half order in convergence, compared with when the truncation error is present along the entire boundary. However, numerical experiments for certain problems show another half order better result. Below we demonstrate how to sharpen the analysis so that a full order gain in convergence is obtained compared with the case when the dominating truncation error is located along an entire boundary.

One difference between the $r^{th}$ error equation (10) and the standard one dimensional error equation (20) is that the coefficient $\tilde{s}^{2}$ in (20) is replaced by $(\tilde{s}^{(r)}_{+})^{2}=\tilde{s}^{2}+h^{2}\lambda^{(r)}$ in (10). We can therefore view $h^{2}\lambda^{(r)}$ as a perturbation to the error equation, which enables us to use the one dimensional results presented in Section 2.2. To obtain a sharp estimate for a two dimensional problem when the large truncation error is located only at a few grid points, we cannot rely on a uniform estimate in $r$ as in (12). Instead, the discrete eigenvalues $\lambda^{(r)}$ are divided into two sets according to indices $r=1,\cdots,r_{\delta}$ and $r=r_{\delta}+1,\cdots,N_{y}$, where $r_{\delta}$ is such that

[TABLE]

for some constant $\delta>0$ independent of $h$. Since we have a consistent spatial discretization the discrete eigenvalues in equation (4), $\lambda^{(r)}$, converge to the eigenvalues, $\lambda_{c}$, of the corresponding continuous problem

[TABLE]

where (38) is closed with the same boundary conditions as the ones approximated by the discrete operator in (4). Therefore, $h^{2}\lambda^{(r)}$ for a particular $r$ decreases when $h$ decreases, and thus $r_{\delta}$ grows with grid refinement.

It is helpful to first consider a specific example. Let

[TABLE]

which corresponds to the standard second order accurate discretization of the Laplace operator in one dimension with Neumann boundary conditions. The operator $Q_{y}/h^{2}$ is diagonalized by the standard unitary cosine transform, for which we have the operator $\Phi$ defined by

[TABLE]

and the discrete eigenvalues are given by

[TABLE]

Note that $\lambda^{(r)}$, converge to the eigenvalues,

[TABLE]

of (38) with Neumann boundary conditions. Furthermore, from the cosine transform of $(1,0,\cdots)^{T}$, we immediately see that

[TABLE]

Here we note the agreement with the uncertainty principle [27], which states that the more concentrated a function is, the more spread out its Fourier transform must be. In other words, locality in $\hat{T}_{0}$ rules out locality in $\hat{\tau}_{0}$.

We shall now consider a more general case, where $Q_{y}/h^{2}$ is a consistent spatial discretization satisfying Assumption 1. An analytical formula of $\lambda^{(r)}$ may not exist, but the eigenvalues of the discrete operator converge to the eigenvalues of the corresponding continuous problem (38) with the same boundary conditions. To ensure a corresponding locality principle in the more general case, we make the following assumption on the spatial discretization in the $y$–direction.

Assumption 2**.**

The diagonalized transformation in Lemma 1 satisfies the uncertainty principle in the sense that when $T_{0}=[\mathcal{O}(1),0,\cdots]^{T}$ then each component of $\hat{\tau}_{0}=\Phi^{-1}T_{0}$ satisfies

[TABLE]

By (35) we know that $\|\hat{\tau}_{0}\|_{1D,y}$ is proportional to $h^{1/2}$. Assumption 2 implies a component–wise bound (41) of $\hat{\tau}_{0}^{(r)}$, and rules out the possibility of a local character of the truncation error also in the spectral representation. Assumption 2 is a reasonable assumption because when the operator $\Phi^{-1}$ acts on $\hat{T}_{0}$, it is an analogue of performing a discrete Fourier transform. It is a discrete Fourier transform if the boundary condition in the $y$ direction is periodic. In the following analysis, we use the component–wise bound for $r\leq r_{\delta}$.

The total error for the two dimensional Laplace–transformed and spectrally decomposed problem is given by

[TABLE]

Note the similarity between (42) and the expression for $\|\hat{\zeta}\|_{2D}^{2}$ in (25). To obtain sharp accuracy results we divide (42) into

[TABLE]

and

[TABLE]

and estimate them separately. For $r\leq r_{\delta}$, by the component–wise bound (41) an additional $1/2$ order is gained in (29). We have

[TABLE]

where $w$ is related to the invertibility of $C(\tilde{s})$ in Lemma 4. In

[TABLE]

the first term can be bounded as

[TABLE]

while the second term

[TABLE]

By Lemma 5, only one admissible root leads to an $h^{-1}$ factor in the error estimate. We therefore only consider $1/(1-|\kappa_{1}^{(r)}|^{2})$ in (47). Together with $(\tilde{s}_{+}^{(r)})^{2}=\tilde{s}^{2}+h^{2}\lambda^{(r)}$, we obtain

[TABLE]

With a common boundary condition like Dirichlet, Neumann or periodic, the continuous eigenvalues $\sqrt{\lambda_{c}}$ in (38) are uniformly distributed. For small $r$ corresponding to the well–resolved components, the discrete eigenvalues converge to the continuous ones. As a consequence, $\sqrt{\lambda^{(r)}}\approx Kr$ for $r<r_{\delta}$. This can also be seen from the formula (39). For $\eta=0$ and $\sqrt{\lambda^{(r)}}\approx Kr$, we have

[TABLE]

where $r_{\delta}h=\delta$ is a constant independent of $h$ and $\gamma\approx 0.58$ is the Euler–Mascheroni constant. This leads to

[TABLE]

Let us now consider $r>r_{\delta}$. In this case, the shift in $\tilde{s}_{+}^{(r)}$ is $h^{2}\lambda^{(r)}>\delta^{2}\sim\mathcal{O}(1)$. Therefore, all $\kappa_{j}^{(r)}$ have absolute values bounded away from 1, which leads to the estimate $1/(1-|\kappa_{j}^{(r)}|^{2})\leq K$ for a constant $K$ independent of $h$. We then have the estimate

[TABLE]

Because of the shift in $\tilde{s}_{+}^{(r)}$, the determinant matrix $C(\tilde{s}_{+}^{(r)})$ is nonsingular by stability. Hence, $\max_{r}\|C^{-1}(\tilde{s}_{+}^{(r)})\|^{2}_{\max}$ is a constant of order $\mathcal{O}(1)$.

Together with (50), we have the combined estimate

[TABLE]

and consequently

[TABLE]

In the above estimate, we use the uniform boundedness of $\|\Phi\|$ in Lemma 1, and $K$ depends on the true solution but not on $h$. The factor $\log(h^{-1})$ grows with mesh refinement, but much slower than $h^{-1}$. Asymptotically, only the exponent of $h$ in (51) is important. Comparing the two dimensional case with a corner truncation error with the one dimensional problem, we can expect a full order gain in convergence from (34) to (51).

3 Summation–by–parts finite difference methods

High order finite difference methods solve wave propagation problems more efficiently than low order methods on smooth domains [12, 15]. Though standard central finite difference stencils can be used in the interior of the domain, it is challenging to derive stable and accurate schemes close to boundaries and interfaces. An approach that has been successfully used to overcome this difficulty is finite difference operators satisfying a summation–by–parts (SBP) property [17, 18] in combination with the simultaneous–approximation–term (SAT) technique [3] to impose boundary and interface conditions. A main advantage of the SBP–SAT method is its unified procedure to construct both provably stable and highly accurate schemes for linear time–dependent problems with boundaries and interfaces. Recent developments also include its use on non–uniform grids [5], and its relation to a form of discontinuous Galerkin method [8]. Introductions of the SBP–SAT finite difference method can be found in [4, 29].

3.1 Semi–discretizations by the SBP operators

To approximate a second derivative, we use SBP operators constructed in [23]. They are central finite difference stencils in the interior of the computational domain and special one–sided stencils at a few grid points near boundaries. The operator $D\approx\partial^{2}/\partial x^{2}$ is called an SBP operator if it can be decomposed as

[TABLE]

where $H$ is a diagonal positive definite operator associated with an L2–equivalent norm, the operator $M$ is symmetric positive semi–definite and $B$ takes the form $\text{diag}(-1,0,0,\cdots,1)$. The first and last row of $S$ approximate the first derivative at the boundary grid points.

The truncation error of the SBP operators in [23] is $\mathcal{O}(h^{2p})$ in the interior and $\mathcal{O}(h^{p})$ near boundaries, $p=1,2,3,4$. Though the accuracy is sacrificed near boundaries in order to satisfy the SBP property, the operators are often termed as $2p^{th}$ order accurate.

In the following, we present the SBP–SAT schemes for the wave equation with Dirichlet and Neumann boundary conditions. For comparison, we also state the result of the accuracy analysis for one dimensional problems, which was derived in [30].

3.1.1 Dirichlet boundary conditions

The semi–discretization of the one dimensional wave equation (18) with a Dirichlet boundary condition is

[TABLE]

where $E_{0}$ picks up $u$ at the boundary with entries zero almost everywhere except $E_{0}(1,1)=1$. The vector $\bar{g}=[g(t),0,\cdots]^{T}$ contains the boundary data and $F_{g}$ is the restriction of the forcing $F(x,t)$ onto the grid. The constant $\iota$ is independent of $h$ but must satisfy $\iota\geq\iota_{0}$ for stability [2], and its lower bound $\iota_{0}$ is presented in [21, 22, 30]. A stability proof of (52) by the energy method can be found in [22]. We also note that if the boundary condition of equation (1) at $y=0,1$ is Dirichlet and the discretization (52) is used also in the $y$ direction, then the operator $Q_{y}$ in (3) satisfies Assumption 1 with $P=h^{-1}H$.

3.1.2 Neumann boundary conditions

For the wave equation (18) with a Neumann boundary condition, it is natural to impose the boundary condition weakly. The corresponding SBP–SAT scheme is

[TABLE]

and a stability proof can be found in [22]. We note that when equation (1) has a Neumann boundary condition at $y=0,1$ and the discretization (53) is used also in the $y$ direction, the operator $Q_{y}$ in (3) satisfies Assumption 1 with $P=h^{-1}H$.

3.1.3 Accuracy analysis for the one dimensional wave equation

The accuracy analysis of both the discretizations (52) and (53) are performed in detail in [30], and we summarize the main result in Table 1.

For the Dirichlet problem, the penalty parameter $\iota$ in (52) has an influence on the convergence rate. If $\iota=\iota_{0}$, the determinant condition is not satisfied, and the gain in convergence is only a half order; if $\iota>\iota_{0}$, the determinant condition is satisfied and $w=0$ in the estimate (34), leading to the optimal two orders gain. In fact, super–convergence is obtained with the sixth order accurate scheme, where the gain in convergence is two and a half orders. The reason of super–convergence is found by a careful analysis of the boundary system (29), which shows that $\sigma_{1}\sim h^{p+3}$ in (33).

For the Neumann problem, the determinant condition is always violated. A careful analysis in [30] shows that the gain in convergence is 1, 2, and 2.5 for the second, fourth and sixth order accurate schemes, respectively. For the fourth and sixth order accurate schemes, the gain in convergence is at least the optimal two orders. This is because the boundary system in the normal mode analysis has a special structure, which leads to $w=0$ in the estimate (34). We refer to a detailed discussion in [30]. We also remark that if the structure of the boundary system is perturbed by modifying the truncation error, the gain in convergence is one order in all three cases.

By Theorem 1, for the SBP–SAT approximation of the Dirichlet and Neumann problem in two space dimensions, the convergence rates are the same as the corresponding one dimensional cases shown in Table 1.

4 Numerical experiments

In this section, we perform numerical experiments to verify our accuracy analysis, with a focus on two different kinds of truncation errors for the wave equation in two space dimensions. The first case is when the large truncation error is located on an entire boundary, which typically occurs in a single–block domain or a multi–block domain with conforming grid interfaces. Such an experiment is presented in Section 4.1. In the second case, the truncation error is only located at a few corner points in a two dimensional domain and the number of such grid points is independent of grid refinement, see Section 4.2. This experiment is constructed as a simplified analogue to a multi–block setting with non–conforming grid interfaces presented in [31]. The dominating truncation error is located on a few corner points, see an illustration in Figure 1(A). More precisely, the $2p^{th}$ order accurate schemes have a dominating truncation error $\mathcal{O}(h^{p-2})$ because of interpolation.

4.1 A standard two dimensional problem

We consider the two dimensional wave equation (1) in the computational domain $[0,1]^{2}$ discretized on a uniform grid. The SBP finite difference operators constructed in [23] are used to approximate the spatial derivatives, and the SAT technique is used to weakly impose all the boundary conditions according to the discretizations (52) and (53). The semi–discretizations (52) and (53) are generalized to two space dimensions using the Kronecker product. We will consider both Dirichlet and Neumann problems, meaning that all boundary conditions are either Dirichlet or Neumann, respectively.

For hyperbolic problems, explicit methods are often used to advance the equation in time. In the following numerical experiments, we employ the fourth order Runge–Kutta method as the time integrator. In [19], it is shown that under reasonable conditions the fully discrete scheme is stable if the semi–discretization is stable and Runge–Kutta methods are used to discretize in time. Since we investigate the convergence rate in space, we choose the time step $\Delta t=0.1h$ small enough so that the temporal error is negligible compared with the spatial error. This step size is well below the stability limit. The final time of the simulations is chosen to be $t=2$.

The numerical solution is compared with an analytic solution constructed using the method of manufactured solutions:

[TABLE]

The L2 errors are computed as the norm of the difference between the exact solution projected onto the grid with grid spacing $h$, $U_{h}$, and the corresponding numerical solution, $u_{h}$, according to

[TABLE]

and the convergence rate is computed by

[TABLE]

In the numerical experiments, we observe that the computed L2 convergence rates are in agreement with the rates in the generalized sense obtained in the analysis.

In the SBP–SAT scheme, the penalty parameter for the Dirichlet problem, $\iota$, is chosen to be $20\%$ larger than its lower bound required by energy stability. It is important to choose $\iota$ larger than the threshold $\iota_{0}$ to obtain the desired rate of convergence, as the choice $\iota=\iota_{0}$ gives suboptimal convergence. This is demonstrated in [30]. We note, however, that an extremely large $\iota$ leads to a very small time step. An increase of $20\%$ of the penalty parameter seems appropriate based on our experiments.

The L2 errors versus the number of grid points in each spatial dimension, $N$, are plotted in Figure 2. The convergence rates computed from the last two mesh refinements for each set of experiments are shown at the end of the corresponding convergence curves in the error plots. According to Theorem 1, Table 1 and the analysis in [30], for both the Dirichlet and Neumann problem, the SBP–SAT approximations have theoretical convergence rates 2, 4 and 5.5 for the second, fourth and sixth order accurate schemes, respectively. As seen in Figure 2, the convergence rates for the corresponding two dimensional problem are 2 and 4 for the second and fourth order accurate methods, which agree well with our analysis. The convergence rate for the sixth order accurate method is 5.75 for the Dirichlet problem and 5.27 for the Neumann problem.

4.2 Truncation error located on a few corner points

The case analyzed in Section 2.3 is a simplified model of the multi–block finite difference discretization with non–conforming grid interfaces. Below we conduct numerical experiments for this simplified model to verify that the analysis in Section 2.3 is sharp.

We use the same setting as in Section 4.1, but in this case we consider the manufactured, smooth solution

[TABLE]

to the two dimensional wave equation (1), and solve until $t=2$ using the methods of interior order $2p=2,4,6$.

To get a truncation error $\mathcal{O}(h^{p-2})$ at a few corner points, as shown in Figure 1(B), we modify the standard SBP–SAT scheme by using erroneous boundary data. For the Dirichlet problem, the true boundary data at the grid point $(0,y_{i})$, denoted by $g_{i}(t)$, can be obtained from (55) as $g_{i}(t)=U(0,y_{i},t)$. However, in the penalty terms in the SAT method, we use $(1+\nu)g_{i}(t)$ on the first and last five grid points. In particular, for the $2p^{th}$ order method we choose $\nu=c_{p}h^{p},\ p=1,2,3$, where the factors $c_{p}$ are chosen to match the coefficients of the truncation error by one–sided difference stencils. Since the boundary data is multiplied by $H^{-1}S^{T}\sim\mathcal{O}(h^{-2})$, this choice results in a truncation error $\mathcal{O}(h^{p-2})$. We note that the stability property of the discretization does not change. We plot the L2 error versus the number of grid points in each space dimension in Figure 3(A). The convergence rates 2.01, 3.07 and 4.21 imply a gain of three orders in all the three cases, and are in agreement with the error estimate (51).

For the Neumann problem we again use erroneous boundary data, this time by perturbing the true data by adding $\nu=c_{p}h^{p-1},\ p=1,2,3$. Here the data is multiplied by $H^{-1}\sim\mathcal{O}(h^{-1})$ in the penalty term, hence the truncation error at those ten grid points is as in the Dirichlet case $\mathcal{O}(h^{p-2})$. As is mentioned towards the end of Section 3.1.3, the perturbation changes the structure of the right–hand side of the boundary system. The analysis for one dimensional problems presented in [30] leads to $w=1$ in (51) for all the three cases, and a corresponding gain of two orders for the two dimensional problem. The L2 errors from the numerical computations are plotted in Figure 3(B). The convergence rates are 1.00, 2.00 and 3.00. In all the three cases, the gain is two orders, and agrees with the accuracy analysis leading to the estimate (51). For neither the Dirichlet nor the Neuman problem do we observe any effect of the logarithmic term in (51), but this is not surprising since asymptotically such an effect would be difficult to detect.

5 Conclusion

In this paper, we extend the accuracy analysis of finite difference methods solving initial–boundary–value problems to two space dimensions. The two dimensional analysis is based on a diagonalization technique to decompose a two dimensional problem into one dimensional problems of the same type. We then continue the analysis by utilizing the results from the one dimensional analysis. We have chosen the second order wave equation as the model problem, but the technique presented in this paper can be used to analyze other equations.

The second contribution of this paper is the analysis of the effect of truncation errors localized at a few grid points in a corner of a two dimensional domain. This kind of truncation errors often occur in multi–block finite difference discretizations with non–conforming grid interfaces. The analysis is performed for a simplified but analogous problem, a single block with large truncation errors at a few grid points close to a corner. The numerical experiments for the simplified problem show that our accuracy analysis is sharp in the limit as the grid spacing approaches zero. Numerical experiments in a multi–block setting presented in [31] also agree well with this analysis.

In addition, we have presented a detailed framework of analyzing the convergence rate of one dimensional problems by the normal mode analysis, and have shown that the critical point is to derive sharp error estimates in the vicinity of $s=0$ in Laplace space. Singularities of the boundary system can also occur on the imaginary axis away from the origin, but they have no influence on the final convergence rate.

Acknowledgements

We would like to thank the help from Stefano Serra–Capizzano on Lemma 1, and Bengt Fornberg on Lemma 5, respectively. This work was performed when the second author was at the University of Bergen and was supported by VISTA (project 6357) in Norway. The first author and the third author are partially supported by the Swedish Research Council (project 106500511).

Appendix A Proof of Lemma 1

Proof.

Let $\tilde{Q}_{y}=PQ_{y}$, then the eigenvalue problem (4) can be written

[TABLE]

Since $P$ is symmetric positive definite, $P^{1/2}$ and $P^{-1/2}$ are also symmetric positive definite. We rewrite (56) as

[TABLE]

With the notation $\hat{Q}_{y}=P^{-1/2}\tilde{Q}_{y}P^{-1/2}$ and $\hat{\varphi}=P^{1/2}\varphi$, we obtain the following new eigenvalue problem

[TABLE]

Because $\tilde{Q}_{y}$ is symmetric negative semi–definite, according to Sylvester’s law of inertia, $\hat{Q}_{y}$ is also symmetric negative semi–definite. Therefore, the eigenvalues of $\hat{Q}_{y}$ are real and non–positive. The matrices $Q_{y}$ in (4) and $\hat{Q}_{y}$ in (57) are similar because $\hat{Q}_{y}=P^{1/2}Q_{y}P^{-1/2}$, so they have the same eigenvalues. This proves that $\lambda\geq 0$ in (4).

Let $\hat{\Phi}$ be the unitary operator $[\hat{\varphi}_{1},\hat{\varphi}_{2},\cdots,\hat{\varphi}_{N_{y}}]$, then its condition number in spectral norm is equal to one, i.e. $\chi(\hat{\Phi})=1$. Because of $\hat{\varphi}=P^{1/2}\varphi$, we have $\hat{\Phi}=P^{1/2}\Phi$. As a consequence, equation (4) can be written as

[TABLE]

where $\Lambda=\text{diag}(\lambda_{1},\lambda_{2},\cdots,\lambda_{N_{y}})$. Because both $P^{-1/2}$ and $\hat{\Phi}$ are invertible, $P^{-1/2}\hat{\Phi}$ is also invertible. Therefore, $Q_{y}$ is diagonalizable as

[TABLE]

Since the spectral norm is unitarily invariant and $\Phi=P^{-1/2}\hat{\Phi}$, we have

[TABLE]

By Assumption 1, $\|\Phi\|$ and $\|\Phi^{-1}\|$ are uniformly bounded. ∎

Appendix B Proof of Lemma 2

Proof.

For the special case $\delta=0$, the lemma is equivalent to Lemma 4 in [30], and its proof is found there. In the following, we only consider the case when $\delta>0$. Let $\tilde{s}=a+bi$ and $\tilde{s}_{+}=c+di$, where $a,b,c,d$ are real numbers and $a\geq\delta>0$. The relation $\tilde{s}_{+}^{2}=\tilde{s}^{2}+\gamma$ gives

[TABLE]

The real and imaginary part of the two sides of the above equation must be equal, which leads to

[TABLE]

Therefore,

[TABLE]

Assume Re$(\tilde{s}_{+})<\delta$, that is $0\leq c<a$. Then we have $c^{2}<a^{2}$, which means that the left hand side of (59) is negative. However, by $c^{2}<a^{2}$ the right hand side of (59)

[TABLE]

This is a contradiction. Therefore, we must have Re$(\tilde{s}_{+})\geq\delta$. ∎

Appendix C Proof of Lemma 3

Proof.

Let

[TABLE]

for some complex number $s$ and $\phi$ is a one dimensional grid function with $\|\phi\|_{1D,x}<\infty$. Substituting (60) to (19), with the notation $\tilde{s}=sh$ we obtain the eigenvalue problem

[TABLE]

We note that (61) is in exactly the same form as (20) with a zero truncation error. As a consequence, the boundary system corresponding to (61) is

[TABLE]

where the left–hand side of (62) is the same as the left–hand side of (22).

If $C(\tilde{s})$ is singular for some $\tilde{s}$ with Re$(\tilde{s})>0$, then (62) has a non–trivial solution $\Sigma\neq\mathbb{0}$. It then follows that (60) with Re$(s)>0$ is a solution of (19). However, this contradicts stability for the reason outlined in Lemma 12.1.1 in [11, pp. 378]. For completeness, we state it below.

We define a sequence of grids indexed by $n$

[TABLE]

The eigenvalue problem reads

[TABLE]

We define a sequence of grid functions $f^{(n)}=\phi^{(1)},\quad n=1,2,\cdots$, which satisfies the eigenvalue problem

[TABLE]

Therefore,

[TABLE]

are solutions of (60) and grow arbitrarily fast, i.e. a contradiction to stability. ∎

Appendix D Estimates of $|\Sigma_{i}|$

D.1 Proof of Lemma 4

Proof.

We decompose $\tilde{s}=i\tilde{\xi}+\eta h$ in the right–half plane into two parts. The first part is when $\tilde{s}$ is in the vicinity of the origin, i.e. $0<|\tilde{s}|\leq\delta$ where $\delta$ is a small constant independent of $h$. Here the range of $\tilde{s}$ includes a small, but $h$–independent interval of the imaginary axis. For example, we can in particular consider Re$(\tilde{s})=\eta h\leq\delta/\sqrt{2}$ and $|\text{Im}(\tilde{s})|\leq\delta/\sqrt{2}$. The second part is when $\tilde{s}$ is away from the origin, Re$(\tilde{s})=\eta h$ and $|\text{Im}(\tilde{s})|>\delta/\sqrt{2}\sim\mathcal{O}(1)$.

We start with the first case when $\tilde{s}$ is in the vicinity of the origin. The entries of $C(\tilde{s})$ are continuous functions in $\tilde{s}$. Every entry can be expanded in terms of $\kappa(\tilde{s})$ or $\tilde{s}$, and $\kappa(\tilde{s})$ can be further expanded to its Taylor series. $C(\tilde{s})$ can be written as

[TABLE]

with the notation $C^{\prime}(\tilde{s})=dC(\tilde{s})/d\tilde{s}$. An estimate of the type (29) can be obtained by using Lemma 3.4 in [24]. For completeness, we state this lemma below (notations of norms are changed to be consistent with the notations used in this paper) and explain how to use it thereafter.

Lemma 6**.**

(Lemma 3.4 in [24]) Consider the $n\times n$ linear system $(A+\delta E)x=b$ where $A$ is singular with rank $n-1$. Let $USV^{*}=A$ be the singular value decomposition of $A$. If $(U^{*}EV)_{nn}\neq 0$ then for all sufficiently small $|\delta|$ we get

[TABLE]

If, in addition, $b$ is in the column space of $A$, then for all sufficiently small $|\delta|$, we have $\|x\|_{max}\leq c\|b\|_{\max}$. Here $c$ is independent of $\delta$.

The first step of using Lemma 6 is to perform the singular value decomposition (SVD) of $C(0)=\mathcal{USV^{*}}$. The value of $w$ in (29) depends on 1) whether $(U^{*}C^{\prime}(0)V)_{nn}$ is equal to zero; 2) whether $\hat{T}_{C}$ is in the column space of $C(0)$.

• •

When $(U^{*}C^{\prime}(0)V)_{nn}\neq 0$: if $T_{C}$ is in the column space of $C(0)$ then $w=0$; otherwise $\|C^{-1}(\tilde{s})\|_{\max}\sim 1/(\eta h)$ and $w=1$.

• •

When $(U^{*}C^{\prime}(0)V)_{nn}=0$: we need to take into account the third term in the Taylor series of $C(\tilde{s})$, i.e. $C^{{}^{\prime\prime}}(0)$, and check whether $(U^{*}C^{{}^{\prime\prime}}(0)V)_{nn}$ is equal to zero. More generally, we obtain the estimate (29) with $m=1,2,\cdots,w-1$ if

[TABLE]

If $w$ in (29) is infinite, then $C(\tilde{s})$ is singular for some $\tilde{s}$ with a positive real part. This contradicts to the stability of the numerical scheme. Therefore, with a stable discretization $w$ is always finite.

Next, we consider the case when $\tilde{s}$ is away from the origin. Assume that $C(\tilde{s})$ is singular at $\tilde{s}=i\tilde{\xi}$ for some $\tilde{\xi}$ where $|\tilde{s}|>\delta\geq 0$ with $\delta$ independent of $h$. Now we may need to use a Puiseux series for $\kappa(\tilde{s})$. This leads to a similar expansion of $C(\tilde{s})$ as above, but with a non–integer leading order exponent. We round the exponent to its nearest integer, denoted by $\alpha$, towards positive infinity. By the same argument as the case $|\tilde{s}|\leq\delta$, at $\tilde{s}=i\tilde{\xi}+\eta h$ we have

[TABLE]

The non–zero components of $\hat{T}_{C}$ are in the form $\partial^{p+2}/\partial x^{p+2}\hat{U}(0,s)$ as given in (21). The key is to realize that

[TABLE]

where $\mathcal{L}$ is the Laplace transform operator in $t$. Here we have used the compatibility condition between the initial and boundary data. At $\tilde{s}=sh=(i\xi+\eta)h$ we have

[TABLE]

Since the singularity is away from the origin, we have

[TABLE]

This proves that no accuracy loss is caused by the singularity of $C(\tilde{s})$ away from the origin. Note that by using this lemma the term $\|\mathcal{L}(\partial^{\alpha}T_{C}/\partial t^{\alpha})\|_{\max}$ is added to the final estimate. ∎

Appendix E Estimates of $1/(1-|\kappa_{j}|^{2})$

E.1 Proof of Lemma 5

Proof.

The characteristic equation takes the general form

[TABLE]

where $a_{j}$ are the coefficients of the standard central finite difference stencils [7]. In the estimate, we only need to consider the $l$ admissible roots and derive estimates for $1/(1-|\kappa_{j}|^{2}),\ j=1,2,\cdots,l$. We start with the case when $\tilde{s}=i\tilde{\xi}+\eta h$ is in the vicinity of the origin, i.e. $|\tilde{s}|\leq\delta$ where $\delta$ and $\eta$ are small constants independent of $h$.

If a root $|\kappa_{j}(0)|<1$ then a perturbation analysis straightforwardly leads to $|\kappa_{j}(\tilde{s})|<1$ and consequently $1/(1-|\kappa_{j}(\tilde{s})|^{2})$ is bounded independently of $h$. However, for a root $|\kappa_{j}(0)|=1$ a careful derivation is needed to obtain the precise dependence of $1/(1-|\kappa_{j}(\tilde{s})|^{2})$ on $h$.

To proceed, we derive another form of the characteristic equation. Since admissible roots are only related to interior stencils, we consider a uniform grid in one space dimension

[TABLE]

with grid spacing $h$. We denote $D^{(2l)}$ the $2l^{th}$ order accurate central finite difference operator approximating $\partial^{2}/\partial x^{2}$, and apply it to the mode $e^{i\omega x_{j}}$, where $-\pi<\omega h\leq\pi$. Similarly to the derivation in [6, pp. 41], the operator $D^{(2l)}$ acting on the mode $e^{i\omega x_{j}}$ gives

[TABLE]

where

[TABLE]

The corresponding characteristic equation is

[TABLE]

A root $|\kappa|=1$ at $\tilde{s}=0$ corresponds to a root $f(l,\omega h)=0$ for $-\pi<\omega h\leq\pi$. It is obvious that $f(l,\omega h)=0$ has a double root $\omega=0$. In addition, $f(l,\omega h)$ is a monotone decreasing function in $l$ so the only roots of $f(l,\omega h)=0$ are the double root $\omega=0$, corresponding to a double root $\kappa=1$ of the characteristic equation (64) with $\tilde{s}=0$. We therefore factorize (64) to

[TABLE]

where $P(1)\neq 0$. When $\tilde{s}=i\tilde{\xi}+\eta h$, the double root is perturbed to two single roots, where precisely one is admissible. Substituting the admissible root $\kappa_{\nu}=1+\nu$ to (66), we obtain

[TABLE]

The value of $\kappa_{\nu}^{l}/P(\kappa_{\nu})$ to the leading order is real and $\mathcal{O}(1)$. Therefore, there exists constant $K_{1}$ and $K_{2}$ independent of $h$ such that

[TABLE]

Note that the admissibility condition leads to Re$(\nu)<0$. We now have the estimate

[TABLE]

for $\tilde{s}$ in a vicinity of the origin.

Next, we consider the case when $\tilde{s}=i\tilde{\xi}+\eta h$ is away from the origin. Assume an admissible root $|\kappa_{a}(i\tilde{\xi}+\eta h)|<1$ satisfies $|\kappa_{a}(i\tilde{\xi})|=1$. The expansion of $\kappa_{a}(i\tilde{\xi}+\eta h)$ around $\tilde{s}=i\tilde{\xi}$ leads to

[TABLE]

where $\beta$ is the leading order exponent in its Puiseux series rounded to the nearest integer towards positive infinity. We note that when estimating the error (25), $\frac{1}{1-|\kappa_{a}(\tilde{s})|^{2}}$ is multiplied by $|\sigma_{a}|^{2}$, which is computed by solving the boundary system (22) and is related to the spatial derivatives of the true solution. We can therefore eliminate the $h$–dependence in (67) in the same manner as in Lemma 4, which increases the order of temporal derivative of the true solution in the final estimate from $\alpha$ to $\alpha+\beta$. ∎