1 Introduction

For a stable and consistent numerical scheme, the solution converges for a fixed time as the grid spacing h approaches zero. However, convergence does not necessarily mean that the error is bounded as the time \(t\rightarrow \infty \). Consequently, the classical definition of stability is not sufficient for an accurate solution after long time integration.

Long time error bounds have previously been studied for hyperbolic problems on first order form [1,2,3]. However, rewriting a problem on first order form makes the problem more computationally demanding [4]. In this note, we consider the long time behavior of the error for the wave equation on second order form.

We discretize using Summation-By-Parts (SBP) operators, and implement the boundary conditions using Simultaneous Approximation Terms (SAT’s) [5]. To minimize reflections, appropriate non-reflecting dissipative boundary conditions are used. We proceed to investigate the scheme and the corresponding error equation, and focus our investigation on the errors after long times.

The paper will proceed as follows. In Sects. 2 and 3, we derive error bounds for the continuous wave equation on second order form. The equation is discretized using the SBP-SAT technique, and it is shown that the continuous analysis carries over to the discrete problem. The theoretical findings are verified by numerical experiments in Sect. 4. Finally, we summarize and draw conclusions in Sect. 5.

2 The Wave Equation on Second Order Form

The wave equation on second order form with non-reflecting boundary conditions is

$$\begin{aligned} \begin{aligned} u_{tt}-c^2u_{xx}&=F(x,t), \quad x\in [0,1], t \ge 0\\ u_t(0,t)-cu_x(0,t)&=g_0(t)\\ u_t(1,t)+cu_x(1,t)&=g_1(t) \end{aligned} \end{aligned}$$
(1)

augmented with the initial conditions \(u(x,0)=f(x)\) and \(u_t(x,0)=h(x)\). In (1), \(c>0\) is a real constant and \(g_{0},g_1,f,h, F\) are given functions. The problem (1) is well-posed, see [6].

2.1 An Error Bound for the Continuous Problem

Let \({\hat{u}}\) be a solution to (1) with a perturbed forcing function \(F+\delta F\). Subtracting (1) from the problem for \({\hat{u}}\) leads to the error equation,

$$\begin{aligned} \begin{aligned} e_{tt}-c^2e_{xx}&=\delta F, \quad x\in [0,1], t \ge 0\\ e_t(0,t)-ce_x(0,t)&=0\\ e_t(1,t)+ce_x(1,t)&=0 \end{aligned} \end{aligned}$$
(2)

for \(e={\hat{u}}-u\) augmented with homogeneous initial conditions.

Multiplying (2) by \(e_t\) and integrating over the spatial domain yields,

$$\begin{aligned} \int _0^1e_te_{tt}dx=c^2\int _0^1e_te_{xx}dx+\int _0^1e_t\delta Fdx. \end{aligned}$$
(3)

Using the fact that \(2e_te_{tt}=\frac{\partial }{\partial t}(e^2_t)\), \(e_te_{xx}=(e_te_x)_x-(e_x^2)_t/2\) and imposing the boundary conditions, we obtain

$$\begin{aligned} \frac{d}{dt}(|||{\bar{e}}|||^2)=-2c(e^2_t(1,t)+e^2_t(0,t))+2\int _0^1e_t\delta F dx. \end{aligned}$$
(4)

In (4), we have defined \(|||{\bar{e}}|||^2=\int _0^1(e_t^2+c^2e_x^2)dx\). Noting that \(\frac{d}{d t}(|||{\bar{e}}|||^2)=2|||{\bar{e}}|||\cdot |||{\bar{e}}|||_t\) and dividing (4) by \(2|||{\bar{e}}|||\) results in

$$\begin{aligned} \frac{d}{d t}(|||{\bar{e}}|||) \le -c\frac{e^2_t(1,t)+e^2_t(0,t)}{|||{\bar{e}}|||^2}|||{\bar{e}}|||+||\delta F|| = -\eta (t)|||{\bar{e}}|||+||\delta F||. \end{aligned}$$
(5)

In (5), we have used that \(2\int _0^1e_t\delta F dx \le 2||e_t||||\delta F||\le 2|||{\bar{e}}|||\cdot ||\delta F||\) and introduced

$$\begin{aligned} \eta (t)=c(e^2_t(1,t)+e^2_t(0,t))/(|||{\bar{e}}|||^2). \end{aligned}$$
(6)

