# A Compact-Type CIP Method for General Korteweg-de Vries Equation

Authors: YuFeng Shi; Biao XU; Yan Guo
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We proposed a hybrid compact-CIP scheme to solve the Korteweg-de Vries equation. The algorithm is based on classical constrained interpolation profile (CIP) method, which is coupled with high-order compact scheme for the third derivatives in Korteweg-de Vries equation. Several numerical examples are presented to confirm the high resolution of the proposed scheme.


---

## Body

## 1. Introduction

The Korteweg-de Vries (KdV) equation, developed by Korteweg and de Vries [1] in 1895 to model weakly nonlinear waves, has been used in many different fields to model various physical phenomena of interest. In recent years, a number of numerical methods are proposed for solving KdV equations. In [2], Zabusky and Kruskal proposed a finite difference method for KdV equations. A local discontinuous Galerkin method was developed for solving KdV type equations containing third derivative terms in one and two space dimensions in [3]. Numerical solutions for general KdV equation with Crank Nicolson method and B-spline FEM were compared with those obtained with Adomian decomposition method (ADM) in [4, 5]. Based on the multisymplectic theory, lots of numerical schemes were proposed for KdV equations [6–8]. Ascher and McLachlan developed a simplified 8-point box scheme [9]. By adding an artificial numerical condition to the periodic boundary, Wang et al. derived some new schemes and proved that they were more efficient than the Preissman scheme in [10]. A completely explicit 6-point multisymplectic scheme is derived in [11]. Recently many other numerical methods have been proposed for solving KdV type equations in [12–18].In recent years, the less diffusive CIP scheme developed by Takewaki et al. [19] for solving hyperbolic equation has become very popular. However, the original CIP method [19–22] utilizing both the point values and its spatial gradients needs auxiliary boundary conditions for the spatial gradient information. Usually, it has to differentiate the equation with spatial variable to get the values of derivation on the node. For the simple case, where the velocity is constant, the procedure is not difficult, but it is not easy for complex equations. In 1992, Lele [23] developed high-order compact (HOC) difference schemes based on implicit interpolations for first and second derivatives. The implicit schemes are very accurate in smooth regions and have spectral-like resolution properties by using the global grid. High-order compact finite difference schemes coupled with high-order low-pass filter are applied to simulate KdV equations in [24].In this paper, a new numerical scheme based on classical CIP and HOC schemes is proposed to solve KdV type equations. The new scheme is based on CIP scheme; as a new ingredient, the classical high-order compact scheme [23] is employed to obtain the derivatives rather than differentiate the equation with spatial variable to construct CIP scheme. By comparing with classical compact scheme for solving KdV equations, no filter is used to formulate the present scheme.The paper is organized as follows. In Section2, we give a brief description of CIP and high-order compact schemes. The numerical arithmetic of the present scheme is also discussed in the last part of this section. The implementation of our present method for KdV type equations is shown in Section 3, and the capability of the method for nonlinear dispersive equations can be observed from the comparison of numerical solutions with exact solutions. A short discussion for the present method is given in Section 4.


## 2. Descriptions of Methods

In this paper, we consider the following generalized KdV equation [25]:
(1)
u
t
+
a
(
u
)
u
x
+
ε
u
x
x
x
=
0
.
The equation can be split into two parts:
(2)
u
t
+
a
(
u
)
u
x
=
0
,
(3)
u
t
=
-
ε
u
x
x
x
,
where ε are real constants. We only consider the advective phase (2) to review CIP method.


### 2.1. The CIP Method

The CIP method in [21] uses cubic-polynomial interpolation to get the values of a function on nodes. The time evolution of spatial derivation is also required (4). We differentiate the advective phase of (2) with the spatial variable, and then the equation for derivatives of u can be obtained [26]. Consider
(4)
∂
g
∂
t
+
a
(
u
)
∂
g
∂
x
=
-
g
∂
a
(
u
)
∂
x
,
where g
=
∂
u
/
∂
x stands for the spatial derivatives of u. The computational domain [
a
,
b
] can be divided into N cells, and the cells are denoted by I
i
=
[
x
i
,
x
i
+
1
]. We only consider a uniform grid and the size of the cell by Δ
x
=
(
b
-
a
)
/
N. The cubic polynomial and the first-order derivative at the nth step can be written as
(5)
U
i
n
(
x
)
=
a
i
X
3
+
b
i
X
2
+
c
i
X
+
d
i
,
U
i
′
n
(
x
)
=
3
a
i
X
2
+
2
b
i
X
+
c
i
,
where X
=
x
-
x
i, and coefficients a
i, b
i, c
i, and d
i will be obtained with the following constrains:
(6)
U
i
n
(
x
i
)
=
u
i
n
,
U
i
n
(
x
i
u
p
)
=
u
i
u
p
n
,
U
i
′
n
(
x
i
)
=
g
i
n
,
U
i
′
n
(
x
i
u
p
)
=
g
i
u
p
n
,
where i
u
p
=
i
-
sgn
⁡
(
a
(
u
i
)
), the sign sgn
⁡
(
a
(
u
i
)
) stands for the sign of a
(
u
i
). Then, the coefficients of the cubic polynomial are given as follows:
(7)
a
i
=
g
i
n
+
g
i
u
p
n
Δ
x
i
2
+
2
(
u
i
n
-
u
i
u
p
n
)
Δ
x
i
3
,
b
i
=
3
(
u
i
u
p
n
-
u
i
n
)
Δ
x
i
2
-
2
g
i
n
+
g
i
u
p
n
Δ
x
i
,
c
i
=
g
i
n
,
d
i
=
u
i
n
,
where Δ
x
i
=
x
i
u
p
-
x
i. Thus, the values of u and g at the (
n
+
1)th step can be obtained as follows:
(8)
u
i
n
+
1
=
U
i
n
(
x
i
-
a
(
u
i
)
Δ
t
)
,
g
i
n
+
1
=
U
i
′
n
(
x
i
-
a
(
u
i
)
Δ
t
)
.We defineξ
i
=
-
a
(
u
i
)
Δ
t, and then the formulas are rewritten as
(9)
u
i
n
+
1
=
a
i
ξ
i
3
+
b
i
ξ
i
2
+
g
i
n
ξ
i
+
u
i
n
,
g
i
n
+
1
=
3
a
i
ξ
i
n
+
2
b
i
ξ
i
n
+
g
i
n
.It can be seen that only two points are used in CIP schemes to getu
i
n
+
1. Then, the advantages of this method can be shown while the computational boundary is complex since less boundary points need to be handled. The CIP method uses only two neighboring stencils but maintains third-order accuracy. In this sense, high order accuracy is obtained though less computational stencils are used. More details for advantages of the CIP schemes can be found in [26].


