A second order accuracy in time, Fourier pseudo-spectral numerical scheme for "Good" Boussinesq equation

The nonlinear stability and convergence of a numerical scheme for the "Good" Boussinesq equation is provided in this article, with second order temporal accuracy and Fourier pseudo-spectral approximation in space. Instead of introducing an intermediate variable \begin{document}$ \psi $\end{document} to approximate the first order temporal derivative, we apply a direct approximation to the second order temporal derivative, which in turn leads to a reduction of the intermediate numerical variable and improvement in computational efficiency. A careful analysis reveals an unconditional stability and convergence for such a temporal discretization. In addition, by making use of the techniques of aliasing error control, we obtain an \begin{document}$ \ell^\infty (0,T^*; H^2) $\end{document} convergence for \begin{document}$ u $\end{document} and \begin{document}$ \ell^\infty (0,T^*; \ell^2) $\end{document} convergence for the discrete time-derivative of the solution in this paper, in comparison with the \begin{document}$ \ell^\infty (0,T^*; \ell^2) $\end{document} convergence for \begin{document}$ u $\end{document} and the \begin{document}$ \ell^\infty (0,T^*; H^{-2}) $\end{document} convergence for the time-derivative, given in [ 19 ].

1. Introduction. The "Good" Boussinesq (GB) equation, similar to Kortewegde Vries (KdV) equation and cubic Schrödinger (CS) equation, provides a balance between dispersion and nonlinearity that leads to the existence of solitons. Many authors have investigated the GB equation and its various extensions and found the unique peculiarities: solitary waves only exist in a finite range of velocities, and they can interact to cause anti-solitons [30], etc.
In detail, the GB equation is formulated as There are many papers concerning the analysis of GB equation. For example, an exact formula which suits the interaction of the solitons is obtained in [29], and the analysis of the solitons interaction mechanism is provided in [30], in which the existence and regularity of solutions of the initial value problem is also proved.
In [32], the authors use the finite difference method to solve this equation and analyze the nonlinear stability and convergence of the method. Other related works could also be found in [1,4,5,15,16,18,19,22,29,30,24], etc.
For simplicity, a periodic boundary condition is assumed for the GB equation (1.1) over an one-dimensional domain Ω = (0, L). With appropriate initial data, u(x, 0) = u 0 (x) and u t (x, 0) = v 0 (x), which are also L-periodic, it is reasonable to assume that the "GB" equation has a unique, periodic and smooth enough solution in the time interval (0, T ). In turn, the Fourier collocation (pseudospectral) differentiation becomes a natural choice to obtain the optimal spatial accuracy. There have been extensive studies about this numerical method; see the fundamental analysis [3,6,20,23]. An application of the pseudospectral methods to the nonlinear equations was originated from the earlier works [27,26,28], in which the steady state Burgers' and Naiver-Stokes equations were analyzed. Moreover, the convergence analysis for the semi-discrete Fourier pseudo-spectral schemes to the time-dependent equations has been provided in [34,33].
For the GB equation (1.1), the pioneering numerical analysis of Fourier pseudospectral method could be found in [19], in which a second order (in time) pseudospectral scheme for (1.1) was analyzed. On the other hand, a severe time step restriction: ∆t ≤ Ch 2 (with C a fixed constant), has to be imposed in the nonlinear stability analysis. As explained in Remark 4.4 of [9], such a stability condition comes more from the technical difficulty in the theoretical analysis than an essential constraint in practical computations. In this article, we propose an alternate numerical scheme, which takes a similar form as the one studied in [19], with slight modifications on the temporal approximation to the linear terms one the right hand side. In more details, the second order approximation 1 4 (u n+1 + 2u n + u n−1 ), as reported in [19], is replaced by an alternate one, 1 2 (u n+1 + u n−1 ). Such a modification leads to an unconditional stability estimate for the linear part, and the nonlinear error terms could be analyzed via a linearized stability estimate. In addition, the theoretical results in [19] were not optimal: only the 2 error estimate for the variable u, combined with the discrete H −2 error estimate for u t , have been reported. A more careful analysis reveals a key difficulty to obtain a higher order error estimate in the energy norm: lack of aliasing error control tool to estimate the numerical error term associated with D 2 N (u 2 ). In this paper, we make use of an aliasing error control technique for the Fourier pseudo-spectral method, developed in [22] and extensively applied in other related works [8,9,14,21,36], so that a convergence estimate in a higher order energy norm in derived: the H 2 error estimate for the variable u, combined with the discrete 2 error estimate for u t , in the full O(∆t 2 + h m ) accuracy order. This is a great improvement over the result reported in [19], in addition to the unconditional convergence estimate.
In addition, the proposed numerical scheme deals with a direct discretization for the solution variable u; no intermediate variable is needed in the numerical implementation. In comparison with a more recent work [9], in which an intermediate variable ψ (to approximate u t ) has to be introduced, this approach simplifies the computational effort and improves the numerical efficiency. In another recent work [36], the second order operator splitting method was analyzed, which in turn leads to a three-stage numerical scheme at each time step. By contrast, our proposed numerical scheme is only involved with a single stage at each time step, and an FFT-based Poisson-like solver is sufficient for the implementation.
The paper is organized as follows. In Section 2, we review Fourier spectral and pseudospectral method, formulate the proposed numerical scheme and state the main theoretical result. The consistency analysis is given by Section 3, and the optimal rate convergence analysis is provided in Section 4. The numerical results are presented in Section 5. Finally, some concluding remarks are made in Section 6.
2. The numerical scheme and the main result.
2.1. Review of Fourier spectral and pseudospectral approximations. For any f (x) ∈ L 2 (Ω), Ω = (0, L), its Fourier series are given by In the Fourier pseudospectral approximation, B N is denoted as the space of the trigonometric polynomials in x of degree is up to N . In turn, we introduce two symbols for the convenience of presentation: P N and I N , the projection and interpolation operators onto B N , respectively. In more details, for any continuous function f with the Fourier expansion (2.1), P N f (x) is defined as Meanwhile, I N is introduced due to the numerical computation at a given set of points. For a given discrete vector function f over a uniform numerical grid with 2N + 1 points, we assume its discrete Fourier series are given by its Fourier interpolation is defined as In other words, the only differences between P N f (x) and I N f (x) are associated with the waves with frequency higher than N . The convergence of the projection and interpolation is given by where · is the standard L 2 norm; see more detailed analysis in [7,23]. For any collocation approximation to the function f (x) at the points x j the discrete differentiation operator D N could be appropriately defined, operating on the vector of grid values f = f (x j ). In more details, one may compute the collocation coefficients (f N c ) via discrete FFT, and then multiply them by the correct values (given by 2 πi) and perform the inverse FFT. At a mathematical level, the differentiation operator D N could be viewed as a matrix, and the above process could be treated as a matrix-vector multiplication. The same process is performed for the second and fourth derivatives ∂ 2 x , ∂ 4 x , where this time the collocation coefficients are multiplied by (−4 2 π 2 /L 2 ) and (16 4 π 4 /L 4 ), respectively. Of course, the differentiation matrix can be applied for multiple times, i.e. the vector f is multiplied by D 2 N and D 4 N , respectively. To facilitate the stability and convergence analysis, we introduce the following discrete L 2 inner product and L 2 norm: for any grid functions f and g. In addition, the summation by parts formulas will play an important role in the later analysis: (2.8) In the nonlinear analysis, some product terms involved with wave frequency greater than N have to be appropriately estimated. In the pseudo-spectral space, its interpolation has to be considered, due to the appearance of aliasing error. The following lemma enables one to obtain an H m bound of the interpolation of the nonlinear term.
2.2. The proposed numerical scheme and the convergence result. The following second-order (in time) scheme is proposed for the "GB" equation (1.1): (2.10) At each time step, only a constant-coefficient linear system needs to be solved: which could be very efficiently implemented via an FFT-based solver, since all the eigenvalues one the left hand side are positive. Before the convergence statement of the numerical scheme, its continuous extension in space has to be introduced, defined by u k ∆t,h = u k N , in which u k N ∈ B N is the continuous version of the discrete grid functions u k , with the interpolation formula given by (2.4). The main convergence result is stated below.
Denote u ∆t,h as the continuous function as the extension of the fully discrete numerical solution given by 2.10. When ∆t, h approaches to 0, the following convergence result is valid: Obviously, the time step ∆t and the space grid size h are bounded by given constants which only depends on the exact solution. The final time T and the exact solution controls the constant C in (2.12). The proof consists of the consistency analysis and convergence analysis, which is shown in later section.
Remark 2.3. The proposed numerical scheme (2.10) takes a very similar form as the one in the earlier work [19]: The primary difference is the replacement of the second order approximation 1 4 (u n+1 + 2u n + u n−1 ) by an alternate stencil, 1 2 (u n+1 + u n−1 ), on the right hand side of the proposed numerical.
However, there are a few essential improvements in terms of the stability property of the proposed numerical scheme over the one reported in [19]. A careful analysis reveals that the numerical stability for (2.13) could only be theoretically justified under a severe time step constraint ∆t ≤ Ch 2 , although an intuitive and linearized stability analysis, as well as the numerical results, indicate that a standard CF L condition ∆t = O(h) is sufficient. In contrast, the special stencil structure of the proposed scheme (2.10) results in an unconditional stability and convergence for a fixed final time, as will be presented in later analysis.
Remark 2.4. There have been other recent works of second order (in time) pseudospectral numerical schemes for the "GB" equation (1.1). For example, an intermediate variable ψ (to approximate u t ) is introduced in [9], which in turn results in more storage memory at each time step. Meanwhile, a second order operator splitting method is proposed and analyzed in [36], which consists of three computational stages at each time step. In contrast, our proposed numerical scheme (2.10) requires the least amount of memory storage and computational time among the existing second order scheme for the "GB" equation, while an unconditional convergence analysis is available.
3.1. Truncation error analysis in space. Given the uniform mesh grid (x i ), 0 ≤ i ≤ 2N , over the domain Ω = (0, L), we set the exact solution as u e , and denote U N as its projection into B N : (3.1) The standard projection estimates indicate that 3) In turn, the following inequalities could be derived: In addition, by making use of estimate (3.16) in [9] (with p = 2), we obtain the following truncation error estimate for U N ; the details are skipped for the sake of brevity.
3.2. Truncation error analysis in time. For simplicity of presentation, we assume T = K∆t with an integer K. The following preliminary estimate is excerpted from an earlier work [2], which will be useful in later consistency analysis.
where C only depends on T .
The following proposition gives the desired consistency result. To simplify the presentation below, for the constructed solution U N , we define its vector grid function U n = IU N as its interpolation: Set U N as the approximation solution constructed by (3.1) and let U as its discrete interpolation. Then we have