Solving (5) for \(|||{\bar{e}}|||\) and assuming that there are constants \(\eta _0\) and \(||\delta F||_{max}\) such that \(\eta (t) \ge \eta _0>0\) and \(||\delta F|| \le ||\delta F||_{max}\) results in

$$\begin{aligned} |||{\bar{e}}||| \le \frac{1- \exp {(-\eta _0t)}}{\eta _0}||\delta F||_{max}, \end{aligned}$$
(7)

after using the homogeneous initial conditions. With \(\eta _0>0\), \(|||{\bar{e}}|||\) is bounded, i.e. it approaches a finite value as \(t\rightarrow \infty \). Since \(|||{\bar{e}}|||^2=\int _0^1(e_t^2+c^2e_x^2)dx\), (7) implies that both \(e_t\) and \(e_x\) are bounded.

Remark 1

The assumption \(\eta (t)\ge \eta _0>0\) is crucial for the result in (7). This is in fact a simplification of a more general condition. It is shown in [2] that the precise requirement is that \(\int _0^t \eta (\tau )d\tau \) is monotonically increasing. The situation is the same here and we clarify this in Appendix B.

The previous analysis does not imply that the error itself is bounded, but rather that its growth is limited. To clarify, we follow [7], and observe that \(\frac{d}{dt}(||e||^2)=2\int _0^1ee_tdx \le 2||e|| ||e_t||\). Using this fact, we get,

$$\begin{aligned} ||e||_t \le ||e_t|| \le |||{\bar{e}}||| \le \frac{1- \exp {(-\eta _0t)}}{\eta _0}||\delta F||_{max}. \end{aligned}$$
(8)

Integrating (8) in time and using the initial conditions yields,

$$\begin{aligned} ||e|| \le \frac{\eta _0t-1+exp(-\eta _0 t)}{\eta _0^2}||\delta F||_{max} \end{aligned}$$
(9)

For large t, the term \(t/\eta _0\) in (9) will dominate, and the error grows linearly in time.

We summarize the results so far in the following proposition.

Proposition 1

In problem (2), both \(e_t\) and \(e_x\) are bounded for long times, and the error e grows at most linearly in time.

Remark 2

The estimate (9) is an upper estimate of ||e||. Consequently, (9) does not necessarily mean that the error is unbounded for long times.

Remark 3

Without the damping term from the boundary conditions, the estimate (5) becomes

$$\begin{aligned} \frac{d}{d t}(|||{\bar{e}}|||) \le ||\delta F||. \end{aligned}$$
(10)

Hence, a linear growth of \(|||{\bar{e}}|||\) and a quadratic growth of ||e|| will be expected.

3 The Semi-discrete Wave Equation on Second Order Form

Next, we discretize (1) using the SBP-SAT technique,

$$\begin{aligned} \begin{aligned} {\bar{v}}_{tt}&= c^2D{\bar{v}}_x-cP^{-1}E_{0}({\bar{v}}_t-c{\bar{v}}_x-{\bar{g}}_0)-cP^{-1}E_{N}({\bar{v}}_t+c{\bar{v}}_x-{\bar{g}}_1)+{\bar{F}} \end{aligned} \end{aligned}$$
(11)

augmented with the initial conditions \({\bar{v}}(0)={\bar{f}}\) and \({\bar{v}}_t(0)={\bar{h}}\). In (11), \(D=P^{-1}Q\) is the SBP finite difference operator that approximates the spatial derivative and \({\bar{v}}_x=D{\bar{v}}\). The matrix \(P=P^T>0\), Q satisfies the SBP property \(Q+Q^T=\)diag\((-1,0,...,0,1)\), \(E_{0}=\)diag(1, 0, ..., 0) and \(E_{N}=\)diag(0, ..., 0, 1). The data related to (11), \({\bar{f}},{\bar{g}}_{0},{\bar{g}}_1,{\bar{h}}, {\bar{F}}\) are grid functions of \(f,g_{0},g_1,h,F\), i.e. the function values are injected at the appropriate grid points. The second and third terms on the right-hand side of (11) are SAT’s that implements the boundary conditions. It can be shown that (11) is a stable scheme and that the solution converges for a fixed time [6].

3.1 An Error Bound for the Semi-discrete Problem

