The Construction and Research of the Modified “Upwind Leapfrog” Difference Scheme with Improved Dispersion Properties for the Korteweg–de Vries Equation

Abstract

This paper covers the construction and research of a scheme to solve the problem with nonlinear dispersion wave equations, described by the model Korteweg–de Vries equation. The article proposes approximating the equation based on improved “Upwind Leapfrog” schemes. Its difference operator is a linear combination of operators of the “Standard Leapfrog” and “Upwind Leapfrog” difference schemes, while the modified scheme is obtained from schemes with optimal weight coefficients. Combining certain values of the weight coefficients mutually compensates for approximation errors. In addition, the modified scheme acquires better properties compared with the original schemes. The results of test calculations of solutions of the nonlinear Korteweg–de Vries equation are presented, illustrating the advantages of the modified scheme.

1. Introduction

Nonlinear wave processes are intensively researched in the various fields of science and technology, including optics, radiophysics and hydrodynamics [1]. Mathematical modeling of issues related to wave processes on the free surface of a heavy liquid leads to hydrodynamic models that consider both nonlinear and dispersion effects [2]. The quasilinear transfer equation is used as a means of testing the methods of theoretical investigation of systems of nonlinear hyperbolic equations, as well as methods of their numerical implementation. One of the generalizations of the quasilinear transport equation for the case of dispersion media is the Korteweg–de Vries equation (KdV) [3]. The equation has soliton solutions localized in space, represented by a classical solitary wave with a single hump and an infinite period, monotonically decaying at infinity [4]. Among the solutions of the KdV equation, both generalized solitary waves and solitary wave packets can be distinguished.

Many physical processes are described by partial differential equations, so the solution of this class of equations is one of the most urgent and complex problems. The complexity is due to the fact that general solution methods have been developed for linear differential equations, such as Fourier, Laplace and others; there are no general solution methods for nonlinear partial differential equations. Therefore, it is necessary to develop and design special solution methods for each type of nonlinear equation or groups of the same type of equation. A wide range of approaches for the numerical integration of the equations of long nonlinear waves is presented in the literature [5,6,7,8].

In [9], a number of problems of the theory of quasilinear equations are considered. Nonstationary solutions of the Cauchy problem are obtained for the model equation; they consider nonlinearity, dispersion and dissipation. The equation under study is used in a wide range of problems of continuum mechanics; in particular, it is used to describe the propagation of nonlinear longitudinal waves in rods. Based on the model equation, complex behavior of traveling waves was revealed, which can be regarded as discontinuity structures in solutions of the same equation without considering dissipation and dispersion. This led to the ambiguity of solutions to standard self-similar problems consisting of a sequence of Riemann waves and shock waves that have a stationary structure. On the basis of the KdV model equation, the interaction of nonlinear waves (when moving towards each other or one after the other) is also studied when the corresponding self-similar problems of collision of discontinuities have a nonunique solution. In the work of S.K. Godunov [10], situations are considered when, as a result of the interaction of nonlinear waves at long times, asymptotics containing discontinuities with an unsteady periodic oscillatory structure are formed.

The classical and modified KdV equations have a significant role in the physics of nonlinear waves in view of their integrability. Other equations from the same family, derived from different physical and technical applications, are less well known. These are the modular KdV equation, the lognormal KdV equation, the Gardner equation, the Shamil equation, the Benjamin–Ono and Kawahara equations, as well as fractional equations. These equations differ in the degree of nonlinearity in the advective term and the order of linear variance (including fractional). The paper [11] is devoted to the discussion of practical applications of the KdV family of equations and the general properties of their solutions. The research methodology is based on the theory of weakly nonlinear waves in weakly dispersed media of various physical nature. It is shown that many KdV equations of the hierarchy have common properties of solutions and do not have explosive instability. For this, the degree of nonlinearity should not be very large.

A generalized scale-invariant analogue of the KdV equation is considered in [12]. The auxiliary equation method allows us to find expressions for a variety of solutions—both bounded and singular—which have explicit closed-form expressions. Generalized solutions of the KdV equation [13] are discussed in [14]. It is not by chance that the problem of waves on the water surface has become central in the development of the theory of nonlinear waves. It belongs to the most significant application of nonlinear hydrodynamics. Waves on the free water surface have always been an important subject for research, since they are a well-known and, at the same time, quite a complex phenomenon that is easily observable, but quite difficult to describe. For example, the original model equations of KdV and Boussinesq were originally derived as an approximation for surface waves. The review [15] provides the derivation of model equations with justification and the use of asymptotic methods where necessary. It is necessary to understand what place the equation in question occupies in the hierarchy of approximations that have a physical meaning, which in turn are consequences of the precise formulation of the chosen problem about waves on the water surface.