### 2.2. High-Order Compact Scheme

A series finite difference scheme to evaluate the spatial derivatives is presented in [23]. The finite difference approximation to the derivative of the function is expressed as a linear combination of the given function values, and then the derivatives of the function are obtained by solving a tridiagonal or pentadiagonal system. Formulas for the first-order and the third-order derivatives are reviewed as below. More results for the approximation to derivatives can be found in [23, 27].


#### 2.2.1. The Derivatives at Interior Nodes

In this paper, we consider the KdV equation on a uniform mesh, the spatial variable at the nodes isx
i
=
i
×
h for 0
≤
i
≤
N and the functions and the derivatives are denoted by u
i, u
i
′. The first-order derivatives for interior nodes are derived by writing approximations of the form [23]. Consider
(10)
u
i
′
+
α
(
u
i
-
1
′
+
u
i
+
1
′
)
+
β
(
u
i
-
2
′
+
u
i
+
2
′
)
=
c
u
i
+
3
-
u
i
-
3
6
h
+
b
u
i
+
2
-
u
i
-
2
4
h
+
a
u
i
+
1
-
u
i
-
1
2
h
.
If the schemes are restricted to β
≥
0 and c
=
0, a one-parameter α-family of fourth-order tridiagonal scheme is obtained. Consider
(11)
β
=
0
,
c
=
0
,
a
=
2
3
(
α
+
2
)
,
b
=
1
3
(
4
α
-
1
)
.A simple sixth-order tridiagonal scheme for first-order derivatives is given withα
=
1
/
3, β
=
0, c
=
0, a
=
14
/
9, and b
=
1
/
9:
(12)
u
i
′
+
1
3
(
u
i
-
1
′
+
u
i
+
1
′
)
=
14
9
u
i
+
1
-
u
i
-
1
2
h
+
1
9
u
i
+
2
-
u
i
-
2
2
h
.For the third derivatives at interior nodes, the following formula is given in [23]:
(13)
α
(
u
i
-
1
′′′
+
u
i
+
1
′′′
)
+
u
i
′′′
=
b
u
i
+
3
-
3
u
i
+
1
+
3
u
i
-
1
-
u
i
-
3
8
h
3
+
a
u
i
+
2
-
2
u
i
+
1
+
2
u
i
-
1
-
u
i
-
2
2
h
3
.The fourth-order tridiagonal schemes can be obtained with the coefficientsa
=
2 and b
=
2
α
-
1. The compact tridiagonal scheme is given with α
=
1
/
2, a
=
2, and b
=
0. And the sixth-order tridiagonal scheme is given with α
=
7
/
16, a
=
2, and b
=
-
1
/
8.


#### 2.2.2. Nonperiodic Boundaries

For those near boundary nodes, approximation formulas for the first-order derivatives of nonperiodic boundary problems are given by one-side formulation as follows [23]:
(14)
u
1
′
+
α
u
2
′
=
1
h
(
a
u
1
+
b
u
2
+
c
u
3
+
d
u
4
)
,
u
N
′
+
α
u
N
-
1
′
=
-
1
h
(
a
u
N
+
b
u
N
-
1
+
c
u
N
-
2
+
d
u
N
-
3
)
.
The coefficients for schemes of the third and fourth order are given by
(15)
Third order:
a
=
-
11
+
2
α
6
,
b
=
6
-
α
2
,
c
=
2
α
-
3
2
,
d
=
2
-
α
6
,
Fourth order:
α
=
3
,
a
=
-
17
6
,
b
=
3
2
,
c
=
3
2
,
d
=
-
1
6
.The sixth-order scheme also be given, for those near boundary nodes, six order approximation formulas for the first-order derivatives can be written as follows:(16)
u
1
′
+
α
u
2
′
=
1
h
(
a
1
u
1
+
a
2
u
2
+
a
3
u
3
+
a
4
u
4
+
a
5
u
5
+
a
6
u
6
)
,
where
(17)
α
=
5
,
a
1
=
-
197
60
,
a
2
=
-
12
5
,
a
3
=
5
,
a
4
=
-
5
3
,
a
5
=
5
12
,
a
6
=
-
1
20
.For the second point, the formula is(18)
α
u
1
′
+
u
2
′
+
α
u
3
′
=
1
h
(
b
1
u
1
+
b
2
u
2
+
b
3
u
3
+
b
4
u
4
+
b
5
u
5
+
b
6
u
6
)
,
where
(19)
α
=
2
11
,
b
1
=
-
20
33
,
b
2
=
-
35
132
,
b
3
=
34
33
,
b
4
=
-
7
33
,
b
5
=
2
33
,
b
6
=
-
1
132
.The dissymmetry condition is used for theNth and (
N
-
1
)th points.


