A High Order Accurate Finite Difference Method for the Drinfel’d–Sokolov–Wilson Equation

We analyse numerically the periodic problem and the initial boundary value problem of the Korteweg-de Vries equation and the Drindfeld–Sokolov–Wilson equation using the summation-by-parts simultaneous-approximation-term method. Two sets of boundary conditions are derived for each equation of which stability is shown using the energy method. Numerical analysis is done when the solution interacts with the boundaries. Results show the benefit of higher order SBP operators.


Introduction
Solitons (or solitary waves) are special kinds of waves with a distinct set of features. Solitons arise due to a balance between nonlinear effects and dispersive effects. Thus, soliton solutions can be found for nonlinear dispersive partial differential equations (PDEs). Solitons are characterized by the following: (i) when propagating at a constant speed, the shape of the soliton remains unchanged; (ii) when a collision occurs between two solitons, both solitons emerge unchanged. In this paper, we study two different PDEs that both exhibit soliton solutions [7].
Originated in the nineteenth century, but it was not until the sixties when the Kortewegde Vries (KdV) equation became a fruitful research object for scientists after they realized its apparent benefits in various applications [25]. The early-known soliton solution of the equation mimics an observed phenomenon of water waves induced by boat engines and initially called by its father the "Waves of Translation". In modern application, the equation is most used to describe the translation of waves with weak dispersion, such as the waves in shallow water which have small amplitude and long wavelength; or the movement of internal flows in some special fluid [13]. It is also helpful in constructing various complex models such as plasma waves, shock waves, tornados, hydrodynamics. Since an exact solution can be derived, the KdV equation is usually useful to verify the accuracy of numerical methods.

The Korteweg-de Vries Equation
We begin by looking at a scalar equation. The initial boundary-value problem to this equation has been analysed in [14]. This section is included to increase readability and to introduce necessary definitions that will be used throughout the paper.
The KdV equation in its most known form can be written as where the subscript notations denote the partial derivatives. The function u depending on the one-dimensional space variable x and the time variable t characterizes the elongation of the wave at point x and time t [25]. The second and third term on the left hand side of (1) is called the nonlinear term and the dispersion term, respectively. An exact solution of (1) is given by (cf, [13]) where c is the wave speed and a is the starting offset. We use the above exact solution to construct the initial condition and calculate the errors for the efficiency analysis. From (2), it is clear that to have a real-valued solution, the wave speed c must be a nonnegative number. Figure 1 illustrates the soliton solution given by (2) for several positive values of c. It can be observed that with those choices of c, the soliton propagates to the right; and its amplitude is proportional to c. In other words, this family of solitary solutions to (1) with higher amplitude move faster than the ones with smaller amplitude.

Continuous Stability
Definition 1 Given two continuous functions u and v, a continuous inner product is defined To analyse continuous stability, we use the so-called "energy method". Given an equation in the form u t = f (u, u x , ...), we multiply both sides of the equation by u T and integrate by parts to obtain (u, u t ). We then add the transpose (u t , u). For the continuous problem to be stable, it requires that there is no energy growth in time, d dt u 2 < 0 for homogeneous boundary conditions. In the context of nonlinear problems, no energy growth in time does not lead to well-posedness as existence, uniqueness and continual data dependence of the solution are not straightforwardly followed. However, considering linear or linearised problems, one could formulate this statement more formally as a sufficient condition for well-posedness, see [5,12]. This property is usually investigated using the relation d dt u 2 = (u t , u) + (u, u t ). To simplify notation, we denote u l = u(x l , t), u r = u(x r , t).
We write (1) in skew-symmetric form, for some constant α. Now we apply the energy method to (3), Adding the transpose (u t , u) yields For the periodic problem, RHS in (4) is identically zero, and (3) is thus stable at the continuous level.

