On compact 4th order finite-difference schemes for the wave equation

We consider compact finite-difference schemes of the 4th approximation order for an initial-boundary value problem (IBVP) for the $n$-dimensional non-homogeneous wave equation, $n\geq 1$. Their construction is accomplished by both the classical Numerov approach and alternative technique based on averaging of the equation, together with further necessary improvements of the arising scheme for $n\geq 2$. The alternative technique is applicable to other types of PDEs including parabolic and time-dependent Schr\"{o}dinger ones. The schemes are implicit and three-point in each spatial direction and time and include a scheme with a splitting operator for $n\geq 2$. For $n=1$ and the mesh on characteristics, the 4th order scheme becomes explicit and close to an exact four-point scheme. We present a conditional stability theorem covering the cases of stability in strong and weak energy norms with respect to both initial functions and free term in the equation. Its corollary ensures the 4th order error bound in the case of smooth solutions to the IBVP. The main schemes are generalized for non-uniform rectangular meshes. We also give results of numerical experiments showing the sensitive dependence of the error orders in three norms on the weak smoothness order of the initial functions and free term and essential advantages over the 2nd approximation order schemes in the non-smooth case as well.


Introduction
Compact higher-order finite-difference schemes for PDEs is a popular subject and a vast literature is devoted to them. The case of such type schemes for the wave equation have recently attracted a lot of interest, in particular, see [2,4,8,12], where much more related references can be found.
We consider compact finite-difference schemes of the 4th approximation order for an initialboundary value problem (IBVP) for the n-dimensional wave equation with constant coefficients, n 1. Their construction on uniform meshes is accomplished by both the classical Numerov approach and alternative technique based simply on averaging of the equation related to the polylinear finite element method (FEM), together with further necessary improvements of the arising scheme for n 2. This alternative technique is applicable to other types of PDEs including parabolic and time-dependent Schrödinger ones. The constructed schemes are implicit and three-point in each spatial direction and time. For n 2, there is a scheme with a splitting operator among them. Notice that we use implicit approximations for the second initial condition in the spirit of the approximations for the equation. Curiously, for n = 1 and the mesh on characteristics of the equation, the 4th order scheme becomes explicit and very close to an exact scheme on a four-point stencil.
We present a conditional stability theorem covering the cases of stability in strong (standard) and weak energy norms with respect to both initial functions and free term in the equation. Note that stability is unconditional for similar compact schemes on uniform meshes for other type PDEs, for example, see [3,11]. Our approach is applied in a unified manner for any n 1 (not separately for n = 1, 2 or 3 as in many papers), the uniform rectangular (not only square) mesh is taken, the stability results are of standard kind in the theory of finite-difference schemes and proved by the energy techniques (not only by getting bounds for harmonics of the numerical solution as in most papers). In particular, the last point allows us to prove rigorously the 4th order error estimate in the strong energy norm for smooth solutions.
Moreover, the main schemes are rather easily generalized for non-uniform rectangular meshes in space and time; we apply averaging technique to this end. Concerning compact schemes on non-uniform meshes for other equations, in particular, see [5,13,14,16].
In our 1D numerical experiments, we concentrate on demonstrating the sensitive dependence of the error orders in the mesh L 2 , uniform and strong energy norms on the weak smoothness order of the both initial functions and the weak dominating mixed smoothness order of the free term. The cases of the delta-shaped, discontinuious or with discontinuos derivatives data are covered. The higher-order practical error behavior is shown compared to standard 2nd approximation order schemes [15,17] thus confirming the essential advantages of 4th order schemes over them in the non-smooth case as well.
The paper is organized as follows. Auxiliary Section 2 contains results on stability of general symmetric three-level method with a weight for hyperbolic equations in the strong and weak energy norms that we need to apply. The main Section 3 is devoted to construction and analysis of the compact 4th order finite-difference schemes. In Section 4, the main compact schemes are generalized to the case of non-uniform rectangular meshes. The results of these sections are received by A. Zlotnik. Section 5 contains results of numerical experiments accomplished by O. Kireeva.