### 2.3. The Present Compact-Type CIP Method

In this section, a new compact-type CIP scheme is proposed for (1). For simplicity, the following equation is considered to introduce the method. Consider
(20)
u
t
+
α
u
u
x
-
δ
u
x
x
x
=
0
,
where α and δ are constants. We split the solution of equation into two phases:
(21)
∂
u
∂
t
+
α
u
∂
u
∂
x
=
0
,
(22)
∂
u
∂
t
=
δ
u
x
x
x
.We consider a 1D mesh, consisting of(
N
+
1
) points: x
0, x
1, x
2
,
…, x
N
-
1, x
N; the values of u on these nodes at the nth step are denoted by u
0
n, u
1
n
…, u
N
-
1
n, u
N
n. At first, (21) is considered and CIP method is applied to the equation. The cubic polynomial at the nth time stage is
(23)
U
i
n
(
X
)
=
a
i
n
X
3
+
b
i
n
X
2
+
c
i
n
X
+
u
i
n
,
where X
=
x
i
-
x and the coefficient a
i
n, b
i
n, and c
i
n are given by (7), where a
(
u
)
=
α
u and u
i
′
n denote the derivative of u
i
n at the ith node. The predictor-corrector scheme is employed to calculate the value u
*.In the present method, the valuesu
i
′
n, 0
≤
i
≤
N are expressed as a linear combination of the given values u
i
n, 0
≤
i
≤
N. On the other hand, the HOC method is employed to evaluate the derivatives u
i
′
n, 0
≤
i
≤
N. A simple sixth-order tridiagonal scheme for interior points and boundary points is used in this paper.Temporal discretization for (22) can be solved by using third-order Runge-Kutta method:
(24)
u
(
1
)
=
u
n
+
δ
t
u
x
x
x
*
,
u
(
2
)
=
3
4
u
n
+
1
4
u
(
1
)
+
1
4
δ
t
u
x
x
x
(
1
)
,
u
n
+
1
=
1
3
u
n
+
2
3
u
(
2
)
+
2
3
δ
t
u
x
x
x
(
2
)
.The HOC scheme (13) is used to solve the third derivatives u
x
x
x in (24). The sixth-order tridiagonal scheme with the periodic boundary condition is used in this paper.The essential ingredients of the computational algorithm for (20) are presented below. Suppose we have got the values u
i
n. The values u
i
n
+
1 are given as follows.(1)
CIP method is used to obtainu
*.
(a)
The values of the first derivative on the all nodes are obtained by using the HOC scheme (12).
(b)
Predictor-corrector CIP scheme is as follows:
(i)
predictor step:(25)
u
i
*
*
=
U
i
n
(
x
i
-
α
u
i
n
Δ
t
)
=
a
i
n
ξ
i
3
+
b
i
n
ξ
i
2
+
c
i
n
ξ
i
+
u
i
n
,
where ξ
i
=
-
α
u
i
n
Δ
t. We also get u
*
*
* at the (n
+
1
/
2)th time stage by using linear interpolation or QUICK scheme based on the value u
i
n;
(ii)
corrector step (CIP method):(26)
u
^
i
*
=
U
i
n
(
x
i
-
α
u
i
⋄
Δ
t
)
=
a
i
n
ξ
i
3
+
b
i
n
ξ
i
2
+
c
i
n
ξ
i
+
u
i
n
,
where u
⋄
=
(
1
/
2
)
(
u
*
*
+
u
*
*
*
);
(iii)
the predictor and corrector step are employed again to getu
*.
(2)
HOC scheme and Runge-Kutta method for (22) are as follows.
(a)
The HOC scheme (13) is used to obtain third-order derivatives.
(b)
Temporal discretization for (22) can be solved by using third-order Runge-Kutta method.The predictor-corrector scheme is an important step in the present method. Periodic boundary condition is applied to (22).


## 2.1. The CIP Method