Let \({\bar{u}}\) be the solution to (1) injected at the grid. Inserting \({\bar{u}}\) into (11) and subtracting (11) with the numerical solution \({\bar{v}}\) results in the discrete error equation,

$$\begin{aligned} \begin{aligned} {\bar{e}}_{tt}&=c^2D{\bar{e}}_x-cP^{-1}E_{0}({\bar{e}}_t-c{\bar{e}}_x)-cP^{-1}E_{N}({\bar{e}}_t+c{\bar{e}}_x)+T_e \end{aligned} \end{aligned}$$
(12)

augmented with homogeneous conditions \({\bar{e}}(0)=0\) and \({\bar{e}}_t(0)=0\). In (12), \({\bar{e}}={\bar{u}}-{\bar{v}}\), \({\bar{e}}_x=D{\bar{e}}\) and \(T_e\) is the truncation error.

Multiplying (12) by \({\bar{e}}_t^TP\) from the left, adding the transpose of the outcome and using the SBP property of Q results in,

$$\begin{aligned} \frac{d}{d t}(|||{\tilde{e}}|||^2_P) =-2c(e_{Nt}^2+e_{0t}^2)+2{\bar{e}}^TPT_e, \end{aligned}$$
(13)

where \(|||{\tilde{e}}|||^2_P={\bar{e}}_t^TP{\bar{e}}_t+c^2{\bar{e}}_x^TP{\bar{e}}_x\) and \(e_{0,N}\) denote the errors at the first and last grid points, respectively.

Remark 4

When applying the energy method to (12) and using the first derivative twice, one arrives directly at (13), with a well-defined first derivative in space. This makes the estimate (13) strikingly similar to (4). With a compact operator in space (which is generally the preferred choice), this would not be the case [8].

Following the path of the continous analysis above, defining

$$\begin{aligned} {\bar{\eta }}(t)=c(e_{Nt}^2+e_{0t}^2)/(|||{\tilde{e}}|||^2), \end{aligned}$$
(14)

and assuming that there is an \({\bar{\eta }}_0\) such that \({\bar{\eta }}\ge {\bar{\eta }}_0>0\), we arrive at the estimate

$$\begin{aligned} |||{\tilde{e}}|||_P \le \frac{1-\exp {(-{\bar{\eta }}_0t)}}{{\bar{\eta }}_0}||T_e||_{max}, \end{aligned}$$
(15)

where \(||T_e||_{max}\) is an upper estimate of \(||T_e||_P\).

Using the fact that \((||{\bar{e}}||_P)_t\le ||{\bar{e}}_t||_P\le |||{\tilde{e}}|||_P\), (15) leads to,

$$\begin{aligned} (||{\bar{e}}||_P)_t \le \frac{1-\exp {(-{\bar{\eta }}_0t)}}{{\bar{\eta }}_0}||T_e||_{max}. \end{aligned}$$
(16)

Integrating (16) in time and using the homogeneous initial conditions yields

$$\begin{aligned} ||{\bar{e}}||_P \le \frac{{\bar{\eta }}_0t-1+\exp {(-{\bar{\eta }}_0t)}}{{\bar{\eta }}_0^2}||T_e||_{max}, \end{aligned}$$
(17)

which is analogous to the continous estimate (9).

Remark 5

The assumption \({\bar{\eta }}(t)\ge {\bar{\eta _0}}>0\) is crucial for obtaining an error bound. This is in fact a simplification of a more general condition. The precise requirement is that \(\int _0^t{\bar{\eta }}(\tau )d\tau \) is monotonically increasing. This situation is exactly the same as in the continous case, see Appendix B for details.

Equation (15) states that \({\bar{e}}_t\) and \({\bar{e}}_x\) remain bounded as \(t \rightarrow \infty \). As in the continous problem, no bound for the actual error \({\bar{e}}\) is obtained. On the contrary, (17) indicates that the error grows linearly in time. We summarize the results so far in the following proposition.

Proposition 2

In problem (12), both \({\bar{e}}_t\) and \({\bar{e}}_x\) are bounded for long times, and the error \({\bar{e}}\) grows at most linearly in time.

Remark 6

Similar to the continuous case discussed above, \(||{\tilde{e}}||\) is predicted to grow at a linear rate and \(||{\bar{e}}||\) at a quadratic rate without the damping from the boundary conditions.

