A hybrid numerical method for the KdV equation by finite difference and sinc collocation method

Applied Mathematics and Computation 355 (2019) 61–72

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

A hybrid numerical method for the KdV equation by finite difference and sinc collocation method Desong Kong a, Yufeng Xu b, Zhoushun Zheng a,∗ a b

Department of Information and Computing Science, Central South University, Changsha 410083, Hunan, PR China Department of Applied Mathematics, Central South University, Changsha 410083, Hunan, PR China

a r t i c l e

i n f o

MSC: 26A33 34A08 Keywords: KdV equation θ -weighted Finite difference method Sinc collocation

a b s t r a c t In this paper, we propose a hybrid numerical method for the KdV equation. More precisely, we discretize the temporal derivative of KdV equation by a θ -weighted scheme and treat the implicitly nonlinear term with the combination of finite difference and sinc collocation method. The stability analysis is presented, and numerical experiments illustrate the efficiency and stability of the proposed hybrid method. The motions and the interactions of solitary waves relying on particular initial boundary conditions are also simulated. © 2019 Elsevier Inc. All rights reserved.

1. Introduction The study of wave equations arising from shallow water dynamics has been a subject of increasing interests in various fields of mathematics and physics. One of the well-known models are named as Korteweg–de Vires (KdV) equation, which is firstly introduced by Boussinesq around 1877 and again derived by Korteweg and de Vires in around 1895 [1]. Nowadays, the development of KdV equation and its generalizations plays a vital role in the description of a wide range of physics phenomena, such as internal waves in coastal waters, waves in plasma physics, magma flow and conduit waves, flow in blood vessels. For more details, we refer to [2–4] the references therein. In mathematics, the classical KdV model has a concise form read as

ut + ε uux + μuxxx = 0, x ∈ R, t ≥ 0,

(1.1)

where the effect of the nonlinear terms uux causing steepening of the waveform is counterbalanced by the dispersion term uxxx representing dispersion effect. This shall make the waveform spread and lead to a solitary wave. The existence and uniqueness of solutions for KdV equation under certain classes of initial functions, smooth enough whilst tending to zeros sufficiently rapidly as x → ∞, have been examined by Lax [5] and Sjöberg [6] in several decades before. Note that there are many KdV equations with particular initial conditions cannot be solved analytically, there are also a large number of excellent monographs available for solving KdV equation numerically, e.g., see [7–9]. In this paper, we shall consider the numerical solution of model (1.1) defined on a bounded domain such that the soliton does not reach the boundary. Therefore, (1.1) can be rewritten as

ut + ε uux + μuxxx = 0, x ∈ (a, b), t ∈ (0, T ], ∗

Corresponding author. E-mail addresses: [email protected] (D. Kong), [email protected] (Y. Xu), [email protected] (Z. Zheng).

https://doi.org/10.1016/j.amc.2019.02.031 0 096-30 03/© 2019 Elsevier Inc. All rights reserved.

(1.2)

62

D. Kong, Y. Xu and Z. Zheng / Applied Mathematics and Computation 355 (2019) 61–72

equipped with initial and boundary conditions



u(x, 0 ) = u0 , x ∈ [a, b], u(a, t ) = u(b, t ) = ux (b, t ) = 0, t ∈ [0, T ].

(1.3)

Due to the easy-to-implement of the finite difference method (FDM) and a sparse system after the discretization in spatial dimension, numerical solution of (1.2) is simulated by many researchers with FDM in previous work such as [10–12]. For instance, in [10], a class of three-level finite difference scheme for solving KdV equation is proposed firstly. The proposed three-level scheme which requires a starting procedure satisfies momentum and energy conservation laws. In [11], a two-level hopscotch algorithm is derived and the stability and dispersion is analyzed. The method is analyzed with respect to stability and dispersion. It is also proven that such scheme is conservative and possesses a minimal phase error. In [12], some three-layer weighted finite-difference schemes for KdV equation are constructed based on conservation laws. In addition to the above earlier endeavor, some methods based on FDM are also developed during the last decades, such as analytical-numerical (ANS) method [13], exponential finite-difference method [14]. Similar with FDM producing a sparse system, Galerkin method using the basis functions vanishing outside a local interval is implemented for KdV equation in [15,16]. In [15], a kind of quadratic B-spline Galerkin finite element method is proposed. A penta-diagonal coefficient matrix is derived by this quadratic B-spline Galerkin method. In [16], a new basis of smooth piecewise cubic polynomials is developed for simulating the motions and the interactions of solitary waves. The coefficient matrix is seven-diagonal. Along with FDM, spectral methods are extensively studied for KdV equation, since it usually achieve high-order accuracy when the ansatz solution is smooth enough. The odd (third) derivative in KdV equation can be efficiently approximated by some Petrov–Galerkin methods [17–19]. What surprises us is that some authors study the sinc collocation–interpolation method for nonlinear initial-boundary value problem [20], which is later applied to some nonlinear wave equations [21] and eventual periodicity of KdV equation [22]. This encourages us to further study the approximation of nonlinear terms in KdV equation by combining finite differences and sinc collocation skill. We note that there are many other numerical methods for KdV equation, such as spline approximation [23] and heat balance integral (HBI) method [24], which will not be discussed for the conciseness of this work. Motivated by the efficient implementation of FDM, the spectral accuracy of sinc collocation method and the fact that both of them are easily defined on equidistant nodes, we shall construct a hybrid method to simulate the numerical solution using the finite difference and sinc collocation method in space for KdV equation. More precisely, we impose a θ -weighted scheme for temporal derivative and expand the implicit nonlinear term by linearization method (more details can be found in [22]). After that a fully discrete scheme for KdV equation is eventually established by combining finite differences and sinc collocation applied to approximate the first-order derivative emerging in nonlinear term and by finite difference scheme with second order accuracy for the third-order derivative term, respectively. According to the matrix analysis technique discussed in several references [22,25], we can further deduce the growth factor for the resulted linear algebraic equation, which directly shows the desired stability conclusion according to von Neumann theorem. Several numerical examples are carefully designed to validate the efficiency of proposed method and detailed comparisons with results from references are also carried out, which again illustrate the efficiency of our method. The rest of this paper is organized as follows. In Section 2, we recall some basic definitions and properties about sinc function for preparation and present numerical scheme for KdV equation. In Section 3, the stability of numerical scheme is derived based on matrix analysis. In Section 4, we present various numerical experiments which validate the accuracy and efficiency of the proposed algorithm. Finally, we conclude with a brief discussion in Section 5. 2. Numerical scheme In this section, some preliminary knowledge about sinc function and its first-order derivative matrix are demonstrated. The semidiscrete and fully discrete schemes of KdV equation (1.2) and (1.3) are also proposed. 2.1. Semidiscrete scheme of KdV equation Let us begin with introducing the semidiscrete scheme of KdV equation on temporal direction. Divide the successive time domain into M subintervals with equivalent time step τ , which yields a monotonically increasing time points, 0 = t0 < t1 < · · · < tM−1 < tM = T . For convenience, denote un as the numerical value at time level t = tn . The semidiscrete scheme on time can be derived by approximating the first-order derivative with θ -weighted scheme