The CIP method in [21] uses cubic-polynomial interpolation to get the values of a function on nodes. The time evolution of spatial derivation is also required (4). We differentiate the advective phase of (2) with the spatial variable, and then the equation for derivatives of u can be obtained [26]. Consider
(4)
∂
g
∂
t
+
a
(
u
)
∂
g
∂
x
=
-
g
∂
a
(
u
)
∂
x
,
where g
=
∂
u
/
∂
x stands for the spatial derivatives of u. The computational domain [
a
,
b
] can be divided into N cells, and the cells are denoted by I
i
=
[
x
i
,
x
i
+
1
]. We only consider a uniform grid and the size of the cell by Δ
x
=
(
b
-
a
)
/
N. The cubic polynomial and the first-order derivative at the nth step can be written as
(5)
U
i
n
(
x
)
=
a
i
X
3
+
b
i
X
2
+
c
i
X
+
d
i
,
U
i
′
n
(
x
)
=
3
a
i
X
2
+
2
b
i
X
+
c
i
,
where X
=
x
-
x
i, and coefficients a
i, b
i, c
i, and d
i will be obtained with the following constrains:
(6)
U
i
n
(
x
i
)
=
u
i
n
,
U
i
n
(
x
i
u
p
)
=
u
i
u
p
n
,
U
i
′
n
(
x
i
)
=
g
i
n
,
U
i
′
n
(
x
i
u
p
)
=
g
i
u
p
n
,
where i
u
p
=
i
-
sgn
⁡
(
a
(
u
i
)
), the sign sgn
⁡
(
a
(
u
i
)
) stands for the sign of a
(
u
i
). Then, the coefficients of the cubic polynomial are given as follows:
(7)
a
i
=
g
i
n
+
g
i
u
p
n
Δ
x
i
2
+
2
(
u
i
n
-
u
i
u
p
n
)
Δ
x
i
3
,
b
i
=
3
(
u
i
u
p
n
-
u
i
n
)
Δ
x
i
2
-
2
g
i
n
+
g
i
u
p
n
Δ
x
i
,
c
i
=
g
i
n
,
d
i
=
u
i
n
,
where Δ
x
i
=
x
i
u
p
-
x
i. Thus, the values of u and g at the (
n
+
1)th step can be obtained as follows:
(8)
u
i
n
+
1
=
U
i
n
(
x
i
-
a
(
u
i
)
Δ
t
)
,
g
i
n
+
1
=
U
i
′
n
(
x
i
-
a
(
u
i
)
Δ
t
)
.We defineξ
i
=
-
a
(
u
i
)
Δ
t, and then the formulas are rewritten as
(9)
u
i
n
+
1
=
a
i
ξ
i
3
+
b
i
ξ
i
2
+
g
i
n
ξ
i
+
u
i
n
,
g
i
n
+
1
=
3
a
i
ξ
i
n
+
2
b
i
ξ
i
n
+
g
i
n
.It can be seen that only two points are used in CIP schemes to getu
i
n
+
1. Then, the advantages of this method can be shown while the computational boundary is complex since less boundary points need to be handled. The CIP method uses only two neighboring stencils but maintains third-order accuracy. In this sense, high order accuracy is obtained though less computational stencils are used. More details for advantages of the CIP schemes can be found in [26].


## 2.2. High-Order Compact Scheme

A series finite difference scheme to evaluate the spatial derivatives is presented in [23]. The finite difference approximation to the derivative of the function is expressed as a linear combination of the given function values, and then the derivatives of the function are obtained by solving a tridiagonal or pentadiagonal system. Formulas for the first-order and the third-order derivatives are reviewed as below. More results for the approximation to derivatives can be found in [23, 27].


### 2.2.1. The Derivatives at Interior Nodes

In this paper, we consider the KdV equation on a uniform mesh, the spatial variable at the nodes isx
i
=
i
×
h for 0
≤
i
≤
N and the functions and the derivatives are denoted by u
i, u
i
′. The first-order derivatives for interior nodes are derived by writing approximations of the form [23]. Consider
(10)
u
i
′
+
α
(
u
i
-
1
′
+
u
i
+
1
′
)
+
β
(
u
i
-
2
′
+
u
i
+
2
′
)
=
c
u
i
+
3
-
u
i
-
3
6
h
+
b
u
i
+
2
-
u
i
-
2
4
h
+
a
u
i
+
1
-
u
i
-
1
2
h
.
If the schemes are restricted to β
≥
0 and c
=
0, a one-parameter α-family of fourth-order tridiagonal scheme is obtained. Consider
(11)
β
=
0
,
c
=
0
,
a
=
2
3
(
α
+
2
)
,
b
=
1
3
(
4
α
-
1
)
.A simple sixth-order tridiagonal scheme for first-order derivatives is given withα
=
1
/
3, β
=
0, c
=
0, a
=
14
/
9, and b
=
1
/
9:
(12)
u
i
′
+
1
3
(
u
i
-
1
′
+
u
i
+
1
′
)
=
14
9
u
i
+
1
-
u
i
-
1
2
h
+
1
9
u
i
+
2
-
u
i
-
2
2
h
.For the third derivatives at interior nodes, the following formula is given in [23]:
(13)
α
(
u
i
-
1
′′′
+
u
i
+
1
′′′
)
+
u
i
′′′
=
b
u
i
+
3
-
3
u
i
+
1
+
3
u
i
-
1
-
u
i
-
3
8
h
3
+
a
u
i
+
2
-
2
u
i
+
1
+
2
u
i
-
1
-
u
i
-
2
2
h
3
.The fourth-order tridiagonal schemes can be obtained with the coefficientsa
=
2 and b
=
2
α
-
1. The compact tridiagonal scheme is given with α
=
1
/
2, a
=
2, and b
=
0. And the sixth-order tridiagonal scheme is given with α
=
7
/
16, a
=
2, and b
=
-
1
/
8.


### 2.2.2. Nonperiodic Boundaries

For those near boundary nodes, approximation formulas for the first-order derivatives of nonperiodic boundary problems are given by one-side formulation as follows [23]:
(14)
u
1
′
+
α
u
2
′
=
1
h
(
a
u
1
+
b
u
2
+
c
u
3
+
d
u
4
)
,
u
N
′
+
α
u
N
-
1
′
=
-
1
h
(
a
u
N
+
b
u
N
-
1
+
c
u
N
-
2
+
d
u
N
-
3
)
.
The coefficients for schemes of the third and fourth order are given by
(15)
Third order:
a
=
-
11
+
2
α
6
,
b
=
6
-
α
2
,
c
=
2
α
-
3
2
,
d
=
2
-
α
6
,
Fourth order:
α
=
3
,
a
=
-
17
6
,
b
=
3
2
,
c
=
3
2
,
d
=
-
1
6
.The sixth-order scheme also be given, for those near boundary nodes, six order approximation formulas for the first-order derivatives can be written as follows:(16)
u
1
′
+
α
u
2
′
=
1
h
(
a
1
u
1
+
a
2
u
2
+
a
3
u
3
+
a
4
u
4
+
a
5
u
5
+
a
6
u
6
)
,
where
(17)
α
=
5
,
a
1
=
-
197
60
,
a
2
=
-
12
5
,
a
3
=
5
,
a
4
=
-
5
3
,
a
5
=
5
12
,
a
6
=
-
1
20
.For the second point, the formula is(18)
α
u
1
′
+
u
2
′
+
α
u
3
′
=
1
h
(
b
1
u
1
+
b
2
u
2
+
b
3
u
3
+
b
4
u
4
+
b
5
u
5
+
b
6
u
6
)
,
where
(19)
α
=
2
11
,
b
1
=
-
20
33
,
b
2
=
-
35
132
,
b
3
=
34
33
,
b
4
=
-
7
33
,
b
5
=
2
33
,
b
6
=
-
1
132
.The dissymmetry condition is used for theNth and (
N
-
1
)th points.


