Classical and Multisymplectic Schemes for Linearized KdV Equation: Numerical Results and Dispersion Analysis

: We construct three ﬁnite difference methods to solve a linearized Korteweg–de-Vries (KdV) equation with advective and dispersive terms and speciﬁed initial and boundary conditions. Two numerical experiments are considered; case 1 is when the coefﬁcient of advection is greater than the coefﬁcient of dispersion, while case 2 is when the coefﬁcient of dispersion is greater than the coefﬁcient of advection. The three ﬁnite difference methods constructed include classical, multisymplectic and a modiﬁed explicit scheme. We obtain the stability region and study the consistency and dispersion properties of the various ﬁnite difference methods for the two cases. This is one of the rare papers that analyse dispersive properties of methods for dispersive partial differential equations. The performance of the schemes are gauged over short and long propagation times. Absolute and relative errors are computed at a given time at the spatial nodes used.


Introduction
In this paper, we solve a linearised Korteweg-de-Vries equation with specified initial and boundary conditions. The three methods include classical, multisymplectic and a modified explicit scheme adapted from Wang et al. [1]. In 1895 [2], two Dutchmen, namely Korteweg and de Vries, derived a nonlinear partial differential equation of the form which describes the long time asymptotic behaviour of a small but finite amplitude of one-dimensional shallow water waves. In the equation above, u = u(x, t) measures the elevation (height of water above equilibrium level) at time t and position x while δ 2 is referred to as the dispersion coefficient. There are two different mechanisms that are present, namely; 1. Nonlinearity (uu x ), which tends to steepen those parts having negative slope.

2.
Dispersion, which makes dispersive wave components of different wave frequencies propagate at different velocities.
The delicate balance between these two effects leads to travelling waves of permanent form, the so-called solitary wave. It is usual to refer to the solitary wave as the single soliton solution, but when more than one of them appears in a solution, they are then termed as solitons. If one of these two competing effects is lost, solitons become unstable and eventually cease to exist. In this respect, solitons are completely different from linear waves.
The KdV equation has also been found to describe a number of important physical phenomena such as magnetohydrodynamic waves in a warm plasma [3], elastic waves in an anharmonic crystal [4], ion-acoustic waves in a plasma [5], the long-lived giant red spot in the highly turbulent Jovian atmosphere and the propagation of short laser pulses in optical fibres [6].
The existence of unique solutions for some classes of initial data has been established in Lax [7]. There exist many integral invariants of the KdV equation, three of these serve as benchmarks to test the efficiency of numerical solvers. If u(x, t) is the solution of the KdV equation given by Equation (1) above, the three invariants are which represent mass, momentum and energy conservation, respectively. Li and vu-Quoc [8] stated that in some areas, the ability to preserve the invariant properties of the original differential equations is a criterion to judge the success of a numerical simulation. The discrete forms of the three conservation laws are There are several good reasons to work with KdV equation as a prototype. The reasons are detailed below:

1.
It is a model of a nonlinear hyperbolic equation with smooth solutions for all times.

2.
It is non-dissipative and therefore is a natural testbed for comparing conservative vs. dissipative discretization.