10)
in which M only depends on the regularity of the exact solution u e .
Proof. To facilitate the consistency estimates, the following quantities are introduced: (3.11) Note that the quantities on the left side are defined on the numerical grid (in space) point-wise, while the ones on the right hand side are continuous functions. We begin with the second order time derivative terms, F 1 and F 1e . With an application of Lemma 3.1, we get (3.12) In turn, an application of discrete summation in Ω leads to due to the fact that U N ∈ B N , and (3.4) was used in the second step. For the terms F 2 and F 2e , we start from a comparison between 1 2 (U n+1 N + U n−1 N ) and U N (· , t n ): (3.14) Meanwhile, an application of Prop. 3.1 gives at each fixed grid point. As a result, we get The terms F 3 and F 3e can be analyzed in the same way. We have For the nonlinear terms F 4 and F 4e , we begin with the following estimate  Therefore, the local truncation error estimate for τ 1 is obtained by combining (3.13), (3.16), (3.17) and (3.19). Obviously, constant M only dependent on the exact solution u e . This completes the consistency analysis. 4. Stability and convergence. Now we begin the analysis of the proposed numerical scheme (2.10). The point-wise numerical error grid function is given bỹ For convenience, we denote u n N ∈ B N and e n N ∈ B N , as the continuous version of numerical solution u n and e n , respectively.
The following preliminary estimate will be used in the later analysis. It is excepted as Lemma 4.1 in [9], and the details are skipped for the sake of brevity. (4.7) The bound for the truncation error term could be obtained by Cauchy inequality: (4.8) The rest work is focused on the nonlinear inner product, and we begin with an application of Cauchy inequality: We notice a W 2,∞ bound for the constructed approximate solution for any n ≥ 0, which comes from the regularity of the constructed solution. In addition, the following a-priori assumption is made for the numerical solution. An a-priori H 2 assumption up to time step t n .
We assume a-priori that the numerical error function (for u) has an H 2 bound at time step t n : e n N H 2 ≤ 1, with e n N = I N e n , (4.11) so that the H 2 and W 1,∞ bound for the numerical solution (up to t n ) is available u n N H 2 = U n N − e n N H 2 ≤ U n N H 2 + e n N H 2 ≤ C * + 1 :=C 0 , u n N W 1,∞ ≤ C|u n N H 2 ≤ CC 0 :=C 1 , (1-D Sobolev embedding),(4.12) As a result of the a-priori assumption, the nonlinear error term could be bounded as follows: 13) with C 1 := C(C * +C 0 ), in which the aliasing error control estimate (2.9) (in Lemma 2.1) is applied in the second step, due to the fact that (U n N + u n N )e n N ∈ B 2N , and the preliminary estimate (4.1) is recalled in the last step. Going back to (4.9), we get (by setting C 2 = CC 1 ) (4.15) In turn, the following energy norm for the numerical error function is defined: Then we arrive at E n+1 − E n ≤ C 2 ∆tE n + (C 2 + 1)∆t(E n+1 + E n ) + 1 2 ∆t τ n Subsequently, an application of discrete Gronwall inequality implies that . It is clear that (4.20) so that we obtain the desired convergence result Recovery of the H 2 a-priori bound (4.11).
With the help of the ∞ (0, T ; H 2 ) error estimate (4.21) for the variable u, we see that the a-priori H 2 bound (4.11) is also valid for the numerical error function e N at time step t n+1 , provided that This completes the convergence analysis. The proof of Theorem 2.2 is finished.