## 2.2.1. The Derivatives at Interior Nodes

In this paper, we consider the KdV equation on a uniform mesh, the spatial variable at the nodes isx
i
=
i
×
h for 0
≤
i
≤
N and the functions and the derivatives are denoted by u
i, u
i
′. The first-order derivatives for interior nodes are derived by writing approximations of the form [23]. Consider
(10)
u
i
′
+
α
(
u
i
-
1
′
+
u
i
+
1
′
)
+
β
(
u
i
-
2
′
+
u
i
+
2
′
)
=
c
u
i
+
3
-
u
i
-
3
6
h
+
b
u
i
+
2
-
u
i
-
2
4
h
+
a
u
i
+
1
-
u
i
-
1
2
h
.
If the schemes are restricted to β
≥
0 and c
=
0, a one-parameter α-family of fourth-order tridiagonal scheme is obtained. Consider
(11)
β
=
0
,
c
=
0
,
a
=
2
3
(
α
+
2
)
,
b
=
1
3
(
4
α
-
1
)
.A simple sixth-order tridiagonal scheme for first-order derivatives is given withα
=
1
/
3, β
=
0, c
=
0, a
=
14
/
9, and b
=
1
/
9:
(12)
u
i
′
+
1
3
(
u
i
-
1
′
+
u
i
+
1
′
)
=
14
9
u
i
+
1
-
u
i
-
1
2
h
+
1
9
u
i
+
2
-
u
i
-
2
2
h
.For the third derivatives at interior nodes, the following formula is given in [23]:
(13)
α
(
u
i
-
1
′′′
+
u
i
+
1
′′′
)
+
u
i
′′′
=
b
u
i
+
3
-
3
u
i
+
1
+
3
u
i
-
1
-
u
i
-
3
8
h
3
+
a
u
i
+
2
-
2
u
i
+
1
+
2
u
i
-
1
-
u
i
-
2
2
h
3
.The fourth-order tridiagonal schemes can be obtained with the coefficientsa
=
2 and b
=
2
α
-
1. The compact tridiagonal scheme is given with α
=
1
/
2, a
=
2, and b
=
0. And the sixth-order tridiagonal scheme is given with α
=
7
/
16, a
=
2, and b
=
-
1
/
8.


## 2.2.2. Nonperiodic Boundaries

For those near boundary nodes, approximation formulas for the first-order derivatives of nonperiodic boundary problems are given by one-side formulation as follows [23]:
(14)
u
1
′
+
α
u
2
′
=
1
h
(
a
u
1
+
b
u
2
+
c
u
3
+
d
u
4
)
,
u
N
′
+
α
u
N
-
1
′
=
-
1
h
(
a
u
N
+
b
u
N
-
1
+
c
u
N
-
2
+
d
u
N
-
3
)
.
The coefficients for schemes of the third and fourth order are given by
(15)
Third order:
a
=
-
11
+
2
α
6
,
b
=
6
-
α
2
,
c
=
2
α
-
3
2
,
d
=
2
-
α
6
,
Fourth order:
α
=
3
,
a
=
-
17
6
,
b
=
3
2
,
c
=
3
2
,
d
=
-
1
6
.The sixth-order scheme also be given, for those near boundary nodes, six order approximation formulas for the first-order derivatives can be written as follows:(16)
u
1
′
+
α
u
2
′
=
1
h
(
a
1
u
1
+
a
2
u
2
+
a
3
u
3
+
a
4
u
4
+
a
5
u
5
+
a
6
u
6
)
,
where
(17)
α
=
5
,
a
1
=
-
197
60
,
a
2
=
-
12
5
,
a
3
=
5
,
a
4
=
-
5
3
,
a
5
=
5
12
,
a
6
=
-
1
20
.For the second point, the formula is(18)
α
u
1
′
+
u
2
′
+
α
u
3
′
=
1
h
(
b
1
u
1
+
b
2
u
2
+
b
3
u
3
+
b
4
u
4
+
b
5
u
5
+
b
6
u
6
)
,
where
(19)
α
=
2
11
,
b
1
=
-
20
33
,
b
2
=
-
35
132
,
b
3
=
34
33
,
b
4
=
-
7
33
,
b
5
=
2
33
,
b
6
=
-
1
132
.The dissymmetry condition is used for theNth and (
N
-
1
)th points.


## 2.3. The Present Compact-Type CIP Method