un+1 − un

τ

+θ

 n+1 n+1  +1 ε u ux + μunxxx + (1 − θ )(ε un unx + μunxxx ) = 0, x ∈ (a, b),

(2.1)

where θ ∈ [0, 1]. When θ = 0, the semidiscrete scheme degenerates to explicit scheme, and when θ = 1, the semidiscrete scheme degenerates to implicit scheme. Whereas, when θ = 12 , the semidiscrete scheme can be regarded as Crank–Nicholson scheme. It is obvious that the (2.1) cannot be solved directly due to the nonlinear term emerging in the implicit time layer, un+1 unx +1 . Inspired by the approach introduced in [22], we deal with the nonlinear term by expanding it on the explicit time layer with Taylor expansion. Namely,

un+1 = u(tn + τ ) = u(tn ) + τ ut (tn ) + O(τ 2 ),

D. Kong, Y. Xu and Z. Zheng / Applied Mathematics and Computation 355 (2019) 61–72

63

and

unx +1 = ux (tn + τ ) = ux (tn ) + τ uxt (tn ) + O(τ 2 ). Then, we get







un+1 unx +1 = u(tn ) + τ ut (tn ) + O(τ 2 ) × ux (tn ) + τ uxt (tn ) + O(τ 2 ) = un unx + τ ux (tn )

 un+1 − un  τ

+ τ u(tn )

 un+1 − un  x

τ

x



+ O ( τ 2 ).

Dropped the high order small quantities with respect to τ 2 , the nonlinear term un+1 unx +1 can be approximated as

un+1 unx +1 ≈ un+1 unx + un unx +1 − un unx .

(2.2)

Substituting the above approximation about un+1 unx +1 into (2.1) yields the the semidiscrete scheme of KdV equation (1.2) as follows:

un+1 + θ τ

 n+1 n 

+1 ε u ux + un unx +1 + μunxxx = un + (2θ − 1 )τ ε un unx − (1 − θ )τ μunxxx .

(2.3)

2.2. Fully discrete scheme of KdV equation In this subsection, we consider the fully discrete scheme of (1.2) with initial and boundary conditions (1.3), the main idea is approximating the first-order derivatives about spatial variable x in (2.2) with the combination of finite difference method and sinc collocation method. To construct the numerical approximation scheme, the spatial interval, [a, b], is discretized by n a set of uniformly distant points {x j }Nj=0 , which satisfy a = x0 < x1 < · · · < xN−1 < xN = b, and h = x j+1 − x j = b−a N . Denote u j as the numerical value of u(x, t) at point (xj , tn ). We approximate the first-order derivative of u(x, t) with respect to x with the help of sinc function which defined as [26–28]

Sh ( x ) =

sin(π x/h ) . π x/h

(2.4)

It can be verified that

Sh (0 ) = 1, Sh ( jh ) = 0, for j = ±1, ±2, . . . . Then the band-limited interpolant of u(x, t) with respect to x can be presented as

C (u, h )(x ) :=

∞ 

u( jh )Sh (x − jh ).

j=−∞

To approximate the first-order derivative of u(x, t) with respect to x over the discrete points, the first-order derivative of sinc function can be calculated as

⎧ i = j, ⎨0, d Sh (xi − jh ) = [Sh (x − jh )]|xi = (−1 )i− j dx ⎩ , i = j. ( i − j )h

(2.5)

For simplicity, let S = (Si j ) = (Sh (xi − jh )). Besides, the first-order derivative of u(x, t) with respect to x can also be approximated by the central difference formula at points xj ,

ux ( x j ) ≈

u j+1 − u j−1 , j = 1, 2, . . . , N − 1. 2h

(2.6)