Numerical accuracy check.
In this subsection we perform a numerical accuracy check for the proposed numerical scheme (2.10). Similar to [9,19,35,36], the exact solitary wave solution of the "GB" equation is given by in which 0 < P ≤ 1. In addition, the amplitude A, the wave speed c 0 and the real parameter P are connected by the following identities: Since the exact profile (5.1) decays exponentially as |x| → ∞, it is natural to apply Fourier pseudospectral approximation on an interval (−L, L), with L large enough. In this numerical experiment, we set the computational domain as Ω = (−40, 40), and a moderate amplitude A = 0.5 is chosen in the test.
To investigate the accuracy in space, the time step size is fixed as ∆t = 10 −4 so that the temporal numerical error is negligible. The numerical scheme is implemented with grid sizes N = 32 to N = 128 in increments of 8, and the computations are carried out up to time T = 4. The discrete ∞ and 2 numerical errors for D 2 N u are displayed in Fig. 1. The spatial spectral accuracy is apparently observed for both the 2 and ∞ norms. Because of the fixed time step ∆t = 10 −4 , a saturation of spectral accuracy appears with an increasing N .  An apparent spatial spectral accuracy is observed for both norms.
To explore the temporal accuracy, we fix the spatial resolution as N = 512 so that the numerical error is dominated by the temporal ones. We compute solutions with a sequence of time step sizes, ∆t = T N K , with N K = 100 to N K = 1000 in increments of 100, and T = 4. Fig. 2 shows the discrete 2 and ∞ norms of the errors between the numerical and exact solutions. A clear second order accuracy is observed for both norms.