In this section, a new compact-type CIP scheme is proposed for (1). For simplicity, the following equation is considered to introduce the method. Consider
(20)
u
t
+
α
u
u
x
-
δ
u
x
x
x
=
0
,
where α and δ are constants. We split the solution of equation into two phases:
(21)
∂
u
∂
t
+
α
u
∂
u
∂
x
=
0
,
(22)
∂
u
∂
t
=
δ
u
x
x
x
.We consider a 1D mesh, consisting of(
N
+
1
) points: x
0, x
1, x
2
,
…, x
N
-
1, x
N; the values of u on these nodes at the nth step are denoted by u
0
n, u
1
n
…, u
N
-
1
n, u
N
n. At first, (21) is considered and CIP method is applied to the equation. The cubic polynomial at the nth time stage is
(23)
U
i
n
(
X
)
=
a
i
n
X
3
+
b
i
n
X
2
+
c
i
n
X
+
u
i
n
,
where X
=
x
i
-
x and the coefficient a
i
n, b
i
n, and c
i
n are given by (7), where a
(
u
)
=
α
u and u
i
′
n denote the derivative of u
i
n at the ith node. The predictor-corrector scheme is employed to calculate the value u
*.In the present method, the valuesu
i
′
n, 0
≤
i
≤
N are expressed as a linear combination of the given values u
i
n, 0
≤
i
≤
N. On the other hand, the HOC method is employed to evaluate the derivatives u
i
′
n, 0
≤
i
≤
N. A simple sixth-order tridiagonal scheme for interior points and boundary points is used in this paper.Temporal discretization for (22) can be solved by using third-order Runge-Kutta method:
(24)
u
(
1
)
=
u
n
+
δ
t
u
x
x
x
*
,
u
(
2
)
=
3
4
u
n
+
1
4
u
(
1
)
+
1
4
δ
t
u
x
x
x
(
1
)
,
u
n
+
1
=
1
3
u
n
+
2
3
u
(
2
)
+
2
3
δ
t
u
x
x
x
(
2
)
.The HOC scheme (13) is used to solve the third derivatives u
x
x
x in (24). The sixth-order tridiagonal scheme with the periodic boundary condition is used in this paper.The essential ingredients of the computational algorithm for (20) are presented below. Suppose we have got the values u
i
n. The values u
i
n
+
1 are given as follows.(1)
CIP method is used to obtainu
*.
(a)
The values of the first derivative on the all nodes are obtained by using the HOC scheme (12).
(b)
Predictor-corrector CIP scheme is as follows:
(i)
predictor step:(25)
u
i
*
*
=
U
i
n
(
x
i
-
α
u
i
n
Δ
t
)
=
a
i
n
ξ
i
3
+
b
i
n
ξ
i
2
+
c
i
n
ξ
i
+
u
i
n
,
where ξ
i
=
-
α
u
i
n
Δ
t. We also get u
*
*
* at the (n
+
1
/
2)th time stage by using linear interpolation or QUICK scheme based on the value u
i
n;
(ii)
corrector step (CIP method):(26)
u
^
i
*
=
U
i
n
(
x
i
-
α
u
i
⋄
Δ
t
)
=
a
i
n
ξ
i
3
+
b
i
n
ξ
i
2
+
c
i
n
ξ
i
+
u
i
n
,
where u
⋄
=
(
1
/
2
)
(
u
*
*
+
u
*
*
*
);
(iii)
the predictor and corrector step are employed again to getu
*.
(2)
HOC scheme and Runge-Kutta method for (22) are as follows.
(a)
The HOC scheme (13) is used to obtain third-order derivatives.
(b)
Temporal discretization for (22) can be solved by using third-order Runge-Kutta method.The predictor-corrector scheme is an important step in the present method. Periodic boundary condition is applied to (22).


## 3. Numerical Results