Remark 7

The results obtained for the second order form hold also for the wave equation rewritten on first order form. See Appendix A for details.

Remark 8

The error bounds obtained in this paper all rely on the terms \(\eta (t)\) in (6) and \({\bar{\eta }}(t)\) in (14). All boundary conditions that produce similar dissipative terms will lead to error bounds.

3.2 The Fully Discrete Numerical Scheme

The numerical scheme used to approximate (1) is,

$$\begin{aligned} \begin{aligned}&D_t{\tilde{D}}_t{\bar{v}}=c^2D_x^2{\bar{v}}-(P_t^{-1}E_{0t}\otimes I_x){\tilde{D}}_t({\bar{v}}-{\bar{h}})\\&\quad -(I_t\otimes P_x^{-1}E_{0x} )({\tilde{D}}_t{\bar{v}}-D_x{\bar{v}}-{\bar{g}}_0) -(I_t\otimes P_x^{-1}E_{Nx} )({\tilde{D}}_t{\bar{v}}+D_x{\bar{v}}-{\bar{g}}_1), \end{aligned} \end{aligned}$$
(18)

where \({\tilde{D}}_t{\bar{v}}=D_t{\bar{v}}+(P_t^{-1}E_{0t}\otimes I_x)({\bar{v}}-{\bar{f}})\), \(D_t=P_t^{-1}Q_t\otimes I_x\) and \(D_x=I_t\otimes P_x^{-1}Q_x\). The subscripts on D, P and Q indicate which derivative that is approximated. The entries of \(E_{0x},E_{0t}\) are zero except entry (1, 1) which is equal to one, and the entries of \(E_{Nx}\) are all zero except entry (N, N) which is equal to one. The vectors \({\bar{g}}_{0,1}\), \({\bar{h}}\) and \({\bar{f}}\) are grid functions where the continuous boundary and initial data are injected at appropriate positions. In (18), the symbol \(\otimes \) denotes the Kronecker product, which is defined as

$$\begin{aligned} A\otimes B= \begin{bmatrix} a_{11}B&\hdots&\hdots&a_{1N}B\\ \vdots&\ddots&\hdots&\vdots \\ a_{N1}B&\hdots&\hdots&a_{NN}B, \end{bmatrix} \end{aligned}$$

for two arbitrary matrices A and B. For further details on the discretization of the wave equation on second order form, the reader is referred to [6].

4 Numerical Results

Consider problem (1) with \(c=1\). The problem is discretized in both time and space using the SBP-SAT technique [6], since the SBP-SAT technique can be applied directly to the second time derivative. In all the calculations below we use the manufactured solution \(u=\sin (2\pi (x-t))\) and choose the boundary and initial data accordingly. We also use \(N=20,40,80\) grid points in space and \(N_t=5000N\) grid points in time for all calculations.

As a quality control, we confirm that the solution converges with the correct rate during mesh refinement for a fixed time. The problem is integrated up to \(T=1\) and the norm of the error at the final time step is given by \(||e||^2_{P}=e^TPe\). The results are shown in Table 1, and the rate of convergence is as expected.

Table 1 The error for different mesh-sizes when using a second (SBP(2,1)), third (SBP(4,2)) and fourth (SBP(6,3)) order SBP scheme when solving (1)
Full size table

Next, the long time performance of the problem for different types of boundary conditions is studied. We integrate up to \(T=500\) and use an SBP scheme with third order overall accuracy (a fourth order accurate inner stencil and second order boundary closures) in both space and time.

First, we apply periodic boundary conditions (implemented using periodic finite difference operators), which lead to \(\eta =0\) in (5) and a linear growth of \(|||{\bar{e}}|||\). See Appendix C for more details. Hence, both \(||e_t||_P\) and \(||e_x||_P\) are expected to grow linearly in time, since the damping effect becomes zero, see Remark 3 and 6. The results are shown in Fig. 1 where \(||e_t||_P\) and \(||e_x||_P\) as well as the error \(||e||_P\) grows at an approximately linear rate.