And the approximation of the third-order derivative is also standard by finite difference method with second order accuracy. It should be noticed that for the approximated values near boundary conditions, we need extra points to approximate them. In these cases, we can interpret the right Neumann boundary condition at the right boundary by 0 = ux (b, t ) ≈ (uN+1 − uN−1 )/(2h ). For the left boundary condition, we set u(x, t ) = 0 when x < a, then the following scheme can be yielded,

uxxx (x1 ) ≈

u3 − 2u2 + 2u0 . 2 h3

(2.7)

Thus, for third-order derivative value at points xj , we have

uxxx (x j ) ≈

u j+2 − 2u j+1 + 2u j−1 − u j−2 , j = 2, . . . , N − 2, 2 h3

(2.8)

and

uxxx (xN−1 ) ≈

−2uN + uN−1 + 2uN−2 − uN−3 . 2 h3

(2.9)

64

D. Kong, Y. Xu and Z. Zheng / Applied Mathematics and Computation 355 (2019) 61–72

Combining (2.5), (2.6), (2.8) and (2.9) for all the discrete points on spatial direction, the approximations in spatial direction about the first and third order derivatives are yielded, respectively, with the matrix form as follows

⎛

0 ⎜−1 ⎜ 1 ⎜0 ux (x, tn ) = ⎜ 2h ⎜

0 0 −1

0 1 0 .. .

⎝

⎞

0 0 1 .. . −1 0

⎟ ⎟ ⎟ ⎟ ⎟ ⎠

..

. 0 0

un := D(1) un ,

(2.10)

1 0 (N+1)×(N+1)

or

ux (x, tn ) = Sun , and

(2.11)

⎛

uxxx (x, tn ) =

0 ⎜2 ⎜−1 ⎜

1 ⎜ ⎜ 2 h3 ⎜

⎜ ⎝



0 0 2 .. .

0 −2 0 .. . −1 0 0

0 1 −2 .. . 2 −1 0

⎞

0 0 1 .. . 0 2 0

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ 1⎟ ⎠ −2

..

. −2 1 0

0

un := D(3) un ,

(2.12)

(N+1 )×(N+1 )

 n

where un = un0 , un1 , . . . , unN−1 , uN . Substituting (2.5), (2.6), (2.8) and (2.9) into (2.3), we get the fully discrete scheme of problem (1.2) and (1.3) as

 

unj +1

+ θτ

ε

unj +1

N−1 

Si j uni

+

unj

+1 +1 unj+1 − unj−1

i=1

= unj + (2θ − 1 )τ ε unj

N−1 

2h



+μ

Si j uni − (1 − θ )τ μ

+1 +1 +1 +1 unj+2 − 2unj+1 + 2unj−1 − unj−2

2 h3

unj+2 − 2unj+1 + 2unj−1 − unj−2

i=1



2 h3

,

(2.13)

j = 1, 2, . . . , N − 2, n = 0, 1, . . . , M − 1, The systems of linearized equations (2.13) can be rewritten in compact form

Aun+1 = bun + Fn+1 , where

A = I + θτ

(2.14)

 

ε diag(Sˆun ) + diag(un )D(1) + μD(3) ,

b = Iˆ + (2θ − 1 )τ ε diag(Sˆun ) − (1 − θ )τ μD(3) , and

Fn+1 = (u(a, tn+1 ), 0, . . . , 0, u(b, tn+1 ) ) . Remark 2.1. In the matrix formula above,



Sˆ = and

Si j , i = 2, . . . , N; j = 1, . . . , N + 1, 0,

i = 1, N + 1; j = 1, . . . , N + 1,

 Iˆ =

1, i = j = 2, . . . , N, 0, others.

Besides, I denotes the identity matrix and diag(Sˆun ) means the (N + 1 ) × (N + 1 ) diagonal matrix with elements Sˆun . In addition, diag(un ) denotes a diagonal matrix where the diagonal entries are elements of un . Meanwhile, the coefficient matrix A is penta-diagonal sparse matrix which can be solved efficiently by any iterative algorithms.

D. Kong, Y. Xu and Z. Zheng / Applied Mathematics and Computation 355 (2019) 61–72

65

3. Stability analysis In this section, we investigate the stability of the numerical scheme obtained from Section 2. Let U nj denote the exact



n , Un value of u(x, t) at point (xj , tn ), and Un = U0n , U1n , . . . , UN−1 N



. Then the evolution of error can be written as

Aen+1 = ben , where ek gives the absolute at t = tk ,





ek = Uk − uk . Now, we consider the spectral radius of C, denoted by ρ (C), where

C :=A−1 b

    −1 ε diag(Sˆun ) + diag(un )D(1) + μD(3)

× Iˆ + (2θ − 1 )τ ε diag(Sˆun ) − (1 − θ )τ μD(3) .



= I + θτ

Denote λ0 , λ3 , λN1 and λN2 as the eigenvalues of I, D(3 ) , diag(Sˆun ) and diag(un )D(1 ) , respectively. It is obvious that diag(Sˆun ) is a diagonal matrix with real numbers. Thus, we yield that λ0 = 1, λN1 is one of the elements of Sˆun , and λN2 , λ3 are complex, generally. The necessary condition of the stability of the proposed numerical scheme is ρ (C ) ≤ 1 + K τ for a positive number K [25]. Denote

λRN1 := τ ελN1 , λRN2 := τ ε Re(λN2 ), λIN2 := τ ε Im(λN2 ), and

λR3 := τ μRe(λ3 ), λI3 := τ μIm(λ3 ), and λIR 3 := τ μλ3 . So for (2.14), a stable scheme must satisfy

   1 + (2θ − 1 )λRN1 + (θ − 1 )λR3 + i(θ − 1 )λI3     1 + θ (λR + λR + λR ) + iθ (λI + λI )  ≤ 1, N1 N2 3 N2 3

