On Long Time Error Bounds for the Wave Equation on Second Order Form

Temporal error bounds for the wave equation expressed on second order form are investigated. We show that, with appropriate choices of boundary conditions, the time and space derivatives of the error are bounded even for long times. No long time bound on the error itself is obtained, although numerical experiments indicate that a bound exists.


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 → ∞. 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.

The Wave Equation on Second Order Form
The wave equation on second order form with non-reflecting boundary conditions is 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].

Remark 1
The assumption η(t) ≥ η 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 t 0 η(τ )dτ 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 d dt (||e|| 2 ) = 2 1 0 ee t dx ≤ 2||e||||e t ||. Using this fact, we get, Integrating (8) in time and using the initial conditions yields, For large t, the term t/η 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 Hence, a linear growth of |||ē||| and a quadratic growth of ||e|| will be expected.

The Semi-discrete Wave Equation on Second Order Form
augmented with the initial conditionsv(0) =f andv t (0) =h. In (11), D = P −1 Q is the SBP finite difference operator that approximates the spatial derivative andv x = Dv. 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),f ,ḡ 0 ,ḡ 1 ,h,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].

An Error Bound for the Semi-discrete Problem
Letū be the solution to (1) injected at the grid. Insertingū into (11) and subtracting (11) with the numerical solutionv results in the discrete error equation, augmented with homogeneous conditionsē(0) = 0 andē t (0) = 0. In (12),ē =ū −v, e x = Dē and T e is the truncation error. Multiplying (12) byē T t P from the left, adding the transpose of the outcome and using the SBP property of Q results in, where |||ẽ||| 2 P =ē T t Pē t + c 2ēT x Pē 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, defininḡ and assuming that there is anη 0 such thatη ≥η 0 > 0, we arrive at the estimate where ||T e || max is an upper estimate of ||T e || P . Using the fact that (||ē|| P ) t ≤ ||ē t || P ≤ |||ẽ||| P , (15) leads to, Integrating (16) in time and using the homogeneous initial conditions yields which is analogous to the continous estimate (9).

Remark 5
The assumptionη(t) ≥η 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 t 0η (τ )dτ is monotonically increasing. This situation is exactly the same as in the continous case, see Appendix B for details. Equation (15) states thatē t andē x remain bounded as t → ∞. As in the continous problem, no bound for the actual errorē 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ē t andē x are bounded for long times, and the errorē grows at most linearly in time.
Remark 6 Similar to the continuous case discussed above, ||ẽ|| is predicted to grow at a linear rate and ||ē|| 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 η(t) in (6) andη(t) in (14). All boundary conditions that produce similar dissipative terms will lead to error bounds.

The Fully Discrete Numerical Scheme
The numerical scheme used to approximate (1) is, 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 N x are all zero except entry (N , N ) which is equal to one. The vectorsḡ 0,1 ,h andf are grid functions where the continuous boundary and initial data are injected at appropriate positions. In (18), the symbol ⊗ denotes the Kronecker product, which is defined as 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].

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π(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 T Pe. The results are shown in Table 1, and the rate of convergence is as expected.
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 η = 0 in (5) and a linear growth of |||ē|||. 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. Table 1 The error for different mesh-sizes when using a second (SBP(2,1)), third (SBP(4,2)) and fourth (SBP (6,3)  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, η(t) becomes zero by inserting the Neumann boundary conditions, yielding a linear growth of |||ē|||. 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.
Finally, we use the non-reflecting dissipative type of boundary conditions that lead to a non-zero η 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.

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, 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π(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.
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).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

A The wave equation on first order form
For comparison, we also consider the wave equation on first order form. By introducing the variables v 1 = u t and v 2 = cu x , the problem (1) can be rewritten as, where ||δF|| max is an upper estimate of the disturbance inF and (cê 2 1 (0, t)+ê 2 1 (1, t))/2||ê|| 2 = η(t) ≥ η 0 > 0. The error in the components v 1 and v 2 areê 1 ,ê 2 respectively, and e = [ê 1 ,ê 2 ] T . The estimate (21) leads to the same conclusions as in the second order case.

C The wave equation with periodic boundary conditions
Consider the continous energy estimate (3), The periodic boundary conditions are implemented using a periodic finite difference operator, D p = −D T p for the spatial derivative. A stable semi-discrete scheme corresponding to (11) is then,v tt = c 2 D pvx +F.