Semi-discrete Stability
Consider a computational domain [x l , x r ] ⊂ R. We discretise the domain using m grid points Definitions of the SBP operators used in this paper are restated in Definition 2 and 3. See [14] for other details on the SBP-SAT approximation technique. Let e 1 = (1, 0, . . . , 0) T , e m = (0, . . . , 0, 1) T be column vectors of length m.
x is said to be a first derivative SBP operator if the matrix H is symmetric positive definite and Q + Q T = 0.
H v is usually utilised. A consistent finite difference approximation of (1), ignoring the boundary treatment, can be written as where the discrete function By Definitions 2 and 3, the SBP operators D 1 , D 3 mimic integration-by-part properties of the corresponding continuous differential operators which leads to the following discrete estimate which exactly mimics (4). Let D 1 and D 3 be periodic operators approximating ∂ ∂ x and ∂ 3 ∂ x 3 , respectively. The stencils of D 1 and D 3 can be constructed in a non-overlapping circular manner. For example, a second order periodic SBP operator approximating the first derivative ∂ ∂ x can be constructed as where the matrices D 1 , H , Q and the identity matrix I are of size (m − 1) × (m − 1), and the solution vector is v = (v 1 , v 2 , . . . , v m−1 ). The first component v 1 approximates v at x 1 = x l and the final component The approximation (7) is stable due to the pure skew-symmetry of D 1 , D 3 . An analysis of runtime-error efficiency is performed on the domain [x l , x r ] = [−20, 20] with m = 51, 101, 201, 301 grid points, wave speed c = 2, time step size t = βh 3 . The periodic scheme (7) is used . Notice that the time step scaling O(h 3 ) comes from the third derivative in (1). The result is shown in Fig. 2. For a fair comparison among different orders of accuracy, sharp estimates of the CFL number β as shown in Table 1 are used. For each accuracy order, starting from a sufficiently large value, the CFL number β is computed by gradually decreasing it until the solution becomes unstable.
A convergence analysis of the periodic scheme (7) is performed on the domain [x l , x r ] = [−20, 20] with m space grid points, the wave speed c = 2 and β = 0.2. The result is shown in Table 2. Let u be the vector containing nodal values of the exact solution. We denote the   error vector as e = u − v, and the l 2 -error as l 2 (e) = e H = e T He. We calculate the convergence rate for two different number of grid points m 1 , m 2 as q = log 10 l 2 (e m 1 ) l 2 (e m 2 ) / log 10 m 2 m 1 .

A Set of Stable 'Dirichlet like' Boundary Conditions
The following set of boundary conditions yields a stable problem to (1) Applying these BCs to (4) leads to, We see that with zero boundary data, g l (t) = g r (t) =g r = 0, the problem (1), (8) is However, since the energy growth is not exclusively bounded by boundary data, (8) is not strongly stable. We refer the readers to [5] for the definition of strong stability.

Lemma 1 A stable SBP-SAT approximation of
Proof The energy method applied to (9) with zero boundary data yields the following energy estimate: We choose the penalty terms which mimics the continuous stability estimate with homogeneous boundary conditions. Snapshots at different time points of the boundary interaction with 'physical' homogeneous data g l (t) = g r (t) =g r = 0 and the total energy over the course of interaction are shown in Fig. 3. The tests are performed with m = 401 grid points using 6th order SBP operators. As can be seen from Fig. 3, after hitting the right 'wall' at t = 10, the solution is reflected and bounces back wave-like without losing energy. Not until t = 15, the left-going waves hit and 'go out of' the left wall. For this test problem, we observe that there is no significant difference in the interaction behavior when using operators of different orders.  Table 3 To verify the accuracy of the approximation (9), time-dependent analytic data is prescribed for g l (t), g r (t),g r (t) in (8), where g l (t), g r (t) are taken directly from the analytic solution and The domain is narrowed to [x l , x r ] = [ −4, 4] so that the solution is initialized non-zero at both boundaries, and the right boundary cuts the soliton in the middle at time t = 2, see Fig. 4. The free SAT penalty parameter τ 4 is set to be -1. Other problem parameters are chosen same to the convergence study in Table 2. The result is shown in Table 3.