which implies that

    1 + (2θ − 1 )λRN1 + (θ − 1 )λR3 + i(θ − 1 )λI3  ≤ 1 + θ (λRN1 + λRN2 + λR3 ) + iθ (λIN2 + λI3 ),

i.e.,



1 + (2θ − 1 )λRN1 + (θ − 1 )λR3



≤ 1 + θ (λRN1 + λRN2 + λR3 )

2

2



+ (θ − 1 )λI3

+



2

2 θ (λIN2 + λI3 ) .

(3.1)

Direct simplification of (3.1) yields

2(θ − 1 )λRN1 − 2θ λRN2 − 2λR3 + 2(θ 2 − 3θ + 1 )λRN1 λR3 − 2θ 2 (λRN1 λRN2 + λRN2 λR3 + λIN2 λI3 ) 2 ≤ (3θ − 1 )(1 − θ )|λRN1 |2 + θ 2 |λN2 |2 + (2θ − 1 )|λIR 3 | ,

(3.2)

which must be held for all the eigenvalues of corresponding matrices. From the fact that if ≤ θ ≤ 1, the right-hand side of (3.2) is nonnegative obviously, whereas the left-hand side may run negative depending on the choice of eigenvalues, we can conclude that the condition 12 ≤ θ ≤ 1 is necessary, but may not sufficient. By numerical calculating, the result shown in Fig. 3.1 illustrates that 12 ≤ θ ≤ 1 is nearly sufficient which satisfies ρ (C ) − 1 < 12 × 10−3 for θ > 12 . For θ = 0, this numerical scheme is conditionally stable with the restriction on time step 1 2

τ≤

2ελN1 + 2μRe(λ3 )

(ελN1 + μRe(λ3 ) )2 + (μIm(λ3 ) )2

.

4. Numerical experiments In this section, several numerical examples are carefully designed to validate the efficiency of the proposed method. In the physical opinion, the motion and interactions of solitary waves will be considered. In addition, some conservation laws that KdV equation satisfies will be examined by numerically calculating the following three invariants corresponding to conservation of mass, momentum and energy, as defined in [16],

 E1 =

b a

u(x )n dx ≈ h

N  j=0

unj ,

66

D. Kong, Y. Xu and Z. Zheng / Applied Mathematics and Computation 355 (2019) 61–72

Fig. 3.1. Numerical demonstration of convergence dependence on θ . Table 4.1 Invariants and error of Example 4.1 at different times with h = τ = 0.001 and θ = 12 .

 E2 =

b a

 E3 =

(u(x )n )2 dx ≈ h

E1

E2

E3

En

0 0.005 0.01

0.144598 0.144598 0.144598

0.086759 0.086759 0.086759

0.046852 0.046852 0.046852

5.60e−5 2.64e−5

(unj )2 ,

j=0



b a

where

( ) = ux nj

N 

Time

( u ( x )n )3 −

 un − un j+1 j−1 2h 0,

  N   3μ (unx (x ))2 dx ≈ h (unj )3 − ((ux )nj )2 , ε ε j=0

3μ

, 1 ≤ j ≤ N − 1, j = 0, N.

To demonstrate the efficiency and accuracy of the presented method for KdV equation, we use the maximum error norm given by





E n = max unj − U jn . 0≤ j≤N

4.1. Numerical simulation for KdV equation with small μ In this subsection, we consider KdV equation with small parameter μ, and choose μ = 4.84 × 10−4 , ε = 1. The motions of single and double solitary waves are simulated in Example 4.1 and Example 4.2, respectively. The numerical solutions obtained from our present method are compared with these values from some other methods reported in [13,15,24]. In addition, we will investigate the accuracy of the present method. Example 4.1. Referring to the numerical experiment in [15,16], the KdV equation has an analytical solution given by

u(x, t ) = 3c sech (κ x − ωt − x0 ), 0 ≤ x ≤ 2, 2

where

κ=

1 2



(4.1)

εc , and ω = ε cκ . μ

In numerical simulation, the initial condition is obtained from (4.1) by taking t = 0, c = 0.3, and x0 = 6. In this case, the analytical values of the three invariants are C1 = 0.144598, C2 = 0.086759, and C3 = 0.046852, respectively. In Table 4.1, we give the numerical values of the three invariants and the maximum error. It can be seen that the present method preserves the invariants well with small maximum error. In Table 4.2, we investigate the convergence rate of the present method in

D. Kong, Y. Xu and Z. Zheng / Applied Mathematics and Computation 355 (2019) 61–72

67

Table 4.2 Error and convergence order of presented method for Example 4.1 at T = 0.1. h (τ = 10−4 )

θ = 1/2

Order

θ =1

Order

τ (h = 1/500 )

θ = 1/2

Order

θ =1

Order

1/50 1/60 1/70 1/80 1/90 1/100

2.15e−2 1.39e−2 1.03e−2 7.62e−3 6.10e−3 4.84e−3

2.40 1.94 2.26 1.90 2.18

1.12e−2 7.59e−3 5.49e−3 4.21e−3 3.38e−3 2.74e−3

2.12 2.10 1.99 1.86 2.01

1/20 1/30 1/40 1/50 1/60 1/70

4.79e−3 2.17e−3 1.24e−3 8.21e−4 5.90e−4 4.48e−4

1.96 1.93 1.86 1.81 1.79