3.
It is notorious. It is well known that unexpected, nonlinear instabilities occasionally arise from reasonable-looking finite difference method. For instance, Zhao and Qin [9] solved u t + ηuu x + δ 2 u xxx = 0 with u(t = 0, x) = u 0 (x) = cos(πx) for x ∈ [0, 2], η = 1, δ = 0.022 using periodic boundary conditions and employing the Zabusky-Kruskal scheme. They tried various ∆x, ∆t combinations satisfying the linear stability bound, yet they always obtained blow-up phenomena for t > 21 π . Setting ∆x = 0.01, ∆t = 0.0001, the solution, using an explicit scheme, is qualitatively correct for a while. Around t = 5, i.e., after 50,000 time steps, error accumulates in the solution so that the linear stability bound is violated. After a few more steps, the solution blows up.
The design and development of symplectic methods for Hamiltonian ODEs has yielded very powerful numerical schemes with beautiful geometric properties. Symplectic and other symmetric methods have been noted for their superior performance, especially for long time integration [10][11][12]. The Korteweg-de-Vries equation has been extensively studied in numerous studies using symplectic and multisymplectic methods [13]. Because of nonlinear instability, all those methods that can remove the phenomenon of unphysical oscillations are entirely implicit or semi-explicit except for the 6-point scheme that was proposed very recently in [14] based on the concept of multisymplectic schemes. However, due to stability constraint on the time step, the 6-point scheme is somewhat slow. Wang et al. [1] attempted to construct a scheme that not only removes oscillation phenomenon for long time propagation but is also faster than the 6-point scheme.
Very recently, the authors in [15] used two existing schemes proposed by Zabusky and Kruskal [4] and Wang et al. [1] to solve u t + 6uu x + βu xxx = 0 with initial conditions u(x, t = 0) = 2µ 2 sech 2 (µx) with µ > 0 [16]. Appadu et al. [15] also constructed two novel methods obtained by modifying the scheme proposed by Zabusky and Kruskal. The performance of the four methods is compared in regard to dispersion and dissipation errors and ability to conserve mass, momentum and energy by using two numerical experiments that involve solitons. It is worthy to note that sine-cosine and canonical transformation methods proposed and employed in [17,18] have proven useful in obtaining exact soliton solutions and solutions of the Navier-Stokes equation, respectively. We next discuss dispersive characteristics of numerical methods.
The term with the lowest even order spatial derivative in the truncation error produces amplitude error in the numerical solution, and this is responsible for numerical dissipation. The leading odd spatial derivative in the truncation error produces small-scale waves as different Fourier components propagate at different phase speeds, and this causes numerical dispersion [19]. The relative phase error (RPE) is a measure of the dispersive character of a scheme. This quantity is a ratio and measures the velocity of the computed waves to that of the physical waves [20,21]. The relative phase error is obtained using the relation [22] where ξ num is the amplification factor of the numerical scheme and ξ exact is the exact amplification factor. The relative phase error of numerical schemes discretizing 1D linear advection equation (u t + βu x = 0) and 1D advection-diffusion equation Plots of relative phase error vs. w at some values of ∆x, ∆t for a few schemes discretizing 1D linear advection and 1D advection-diffusion equation were obtained in [21] and [24], respectively. Plots of RPE vs. w x vs. w y were obtained in [23,25] for numerical schemes discretizing 2D advection-diffusion and 3D advection-diffusion equations. Ascher and McLachlan [26] obtained the dispersion relation of the partial differential equation On considering the linearized version of the PDE Appadu et al. [15] considered the linearized form u t + βu xxx = 0 of the PDE u t + γuu x + βu xxx = 0 to study dispersion analysis of the KdV equation. Some work on computing the optimal temporal step size by minimizing the integrated error at a given ∆x was done by Appadu et al. [27]. The integrated error was obtained as the square of dispersion error.
The paper is organized as follows. The two numerical experiments considered are described in Section 2. Section 3 is devoted to Scheme 1 applied to numerical experiment 1 and details the following: derivation, stability, consistency, presentation of numerical results and dispersion analysis. Sections 4 and 5 give information on work done on Schemes 2 and 3 when used to solve numerical experiment 1. Sections 6-8 concerns the derivation, stability, consistency, presentation of numerical results and dispersion analysis of Schemes 1, 2, 3 when solving numerical experiment 2. We highlight salient features of the paper in Section 9.

Numerical Experiment
This paper is dedicated to the analysis of three numerical methods for the approximation of two different forms of linear KdV equations.

Case 2 (Numerical Experiment 2)
We considered a new case when the coefficient of dispersion was greater than that of advection. We solve The initial condition is u(x, 0) = sin(x), and the boundary conditions are u(0, t) = sin(3t) (16) and u(2π, t) = sin(2π + 3t).

Stability of Scheme 1 for Numerical Experiment 1
In this section, the stability region of Scheme 1 is obtained. We substitute v n i by ξ n e Iθih in Equation (18). This gives After some rearrangements, we get where Therefore, We denote the amplification factor of the physical and the computational modes by ξ 1 and ξ 2 , respectively. and If we seek amplification factor ξ that is not purely imaginary, we require 4 − A 2 ≥ 0, which implies that |A| ≤ 2. This implies that we need k such that Note that the inequality is dominated by 2k . We obtain the stability region as If we choose h = π/10, (22) gives k ≤ 0.012775. We can also obtain the stability region by use of 2D plots, as shown in Figures 1 and 2. The stability region when h = π/10 is 0 < k ≤ 0.012, as illustrated in Figures 1b and 2b.
The procedure used in this section is the von Neumann approach. This procedure will be employed to determine the expression for the amplification factor of the schemes in the later sections. The region is obtained by imposing |ξ ≤ 1|.

Consistency of Scheme 1 for Numerical Experiment 1
The scheme in Equation (18) is given by We use the Taylor series about the point (n, i) in order to determine the order of accuracy of the scheme. By the Taylor's expansion, The accuracy of the scheme is quadratic both in space and in time.
The detailed consistency analysis is discussed above. The same procedure will be followed throughout the entire paper. We may not discuss it in detail for other methods.