Research on the convective flows of a viscous incompressible fluid is one of the most common problems in a variety of theoretical and applied scientific disciplines. In recent decades, interest in research on solutions describing natural convection has increased significantly. Convection is an example of the self-organization of dynamic systems observed during experiments. The first named structure of self-organization also refers to convection—the Benard cells called hexagonal prisms with a liquid rise in the center. During hydrodynamic experiments, Bernard suggested that an important role in the occurrence of convection is played not only by gravity, but also by the thermocapillary Marangoni effect. It not only explains the reasons for the generation of thermocapillary convection, but also the mechanisms of the occurrence of bioconvection and concentration convection. The paper of Velarde M., Nepomnyashchy A. and Hennenberg M. is devoted to research on the mechanism of the occurrence of oscillatory interphase instability and wave movements in the Benard layers [16]. Research on the fourth-order nonlinear Schrödinger equation (NSEIV) was performed by Sedletsky Y. [17] for the envelope of Stokes waves on the water surface of finite depth. The NSEIV describes the amplitude modulations of the fundamental harmonic of Stokes waves on the surface of the medium and deep (compared to the wavelength) liquid layer. It is derived based on the multi-scale method used to solve equations of mathematical physics. The added terms of the equation describe the effect of third-order linear dispersion and the effects of nonlinearity dispersion. As the nonlinearity and dispersion decrease, the equation is uniformly transformed into the NSE for Stokes waves on the surface of a liquid of finite depth, first obtained by Hashimoto and Ono. The coefficients of the resulting equation are given explicitly as functions of $k\cdot h$, where $k$ is the wave number and $h$ is the water depth. If $k\cdot h$ tends to infinity, then these coefficients are converted into the NSEIV coefficients, first obtained for infinite depth.

The research of Popov S.P. is devoted to investigating the properties of two-dimensional soliton solutions of the Lidke–Spachek type evolutionary equation [18]. This equation is a generalization of the two-dimensional case of the KdV equation. An algorithm for the numerical solution of an evolutionary equation of the Lidke–Spachek type is described. The calculation data of the structure of localized soliton solutions and the process of their interaction with plane solitons are discussed. The results of numerical experiments are compared with the solutions of the Kadomtsev–Petviashvili equations of the first kind and the Zakharov–Kuznetsov equation, which is similar in structure to the considered evolutionary equation.

In the paper of S.J.D. D’Alessio, J.P. Pascal, E. Ellaban, and C. Ruyer-Quil, the liquid flow stability under the influence of gravity, loaded with a soluble surfactant down a heated slope, is investigated [19]. The Marangoni instability associated with settling films containing heated surfactants is researched. A mathematical model is proposed that considers the variations of the surface tension depending on the surfactant concentration and temperature. The linear stability is analyzed both asymptotically for small wave numbers and numerically for arbitrary wave numbers. An expression for the critical Reynolds number is obtained that considers thermocapillary and salt capillary effects and reduces to known published results for special cases. The nonlinear model was solved numerically to research the instability of the equilibrium flow and the development of the emerging constant surface waves. Studies have shown that the results of the numerical implementation of the nonlinear model are in good agreement with the analysis of linear stability.

Analysis of the works of other authors has shown that among the main methods, methods based on the pseudo-spectral approach can be distinguished, as well as those based on the use of explicit and implicit finite-difference schemes.

During the numerical solution of hydrodynamic equations, discontinuous solutions may arise and propagate. This imposes certain requirements on the design of optimal difference schemes that can be used in a wide range. To study waves of small but finite amplitude in dispersive media, it becomes necessary to develop numerical methods with good dissipative and dispersive properties; when implementing the study, nonstrong, nonphysical oscillations in the solution are required to achieve at least the second order of accuracy.

In this regard, this article presents a solution to the problem of nonlinear wave dynamics of solitons, including the Korteweg–de Vries equation, using the example of a modified “Upwind Leapfrog” scheme. A detailed analysis of this scheme was performed in [20,21].

2. Materials and Methods

Test problems of a special type are used when developing difference schemes for nonlinear problems. The KdV equation is considered as a test problem to research the waves of small but finite amplitude in a dispersion environment.

2.1. The KdV Mathematical Model

In the rectangle ${\overline{Q}}_T=\left\{\left(x,t\right):0\le x\le l,0\le t\le T\right\}$ for the KdV equation:

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+\beta \frac{{\partial }^{3}u}{\partial x^3}=0$$

(1)

or

$$u_t+u{u}_x+\beta u_{xxx}=0, \beta =const>0,$$

It is necessary to define a solution $u\left(x,t\right)$ for the Cauchy problem under the following conditions:

at the initial time moment $t=0$:

$$u\left(x,0\right)=u_0\left(x\right), 0\le x\le l$$

(2)

on the boundaries:

$$u\left(0,t\right)=u\left(l,t\right)=u_{xx}\left(0,t\right)=0, t\ge 0$$

(3)

or periodic conditions in the space:

$$u\left(0,t\right)=u\left(l,t\right),\hspace{0.17em}\hspace{0.17em}{u}_x\left(0,t\right)=u_x\left(l,t\right),\hspace{0.17em}\hspace{0.17em}{u}_{xx}\left(0,t\right)=u_{xx}\left(l,t\right),t\ge 0$$

(3a)

where $u\left(x,t\right)=v+\delta v\left(x,t\right)$; $v=const$ is the average velocity; $\delta v\left(x,t\right)$ is the oscillatory velocity; $\beta$ is a dispersion coefficient.

2.2. Numerical Solution of the KdV Equation Based on the “Left Corner” Difference Scheme

Let us define a numerical solution to the model problem (1)–(3) using the “left corner” finite-difference scheme [22]. We will use this problem as a test one.

A uniform computational grid on the domain of the problem definition $\omega ={\omega }_{\tau }\times {\overline{\omega }}_h$ is introduced, where ${\omega }_{\tau }=\left\{t^n|n=0,1,\hspace{0.17em}\dots ,T\right\}$, with the constant time step $\tau =t_{n+1}-t_n$, $T$ is the time steps number, ${\overline{\omega }}_h=\left\{x_i|x_i=ih;i=0,\hspace{0.17em}1,\hspace{0.17em}\dots ,\hspace{0.17em}N;Nh=l\right\}$ with the step $h$ by a spatial variable, $N$ is the number of steps in the space.