4.18e−2 3.03e−2 2.37e−2 1.94e−2 1.64e−2 1.42e−2

0.80 0.86 0.89 0.92 0.93

Fig. 4.1. Numerical solution and distribution of absolute error of Example 4.1 with h = τ = 10−3 and θ = 1/2. Table 4.3 Comparison of numerical and exact solution of Example 4.1 at T = 0.005 with h = 0.0125, τ = 0.001 and θ = 1. x

Numerical solution

Exact solution

Present

B-spline [15]

ANS [13]

HBI [24]

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4  2.0

0.0 0 0 0 0 0 0 0 0.0 0 025529 0.0 03090 04 0.03657255 0.35619410 0.86353642 0.17776094 0.01628422 0.00136224 0.0 0 011306 0.0 0 0 0 0938 0.0 0 0 0 0 078 0.0 0 0 0 0 0 06 0.0 0 0 0 0 0 01 0.0 0 0 0 0 0 0 0  0.0 0 0 0 0 0 0 0

0.0 0 0 0 0 0 0 0 0.0 0 025692 0.00309265 0.03658902 0.35578320 0.86299040 0.17789270 0.01627284 0.00136096 0.0 0 011295 0.0 0 0 0 0937 0.0 0 0 0 0 078 0.0 0 0 0 0 0 06 0.0 0 0 0 0 0 01 0.0 0 0 0 0 0 0 0  0.0 0 0 0 0 0 0 0

0.0 0 0 0 0 0 0 0 0.0 0 025686 0.00309148 0.03656208 0.35613919 0.86253589 0.17792331 0.01626755 0.00136019 0.0 0 011289 0.0 0 0 0 0936 0.0 0 0 0 0 078 0.0 0 0 0 0 0 06 0.0 0 0 0 0 0 01 0.0 0 0 0 0 0 0 0  0.0 0 0 0 0 0 0 0

0.0 0 0 0 0 0 0 0 0.0 0 026665 0.00320978 0.03794787 0.36617911 0.85626808 0.17201141 0.01568003 0.00131099 0.0 0 010880 0.0 0 0 0 0902 0.0 0 0 0 0 075 0.0 0 0 0 0 0 06 0.0 0 0 0 0 0 01 0.0 0 0 0 0 0 0 0  0.0 0 0 0 0 0 0 0

0.0 0 0 0 0 0 0 0 0.0 0 025688 0.00309232 0.03658531 0.35574620 0.86305610 0.17787410 0.01627119 0.00136083 0.0 0 011294 0.0 0 0 0 0937 0.0 0 0 0 0 078 0.0 0 0 0 0 0 06 0.0 0 0 0 0 0 01 0.0 0 0 0 0 0 0 0  0.0 0 0 0 0 0 0 0

maximum error norm which shows that convergence order on spatial direction is about two independent of the choice of θ , meanwhile, the convergence order on temporal direction is around one for θ = 1 and two for θ = 12 , respectively. To study the stability of the present method for long time behavior of a single solitary wave, the numerical solution and distribution of absolute error are shown in Fig. 4.1. From the results shown in Fig. 4.1, we know that the present method can simulate the exact solution well and the maximum error of the present method just occurs around the peak position of wave amplitude. In Tables 4.3 and 4.4, we show the numerical values by the present method, Galerkin B-spline finite element method, ANS, HBI and the exact values at the the selected notes for T = 0.005 and T = 0.01, respectively. On the other hand, Fig. 4.2

68

D. Kong, Y. Xu and Z. Zheng / Applied Mathematics and Computation 355 (2019) 61–72 Table 4.4 Comparison of numerical and exact solution of Example 4.1 at T = 0.01 with h = 0.0125, τ = 0.001 and θ = 1. x

Numerical solution

Exact solution

Present

B-spline [15]

ANS [13]

HBI [24]

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4  2.0

0.0 0 0 0 0 0 0 0 0.0 0 024031 0.00297413 0.03524336 0.34647976 0.87021124 0.18367914 0.01691035 0.00141549 0.0 0 011749 0.0 0 0 0 0975 0.0 0 0 0 0 081 0.0 0 0 0 0 0 07 0.0 0 0 0 0 0 01 0.0 0 0 0 0 0 0 0  0.0 0 0 0 0 0 0 0

0.0 0 0 0 0 0 0 0 0.0 0 024750 0.00297987 0.03527915 0.34553080 0.86920580 0.18394800 0.01688850 0.00141291 0.0 0 011727 0.0 0 0 0 0972 0.0 0 0 0 0 080 0.0 0 0 0 0 0 07 0.0 0 0 0 0 0 01 0.0 0 0 0 0 0 0 0  0.0 0 0 0 0 0 0 0

0.0 0 0 0 0 0 0 0 0.0 0 024746 0.00297741 0.03522687 0.34632607 0.86833708 0.18396796 0.01687704 0.00141118 0.0 0 011712 0.0 0 0 0 0971 0.0 0 0 0 0 081 0.0 0 0 0 0 0 07 0.0 0 0 0 0 0 01 0.0 0 0 0 0 0 0 0  0.0 0 0 0 0 0 0 0

0.0 0 0 0 0 0 0 0 0.0 0 026665 0.00320979 0.03794801 0.36618044 0.85627119 0.17201203 0.01568009 0.00131099 0.0 0 010881 0.0 0 0 0 0902 0.0 0 0 0 0 075 0.0 0 0 0 0 0 06 0.0 0 0 0 0 0 01 0.0 0 0 0 0 0 0 0  0.0 0 0 0 0 0 0 0