Numerical Results
In this section, the approximation of the linear dispersive Equation (12) using Scheme 1 given in Equation (18) is presented. Within the stability region, the spatial and the temporal step sizes employed are h = π/10 and k = 0.001, respectively. The solution profiles are shown in Figure 3. These are compared with the exact solution u(x, t) = sin(x − t). The profile of absolute errors vs. x is also presented at T = 2.0 and T = 4.0. Figure 4 displays the absolute error profiles for the scheme. Table 1

Error
Step Sizes k

Dispersion Analysis
where α is the dispersion relation. We consider Substituting these values in Equation (12), we obtain This implies that the perturbation for u(x, t) can be written as a function of θ alone as e Iθ[(θ 2 −2)t+x] . The amplification factor is determined from the relation Relative phase error (RPE) is defined [27] as We obtain plots of the arguments of ξ exact , ξ 1 , ξ 2 vs. w ∈ [0, π] in Figure 5. For values of w ∈ [0, 1], the graphs of the arguments of both ξ exact and ξ 1 are very close to each other, as illustrated in Figure 5a. For values of w ∈ [0, 0.5], the graphs for the arguments of ξ exact and ξ 2 are very close to each other, as depicted in Figure 5b, and we note that arg(ξ 2 ) and arg(ξ exact ) are of opposite signs for 0.5 < w < π.

Scheme 2 for Numerical Experiment 1
In 2007, Wang et al. [14] constructed the following scheme to solve where η is the coefficient of the nonlinear term uu x . The scheme is multisymplectic and can remove the dispersive oscillations and preserves approximately several conservation laws of the KdV equation [14].

Stability of Scheme 2 for Numerical Experiment 1
Substituting v n i by ξ n e Iθih and simplifying gives After some rearrangements, we get where and therefore, we get The range of values of the temporal step size k for which |ξ| ≤ 1 is shown in the profiles of the amplification factor ξ (see Figures 6 and 7) given that w = θh ∈ [−π, π] and h = π 10 . It is observed that

Consistency of Scheme 2 for Numerical Experiment 1
We recall that the scheme is given by Taylor series expansion about (n, i) after simplification gives The last equation shows that the scheme is second order accurate both spatially and temporally.

Numerical Results
In this section, the approximation of the linear dispersive Equation (12) by Scheme 2 given in (25) is presented. Within the stability region, the spatial and the temporal step size employed are h = π/10 and k = 0.0001, respectively. The solution profiles are shown in Figure 8. These are compared with exact solution u(x, t) = sin(x − t). The profile of absolute errors vs. x is also presented at T = 2.0 and T = 4.0. Figure 9 displays the absolute error profiles for the multisymplectic scheme. We display L 1 and L ∞ errors in Table 2 (a) T = 2.0 (b) T = 4.0

Error
Step Sizes k T = 2 T = 4

Dispersion Analysis
Arguments of numerical amplification factors ξ 1,2 derived in Section 4.1 are compared with the exact amplification factor within the stability region. The comparison is shown in Figure 10. We obtain plots of the arguments of ξ exact , ξ 1 , ξ 2 vs. w ∈ [0, π] in Figure 10. For values of w ∈ [0, 1], the graphs for the arguments of both ξ exact and ξ 1 are very close to each other, as illustrated in Figure 10a. Figure 10b shows that the arguments of ξ exact and ξ 1 are not close for even small values of w.
(a) ξ 1 (b) ξ 2 Figure 10. Plot of the arguments of amplification factors of Scheme 2 for Numerical Experiment 1 with spatial and temporal step sizes h = π/10 and k = 0.001, respectively

Stability of Scheme 3 for Numerical Experiment 1
Substituting v n i by ξ n e Iθih and simplifying gives After some rearrangements, we get where Therefore, We employ the condition |ξ| ≤ 1 to determine the appropriate range for the temporal step size when the spatial step size is π 10 with w ∈ [−π, π]. The temporal step size is based on the behaviour of both . The region of stability, as shown in Figures 11 and 12, implies that the scheme is stable for all 0 < k ≤ 0.001.

Consistency of Scheme 3 for Numerical Experiment 1
In this section, the pointwise accuracy of Scheme (26) will be discussed. The scheme is written as Using Taylor's series expansion about (n, i) gives After some simplifications, we get, We conclude here that the scheme is second order accurate both in space and time.