General symmetric three-level method for second order hyperbolic equations and its stability theorem
Let H h be a family of Euclidean spaces endowed with an inner product (·, ·) h and the corresponding norm · h , where h is the parameter (related to a spatial discretization). Let linear operators B h and A h act in H h and have the properties We assume that they are related by the following inequality For methods of numerical solving 2nd order elliptic equations, usually α h = c 0 /h min , where h min is a minimal size of the spatial discretization. We introduce the uniform mesh . We introduce the mesh averages and difference operators with y m = y(t m ),y m = y m−1 andŷ m = y m+1 , as well as the summation operator with the variable upper limit I m ht y = h t m l=1 y l for 1 m M and I 0 ht y = 0. We consider a general symmetric three-level in t method with a weight σ: where v: ω ht → H h is the sought function and the functions v 0 , u 1 ∈ H h and f : {t m } M −1 m=0 → H h are given; we omit their dependence on h for brevity. Note that the parameter σ can depend on h := (h, h t ). Recall that linear algebraic systems with one and the same operator B h + σh 2 t A h has to be solved at time levels t m to find the solution v m+1 , 0 m M − 1.
Let the following conditions related to σ hold: either σ 1 4 and ε 0 = 1, or Then one can introduce the following σ-and h t -dependent norm in H h and bound it from below: Obviously, for σ 1 4 , one also has w 0,h w B h , and then the norms · 0,h and · B h are equivalent uniformly in h.
We present the stability theorem for method (2.2)-(2.3) with respect to the initial data v 0 and u 1 and the free term f in the strong (standard) and weak energy mesh norms.
Define the norm y L 1 one can replace the f -term with 2I M −1 Proof. Similar bounds have recently been proved in [18] for the method of a more general form, with the parameter τ > 0 and an operator B 1h = B * 1h > 0 acting in H h . In these bounds, one can take τ = 1 and easily see from their proofs that the bounds mainly remain valid for B 1h = B * 1h 0, in particular, B 1h = 0 (the case considered here), up to the norm of f standing in (2.6) and the norm of g − s t g 0 mentioned in Item 2.
To verify the validity of the bounds precisely with the norms of f andg := g −s t g 0 indicated in this theorem, it suffices to modify bounds for the following summands with f in the strong energy equality in [18, Theorem 1] for 1 m M , and the relationsδ t = 1 2 (δ t +δ t ) and (2.5) have been applied. Clearly in fact the norm · 0,h stands on the left in (2.6) and on both sides in (2.7). Bounds of type (2.6) with a stronger norm of f can be found in [11].
Below we also refer to the following stability result.
whereas bound (2.7) remains valid (its proof does not change for ε 0 0). To be convinced of the latter bound, it is necessary to transform and bound differently the terms with v 0 and u 1 in the case f = 0 in the strong energy equality in [18]. Namely, using the formulas t v 1 = v 0 + 1 2 h tδt v 1 and equation (2. 3) with f 0 = 0, we can set C h := (B h + σh 2 t A h ) −1 and obtain A h under the assumptions either σ 1 4 and ε 0 = 1, or (2.4) with 0 ε 0 < 1 and, as a corollary, h . But, for ε 0 = 0, the quantity w 0,h could be (in general) only a semi-norm in H h , and its lower bound by w B h uniformly in h is not valid any more.
It is well-known that each of bounds (2.6)-(2.7) implies existence and uniqueness of the solution to method (2.2)-(2.3) for any given v 0 , u 1 ∈ H h and f : {t m } M −1 m=0 → H h . The same applies to finite-difference schemes below.
3 Construction and properties of compact finite-difference schemes of the 4th order of approximation We consider the following IBVP with the nonhomogeneous Dirichlet boundary condition for the slightly generalized wave equation Here a 1 > 0, . . . , a n > 0, x = (x 1 , . . . , x n ), Ω = (0, X 1 )×. . .×(0, X n ), n 1, ∂Ω is the boundary of Ω and Γ T = ∂Ω × (0, T ) is the lateral surface of Q T . Hereafter the summation from 1 to n over the repeated indices i, j (and only over them) is assumed. Below δ (ij) is the Kronecker symbol.
Let below H h be the space of functions defined onω h , equal 0 on ∂ω h and endowed with the inner product (v, w) h = h 1 . . . h n x k ∈ω h v k w k and the norm w h = (w, w) Lemma 3.1. For sufficiently smooth inQ T solution u to equation (3.1), the following formula holds where and I is the identity operator. Note that s Nĵ = I for n = 1.
Proof. We give two different proofs.
1. The first one follows to the classical Numerov approach. We take the simplest explicit three-level discretization of equation (3.1) having the form and, under the assumption of sufficient smoothness of u, select the leading term of its approximation error We express the derivatives ∂ 4 t u and ∂ 4 k u in terms of mixed derivatives by differentiating equation (3.1): Then formula (3.4) takes the form Here all the 2nd order derivatives can be replaced by the corresponding symmetric three-point difference discretizations preserving the order of the remainder: Recalling the definition of ψ e in (3.4), we can rewrite the last formula as (3.3).
2. The second proof is shorter and is based simply on averaging of equation (3.1) related to the polylinear finite elements. We define the well-known average in the variable x k related to the linear finite elements For a function w(x k ) smooth on [0, X k ], the following formulas hold and The first formula is checked by integrating by parts and other formulas hold owing to the Taylor formula at x kl with the residual in the integral form The respective formulas hold for the averaging operator q t in the variable t = x n+1 as well (since one can set X n+1 = T and h n+1 = h t ).
We apply the operatorqq t withq := q 1 . . . q n to equation (3.1) at the nodes of ω h and get The multiple application of the above formulas for the averages leads to and thus formula (3.3) is derived once again.