The “left corner” finite-difference scheme for the KdV equation has the form:

$$\frac{u_i^{n+1}-u_i^n}{\tau }+u_{i-1/2}\frac{u_i^n-u_{i-1}^n}{h}+\beta \frac{u_{i+1}^n-3u_i^n+3u_{i-1}^n-u_{i-2}^n}{h^3}=0$$

(4)

where $u_{i-1/2}=\left(u_i+u_{i-1}\right)/2$.

2.3. Research the “Left Corner” Scheme Stability

We investigate the stability of the “left corner” scheme (4) for the KdV equation. Define $u_i^n={\phi }^n\cdot e^{jki}$, where $j=\sqrt{-1}$. We substitute this expression inequality (4), and divide the equality by ${\phi }^n\cdot e^{jki}$:

$$\frac{\phi -1}{\tau }+u\cdot \frac{1-e^{-jk}}{h}+\beta \cdot \frac{e^{jk}-3+3\cdot e^{-jk}-e^{-2jk}}{h^3}=0$$

After the transformations, this expression will take the form:

$$\phi =1-\frac{u\tau }{h}\left(1-\cos k\right)-\frac{u\tau j}{h}\sin k-\beta \tau \frac{4\cos k-3-\cos 2k-2j\sin k+j\sin 2k}{h^3}=0$$

If we introduce the notations, we get: $a=\frac{u\tau }{h}, b=\frac{\beta \tau }{h^3}$.

$$\phi =1-a\left(1-\cos k\right)-a\hspace{0.17em}j\sin k-b\left(4\cos k-3-\cos 2k\right)+b\hspace{0.17em}j\left(2\sin k-\sin 2k\right)$$

We select the real and imaginary parts:

$$\phi =1+\left(\cos k-1\right)\left(a+2b\left(\cos k-1\right)\right)-j\left(a+2b\left(\cos k-1\right)\right)\sin k$$

If we introduce the notation $c=a+2b\left(\cos k-1\right)$, we get:

$$\phi =1+\left(\cos k-1\right)c-j\hspace{0.17em}c\sin k$$

A necessary and sufficient condition for the stability of the equation is the fulfillment of inequality $\left|\phi \right|\le 1$.

$$\begin{gathered} {\left|\phi \right|}^2={\left(1-c+c\cos k\right)}^2+c^2{\sin }^{2}k={\left(1-c\right)}^2+2c\cos k\left(1-c\right)+c^2 \\ \le {\left(1-c\right)}^2+2c\left(1-c\right)+c^2=1 \mathrm{at} 0\le c\le 1 \mathrm{or} 0\le \frac{u\tau }{h}+2\frac{\beta \tau }{h^3}\left(\cos k-1\right)\le 1 \end{gathered}$$

Therefore, the scheme is stable under the conditions $\beta \le u{h}^2/4, \frac{u\tau }{h}\le 1$. The disadvantage of the “left corner” scheme is large dissipation. To eliminate large grid viscosity, it is necessary to reduce the step $h$ of the spatial variable, which leads to restrictions that are more stringent on the parameter $\beta .$ As a result, the class of problems for which the scheme of the form (4) is applicable is narrowed.

2.4. “Upwind Leapfrog” and “Standard Leapfrog” Difference Schemes for the KdV Equation

The basic constructive aspects of developing effective grid schemes of the “Upwind Leapfrog” type with improved dispersion properties [9,10] for the nonlinear KdV equations are illustrated by solving the test problem (1)–(3).

A uniform computational grid on the domain of the problem definition $\omega ={\omega }_{\tau }\times {\overline{\omega }}_h,$ is introduced, where ${\omega }_{\tau }=\left\{t^n|n=0,1,\hspace{0.17em}\dots ,T\right\}$, with the constant time step $\tau =t_{n+1}-t_n$, $T$ is the time steps number, ${\overline{\omega }}_h=\left\{x_i|x_i=ih;i=0,\hspace{0.17em}1,\hspace{0.17em}\dots ,\hspace{0.17em}N;Nh=l\right\}$ with the step $h$ by the spatial variable, and $N$ is the number of steps in the space.

The numerical solution of the problem (1)–(3) can be carried out using finite-difference schemes [23,24]:

the “Upwind Leapfrog” scheme:

$$\frac{u_i^{n+1}-u_i^n}{2\tau }+\frac{u_{i-1}^n-u_{i-1}^{n-1}}{2\tau }+u_{i-1/2}^n\frac{u_i^n-u_{i-1}^n}{h}+\beta \frac{u_{i+1}^n-3u_i^n+3u_{i-1}^n-u_{i-2}^n}{h^3}=0$$

(5)

the “Standard Leapfrog” scheme:

$$\frac{u_i^{n+1}-u_i^{n-1}}{2\tau }+u_{i+1/2}^n\frac{u_{i+1}^n-u_i^n}{2h}+u_{i-1/2}^n\frac{u_i^n-u_{i-1}^n}{2h}+\beta \frac{u_{i+2}^n-2u_{i+1}^n+2u_{i-1}^n-u_{i-2}^n}{2h^3}=0$$

(6)

In the “Standard Leapfrog” scheme for a finite-difference approximation of a third-order derivative, a symmetric difference derivative is used, obtained as a half-sum of the right and left difference derivatives, where the approximation of the second derivative is used as a function.