Numerical Results
In this section, the approximation of the linear dispersive Equation (12) using Scheme 3 is presented. Within the stability region, the spatial and the temporal step size employed are h = π/10 and k = 0.001, respectively. The solution profiles are shown in Figure 13. These are compared with exact solution u(x, t) = sin(x − t). The profile of absolute errors is also presented at T = 2.0 and T = 4.0. Figure 14 displays the absolute error profiles for the multisymplectic scheme. We compare L 1 and L ∞ errors in Table 3.  Table 3. L 1 and L ∞ errors for various k when h = π/10 when Scheme 3 is employed to approximate Numerical Experiment 1 at T = 2 and T = 4.

Error
Step Sizes k T = 2 T = 4

Dispersion Analysis
The arguments of the numerical amplification factors ξ 1,2 derived in Section 7.1 are compared with the exact amplification factor within the stability region. The profiles of the arguments of the amplification factors are shown in Figure 15. with the spatial and temporal step sizes h = π/10 and k = 0.001, respectively.

Scheme 1 for Numerical Experiment 2
A classical finite difference scheme given by v n+1 was proposed by Zabusky and Kruskal [4] for the solution of the nonlinear equation u t + uu x + δ 2 u xxx = 0. The scheme will be adapted to solve u t + 2u x + 5u xxx = 0 and will be referred to as Scheme 1. Scheme 1 is given by where i = 1, 2, 3, . . . , NP − 3 and n = 2, 3, . . . , itmax − 1.

Stability of Scheme 1 for Numerical Experiment 2
The amplification factor of this scheme satisfies the equation Therefore, If we seek the amplification factor that is not purely imaginary, we require 4 − A 2 ≥ 0 which implies that |A| ≤ 2. This implies that we need k such that Note that the inequality is dominated by 10k h 3 (sin(2w) − 2 sin(w)) as h → 0. The maximum of this is at the point cos(w) = − 1 2 . The stability region is If we choose h = π/10, we obtain k ≤ 0.00241869. The range of values of the temporal step size k for which |ξ| ≤ 1 is also shown in the profile of the amplification factors ξ shown in Figures 16 and 17 given that w = θh ∈ [−π, π] and h = π 10 . The profile is shown both . It can be deduced from the two profiles in Figures 16 and 17 that 0 < k ≤ 0.0025.

Consistency of Scheme 1 for Numerical Experiment 2
We use Taylor's series expansion about (n, i) in order to determine its order of accuracy. We obtain and we conclude that the accuracy of the scheme is quadratic both in space and in time.

Numerical Results
In this section, the approximation of the linear dispersive Equation (15) by Scheme 1 (29) is presented. Within the stability region, the spatial and the temporal step sizes employed are h = π/10 and k = 0.001, respectively. The solution profiles are shown in Figure 18. These are compared with exact solution u(x, t) = sin(x + 3t). The profile of absolute errors is also presented at T = 2.0 and T = 4.0. Figure 19 displays the absolute error profiles of Scheme 1 for Numerical Experiment 2. We compare L 1 and L ∞ errors in Table 4.  Table 4. L 1 and L ∞ errors for various k when h = π/10 when Scheme 1 is employed to approximate Numerical Experiment 2 at T = 2 and T = 4.

Error
Step Sizes k T = 2 T = 4

Dispersion Analysis
We consider u t + 2u x + 5u xxx = 0. We use same techniques as in Section 3.4 The amplification factor is determined from the relation The relative phase error (RPE) is defined [27] We obtain plots of the arguments of ξ exact , ξ 1 , ξ 2 vs. w ∈ [0, π] in Figure 20. For values of w ∈ [0, 1] the graphs of the arguments of both ξ exact and ξ 1 are very close to each other, as illustrated in Figure 20a.