Formula (3.3) means that the discretization of equation (3.1) of the form
has the approximation error of the order O(|h| 4 ).
Notice that the coefficients of formulas respectively on ω ht and ω h are independent of h. For discretization (3.10), we consider the corresponding equation at t 0 = 0 3), and find out for which u 1N and f 0 N its approximation error also has the order O(|h| 4 ). Let 0 <h t T and h t h t .
, the approximation error of equation (3.11) satisfies the following formula Proof. Let 0 t h t . Once again we give two proofs. 1. Using Taylor's formula in t and grouping separately terms with the time derivatives of odd and even orders, we obtain In virtue of equation (3.1) we have Moreover, (∂ t u) 0 = u 1 , therefore we find Using (3.1) for t = 0 and the formula s N = s Nk + 1 12 h 2 k Λ k , we also have Therefore we have proved the formula This formula and (3.15) under choice (3.12)-(3.13) lead to formula (3.14). 2. Again the second proof is based on averaging of equation (3.1). We define the related one-sided average in t over (0, h t ) and apply ht Using Taylor's formula at t = 0 and calculating the arising integrals, we find Here we omit the integral representations for O(h 4 t )-terms for brevity. As in the proof of Lemma 3.1 and owing to the last expansion, we haveq(δ t u) 0 = s N (δ t u) 0 + O(|h| 4 ) and Also owing to Taylor's formula in t at t = 0 we can write down Thus similarly first to (3.18) and second to (3.19) we obtain Inserting all the derived formulas into (3.17), we again obtain the desired result. (3.13)) holds for the following three-and two-level approximations One can easily check this using the Taylor formula in t at t = 0.
, then clearly the same property holds for the one more three-level approximation Remark 3.3. Below we consider the case of non-smooth f . Namely the above second proofs of Lemmas 3.1-3.2 clarify that then f m N should be replaced withqq t f m , 0 m M −1, according to (3.9) and (3.17) and identically to the polylinear FEM with the weight [17], or with some its suitable approximation.
In the simplest case n = 1, equations (3.10)-(3.11) supplemented with the boundary condition take the following form where equations are valid respectively on ω h and ω h . Hereafter we assume that the function v 0 is given onω h and take the general nonhomogeneous Dirichlet boundary condition. This scheme can be interpreted as the particular case of scheme (2.2)-(2.3) with the operators B h = I and A h = −a 2 1 Λ 1 and the weight σ = σ(h) = 1 (a similar choice of σ was used in [11] in the 1D parabolic case) or the bilinear finite element method [17] (though the right-hand sides of the equations are not the same; but see also Remark 3.3).
But for n 2 the above constructed equations (3.10)-(3.11) are not of type (2.2)-(2.3). Therefore we replace them with the following one for a function u sufficiently smooth inQ T , and thus the approximation errors of the both equations of this scheme are also of the order O(|h| 4 ).
But the latter scheme fails for n 3 similarly to [3] in the case of the time-dependent Schrödinger equation. The point is that s N should approximate I adequately, but for the minimal and maximal eigenvalues of s N < I as the operator in H h we have Therefore λ min (s N ) > 1 − n 3 and λ min ( that is suitable for n = 1, 2, but s N becomes almost singular for n = 3 and even λ min (s N ) < 0 (i.e., s N is not positive definite any more) for n 4, for small |h|.
Thus for n = 3 it is of sense to replace the last scheme with the scheme Moreover, for any n 1 we can use the following scheme Herewith for the minimal and maximal eigenvalues ofs N < I as the operator in H h we have Moreover, the following relation betweens N and s N holds   For n 2, all the constructed schemes can be effectively implemented by means of solving the systems of linear algebraic equations with the same matrix arising at each time level using FFT with respect to sines in all (or n − 1) spatial directions (after excluding the given valueŝ v| ∂ω h =ĝ in the equations at the nodes closest to ∂ω h ). The matrices are non-singular (more exactly, symmetric and positive definite after the mentioned excluding) is definitely guaranteed under the hypotheses of Theorem 3.1 below.
For n 2, we also write down the schemē with the following splitting operator at the upper time level Splitting of such type is well-known and widely used, in particular, see [11,17], and the implementation of this scheme is most simple and comes down to sequential solving of systems with tridiagonal matrices in all n spatial directions which are definitely non-singular under the hypotheses of Theorem 3.1 below. The following relation betweenB N ands N holdsB N =s N + 1 12 h 2 tĀ N + R with the "residual" operator Clearly R as the operator in H h satisfies R = R * > 0. In particular, one has for a function u sufficiently smooth inQ T , scheme (3.29)-(3.30) has the approximation error O(|h| 4 ) as scheme (3.26)-(3.27). Note that some other known methods of splitting are able to deteriorate this order of approximation. Now we study the operator inequality in (2.1) for the above arisen operators.
where C 0 = 4 3 in the first case of (B h , A h ) or C 0 = 1 in other cases.
Proof. Let 1 k n and {λ be the collection of eigenvalues of the operator −Λ k in H h , with the maximal of them λ . The inequality −Λ k α 2 1h s kN in H h is equivalent to the following inequality between the eigenvalues of these operators λ Consequently the sharp constant is Herewith α 2 1h = 6 1 , thus the last bound is asymptotically sharp. Similarly for n = 2 the inequality A N α 2 h s N in H h holds with It is not difficult to check that the function under the max sign has the positive partial derivatives with respect to arguments λ This implies (3.33) in the first case. The last bound is asymptotically sharp too.
Next, in virtue of the inequalities s Nî <s Nî for n = 3 (see formula (3.28) for n = 2) and −Λ k < 3 2 λ (k) max s kN in H h , the following inequalities in H h hold: for n 2. Therefore inequality (3.33) has been proved in all the cases.
Now we state a result on conditional stability in two norms for the constructed schemes. Then the solutions to all the listed schemes satisfy the following bounds Proof. The theorem follows immediately from the above general stability Theorem 2.1 applying assumption (2.4) for σ = 1/12, in virtue of inequality (2.5) and Lemma 3.3.
The proof is standard (for example, see [11]) and follows from the stability bound (3.35) applied to the error r := u − v (herewith r| ∂ω h = 0, r 0 = 0). The approximation errors play the role of f m N , 1 m M − 1, and u 1N in the equations of the schemes, and the above checked conclusion that they have the order O(|h| 4 ) for all the listed schemes, as well as h t = O(|h|) in Theorem 3.1.
Notice that, in the very particular case h 1 a 1 = . . . hn an = h t , schemes (3.20)-(3.21) and (3.29)-(3.30) become explicit (since thenB N = I, see (3.31)) and, moreover, the latter scheme differs from the simplest explicit scheme only by the above derived approximations of the free terms in its equations. Herewith, for scheme (3.20)-(3.21), condition (3.34) is valid with C 0 = 1 and only ε 0 = 0 (actually, one can check that with some 0 < ε 0 = ε 0 (h) < 1). But, for scheme (3.29)-(3.30) and n 2, the condition even with ε 0 = 0 fails; more careful analysis of inequality (3.33) for this scheme still allows to improve the bound for α 2 h but not the drawn conclusion itself. According to Remark 2.1 even in this particular case some stability bounds still hold for scheme (3.20)-(3.21). The bounds contain terms of the following type Thus w 0,h remains a norm in H h but clearly is no longer bounded from below by w h uniformly in h (since the constant in the last inequality is sharp and has the order O 1 ). The explicit scheme for n = 1 is very specific. Its equations are rewritten using a 4-point stencil For clarity, let us pass to the related Cauchy problem with any k ∈ Z, x k = kh, h = h 1 , a = a 1 and the omitted boundary condition. Then the following explicit formula holds It can be verified most simply by induction with respect to m. Notice that all the mesh nodes lie on the characteristics x − x k = ±a 1 t of the equation. Of course, the stability of the scheme can be directly proved applying this formula.
Let us take v 0 = u 0 and reset u 1N k = 1 2h x k+1 x k−1 u 1 (x) dx and hereafter T m k is the triangle with the vertices {(x k±m , 0), (x k , t m )} and R m k is the rhomb with the vertices {(x k±1 , t m ), (x k , t m±1 )}. Then the above formula for v m k takes the form v m k = 1 2 (u 0 (x k−m ) + u 0 (x k+m )) + 1 2a 1 x k+m thus at the mesh nodes it reproduces the classical d'Alembert formula for the solution u to the Cauchy problem for the 1D wave equation, where the approximate and exact solution coincide: v m k ≡ u(x k , t m ) for any k ∈ Z and 0 m M . In this respect, see also [7].

The case of non-uniform rectangular meshes
This section is devoted to a generalization to the case of non-uniform rectangular meshes. Let 1 k n. Define the general non-uniform meshes 0 = x k0 < x k1 < . . . < x kN k = X k in x k with the steps h kl = x kl − x k(l−1) and ω ht with the nodes 0 = t 0 < t 1 < . . . < t M = T and steps as well as h k max = max 1 l N k h kl and h t max = max 1 m M h tm . Define the difference operators where w l = w(x kl ) and y m = y(t m ). The last four operators are the generalizations of those defined above so that their notation is the same. We extend the above technique based on averaging equation (3.1) and generalize the above average in x k : For a function w(x k ) smooth on [0, X k ], formula (3.6) remains valid and on ω hk , and the first bound (3.7) remains valid for s = 1, 3 with h k replaced with h * k , see also (3.8), that follows from Taylor's formula after calculating the arising integrals over [x k(l−1) , x k(l+1) ]. Due to Taylor's formula we also have thus the second expansion for q k w implies that i.e., s kN = I + 1 12 (h k+ β k δ k − h k α kδk ) or, in the averaging form, all the presented formulas are valid on ω hk . The operator s kN generalizes one defined above. Its another derivation was originally given in [5], see also [10,13]. Recall that the natural (though not necessary) property α kl 0 and β kl 0 is equivalent to the rather restrictive condition on the ratio of the adjacent mesh steps 2 On ω ht , the average q t w = q n+1 w is defined similarly, and thus with s tN = I + 1 12 (h t+ β t δ t − h t α tδt ) or, in the averaging form, Let ω h = ω h1 × . . . × ω hn . Formula (3.9) for u remains valid and implies where h max = max{h 1 max , . . . , h n max , h t max }. Formula (3.17) for u remains valid as well. It involves only two first time levels thus easily covers the case of the non-uniform mesh in t and implies now on ω h , where q t y 0 is given by formula (3.16) with h t1 in the role of h t . Owing to the above formulas, see also Remark 3.4, the last two formulas with u lead us to the generalized scheme (3.26)-(3.27) where equations are valid respectively on ω h and ω h and have the approximation errors of the order O(h 3 max ). For the uniform mesh in t, the left-hand side of (4.1) takes the previous form whereas the terms N s tN f can be simplified keeping the same order of the approximation error: The splitting version of equation (4.1) can be got by replacing the operators in front of δ t v andδ t v by the operators of the form B N = (s 1 + 1 12 h * tht σ t a 2 1 Λ 1 ) . . . (s n + 1 12 h * tht σ t a 2 n Λ n ), where respectivelyh t = h t+ and σ t = β t , orh t = h t and σ t = α t . Sincē where the operator R satisfies formula (3.32) with h 2 t replaced with h * tht σ t , this replacement conserves the approximation error of the order O(h 3 max ). The splitting version of equation (4.2) is got simply by replacings N + h 2 t1 12Ā N with the above operator (3.31) with h t1 in the role of h t . One can check also that the approximation errors still has the 4th order O(h 4 max ) for nonuniform meshes close in a sense to the uniform one, cp. [14], provided that f 0 N is taken as in (3.13) with h t1 in the role of h t .
Here we do not touch the stability study in the case of the non-uniform mesh (even only in space) but this is noticeably more cumbersome like in [14] (since the operator s kN is not self-adjoint any more) and, moreover, imposes stronger conditions on h t , see also [16].

Numerical experiments for non-smooth data
In the IBVP (3.1)-(3.2) in the 1D case, we now take Ω := (−X/2, X/2) and rewrite the boundary condition as u| x=−X/2 = g 0 (t) and u| x=X/2 = g 1 (t), t ∈ (0, T ). We intend to analyze the practical error orders γ pr of r = u − v in three uniform in time mesh norms which below are denoted respectively as L 2 h , C h and E h (the 2nd and 3rd norms are the uniform and strong energy-type ones). Here w h = h N k=1 w 2 k 1/2 and N = N 1 . The respective expected theoretical error orders γ th are min 4 5 α, 4 , α 0; 4 5 (α − 1 2 ), 1 2 < α 11 2 ; 4 5 (α − 1), 1 α 6 (5.2) (in the spirit of [1]), where α is the parameter defining the weak smoothness of the data, see details below (concerning the first order, for α 1, it should refer to the continuous L 2 norm rather than the mesh one but that we will ignore). The proof of the first order in the case u 1 = f = 0 see in [6]. For comparison, recall that for the 2nd approximation order methods the corresponding theoretical error orders γ (2) th are according to [17]; recall that the middle error order is derived from two other ones. These orders also have recently been confirmed practically in [15]. Let P 0 (x) = (sgn x + 1)/2 be the Heaviside-type function, P 1 (x) = 1 − 2|x|, P k (x) = (sgn x)(2x) k (k 2) and Q k (t) = P 0 (t − t * )(t − t * ) k (k 0 and 0 < t * < T ) be piecewisepolynomial functions. For uniformity, we also set P −1 (x) = δ(x) and Q −1 (t) = δ(t − t * ) as the Dirac delta-functions concentrated at x = 0 and t = t * . We put X = T = 1.
To identify error orders more reliably, we compute the errors for respectively N = 200, 400, . . ., N max , where N max = 3200, 2000, 800 respectively for 1 2 α 5 2 , α = 7 2 , 9 2 ; also N = 200, 300, . . ., N max with N max = 600 for α = 11 2 (N max is lesser for α 7 2 to avoid an impact of the round-off errors on γ pr ). We plot graphs of log 10 r versus log 10 N , where r is each of the three norms (5.1), and seek the almost linear dependence between them by the least square method. Thus we calculate the dependence r ≈ c 0 h γpr = c 0 ( X N ) γpr . For α = 3 2 , 7 2 and 11 2 and the extended set N = 200, 400, . . . , 3200, we present E, C h , L 2 h -norms of the error denoted respectively by , , ♦ on Figs. 1-2. Notice the abrupt decrease of the error range as α grows. We also observe the slight oscillation of the data for α = 3 2 that is an exception (they also present for α = 1 2 ); instead, the linear behavior is typical for other α and the values of N mentioned at the beginning of the paragraph. The growth of the data for α = 7 2 and much more significant for α = 11 2 (as N increases) reflects the impact of the round-off errors.
a) α = 3/2 b) α = 7/2 The computed c 0 and γ pr together with the respective theoretical orders γ th and γ (2) th , see (5.2)-(5.3), and the error norms r (N ) for N = 200, N max are collected in Table 1. The main observation is the nice agreement between γ pr and γ th for all three norms in all Examples E α , thus the sensitive dependence of γ pr on the data smoothness order α becomes quite clear. This agreement is mainly better for the first and second norms (5.1) (similarly to [15]). Notice that γ (2) th /γ pr and the error r (200) in each norm decrease rapidly as α grows; the former demonstrates the essential advantages of the 4th approximation order scheme over 2nd order ones in the important case of non-smooth data as well.
Finally we remind the explicit scheme (3.36)-(3.37). For the same X and a but h t = h/a and T = M h t > 1, for example, the C h -norm of the error equals 0.311E−14 even for N = 20 and M = 10 already in Example E 3/2 ; thus clearly it is caused purely by the round-off errors.