A Set of Strongly Stable Boundary Conditions
We show that the following set of boundary conditions yields a strongly stable problem subjected to (1), Denote (11) can then be rewritten as . Inserting (11) and completing the squares leads to It is clear that all the terms on the right hand side are either bounded or less than or equal to zero when |u l | 0, |u r | 0. However, even if |u l | −→ 0 we can still control 4|u l | = 0. Similar argument can be applied for u r to deduce that the right hand side is still bounded. Therefore, it is sufficient to conclude that (11) is a set of strongly stable boundary conditions for (1). We show that there exists real numbers τ l , τ r , σ r such that the following SBP-SAT approximation of (1) subjected to (11) (12) is strongly stable if

Lemma 2 Equation
Proof The energy method applied to (12) leads to, To cancel the cubic terms and the terms involving second derivatives, we choose τ l = −1, τ r = 1. The equation becomes If 2σ r + 1 = 0, a desirable estimate cannot be obtained becauseg r d 1;m v is not bounded. When 2σ r + 1 = 0, we have It can be seen that with 2σ r + 1 < 0, all terms of the above estimate are either bounded or non-positive.
Snapshots at different time points of the boundary interaction with homogeneous data g l (t) = g r (t) =g r = 0 and the total energy over the course of interaction are shown in Fig. 5. The tests are performed with m = 401 grid points using 6th order SBP operators. The behavior is fundamentally different from the previously derived boundary conditions (8).
Most of the solution is absorbed by the right wall, while the small remaining part is reflected. The energy decreases smoothly in time in contrast to the fluctuating downtrends in Fig. 3  There is also no difference in the interaction behaviour between operators of different orders. We observe same interior and boundary accuracy properties for this set of boundary conditions as reported in Table 3.

The Coupled Drinfeld-Sokolov-Wilson Equation
We consider the following system of nonlinear equations, An exact solution of (14) is given as follows [26], where c is the wave speed and a is the starting offset. An illustration is given in Fig. 6 showing the soliton solution at two different time points. Unlike the Korteweg-de Vries equation, a similar technique to perform a skew-symmetric splitting does not straightforwardly give a nonlinear estimate. Instead, we define a linearised problem to (14) as whereū andv are constant fields, or frozen coefficients [20]. Applying the energy method leads to We thus can obtain the energy estimate for each component.
It can be shown that in the standard norm u 2 + v 2 , we cannot obtain an upper bound in the energy since besides the boundary terms, the remainder terms cannot be controlled.
However, in the modified norm 1 3 u 2 + v 2 , it is possible to obtain an estimate, d dt At continuous level, (16) shows that the periodic linearised problem (15) is well-posed. In this paper, we refer to this property of the original nonlinear problem (14) as it is "linearly stable". Considering zero boundary data, the four-by-four matrix in (16) has three eigenvalues. Therefore, the correct number of boundary conditions is three. A consistent SBP approximation to the linearised problem (15) reads where u, v are the discrete variables (here we use the same notation as the corresponding continuous variables) andŪ ,V denote the square diagonal constant matrices of which diagonals contain the discrete analogues of frozen coefficientsū andv, respectively. In that case, we call (17) a "linearised scheme" which we utilise to derive a discrete l 2 estimate. The corresponding "nonlinear scheme" to the original problem (14) is then obtained by replacinḡ U ,V with the discrete variables u, v in the diagonals. We discuss stability of the nonlinear problem (14) and the linearised scheme (17) in the following remark.
Remark Assuming smoothness of u, v, [20] shows that discrete stability of the linearised scheme and the nonlinear scheme are equivalent. Because our main focus is to derive provably stable discretisations, we omit the details on nonlinear well-posedness and nonlinear stability and refer the interested readers to [18] and [11]. The main argument to include here is that under the smoothness assumption, studying the linearised scheme can help guide the discrete boundary treatment of the nonlinear scheme (17). It is important to highlight that the approach of deriving nonlinear schemes through analysing the corresponding linearised problem has been successful for systems such as the Navier-Stokes equations, compressible Euler equations, and shallow water equations [12,16,19].
Applying the energy method to (17) gives where the sum of the third derivatives can be simplified as  Analogously to the continuous estimate (16), we have the discrete energy estimate d dt This estimate shows that in the modified norm, the linearised scheme (17) approximating the periodic problem (15) is stable.
The discrete error is calculated as e =  Table 4. An analysis of runtime-error efficiency is performed on the same domain and problem settings with different number of grid points. The result is shown in Fig. 7. For a fair comparison among different orders of accuracy, sharp estimates of the CFL number β, as shown in Table 5, are used. For each accuracy order, starting from a sufficiently large value, the CFL number β is computed by gradually decreasing it until the solution becomes unstable.