The extended templates of these schemes, considering the approximation of the third derivative by a spatial variable, are given in Figure 1 and Figure 2. Both schemes have three layers.

To solve the KdV model equation, we consider the application of a modified scheme, which is a linear combination of the difference scheme “Upwind Leapfrog” with weight coefficients 2/3 and “Standard Leapfrog” with weight coefficients 1/3 [21]. The approximation of Equation (1) has the form:

$$\begin{gathered} \frac{u_i^{n+1}-u_i^n}{\tau }+2\frac{u_{i-1}^n-u_{i-1}^{n-1}}{3\tau }+\frac{u_i^n-u_i^{n-1}}{3\tau }+u_{i+1/2}^n\frac{u_{i+1}^n-u_i^n}{3h}+5u_{i-1/2}^n\frac{u_i^n-u_{i-1}^n}{3h}+ \\ \frac{1}{3}\beta \frac{u_{i+2}^n+2u_{i+1}^n-12u_i^n+14u_{i-1}^n-5u_{i-2}^n}{h^3}=0 \end{gathered}$$

(7)

2.5. Research the Scheme Stability, Obtained on the Basis of a Linear Combination of the “Upwind Leapfrog” and “Standard Leapfrog” Schemes

We investigate the stability of the scheme (7) using the harmonic method [25,26]. Let $u_i^n={\phi }^n\cdot e^{jki}$, where $j=\sqrt{-1}$; we substitute this expression in the equality (7):

$$\begin{gathered} \frac{{\phi }^{n+1}\cdot e^{jki}-{\phi }^n\cdot e^{jki}}{\tau }+2\frac{{\phi }^n\cdot e^{jk\left(i-1\right)}-{\phi }^{n-1}\cdot e^{jk\left(i-1\right)}}{3\tau }+\frac{{\phi }^n\cdot e^{jki}-{\phi }^{n-1}\cdot e^{jki}}{3\tau }+ \\ +u\frac{{\phi }^n\cdot e^{jk\left(i+1\right)}+4{\phi }^n\cdot e^{jki}-5{\phi }^n\cdot e^{jk\left(i-1\right)}}{3h}+ \\ \frac{1}{3}\beta \frac{{\phi }^n\cdot e^{jk\left(i+2\right)}+2{\phi }^n\cdot e^{jk\left(i+1\right)}-12{\phi }^n\cdot e^{jki}+14{\phi }^n\cdot e^{jk\left(i-1\right)}-5{\phi }^n\cdot e^{jk\left(i-2\right)}}{h^3}=0 \end{gathered}$$

Divide the equality by ${\phi }^n\cdot e^{jki}$; thus, we have:

$$\begin{gathered} \frac{\phi -1}{\tau }+2\cdot \frac{e^{-jk}-{\phi }^{-1}\cdot e^{-jk}}{3\tau }+\frac{1-{\phi }^{-1}}{3\tau }+u\cdot \frac{e^{jk}+4-5\cdot e^{-jk}}{3h}+ \\ \beta \cdot \frac{e^{2jk}+2e^{jk}-12+14\cdot e^{-jk}-5e^{-2jk}}{3h^3}=0 \end{gathered}$$

multiply the equation by $\phi \cdot \tau$:

$${\phi }^2+\left(\frac{2}{3}\cdot e^{-jk}-\frac{2}{3}+u\cdot \tau \cdot \frac{e^{jk}+4-5\cdot e^{-jk}}{3h}+\beta \cdot \tau \cdot \frac{e^{2jk}+2e^{jk}-12+14\cdot e^{-jk}-5e^{-2jk}}{3h^3}\right)\phi -\frac{2e^{-jk}+1}{3}=0.$$

Let $a=u\tau /h$, $b=\beta \tau /h^3$, then

$${\phi }^2+\left(\frac{2}{3}\cdot e^{-jk}-\frac{2}{3}+a\cdot \frac{e^{jk}+4-5\cdot e^{-jk}}{3}+b\cdot \frac{e^{2jk}+2e^{jk}-12+14\cdot e^{-jk}-5e^{-2jk}}{3}\right)\phi -\frac{2e^{-jk}+1}{3}=0.$$

We use the Euler formula $e^{jk}=\cos k+j\sin k$:

$${\phi }^2+(\frac{2}{3}\cdot \left(\cos k-1\right)-j\frac{2}{3}\cdot \sin k+\left(a+2b\cdot \left(\cos k-1\right)\right)\frac{-4\cos k+6j\sin k+4}{3})\phi -\frac{2e^{-jk}+1}{3}=0.$$

We denote $c=a+2b\left(\cos k-1\right)$; then, the equation has the form:

$${\phi }^2+(\frac{2}{3}\cdot \left(\cos k-1\right)-j\frac{2}{3}\cdot \sin k+c\cdot \frac{-4\left(\cos k-1\right)+6j\sin k}{3})\phi -\frac{2e^{-jk}+1}{3}=0$$

$${\phi }^2+\left(\frac{2}{3}\cdot \left(\cos k-1\right)\left(1-2c\right)-j\cdot \sin k\left(\frac{2}{3}-2c\right)\right)\phi -\frac{2e^{-jk}+1}{3}=0$$

The solution of the quadratic equation with respect to $\phi$:

$$\begin{array}{l}{\phi }_{1,2}=-\left(\frac{1}{3}\cdot \left(\cos k-1\right)\left(1-2c\right)-j\cdot \sin k\left(\frac{1}{3}-c\right)\right)\pm \\ \sqrt{{\left(\frac{1}{3}\cdot \left(\cos k-1\right)\left(1-2c\right)-j\cdot \sin k\left(\frac{1}{3}-c\right)\right)}^2+\frac{2e^{-jk}+1}{3}}.\end{array}$$

We consider the case $c=0$:

$${\phi }_{1,2}=-\left(\frac{1}{3}\cdot \left(\cos k-1\right)-\frac{1}{3}j\cdot \sin k\right)\pm \sqrt{{\left(\frac{1}{3}\cdot \left(\cos k-1\right)-j\cdot \frac{1}{3}\sin k\right)}^2+\frac{2e^{-jk}+1}{3}}$$

or

$${\phi }_{1,2}=-\frac{1}{3}\cdot e^{-jk}+\frac{1}{3}\pm \left(\frac{1}{3}\cdot e^{-jk}+\frac{2}{3}\right)$$

As a result, we have ${\phi }_1=1$, ${\phi }_2=-\frac{2}{3}\cdot e^{-jk}-\frac{1}{3}$; therefore, ${\phi }_2$ is not a solution.

We consider the case $c=1$:

$$\begin{gathered} {\phi }_{1,2}=-\left(-\frac{1}{3}\cdot \left(\cos k-1\right)+\frac{2}{3}j\cdot \sin k\right)\pm \\ \sqrt{{\left(-\frac{1}{3}\cdot \left(\cos k-1\right)+\frac{2}{3}j\cdot \sin k\right)}^2+\frac{2e^{-jk}+1}{3}} \end{gathered}$$

or ${\phi }_{1,2}=\frac{1}{3}\cdot \left(\cos k-1\right)-\frac{2}{3}j\cdot \sin k\pm \left(\frac{2}{3}\cos k-\frac{1}{3}j\sin k+\frac{1}{3}\right)$.

As a result, we obtain ${\phi }_1=e^{-jk}$, i.e., ${\phi }_1=\cos k-j\sin k$, so ${\left|{\phi }_1\right|}^2={\cos }^{2}k+{\sin }^{2}k=1$, $\left|{\phi }_1\right|=1,$${\phi }_2=-\frac{1}{3}{e}^{-jk}-\frac{2}{3}$; therefore, ${\phi }_2$ is not a solution.

Denote $\psi \left(k,c\right)$ absolute values of the function ${\phi }_{1,2}\left(k,c\right)$:

$$\begin{gathered} \psi =|-\left(\frac{1}{3}\cdot \left(\cos k-1\right)\left(1-2c\right)-j\cdot \sin k\left(\frac{1}{3}-c\right)\right)\pm \\ \sqrt{{\left(\frac{1}{3}\cdot \left(\cos k-1\right)\left(1-2c\right)-j\cdot \sin k\left(\frac{1}{3}-c\right)\right)}^2+\frac{2e^{-jk}+1}{3}}|. \end{gathered}$$

The behavior of the function values $\psi \left(k,c\right)$ can be checked numerically. As a result, for the values $k\in \left[0,2\pi \right]$ and $c\in \left[0,1\right]$ we can make sure that the inequality $\psi \le 1$ is observed, i.e., $\left|\phi \right|\le 1$, indicating the stability of the difference scheme.

2.6. Research the Accuracy of the Modified Scheme

The function $u\left(x,t\right)$ was written as a finite trigonometric Fourier series [17] in the complex form:

$$u\left(x,t\right)={\displaystyle \sum _{m=-N}^N{C}_m\left(t\right)e^{j\omega mx}},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{u}^2\left(x,t\right)/2={\displaystyle \sum _{m=-N}^N{D}_m\left(t\right)e^{j\omega mx}}$$

(8)

where $\omega =\frac{\pi }{l}$; $m$ is the harmonic number; $C_m\left(t\right)=\frac{2}{l}{\displaystyle \underset{0}{\overset{l}{\int }}u\left(x,t\right)e^{-j\omega mx}dx}$ is the complex amplitude of the m-the harmonic; $j=\sqrt{-1}$. After substituting (8) into (1), we get:

$${\left({\displaystyle \sum _{m=-N}^N{C}_m\left(t\right)e^{j\omega mx}}\right)}_t^{\prime }+{\left({\displaystyle \sum _{m=-N}^N{D}_m\left(t\right)e^{j\omega mx}}\right)}_x^{\prime }+\beta {\left({\displaystyle \sum _{m=-N}^N{C}_m\left(t\right)e^{j\omega mx}}\right)}_{xxx}^{‴}=0$$

We can change the sequence of differentiation and summation operations in the partial sum of the series and calculate the derivative by the space:

$${\displaystyle \sum _{m=-N}^N{\left(C_m\left(t\right)\right)}_t^{\prime }{e}^{j\omega mx}}+j\omega m{\displaystyle \sum _{m=-N}^N{D}_m\left(t\right)e^{j\omega mx}}-j\beta {\omega }^3{m}^3{\displaystyle \sum _{m=-N}^N{C}_m\left(t\right)e^{j\omega mx}}=0$$

Considering that the functions $e^{j\omega mx}$ are linearly independent for different values of m, we find:

$${\left(C_m\left(t\right)\right)}_t^{\prime }=-j\omega m{D}_m\left(t\right)+j\beta {\omega }^3{m}^3{C}_m\left(t\right)$$