Scheme 2 for Numerical Experiment 2
Wang et al. [14] proposed 1 2k ( For the approximation of the nonlinear KdV equation u t + ηuu x + δ 2 u xxx = 0, we adapt the method in Equation (30) to approximate Equation (15). The proposed method is 1 2k ( Analysis the stability, consistency and the dispersion properties of the Scheme given by ((31)) are discussed below.

Consistency of Scheme 2 for Numerical Experiment 2
In this section, the consistency of the numerical scheme (31) is discussed. We consider Equation (31) and use Taylor's series expansion about (n, i) to get We note that v tx + 2v xx + 5v xxxx ≈ 0 and we therefore obtain We deduce that the scheme is of second order accuracy in both time and space.

Numerical Result
In this section, the approximation of the linear dispersive Equation (15) by Scheme 2, given in Equation (31), is presented. Within the stability region, the spatial and the temporal step sizes employed are h = π/10 and k = 0.001, respectively. The solution profiles are shown in Figure 23. These are compared with exact solution u(x, t) = sin(x + 3t). The profile of absolute errors are presented at T = 2.0 and T = 4.0. Figure 24 displays the absolute error profiles. We display L 1 and L ∞ errors in Table 5.  Table 5. L 1 and L ∞ errors for various k when h = π/10 when Scheme 2 is employed to approximate Numerical Experiment 2 at T = 2 and T = 4.

Error
Step Sizes k T = 2 T = 4

Dispersion Analysis
We plot the arguments of the numerical amplification factors as derived in Section with the exact amplification factor in Figure 25. The argument of ξ 1 and the argument of exact amplification factor are quite close to each other for w ∈ [0, 1].
(a) ξ 1 (b) ξ 2 Figure 25. Plot of the arguments of amplification factors of Scheme 2 for Numerical Experiment 2 with the spatial and temporal step sizes h = π/10 and k = 0.001, respectively.

Scheme 3 for Numerical Experiment 2
Wang et al. [1] proposed the scheme for the numerical approximation of the nonlinear equation. We adapt this scheme to solve the linear dispersive equation u t + 2u x + 5u xxx = 0. The scheme for Equation (15)  ( In explicit form, after some index shift, this equation is written as

Stability of Scheme 3 for Numerical Experiment 2
We substitute v n i by ξ n e Iθih , and after simplification we get This is a quadratic equation in ξ, written as The amplification factors ξ 1,2 are, therefore, obtained as The stability constraint for the scheme is computed from each of these amplification factors using h = π 10 for θh = w ∈ [−π, π]. The stability region is deduced from Figures 26 and 27 as 0 < k ≤ 0.001.

Consistency of Scheme 3 for Numerical Experiment 2
To discuss the consistency of the scheme in Equation (35), consider the Taylor expansion about (n, i). After simplification, we get hence, the scheme is spatially and temporally accurate of order 2.

Numerical Results
In this section, the approximate solution of the linear dispersive Equation (15) by Scheme 3, given in Equation (35), is presented. Within the stability region, the spatial and the temporal step sizes employed are h = π/10 and k = 0.001, respectively. The solution profiles are shown in Figure 28. These are compared with exact solution u(x, t) = sin(x + 3t). The profile of absolute errors is also presented at T = 2.0 and T = 4.0. Figure 29 displays the absolute error profiles for the multisymplectic scheme. We display L 1 and L ∞ errors in Table 6.  Table 6. L 1 and L ∞ errors for various k when h = π/10 when Scheme 3 is employed to approximate Numerical Experiment 2 at T = 2 and T = 4.

Error
Step Sizes k T = 2 T = 4

Dispersion Analysis
The arguments of the numerical amplification factors ξ 1,2 are compared with the argument of the exact solution of the linear dispersive Equation (15). The amplification factors are derived in Section 8.1. The comparisons are shown in Figure 30.
(a) ξ 1 (b) ξ 2 Figure 30. Plot of the arguments of amplification factors of Scheme 3 for Numerical Experiment 2 with the spatial and temporal step sizes are h = π/10 and k = 0.001, respectively.

Conclusions
This paper considers, proposes and analyzes three finite difference schemes for two linear dispersive KdV equations. One of the cases is such that the advective term dominates the dispersive term, while in the second case, the dispersive term dominates the advective term. Stability regions were derived for the schemes, their accuracies were discussed and their dispersion analyses were conducted and presented.
It is observed that for Numerical Experiment 1, the stability region for Scheme 3, 0 < k ≤ 0.001, is narrower compared to the stability region 0 < k ≤ 0.01 for Scheme 1 and Scheme 2. The maximum absolute error of Scheme 1 is relatively less than the maximum absolute error for Scheme 2 and Scheme 3 both for short time and long time integration. In addition, the stability restraint is relaxed when Scheme 1 is employed to approximate Numerical Experiment 2 than Scheme 2 and Scheme 3. In this case, the maximum absolute error, when Scheme 2 is employed to solve Numerical Experiment 2 is less compared to Scheme 1 and Scheme 3. In fact, the highest error is recorded for Scheme 1.
A comparison of results of classical finite difference methods with Laplace Adomian Decomposition Method for Numerical Experiment 1 was done in [30], and it was observed that at a low propagation time of 0.1 s, LADM is better than the classical scheme but that at long time propagation (say when time = 4.0), the performance of LADM deteriorates significantly.
We intend to use the knowledge acquired in this paper to construct methods for linearized KdV equations in both 1 − D and 2 − D with more challenging initial conditions and check which of the three schemes is the best performing. Moreover, we can also construct similar methods for PDEs quite related to the KdV equations, especially those with uu x and u xxx terms present and these examples of PDEs are fractional KdV, modified KdV, KdV-Burgers-Kuramoto, KdV-Burgers and stochastic KdV equations.