0.0 0 0 0 0 0 0 0 0.0 0 024746 0.00297915 0.03527078 0.34551640 0.86931950 0.18391210 0.01688449 0.00141257 0.0 0 011724 0.0 0 0 0 0972 0.0 0 0 0 0 081 0.0 0 0 0 0 0 07 0.0 0 0 0 0 0 01 0.0 0 0 0 0 0 0 0  0.0 0 0 0 0 0 0 0

Fig. 4.2. Absolute error distribution of Example 4.1.

shows a semi-logarithmic axis figure about the absolute errors between the numerical solution and exact solution from which we can derive that the present method have higher accuracy than HBI method and as well as the ANS method. Example 4.2. Next, we consider the interaction of two solitary waves where the initial condition is given by the analytical solution

u(x, t ) = 12μ(log F )xx , 0 ≤ x ≤ 4, where

F = 1 + eη1 + eη2 +

(4.2)

 α − α 2 1 2 eη1 +η2 , α1 + α2

ηi = αi x − αi3 μt + bi , i = 1, 2,   0.3 0.1 α1 = , α2 = , b1 = −0.48α1 , b2 = −1.07α2 . μ μ We calculate the three invariants where the analytic values are C1 = 0.099055, C2 = 0.019513 and C3 = 0.002380, the numerical values of three invariants and the maximum error on T = 0.005 and T = 0.01 are shown in Table 4.5. Fig. 4.3 shows the numerical solution and distribution of absolute error at time t = 0.01 with h = 0.0125, τ = 0.001 and θ = 12 . In addition, we simulate the numerical solution shown in Fig. 4.4 for T = 0.01. All of the results can illustrate that the numerical solution highly agrees with the physical phenomena. In Tables 4.6 and 4.7, we show the numerical values and exact ones at different points for h = 0.0125, τ = 0.001 and θ = 12 at T = 0.005 and T = 0.01, respectively. Beside, a semi-logarithmic axis figure shows the absolute error between the

D. Kong, Y. Xu and Z. Zheng / Applied Mathematics and Computation 355 (2019) 61–72 Table 4.5 Invariants and error of Example 4.2 with h = 0.0125, τ = 0.001 and θ = 1/2. Time

E1

E2

E3

En

0 0.005 0.01

0.099055 0.099055 0.099055

0.019513 0.019513 0.019513

0.002380 0.002380 0.002379

6.90e−3 1.41e−2

Fig. 4.3. Numerical solution and distribution of absolute error of Example 4.2.

Fig. 4.4. Numerical solution of Example 4.2. Table 4.6 Comparison of numerical and exact solution of Example 4.2 at T = 0.005 with h = 0.0125, τ = 0.001 and θ = 1/2. x

Numerical solution

Exact solution

Present

B-spline [15]

ANS [13]

HBI [24]

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8  4.0

0.0 0 0 0 0 0 0 0 0.00141224 0.16405407 0.07305194 0.00132886 0.01291987 0.11306093 0.05037215 0.00354271 0.0 0 020258 0.0 0 0 01144 0.0 0 0 0 0 065 0.0 0 0 0 0 0 04 0.0 0 0 0 0 0 0 0 0.0 0 0 0 0 0 0 0  0.0 0 0 0 0 0 0 0

0.0 0 0 0 0 0 0 0 0.00141361 0.16388340 0.07309954 0.00132846 0.01291965 0.11307170 0.05037538 0.00354259 0.0 0 020257 0.0 0 0 01144 0.0 0 0 0 0 065 0.0 0 0 0 0 0 04 0.0 0 0 0 0 0 0 0 0.0 0 0 0 0 0 0 0  0.0 0 0 0 0 0 0 0

0.0 0 0 0 0 0 0 0 0.00141571 0.16502401 0.07308242 0.00132717 0.01291787 0.11306786 0.05036473 0.00354275 0.0 0 020256 0.0 0 0 01144 0.0 0 0 0 0 065 0.0 0 0 0 0 0 04 0.0 0 0 0 0 0 0 0 0.0 0 0 0 0 0 0 0  0.0 0 0 0 0 0 0 0

0.0 0 0 0 0 0 0 0 0.00146231 0.16469787 0.07126321 0.00130927 0.01295362 0.11258566 0.05012470 0.00350801 0.0 0 020 051 0.0 0 0 01132 0.0 0 0 0 0 064 0.0 0 0 0 0 0 04 0.0 0 0 0 0 0 01 0.0 0 0 0 0 0 0 0  0.0 0 0 0 0 0 0 0

0.0 0 0 0 0 0 0 0 0.00141318 0.16055730 0.07392883 0.00132844 0.01290511 0.11263310 0.05056182 0.00354377 0.0 0 020258 0.0 0 0 01144 0.0 0 0 0 0 065 0.0 0 0 0 0 0 04 0.0 0 0 0 0 0 0 0 0.0 0 0 0 0 0 0 0  0.0 0 0 0 0 0 0 0

69

70

D. Kong, Y. Xu and Z. Zheng / Applied Mathematics and Computation 355 (2019) 61–72 Table 4.7 Comparison of numerical and exact solution of Example 4.2 at T = 0.01 with h = 0.0125, τ = 0.001 and θ = 1/2. x

Numerical solution

Exact solution

Present

B-spline [15]

ANS [13]

HBI [24]

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8  4.0