(9)

Let us consider the application of a modified scheme, which is a combination of the “Upwind Leapfrog” and “Standard Leapfrog” difference scheme (7) for solving the KdV model equation.

Considering (8) and $x_i=ih$, the Equation (7) has the form:

$$\begin{gathered} {\displaystyle \sum _{m=-N}^N\frac{C_m^{n+1}{e}^{j\omega mhi}-C_m^n{e}^{j\omega mhi}}{\tau }}+2{\displaystyle \sum _{m=-N}^N\frac{C_m^n{e}^{j\omega mh\left(i-1\right)}-C_m^{n-1}{e}^{j\omega mh\left(i-1\right)}}{3\tau }}+ \\ {\displaystyle \sum _{m=-N}^N\frac{C_m^n{e}^{j\omega mhi}-C_m^{n-1}{e}^{j\omega mhi}}{3\tau }}+{\displaystyle \sum _{m=-N}^N\frac{D_m^n{e}^{j\omega mh\left(i+1\right)}+4D_m^n{e}^{j\omega mhi}-5D_m^n{e}^{j\omega mh\left(i-1\right)}}{3h}}+ \\ \frac{1}{3}\beta {\displaystyle \sum _{m=-N}^N\frac{C_m^n{e}^{j\omega mh\left(i+2\right)}+2C_m^n{e}^{j\omega mh\left(i+1\right)}-12C_m^n{e}^{j\omega mhi}+14C_m^n{e}^{j\omega mh\left(i-1\right)}-5C_m^n{e}^{j\omega mh\left(i-2\right)}}{h^3}}=0. \end{gathered}$$

After the transformations, we have:

$$\begin{gathered} {\displaystyle \sum _{m=-N}^N\frac{C_m^{n+1}-C_m^n}{\tau }{e}^{j\omega mhi}}+2{\displaystyle \sum _{m=-N}^N\frac{C_m^n-C_m^{n-1}}{3\tau }{e}^{-j\omega mh}{e}^{j\omega mhi}}+ \\ {\displaystyle \sum _{m=-N}^N\frac{C_m^n-C_m^{n-1}}{3\tau }{e}^{j\omega mhi}}+{\displaystyle \sum _{m=-N}^N\frac{D_m^n{e}^{j\omega mh}+4D_m^n-5D_m^n{e}^{-j\omega mh}}{3h}{e}^{j\omega mhi}}+ \\ \frac{1}{3}\beta {\displaystyle \sum _{m=-N}^N\frac{C_m^n{e}^{2j\omega mh}+2C_m^n{e}^{j\omega mh}-12C_m^n+14C_m^n{e}^{-j\omega mh}-5C_m^n{e}^{-2j\omega mh}}{h^3}{e}^{j\omega mhi}}=0 \end{gathered}$$

Due to linear independence of $e^{j\omega mi}$, we can rewrite the last expression in the form:

$$\begin{gathered} \frac{C_m^{n+1}-C_m^n}{\tau }+2\frac{C_m^n-C_m^{n-1}}{3\tau }{e}^{-j\omega mh}+\frac{C_m^n-C_m^{n-1}}{3\tau }+D_m^n\frac{e^{j\omega mh}+4-5e^{-j\omega mh}}{3h}+ \\ \frac{1}{3}\beta C_m^n\frac{e^{2j\omega mh}+2e^{j\omega mh}-12+14e^{-j\omega mh}-5e^{-2j\omega mh}}{h^3}=0 \end{gathered}$$

(10)

At $\tau \to 0$ from (10), it is followed:

$$\begin{gathered} {\left(C_m\left(t\right)\right)}_t^{\prime }\left(\frac{4}{3}+2\frac{e^{-j\omega mh}}{3}\right)+D_m^n\frac{e^{j\omega mh}+4-5e^{-j\omega mh}}{3h}+ \\ \frac{1}{3}\beta C_m^n\frac{e^{2j\omega mh}+2e^{j\omega mh}-12+14e^{-j\omega mh}-5e^{-2j\omega mh}}{h^3}=0 \end{gathered}$$

or

$${\left(C_m\left(t\right)\right)}_t^{\prime }=-\frac{e^{j\omega mh}+4-5e^{-j\omega mh}}{2h\left(2+e^{-j\omega mh}\right)}{D}_m^n-\beta \frac{e^{2j\omega mh}+2e^{j\omega mh}-12+14e^{-j\omega mh}-5e^{-2j\omega mh}}{2h^3\left(2+e^{-j\omega mh}\right)}{C}_m^n.$$

(11)

Lemma 1.

At approximation of the problem (1)–(3) by the difference scheme (7), for each harmonic, the solutions to the problem of the propagation velocity of the wave and the dispersion term are less than the real values and differ, respectively, by magnitudes:

$${\alpha }_1=1-\frac{e^{j\omega mh}+4-5e^{-j\omega mh}}{2j\omega mh\left(2+e^{-j\omega mh}\right)},\to {\alpha }_2=1+\frac{e^{2j\omega mh}+2e^{j\omega mh}-12+14e^{-j\omega mh}-5e^{-2j\omega mh}}{2j{\omega }^3{m}^3{h}^3\left(2+e^{-j\omega mh}\right)}$$

Proof of Lemma 1.

At $\tau \to 0$ from (11) it follows that:

$$\begin{array}{l}{\left(C_m\left(t\right)\right)}^{\prime }{}_t=-j\omega m{D}_m^n\frac{e^{j\omega mh}+4-5e^{-j\omega mh}}{2j\omega mh\left(2+e^{-j\omega mh}\right)}-\\ -j{\omega }^3{m}^3\beta C_m^n\frac{e^{2j\omega mh}+2e^{j\omega mh}-12+14e^{-j\omega mh}-5e^{-2j\omega mh}}{2j{\omega }^3{m}^3{h}^3\left(2+e^{-j\omega mh}\right)}.\end{array}$$

Given equality (9), the solution using scheme (7) corresponds to the solution of the equation:

$$u_t^{\prime }=-{\left(u^2/2\right)}_x\left(1-{\alpha }_1\right)-\beta u_{xxx}^{‴}\left(1-{\alpha }_2\right)$$

$${\alpha }_1=1-\frac{e^{j\omega mh}+4-5e^{-j\omega mh}}{2j\omega mh\left(2+e^{-j\omega mh}\right)},\to {\alpha }_2=1+\frac{e^{2j\omega mh}+2e^{j\omega mh}-12+14e^{-j\omega mh}-5e^{-2j\omega mh}}{2j{\omega }^3{m}^3{h}^3\left(2+e^{-j\omega mh}\right)}$$

□

Lemma 2.

The modified difference scheme for the problem (1)–(3) has an approximation error equal to the$O\left({\tau }^2+h^3\right)$.

Proof of Lemma 2.

At $h\to 0$ from (10), it follows that:

$$\begin{gathered} \frac{C_m^{n+1}-C_m^n}{\tau }+2\frac{C_m^n-C_m^{n-1}}{3\tau }{e}^{-j\omega mh}+\frac{C_m^n-C_m^{n-1}}{3\tau }+D_m^n\frac{e^{j\omega mh}+4-5e^{-j\omega mh}}{3h}+ \\ +\frac{1}{3}\beta C_m^n\frac{e^{2j\omega mh}+2e^{j\omega mh}-12+14e^{-j\omega mh}-5e^{-2j\omega mh}}{h^3}=\frac{C_m^{n+1}-C_m^{n-1}}{\tau }+2j\omega m{D}_m^n- \\ 2j\beta {\omega }^3{m}^3{C}_m^n=2{\left(C_m\left(t\right)\right)}_t^{\prime }+2j\omega m{D}_m^n-2j\beta {\omega }^3{m}^3{C}_m^n+O\left({\tau }^2\right) \end{gathered}$$

We can investigate the order of the approximation error using a modified difference scheme of the convective term in the space. We made the replacement $j\omega mh=s$:

$$\begin{gathered} {\alpha }_1=1-\frac{e^s+4-5e^{-s}}{2s\left(2+e^{-s}\right)}=-\frac{s^3}{36}+O\left(s^4\right) \\ \mathrm{or} {\alpha }_1=j\frac{{\left(\omega mh\right)}^3}{36}+O\left(h^4\right) \end{gathered}$$

According to the obtained expression, the modified difference scheme approximates the convective term with the third order of accuracy in the space.

We can estimate the order of the approximation error of the dispersion term in the space:

$${\alpha }_2=1-\frac{e^{2s}+2e^s-12+14e^{-s}-5e^{-2s}}{2s^3\left(2+e^{-s}\right)}=-\frac{s^2}{12}+O\left(s^3\right)$$

The research on the order of approximation error by a modified difference scheme for solving problems described by the Korteweg–de Vries model equation showed that the difference scheme (9) has an approximation error equal to the $O\left({\tau }^2+h^3\right)$. □

3. Results and Discussion

A numerical implementation of a nonlinear dispersion model used to solve the KdV equation in the form of a software module is carried out. Comparison of numerical (based on schemes (4)–(7)) and analytical solutions of the KdV problem with the initial conditions $\delta v_0\left(x\right)=0.1\sin \left(\pi \left(x-10\right)/20\right)\left(h(10-x)-h(30-x)\right)$, where $h(x)$ is the Heaviside function, is shown in Figure 3.

Input data for the presented numerical experiments: average velocity $v=0.5$ m/s, time step $\tau =0.1$ s, $\beta =0.1$ m³/s, space step h = 1 m, the length of the time interval T is 100 s. A detailed idea of the behavior of the wave during its propagation, described by the KdV equation using modified schemes, can be obtained from the analysis of the wave profile (along the abscissa axis, the distance is postponed; along the ordinate axis, we see the oscillatory component of the velocity $\delta v$).

Figure 3 demonstrates the solution to the linearized KdV problem. The red line is the numerical solution of the linearized KdV problem; the blue line is its analytical (accurate) solution. The numerical implementation of the problem obtained on the basis of the schemes is: the “left corner” shown in Figure 3a; the “Standard Leapfrog” shown in Figure 3b; the “Upwind Leapfrog” shown in Figure 3c; and the proposed scheme of the form (7) (a linear combination of the schemes “Upwind Leapfrog” and “Standard Leapfrog”) shown in Figure 3d.

The results of numerical calculations of the linearized KdV problem based on a modified difference scheme and a linear combination of cabaret and cross schemes practically coincide with the analytical solution. The use of the “left corner” scheme is impractical for solving the problem because this scheme has a large error caused by dissipation, which is evident from the calculation results.

Table 1. The dependence of the error of solving the problem on time.

	50 s	100 s	200 s	500 s	1000 s