In this section, some numerical tests for KdV and general KdV equations are carried out. The discreteL
2 and L
∞ error norms are defined as follows:
(27)
∥
e
n
∥
L
∞
=
max
⁡
0
≤
j
≤
N
|
u
j
-
u
~
j
|
,
∥
e
n
∥
L
2
=
(
∑
j
=
0
N
‍
|
u
j
-
u
~
j
|
2
Δ
x
)
1
/
2
,
where u and u
~ are exact and numerical solution, respectively. For KdV equations, there are an infinite number of conservation laws [28]. We will focus our analysis on the following three conservation laws:
(28)
I
1
=
∫
a
b
‍
u
(
x
,
t
)
d
x
,
I
2
=
∫
a
b
‍
u
(
x
,
t
)
2
d
x
,
I
3
=
∫
a
b
‍
(
u
(
x
,
t
)
2
-
1
3
u
(
x
,
t
)
3
)
d
x
,
where I
1, I
2, and I
3 represent mass, momentum, and energy.The nonperiodic boundary formulation is applied to (21) (HOC approximation formulas for the third-order derivatives are used) and periodic boundary conditions for third-order derivatives are used in the following examples.Example 1.
In this example, we consider the following classical KdV equation:(29)
u
t
+
6
u
u
x
+
u
x
x
x
=
0
,
-
15
≤
x
≤
15
.
The analytical solution for (29) is
(30)
u
(
x
,
t
)
=
0.5
sec
h
2
(
0.5
(
x
-
t
)
)
.
The time-steps are set by the relationΔ
t
=
C
(
Δ
x
)
3. The L
2 and L
∞ errors, orders, invariants, and time costs at time T
=
0.1 are illustrated in Table 2. It can be observed that the proposed scheme is third-order accurate in the spatial dimension. It is well known that high-order TVD Runge-Kutta methods suffer from small time-step restrictions. In this case, we observe that numerical errors are still not dominated by the spatial discretization with the relation Δ
t
=
C
(
Δ
x
)
3.
Table1 indicates L
2 and L
∞ errors and invariants with N
=
100 at time T
=
0.1
,
0.2
,
…
,
0.9. The present method can also be shown to have the conservative property.Table 1
Numerical errors and orders of CIP-HOC method for the KdV equation (29).
N
L
∞ error
L
∞ order
L
2 error
L
2 order
I
1
I
2
I
3
Cost(s)
20
5.72
E
-
03
9.94
E
-
03
1.9968246
0.6627258
0.5745942
6.25
E
-
05
40
1.32
E
-
03
2.12
1.89
E
-
03
2.39
1.9996543
0.6663460
0.5774709
2.50
E
-
04
80
2.11
E
-
04
2.65
2.85
E
-
04
2.73
1.9999840
0.6666518
0.5777566
2.50
E
-
03
160
2.77
E
-
05
2.93
3.61
E
-
05
2.98
1.9999984
0.6666662
0.5777761
3.13
E
-
02
320
3.49
E
-
06
2.99
4.54
E
-
06
2.99
1.9999988
0.6666667
0.5777776
6.09
E
-
01Table 2
Numerical errors and invariants withN
=
100 at different times.
Time
L
∞ error
L
2 error
I
1
I
2
I
3
0.10
1.09
E
-
04
1.47
E
-
04
1.9999939
0.6666616
0.5777688
0.20
1.65
E
-
04
2.47
E
-
04
1.9999888
0.6666565
0.5777615
0.30
1.98
E
-
04
3.14
E
-
04
1.9999830
0.6666515
0.5777556
0.40
2.14
E
-
04
3.61
E
-
04
1.9999775
0.6666464
0.5777505
0.50
2.28
E
-
04
3.94
E
-
04
1.9999760
0.6666413
0.5777461
0.60
2.41
E
-
04
4.18
E
-
04
1.9999660
0.6666363
0.5777421
0.70
2.46
E
-
04
4.37
E
-
04
1.9999590
0.6666312
0.5777383
0.80
2.42
E
-
04
4.52
E
-
04
1.9999651
0.6666262
0.5777347
0.90
2.51
E
-
04
4.63
E
-
04
1.9999605
0.6666211
0.5777312
1.00
2.59
E
-
04
4.73
E
-
04
1.9999441
0.6666161
0.5777277Example 2.
In this example, we consider the following classical KdV equation:(31)
u
t
+
u
u
x
+
u
x
x
x
=
0
,
0
≤
x
≤
40
,
with initial condition
(32)
u
(
x
,
0
)
=
12
cosh
⁡
(
x
-
15
)
2
,
0
≤
x
≤
40
.
The numerical solutions are obtained withd
x
=
0.2 and d
t
=
0.0001; numerical and exact solutions at times t
=
1
,
3
,
5 are presented in Figure 1; the figure shows that the numerical dissipation of soliton is very small.Figure 1
Numerical and analytical solutions for equationu
t
+
u
u
x
+
u
x
x
x
=
0 of Example 2 at various time stages.Example 3.
The general KdV equation is presented by the following equation [29]:
(33)
u
t
+
6
u
p
u
x
+
u
x
x
x
=
0
.
Consider the initial value problem associated with (33) using the initial condition for p
=
2:
(34)
u
(
x
,
0
)
=
c
sech
(
c
x
)
,
-
10
≤
x
≤
12
,
where c
=
4. The numerical solutions for p
=
2 are obtained with d
x
=
0.1 and d
t
=
0.05
×
d
x
3 for p
=
2. The progress of the numerical and analytical solutions at times t
=
1,2 and the absolute error at t
=
1 are shown in Figure 2. The L
∞ and L
2 error estimates for the case p
=
2 are given in Table 3.
For the casep
=
1 with the initial condition [14]
(35)
u
(
x
,
0
)
=
c
2
sec
h
2
(
c
2
x
-
7
)
,
0
≤
x
≤
40
,
and the analytical solution is
(36)
u
(
x
,
t
)
=
c
2
sec
h
2
(
c
2
(
x
-
c
t
)
-
7
)
.
TheL
∞ and L
2 error estimates with d
x
=
0.2 are given in Table 4 for t
=
1,2
,
3,4
,
5, from which it is not difficult to see that the present results are comparable with those present in [14]. The numerical and the analytical solutions at times t
=
1
,
3
,
5 and the absolute error for t
=
3 are shown in Figure 3.Table 3
Error norms at different time stages for equationu
t
+
6
u
2
u
x
+
u
x
x
x
=
0, with u
(
x
,
0
)
=
c
sech
⁡
(
c
x
)
⁡
,
-
10
≤
x
≤
12
,
c
=
0.5.  d
x
=
0.01, and d
t
=
0.05
×
d
x
3.
T
L
∞
L
2
0.1
5.354322207020701
E
-
004
4.972071746712235
E
-
004
0.5
6.080056468253936
E
-
004
7.277368149932270
E
-
004
1.0
2.970167786819244
E
-
003
3.320246422392297
E
-
003Table 4
Error norms at different time stages for equationu
t
+
6
u
u
x
+
u
x
x
x
=
0, with u
(
x
,
0
)
=
(
c
/
2
)
sech
2
⁡
(
(
c
/
2
)
x
-
7
)
,
c
=
0.5.  d
x
=
0.2, and d
t
=
0.1
×
d
x
3.
T
Present
Dehghan and Shokri [14]
L
∞
L
2
L
∞
L
2
1.0
1.9363
E
-
005
2.7752
E
-
005
1.8048
E
-
005
6.2366
E
-
005
2.0
2.8040
E
-
005
4.4840
E
-
005
3.0373
E
-
005
1.1264
E
-
005
3.0
3.5184
E
-
005
6.2377
E
-
005
4.0088
E
-
005
1.5537
E
-
005
4.0
4.2136
E
-
005
8.2619
E
-
005
4.8347
E
-
005
1.9400
E
-
005
5.0
5.3128
E
-
005
1.0617
E
-
004
5.6090
E
-
005
2.2943
E
-
004Figure 2
Numerical and analytical solutions (line) for equationu
t
+
6
u
2
u
x
+
u
x
x
x
=
0 of Example 3 at various time stages and the absolute error for t
=
1.Figure 3
Numerical and analytical solutions for equationu
t
+
6
u
u
x
+
u
x
x
x
=
0 with initial condition u
(
x
,
0
)
=
(
c
/
2
)
sec
h
2
(
(
c
/
2
)
x
-
7
)
,
0
≤
x
≤
40
,
c
=
0.5 of Example 3 at various time stages and the absolute error for t
=
3.Example 4.
In this example, another type of general KdV equation is considered:(37)
u
t
+
u
n
u
x
+
ε
u
x
x
x
=
0
,
with the initial value problem [30]
(38)
u
(
x
,
0
)
=
(
c
(
n
+
1
)
(
n
+
2
)
2
)
1
/
n
sec
h
2
/
n
(
n
2
c
ε
(
x
-
x
0
)
)
,
where c
=
0.3, ε
=
0.000484. The single soliton solutions for n
=
1 are computed in x
∈
[
0,2
] with space step d
x
=
0.02 and time step d
t
=
0.25
×
d
x
2 and are shown in Figure 4. The numerical and analytical solutions for n
=
2 with d
x
=
0.01 and d
t
=
0.5
×
d
x
2 are shown in Figure 4. The L
∞ and L
2 error estimates for n
=
1 at times t
=
0.5
,
1
,
1.5
,
2 are given in Table 5. We can observe that the present results are slightly more accurate than those present in [14]. Table 6 shows the L
∞ and L
2 error estimates for the case of n
=
2.Table 5
Error norms at different time stages for equationu
t
+
u
u
x
+
ε
u
x
x
x
=
0, with initial condition (38) for n
=
1, d
x
=
0.02, and d
t
=
0.25
×
d
x
2.
T
Present
Dehghan and Shokri [14]
L
∞
L
2
L
∞
L
2
0.5
1.7139
E
-
004
5.4536
E
-
005
7.2329
E
-
004
2.5286
E
-
003
1.0
1.6824
E
-
004
5.3897
E
-
005
1.7957
E
-
004
6.2172
E
-
003
1.5
1.9094
E
-
004
6.9749
E
-
005
3.8906
E
-
003
1.3010
E
-
002
2.0
3.8605
E
-
004
1.3214
E
-
004
6.6701
E
-
003
2.2965
E
-
002Table 6
Error norms at different time stages for equationu
t
+
u
2
u
x
+
ε
u
x
x
x
=
0, with initial condition (38) for n
=
2, d
x
=
0.01, and d
t
=
0.5
×
d
x
2.
T
L
∞
L
2
0.5
3.252960639366087
E
-
004
9.332884478549439
E
-
005
1.0
3.962848426597443
E
-
004
1.333167429458641
E
-
004
1.5
1.160144038488964
E
-
003
3.754080717733928
E
-
004
2.0
2.442770975886521
E
-
003
7.776565570880883
E
-
004Numerical and analytical solutions for equationu
t
+
u
2
u
x
+
ε
u
x
x
x
=
0 and u
t
+
u
2
u
x
+
ε
u
x
x
x
=
0 of Example 4 at various time stages.
(a)
u
t
+
u
u
x
+
ε
u
x
x
x
=
0
(b)
u
t
+
u
2
u
x
+
ε
u
x
x
x
=
0Example 5.
We also consider the equation [3]
(39)
u
t
+
u
u
x
+
ε
u
x
x
x
=
0
,
0
≤
x
≤
2
.
The double soliton collision case has the initial condition(40)
u
(
x
,
0
)
=
3
c
1
sec
h
2
(
k
1
(
x
-
x
1
)
)
+
3
c
2
sec
h
2
(
k
2
(
x
-
x
2
)
)
,
where c
1
=
0.3, c
2
=
0.1, x
1
=
0.4
,
x
2
=
0.8, k
i
=
0.5
c
i
/
ε, and ε
=
4.84
×
1
0
-
4. The solution is computed in x
∈
[
0,2
] and is shown in Figure 5. We can observe that nonoscillate numerical solutions can be obtained by using the present method.
The triple soliton collision case has the initial condition(41)
u
(
x
,
0
)
=
2
3
sec
h
2
(
x
-
1
108
ε
)
,
with ε
=
1
0
-
4 [3]. The numerical solution is computed in x
∈
[
0,3
] and is shown in Figure 6. From the figure we can conclude that the numerical algorithm captures the numerical solutions without oscillations.Figure 5
Double soliton solutions for equationu
t
+
u
u
x
+
ε
u
x
x
x
=
0 of Example 5 at various time stages.Figure 6
Triple soliton solutions for equationu
t
+
u
u
x
+
ε
u
x
x
x
=
0 of Example 3 at various time stages.


## 4. Conclusions

In this paper, we have presented a new scheme based on the traditional CIP and HOC scheme. A conclusion can be drawn from the comparison between the numerical and the exact solutions that the present compact-type CIP method provides highly accurate numerical solutions of KdV type equations. The numerical results also show that the present method works well for some nonlinear problems.



---

*Source: 640194-2014-08-03.xml*