0.0 0 0 0 0 0 0 0 0.00135957 0.16302807 0.07456142 0.00134528 0.01284622 0.11319380 0.05046562 0.00356678 0.0 0 020405 0.0 0 0 01152 0.0 0 0 0 0 065 0.0 0 0 0 0 0 04 0.0 0 0 0 0 0 0 0 0.0 0 0 0 0 0 0 0  0.0 0 0 0 0 0 0 0

0.0 0 0 0 0 0 0 0 0.00136239 0.16276950 0.07465229 0.00134446 0.01284584 0.11321610 0.05047162 0.00356654 0.0 0 020403 0.0 0 0 01152 0.0 0 0 0 0 065 0.0 0 0 0 0 0 04 0.0 0 0 0 0 0 0 0 0.0 0 0 0 0 0 0 0  0.0 0 0 0 0 0 0 0

0.0 0 0 0 0 0 0 0 0.00136702 0.16495009 0.07467392 0.00134190 0.01284207 0.11320585 0.05045165 0.00356687 0.0 0 020401 0.0 0 0 01152 0.0 0 0 0 0 065 0.0 0 0 0 0 0 04 0.0 0 0 0 0 0 0 0 0.0 0 0 0 0 0 0 0  0.0 0 0 0 0 0 0 0

0.0 0 0 0 0 0 0 0 0.00146231 0.16469794 0.07126324 0.00130927 0.01295362 0.11258571 0.05012472 0.00350801 0.0 0 020 051 0.0 0 0 01132 0.0 0 0 0 0 064 0.0 0 0 0 0 0 04 0.0 0 0 0 0 0 01 0.0 0 0 0 0 0 0 0  0.0 0 0 0 0 0 0 0

0.0 0 0 0 0 0 0 0 0.00136152 0.15599870 0.07645022 0.00134442 0.01281754 0.11233430 0.05084663 0.00356898 0.0 0 020404 0.0 0 0 01152 0.0 0 0 0 0 065 0.0 0 0 0 0 0 04 0.0 0 0 0 0 0 0 0 0.0 0 0 0 0 0 0 0  0.0 0 0 0 0 0 0 0

Fig. 4.5. Absolute error distribution of Example 4.2.

numerical solutions and exact ones in Fig. 4.5. One can see from these tables and figures that the numerical solutions obtained by the present method are more accurate than the ones shown in [13,24] and are as well as the ones in [15]. 4.2. Numerical solution for KdV equation for long time evolution In this subsection, we will investigate the efficiency of the present method for long time evolution of KdV equation. The motions and interactions of solitary waves with initial condition given by [17,19] are simulated. The parameters in the following numerical experiment is chosen as ε = μ = 1. Example 4.3. Let us consider the motions and interactions of KdV equation with three and four solitons. Their initial conditions are given by

u(x, t ) =

3 

12κi2 sech (κi (x − 4κi2 t − xi )), −100 ≤ x ≤ 100, t = 0, 2

(4.3)

i=1

where

κ1 = 0.3, κ2 = 0.25, κ3 = 0.2, x1 = −60, x2 = −44, x3 = −26. for three solitary waves and

u(x, t ) =

4  i=1

12κi2 sech (κi (x − 4κi2 t − xi )), −150 ≤ x ≤ 150, t = 0, 2

(4.4)

D. Kong, Y. Xu and Z. Zheng / Applied Mathematics and Computation 355 (2019) 61–72

71

Table 4.8 Invariants of three solitary waves (left) with h = τ = 0.1 and θ = 1 and of four solitary waves (right) with h = 0.2, τ = 0.1 and θ = 1. Time

E1

E2

E3

Time

E1

E2

E3

0 56 112 168 224 280

18.0 0 0 0 18.0018 17.9974 17.9971 17.9985 17.9995

9.8274 9.5936 9.5138 9.3228 9.0697 8.8327

5.2636 5.0328 4.9651 4.7808 4.5362 4.3141

0 80 160 240 320 400

21.60 0 0 21.6028 21.6049 21.6011 21.6007 21.6074

10.3887 9.9723 9.7448 9.7023 9.4774 9.1922

5.2740 4.8594 4.6426 4.6035 4.3943 4.1368

Fig. 4.6. Numerical simulation for the interaction of three (left) and four (right) solitary waves.

where

κ1 = 0.3, κ2 = 0.25, κ3 = 0.2, κ4 = 0.15, x1 = −85, x2 = −60, x3 = −35, x3 = −10. for four solitary waves, respectively. By taking T = 280, h = τ = 0.1 and θ = 1, the numerical values of the three invariants of three solitary waves are calculated shown in Table 4.8 (left) and the three invariants of four solitary waves are calculated with T = 400, h = 0.2, τ = 0.1 and θ = 1 shown in Table 4.8 (right). The results illustrate that the present method is efficient to preserve the three invariants with exact values. The interactions of three and four solitary waves at different time and its global view of their interactions are shown in Fig. 4.6 (the left for three solitary waves and the right for four solitary waves). The numerical

72

D. Kong, Y. Xu and Z. Zheng / Applied Mathematics and Computation 355 (2019) 61–72