“left corner”	0.373	0.524	0.719	0.846	0.908
“Standard Leapfrog”	0.146	0.224	0.324	0.505	0.713
“Upwind Leapfrog”	0.103	0.146	0.197	0.283	0.391
a linear combination of the schemes “Upwind Leapfrog” and “Standard Leapfrog”	0.048	0.059	0.069	0.085	0.101

shows the dependence on the error of solving the problem normally $L_1$ on time. The relative error was estimated using the formula: ${\psi }^n={\displaystyle \sum _i{\psi }_i^n}/{\displaystyle \sum _i\left|q\left(x_i,t^n\right)\right|},$ ${\psi }_i^n=\left|q_i^n-q\left(x_i,t^n\right)\right|$, where $q\left(x_i,t^n\right)$—the accurate solution in the node $\left(x_i,t^n\right)$, $q_i^n$—numerical solution on the nth time layer. Calculation data for the table: time step $\tau =0.1$ s, $\beta =0.1$ m/s, space step h = 1 m.

The thickening of the computational grid leads to a quadratic increase in labor intensity, and is possible only with small values of the parameter $\beta$. From the restriction $\beta \le u{h}^2/4$, it follows that the grinding of the grid step $h$ by spatial coordinate leads to a narrowing of the range of input parameters set for the initial KdV problem, in which the difference schemes described above are stable; namely, the upper value of the parameter $\beta$. For large parameter values $\beta$, a stable solution can be obtained only on coarse grids; therefore, it is necessary to increase the accuracy of the difference scheme to reduce the error of the solution. We set the following input data for calculations when solving the linearized test problem KdV: time step $\tau =0.05$ s, $\beta =0.025$ m³/s, space step h = 0.5 m.

Table 2. Dependence of the error on solving the problem on time.

	50 s	100 s	200 s	500 s	1000 s
“Left corner”	0.229	0.368	0.541	0.772	0.934
“Standard Leapfrog”	0.069	0.1	0.144	0.237	0.34
“Upwind Leapfrog”	0.046	0.068	0.101	0.153	0.223
a linear combination of the schemes “Upwind Leapfrog” and “Standard Leapfrog”	0.022	0.027	0.034	0.046	0.057

shows the dependence of the error on solving the problem on time for a more detailed calculation grid.

Figure 4 demonstrates the solution to the nonlinear KdV problem. The numerical solution of the problem obtained on the basis of the schemes is: the “left corner” shown in Figure 4a; the “Standard Leapfrog” shown in Figure 4b; the “Upwind Leapfrog” shown in Figure 4c; and the modified scheme of the form (7) (a linear combination of schemes “Upwind Leapfrog” and “Standard Leapfrog”) shown in Figure 4d.

Comparing the data in Table 1 and Table 2 (for the initial and smaller grids), we conclude that when the grid thickens, there is no increase in the accuracy of the numerical solution of the problem. Therefore, it is advisable to use a more accurate scheme to solve the problem, which, according to the experiments, is a modified scheme of the form (7).

Figure 5 shows the solution to the KdV problem in the dynamics. The solution of the problem based on the proposed scheme (7) is given. Figure 5a shows the graphs of the solution, where initial condition is bell solitary wave solution $\delta v_0\left(x\right)=0.1\sin \left(\pi \left(x-10\right)/20\right)\left(h(10-x)-h(30-x)\right)$. Figure 5b shows the graphs of the solution, where the initial distribution is kink solitary wave solution $\delta v_1\left(x\right)=0.1\left(h(10-x)-h(30-x)\right)$.

According to the graphs shown in Figure 5, we can follow the changing shape of the curve to solve the KdV problem and its shift along the Ox axis at $t\in \left[0,100\right]$, with a 25 s time step. The developed difference schemes for the numerical implementation of the nonlinear and linearized KdV problem do not differ much from each other in their properties (accuracy and stability). The behavior of the numerical solution, in this case, depends more on the choice of the difference scheme. The executed numerical experiments based on the proposed scheme of the form (7) are in good agreement with the analytical solution of the linearized KdV problem. The obtained numerical solutions of the KdV problem are used in the study of the dynamics of changes in the wave profile when studying the propagation of wave processes over long distances. It is important to develop and use high-precision difference schemes in the study of these processes since error tends to accumulate. The proposed scheme allows us to obtain a numerical solution to the KdV problem in the nonlinear case, which is relevant in solving problems of wave hydrodynamics.

4. Conclusions

The paper proposes an approximation of the KdV equation of a modified difference scheme with improved dispersion properties. It is the result of a linear combination of the difference operators of the “Upwind Leapfrog” and “Standard Leapfrog” schemes, for which the weight coefficients were obtained as a result of minimizing the order of approximation error. A study of its behavior on uniform rectangular grids was carried out.

A comparison of the properties of the schemes “left corner”, “Standard Leapfrog”, and “Upwind Leapfrog” and the proposed scheme showed that the modified scheme has great advantages over the others. The disadvantage of the “left corner” scheme is a large dissipation. To eliminate large grid viscosity, it becomes necessary to reduce the step $h$ by the spatial variable; however, the thickening of the computational grid leads to a quadratic increase in labor intensity and more stringent restrictions on the parameters $\beta$. The stability and accuracy research results showed that the proposed difference scheme has the same limitations as the “left corner” scheme with an approximation error equal to the $O\left({\tau }^2+h^3\right)$. According to the results of numerical calculations, it follows that it is expedient to use a modified difference scheme to solve the KdV problem.