Secondly, we apply Neumann boundary conditions \(u_x(0,t)=g_0(t)\) and \(u_x(1,t)=g_1(t)\), which also should result in a linear growth of \(||e_t||_P\) or \(||e_x||_P\) (or both). As in the periodic case, \(\eta (t)\) becomes zero by inserting the Neumann boundary conditions, yielding a linear growth of \(|||{\bar{e}}|||\). More details on the energy estimates when using Neumann boundary conditions can be found in [9]. Figure 2 shows that \(||e_x||_P\) grows at a linear rate, which is expected from (10). The actual error, \(||e||_P\), also grows at a linear rate which indicates that the estimates (7) and (17) may be too pessimistic.

Fig. 1
figure 1

The P-norm of e (upper figure), \(e_x\) (middle figure) and \(e_t\) (lower figure) for different mesh sizes. Blue line: \(N=20\), red line: \(N=40\) and green line: \(N=80\). The dashed line indicate a linear growth, and periodic boundary conditions are used (Color figure online)

Full size image
Fig. 2
figure 2

The P-norm of e (upper figure), \(e_x\) (middle figure) and \(e_t\) (lower figure) for different mesh sizes. Blue line: \(N=20\), red line: \(N=40\) and green line: \(N=80\). The dashed line indicate a linear growth, and the Neumann boundary conditions are used (Color figure online)

Full size image
Fig. 3
figure 3

The P-norm of e (upper figure), \(e_x\) (middle figure) and \(e_t\) (lower figure) for different mesh sizes. Blue line: \(N=20\), red line: \(N=40\) and green line: \(N=80\). Non-reflecting boundary conditions are used (Color figure online)

Full size image

Finally, we use the non-reflecting dissipative type of boundary conditions that lead to a non-zero \(\eta _0>0\). Figure 3 shows the new theoretical result that \(||e_t||_P\) and \(||e_x||_P\) do not grow but rather stay constant in time. Also \(||e||_P\) is bounded, which is a better result than predicted by the non-sharp estimates.

For brevity, we have restricted the numerical experiments to the cases discussed above. However, error boundedness for dissipative boundary conditions and linear growth otherwise is valid in general.

4.1 The Wave Equation in Two Space Dimensions

In this section, we consider the wave equation in two space dimensions with first order absorbing boundary conditions,

$$\begin{aligned} \begin{aligned} u_{tt}-c^2(u_{xx}+u_{yy})&=F(x,y,t)\\ u_t(0,y,t)-cu_x(0,y,t)&=g_{x0}(y,t)\\ u_t(1,y,t)-cu_x(1,y,t)&=g_{x1}(y,t)\\ u_t(x,0,t)-cu_y(x,0,t)&=g_{y0}(x,t)\\ u_t(x,1,t)-cu_y(x,1,t)&=g_{y1}(x,t). \end{aligned} \end{aligned}$$
(19)

Equations (19) is discretized in space and time using an SBP scheme of third order accuracy and with boundary and initial conditions enforced using SAT’s. We use \(N=20,40,80\) grid points in either space direction and \(N_t=500N\) grid points in time. Furthermore, we choose \(c=1\), let the exact solution be \(u=\sin (2\pi (x+y-2t))\) and choose the data accordingly. As shown in Figure 4, the results observed in the one-dimensional cases are also observed in the two-dimensional setting.

Fig. 4
figure 4

The P-norm of e (upper figure), \(e_x\) (middle figure) and \(e_t\) (lower figure) for the wave equation in two space dimensions and different mesh sizes. Blue line: \(N=20\), red line: \(N=40\) and green line: \(N=80\). First order non-reflecting boundary conditions are used (Color figure online)

Full size image

5 Summary and Conclusions

Long time error development for the wave equation on second order form has been investigated. We have theoretically shown that both the spatial and temporal derivatives of the errors are bounded if dissipative boundary conditions are given. It is also shown that the actual error grows at most linearly in time.

For non-dissipative boundary conditions, a linear growth for the spatial and temporal derivatives of the error is predicted and at most a quadratic one for the actual error.

Numerical experiments confirm that both the time and space derivatives of the error are bounded for dissipative boundary conditions. It is found that the actual error is also bounded. For non-dissipative boundary conditions, the time and space derivatives of the error grow linearly as well as the actual error.

The theoretical predictions of the growth for the time and space derivatives are validated by the numerical experiments, while the (non-sharp) theoretical estimate for the error seems a bit too pessimistic.

The results of this paper can be generalized to any scheme that leads to an energy estimate that mimics the continuous one in (7).