results keep the agreement with the physical phenomena that the solitary waves are stable against mutual collision and retain their identities, the soliton with higher amplitude travels with faster speed and so on. 5. Conclusion The numerical scheme for the KdV equation was presented with a finite difference and sinc collocation method in this paper. A detail stability analysis was provided for the fully discrete scheme of KdV equation. The result illustrated that the stability of this method for nonlinear problem was associated with the initial condition and the selection of the values of coefficients when the parameter θ took different values. The performance of our present method was evaluated in terms of the maximum error norm and the three invariants values of KdV equation. Numerical experiments illustrated the convergence rates and the stability of the present method. The results by comparing with the proposed results illustrated the present method can also simulate the exact solution well. Through the numerical examples that the motions and interactions of solitary waves were simulated in the cases of small parameter and large time domain, respectively, we conclude that the numerical results keep the agreement with the physical phenomena with the proposed method. Acknowledgments This work was supported by the National Key Research and Development Program of China (Nos. 2017YFB0305601 and 2017YFB0701700). The second author is partly supported by NSFC (No. 11501581). References [1] D.J. Korteweg, G. de Vries, On the change of form of long waves advancing in a rectangular canal and on a new type of long stationary waves, Philos. Mag. J. Sci. 39 (240) (1895) 422–443. [2] R.C. Cascaval, Variable coefficient KdV equations and waves in elastic tubes, Lect. Notes Pure Appl. Math. 234 (2003) 57–70. [3] D.G. Crighton, Applications of KdV, Acta Appl. Math. 39 (1995) 39–67. [4] A. Jeffrey, T. Kakutani, Weak nonlinear dispersive waves: a discussion centered around the Korteweg–de Vries equation, SIAM Rev. 14 (4) (1972) 582–643. [5] P.D. Lax, Integrals of nonlinear equations of evolution and solitary waves, Commun. Pure Appl. Math. 21 (5) (1968) 467–490. [6] A. Sjöberg, On the Korteweg–de Vries equation: existence and uniqueness, J. Math. Anal. Appl. 29 (3) (1970) 569–579. [7] J. Yan, C.W. Shu, A local discontinuous Galerkin method for KdV type equations, SIAM J. Numer. Anal. 40 (2) (2002) 769–791. [8] J. Jackaman, G. Papamikos, T. Pryer, The design of conservative finite element discretisations for the vectorial modified KdV equation, Appl. Numer. Math. 137 (2019) 230–251. [9] L. Brugnano, G. Gurioli, Y. Sun, Energy-conserving Hamiltonian boundary value methods for the numerical solution of the Korteweg–de Vries equation, J. Comput. Appl. Math. 351 (2019) 117–135. [10] N.J. Zabusky, M.D. Kruskal, Interaction of “solitons” in a collisionless plasma and the recurrence of initial states, Phys. Rev. Lett. 15 (6) (1965) 240–243. [11] I.S. Greig, J.L. Morris, A hopscotch method for the Korteweg–de Vries equation, J. Comput. Phys. 20 (1) (1976) 64–80. [12] V.I. Mazhukin, P.P. Matus, I.A. Mikhailyuk, Finite-difference schemes for the Korteweg–de Vries equation, Differ. Equ. 36 (5) (20 0 0) 789–797. [13] S. Özer, S. Kutluay, An analytical-numerical method for solving the Korteweg–de Vries equation, Appl. Math. Comput. 164 (3) (2005) 789–797. [14] A.R. Bahadir, Exponential finite-difference method applied to Korteweg–de Vries equation for small times, Appl. Math. Comput. 160 (3) (2005) 675–682. [15] E.N. Aksan, A. Özdes¸ , Numerical solution of Korteweg–de Vries equation by Galerkin b-spline finite element method, Appl. Math. Comput. 175 (2) (2006) 1256–1265. [16] S.Y. Hao, S.S. Xie, S.C. Yi, The Galerkin method for the KdV equation using a new basis of smooth piecewise cubic polynomials, Appl. Math. Comput. 218 (17) (2012) 8659–8671. [17] H.P. Ma, W.W. Sun, A Legendre–Petrov–Galerkin and Chebyshev collocation method for third-order differential equations, SIAM J. Numer. Anal. 38 (5) (20 0 0) 1425–1438. [18] H.P. Ma, W.W. Sun, Optimal error estimates of the Legendre–Petrov–Galerkin method for the Korteweg–de Vries equation, SIAM J. Numer. Anal. 39 (4) (2001) 1380–1394. [19] J. Shen, A new dual-Petrov–Galerkin method for third and higher odd-order differential equations: application to the KdV equation, SIAM J. Numer. Anal. 41 (5) (2003) 1595–1619. [20] N. Bellomo, L. Ridolfi, Solution of nonlinear initial-boundary value problems by sinc collocation–interpolation methods, Comput. Math. Appl. 29 (4) (1995) 15–28. [21] R. Revelli, L. Ridolfi, Sinc collocation–interpolation method for the simulation of nonlinear waves, Comput. Math. Appl. 46 (8) (2003) 1443–1453. [22] K. Al-Khaled, N. Haynes, W. Schiesser, M. Usman, Eventual periodicity of the forced oscillations for a Korteweg–de Vries type equation on a bounded domain using a sinc collocation method, J. Comput. Appl. Math. 330 (2018) 417–428. ˙ [23] D. Irk, D. Idris , S. Bülent, A small time solutions for the Korteweg–de Vries equation using spline approximation, Appl. Math. Comput. 173 (2) (2006) 834–846. [24] S. Kutluay, A.R. Bahadír, A. Özdes¸ , A small time solutions for the Korteweg–de Vries equation, Appl. Math. Comput. 107 (2–3) (20 0 0) 203–210. [25] J.C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, Vol. 88, SIAM, 2004. [26] L.N. Trefethen, Spectral Methods in Matlab, SIAM, 20 0 0. [27] J. Lund, K.L. Bowers, in: Sinc Methods for Quadrature and Differential Equations, SIAM, 1992. [28] F. Stenger, Numerical Methods Based on Sinc and Analytic Functions, Springer Science & Business Media, 2012.