A Set of Stable 'Dirichlet like' Boundary Conditions
A set of linearly well-posed boundary conditions is given by A consistent SBP-SAT approximation of (14) subjected to the boundary conditions (19) reads We discuss stability of (20) in the following lemma.
Snapshots at t = 20 of the boundary interaction and the total energy over the course of interaction are shown in Fig. 8. The numerical experiment is performed on the corresponding nonlinear scheme to (20) with SBP operators of 2nd, 4th and 6th order using the same number of grid points m = 201. Due to the boundary condition v l = v r = (v x ) r = 0, a stationary profile after the interaction appears near the right boundary. In Fig. 8, we cut off showing this  Table 6  stationary pulse because the difference in the interior behavior among different orders is more interesting. The less accurate second order operators generate spurious oscillations which are amplified from the error in capturing the boundary interaction. This so-called "π-mode" pollutes the solution and does not disappear with grid refinement. Though the energy grows after the interaction at around t = 15, it is still bounded in long time simulation as seen in the energy plot (Fig. 8b). The energy plots (Fig. 8b, d, f) numerically verify that our boundary treatment is stable. Notice that the energy jump for 4th and 6th order operators are of lower magnitude. Similar convergence study after boundary interaction is performed in [x l , x r ] = [−4, 4]. The boundary data g l (t), g r (t) in (19) (20) is derived from the exact solution v exact , and Note again that at the final time t = 2, the right boundary cuts the soliton profile in the middle, see Fig. 9. The free parameter σ 2 is set to be −2. The convergence rates are reported in Table 6. Because the boundary conditions (19) only specify the analytic solution to v, optimal rates can be seen for the v component, but degraded rates are obtained on the u component.

Characteristic Boundary Conditions
We denote BT the boundary terms given by (16). Recall that BT = w T Mw Since M is symmetric, it can be diagonalized as where g r (t), g l (t) are 4 × 1 vectors of time-dependent real functions. We apply (21) to the boundary term given by (18) to obtain the continuous estimate which contains the terms that are either bounded by given data or non-positive. Therefore,   Fig. 10a, c, e. The observed behaviour differs from the previous set of boundary conditions (19). As can be compared to the result in Fig. 8, in this case, the left-going reflections are not visible. Although a stationary pulse in variable u is still left on the right boundary, from the energy plots (Fig. 10b, d, f), it can be seen that there is almost no energy growth over time except for a small increase at t = 15.

Conclusion
In this paper, we investigate the numerical properties of solitons by studying two PDEmodels, the Korteweg-de Vries equation and the coupled Drinfeld-Sokolov-Wilson equation, that can be used to describe such phenomenon. The SBP operators for finite difference method are applied to solve the periodic problem. We analyse the involvement of boundary conditions by deriving two different sets of stable boundary conditions for each equation. The implementation of those boundary conditions utilizes the SAT technique, of which accuracy and stability is shown in the numerical result by different orders of accuracy. For the KdV equation, we show nonlinear stability using the energy method and skew-symmetric splitting.
For the DSW equation, we show linear stability in a modified norm. The numerical solution converges to the exact solution with the expected convergence rate. The boundary interaction is illustrated in figures. The numerical study also strengthens the motivation for higher order accurate SBP operators. Future work may consist of applying optimized time-stepping to the discrete approximations. One could also investigate the possibility of obtaining a nonlinear energy estimate for the DSW equation.
Funding Open access funding provided by Uppsala University.
Data Availibility Statement Data sharing not applicable to this article as no datasets were generated or analysed during the current study

Conflict of interest Not applicable.
Code availability The code that support the findings of this study is available from the corresponding author upon reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.