5.2.
Long time simulation results. In this subsection we present the longer time simulation results for the "Good" Boussinesq equation (1.1), using the proposed numerical scheme (2.10). The initial data is taken the same as the exact solution (5.1) at t = 0: The physical domain is set as Ω = (−40, 40), and we take the temporal step size as ∆t = 10 −3 . The time snapshots of the evolution computed by the proposed 6. Concluding remarks. In this paper, we propose a Fourier pseudo-spectral method for the "Good" Boussinesq equation, with a second-order temporal accuracy, and provide a stability and convergence analysis. In comparison with a few recent related works, we apply a direct approximation to the second order temporal derivative, instead of introducing an intermediate variable ψ to approximate the first order temporal derivative. Only a single stage is computed at each time step, and no splitting is needed in this approach. In turn, the proposed numerical algorithm results in an improvement in the computational efficiency. In addition, an alternate temporal stencil is applied , which replaces the approximation 1 4 (u n+1 +2u n +u n−1 ) by a simpler one, 1 2 (u n+1 + u n−1 ), in the linear energy parts. A careful linear analysis reveals an unconditional energy stability for such an alternate temporal stencil, and a linearized stability analysis yields an unconditional convergence estimate for its application to the nonlinear problem. With the help of the aliasing error control technique, an ∞ (0, T * ; H 2 ) convergence for u and ∞ (0, T * ; 2 ) convergence for the discrete time-derivative of the solution are obtained. Various numerical experiments have validated such a theoretical analysis.
The idea of such an alternate temporal stencil could be applied to many other nonlinear hyperbolic equations, in which the wave equation structures are preserved in the linear part.