# A Dufort-Frankel Difference Scheme for Two-Dimensional Sine-Gordon Equation

Authors: Zongqi Liang; Yubin Yan; Guorong Cai
Journal: Discrete Dynamics in Nature and Society
 (2014)

---

## Abstract

A standard Crank-Nicolson finite-difference scheme and a Dufort-Frankel finite-difference scheme are introduced to solve two-dimensional damped and undamped sine-Gordon equations. The stability and convergence of the numerical methods are considered. To avoid solving the nonlinear system, the predictor-corrector techniques are applied in the numerical methods. Numerical examples are given to show that the numerical results are consistent with the theoretical results.


---

## Body

## 1. Introduction

The sine-Gordon equation arises in extended rectangular Josephson junctions, which consist of two layers of superconducting materials separated by an isolating barrier. A typical arrangement is a layer of lead and a layer of niobium separated by a layer of niobium oxide. A quantum particle has a nonzero significant probability of being able to penetrate in the other side of a potential barrier that would be impenetrable to the corresponding classical particle. Mathematical model of light bullets in Maxwell-Bloch system can also be described by the two-dimensional undamped sine-Gordon equation [1]. This equation is also applied in a large number of areas of physics, for example, crystal dislocation theory [2], self-induced transparency [2], laser physics [2], and particle physics [3, 4]. It is also a special case of the Bady Skyrme model which describes baryons in a nonlinear manner [5].Like some well-known partial differential equations such as KdV, mKdV, and nonlinear Schrödinger equations [6], the two-dimensional sine-Gordon equation possesses various types of soliton solutions (such as line solitons, elliptic ring soliton, and ring solitons). For the undamped (ρ
=
0) sine-Gordon equation in higher dimensions exact soliton solutions have been obtained in Hirota [7], in Zagrodzinsky [8] using Lamb’s method, in Leibbrandt [9] by Bläcklund transformation, and in Kaliappan and Lakshmanan [10] by Painlevé transcendents.Numerical solutions for two-dimensional undamped sine-Gordon equation have been given among others by Guo et al. [11] using two finite difference schemes, Xin [1] studying sine-Gordon equation as an asymptotic reduction of the two level dissipationless Maxwell-Bloch system, Christiansen and Lomdahl [12] using a generalized leapfrog method, Argyris et al. [13] using the finite element method. Sheng et al. [14] introduced a split cosine scheme, Bratsos [15] used a three-time level fourth-order explicit finite-difference scheme, Mirzaei and Dehghan [16] applied the continuous linear boundary elements method, Chen et al. [17] applied the multilevel augmentation method for solving the sine-Gordon equations, and so forth. Numerical approaches to the damped sine-Gordon equation can also be found in Nakajima et al. [18] who considered dimensionless loss factors and unitless normalized bias and Gorria et al. [19] who studied the nonlinear wave propagation in a planar wave guide consisting of two rectangular regions joined by a bent of constant curvature. Bratsos [15] and Djidjeli et al. [20] used a two-step one-parameter leap-frog scheme, which is a generalization to that used by Christiansen and Lomdahl [12]. Dehghan and Shokri [21] used the radial basis functions as a truly meshfree method, to solve the two-dimensional damped/undamped sine-Gordon equation.Consider the two-dimensional sine-Gordon equation as follows:(1)
∂
2
u
∂
t
2
+
ρ
∂
u
∂
t
=
∂
2
u
∂
x
2
+
∂
2
u
∂
y
2
-
ϕ
(
x
,
y
)
sin
u
,
with u
=
u
(
x
,
y
,
t
) in the region Ω
=
{
(
x
,
y
)
∈
[
L
x
0
,
L
x
1
]
×
[
L
y
0
,
L
y
1
]
} for t
>
0 and the parameter ρ being the so-called dissipative term, which is assumed to be a real number with ρ
≥
0. When ρ
=
0, (1) reduces to the undamped sine-Gordon equation in two space variables, whereas, when ρ
>
0 to the damped one. The function ϕ
(
x
,
y
) can be explained as a Josephson current density.Initial conditions associated with (1) will be assumed to be of the form
(2)
u
(
x
,
y
,
0
)
=
f
(
x
,
y
)
,
(
x
,
y
)
∈
Ω
,
with initial velocity
(3)
∂
u
∂
t
(
x
,
y
,
0
)
=
g
(
x
,
y
)
,
(
x
,
y
)
∈
Ω
.In (2) and (3), the functions f
(
x
,
y
) and g
(
x
,
y
) represent wave modes or kinks and velocity, respectively. Boundary conditions will be assumed to be of the form
(4)
∂
u
∂
x
(
L
x
0
,
y
,
t
)
=
∂
u
∂
x
(
L
x
1
,
y
,
t
)
=
p
(
L
x
1
,
y
,
t
)
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
L
y
0
<
y
<
L
y
1
,
t
>
0
,
∂
u
∂
y
(
x
,
L
y
0
,
t
)
=
∂
u
∂
y
(
x
,
L
y
1
,
t
)
=
q
(
x
,
L
y
1
,
t
)
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
L
x
0
<
x
<
L
x
1
,
t
>
0
,
where p
(
x
,
y
,
t
) and q
(
x
,
y
,
t
) are normal gradients along the boundary of the region Ω.If partial differential equation contains the second order term∂
2
u
/
∂
x
2, the equation is discredited by using the standard Crank-Nicolson scheme. However, it is verified that the scheme is unstable unconditionally for the heat equation ∂
u
(
x
,
t
)
/
∂
t
=
∂
2
u
(
x
,
t
)
/
∂
x
2 when it is discretized by the leap-frog scheme in time direction and the Crank-Nicolson scheme in space direction, see [22]. The Dufort-Frankel method is very similar to the leap-frog scheme but it has better numerical stability, and sometimes it will bring unpredicted effect (if Dufort-Frankel approximate scheme is applied in spatial direction, the heat equation is unconditionally stable [22]). Wu [23] applied the Dufort-Frankel scheme for solving the linear and nonlinear one-dimensional Schrödinger equations and obtained an unconditionally stable scheme. Markowich et al. [24] applied the Wigner-measure analysis for investigating the convergence of the Dufort-Frankel scheme for the Schrödinger equation in semiclassical regime. Lai et al. [25] used a simple Dufort-Frankel type scheme for solving the time-dependent Gross-Pitaevskii equation. In this paper, we will use the Dufort-Frankel scheme to solve two-dimensional sine-Gordon equation.The organization of the paper is as follows. In Section2, a Crank-Nicolson finite-difference (CNFD) scheme and a Dufort-Frankel finite-difference (DFFD) scheme for solving two-dimensional sine-Gordon equation are introduced. In Section 3, the stabilities of the two schemes are discussed. In Section 4, the error estimates of CNFD and DFFD schemes are proved. To avoid solving nonlinear equations, the predictor-corrector methods of the two schemes are proposed in Section 5. In Section 6, numerical results are investigated.


## 2. Numerical Methods

### 2.1. The CNFD Scheme and the DFFD Scheme

For the numerical solution the regionΩ
×
[
t
>
0
] with its boundary ∂
Ω consisting of the lines and t
=
0 being covered with a rectangular mesh, the points with coordinates (
x
,
y
,
t
)
=
(
x
k
,
y
m
,
t
n
)
=
(
L
x
0
+
k
h
x
,
L
y
0
+
m
h
y
,
n
τ
), Ω
h
=
{
(
x
k
,
y
m
)
} with k
,
m
=
0,1
,
…
,
N, and n
=
0,1
,
…. The h
x
=
(
L
x
1
-
L
x
0
)
/
N and h
y
=
(
L
y
1
-
L
y
0
)
/
N represent the space steps along x direction and y direction, while τ represents the time step. The solution of an approximating difference scheme at the same point will be denoted by u
k
,
m
n (k
,
m
=
0,1
,
…
,
N) for the purpose of analyzing stability, the numerical value actually obtained will be denoted by u
~
k
,
m
n.The CNFD scheme of (1) is
(5)
u
k
,
m
n
+
1
-
2
u
k
,
m
n
+
u
k
,
m
n
-
1
τ
2
+
ρ
u
k
,
m
n
+
1
-
u
k
,
m
n
-
1
2
τ
=
u
k
+
1
,
m
n
-
2
u
k
,
m
n
+
u
k
-
1
,
m
n
h
x
2
+
u
k
,
m
+
1
n
-
2
u
k
,
m
n
+
u
k
,
m
-
1
n
h
y
2
-
1
4
ϕ
k
,
m
(
sin
u
k
,
m
n
+
1
+
2
sin
u
k
,
m
n
+
sin
u
k
,
m
n
-
1
)
,
where ϕ
k
,
m
=
ϕ
(
x
k
,
y
m
), k
,
m
=
0,1
,
…
,
N.The DFFD scheme for (5) is obtained by replacing the term 2
u
k
,
m
n by u
k
,
m
n
+
1
+
u
k
,
m
n
-
1 in the CNFD scheme in the discretizations of ∂
u
(
t
,
x
,
y
)
/
∂
x
2 and ∂
u
(
t
,
x
,
y
)
/
∂
y
2, that is,
(6)
u
k
,
m
n
+
1
-
2
u
k
,
m
n
+
u
k
,
m
n
-
1
τ
2
+
ρ
u
k
,
m
n
+
1
-
u
k
,
m
n
-
1
2
τ
=
u
k
+
1
,
m
n
-
(
u
k
,
m
n
+
1
+
u
k
,
m
n
-
1
)
+
u
k
-
1
,
m
n
h
x
2
+
u
k
,
m
+
1
n
-
(
u
k
,
m
n
+
1
+
u
k
,
m
n
-
1
)
+
u
k
,
m
-
1
n
h
y
2
-
1
4
ϕ
k
,
m
(
sin
u
k
,
m
n
+
1
+
2
sin
u
k
,
m
n
+
sin
u
k
,
m
n
-
1
)
.
Note that
(7)
u
k
+
1
,
m
n
-
(
u
k
,
m
n
+
1
+
u
k
,
m
n
-
1
)
+
u
k
-
1
,
m
n
h
x
2
=
u
k
+
1
,
m
n
-
2
u
k
,
m
n
+
u
k
-
1
,
m
n
h
x
2
-
μ
x
2
u
k
,
m
n
+
1
-
2
u
k
,
m
n
+
u
k
,
m
n
-
1
τ
2
,
where μ
x
=
τ
/
h
x. The DFFD scheme can be expressed as follows, with μ
y
=
τ
/
h
y:
(8)
(
1
+
μ
x
2
+
μ
y
2
)
(
u
k
,
m
n
+
1
-
2
u
k
,
m
n
+
u
k
,
m
n
-
1
)
+
ρ
τ
2
(
u
k
,
m
n
+
1
-
u
k
,
m
n
-
1
)
=
μ
x
2
(
u
k
+
1
,
m
n
-
2
u
k
,
m
n
+
u
k
-
1
,
m
n
)
+
μ
y
2
(
u
k
,
m
+
1
n
-
2
u
k
,
m
n
+
u
k
,
m
-
1
n
)
-
τ
2
4
ϕ
k
,
m
(
sin
u
k
,
m
n
+
1
+
2
sin
u
k
,
m
n
+
sin
u
k
,
m
n
-
1
)
=
0
.


### 2.2. Local Truncation Error

By Tayor expansion, it can be easily obtained that for the smooth functionu
(
x
,
y
,
t
), the principal part of the local truncation error of the DFFD (8) is
(9)
(
1
+
μ
x
2
+
μ
y
2
)
×
u
(
x
k
,
y
m
,
t
n
+
1
)
-
2
u
(
x
k
,
y
m
,
t
n
)
+
u
(
x
k
,
y
m
,
t
n
-
1
)
τ
2
+
ρ
u
(
x
k
,
y
m
,
t
n
+
1
)
-
u
(
x
k
,
y
m
,
t
n
-
1
)
2
τ
-
u
(
x
k
+
1
,
y
m
,
t
n
)
-
2
u
(
x
k
,
y
m
,
t
n
)
+
u
(
x
k
-
1
,
y
m
,
t
n
)
h
x
2
-
u
(
x
k
,
y
m
+
1
,
t
n
)
-
2
u
(
x
k
,
y
m
,
t
n
)
+
u
(
x
k
,
y
m
-
1
,
t
n
)
h
y
2
=
∂
2
u
∂
t
2
(
x
k
,
y
m
,
t
n
)
+
ρ
∂
u
∂
t
(
x
k
,
y
m
,
t
n
)
-
∂
2
u
∂
x
2
(
x
k
,
y
m
,
t
n
)
-
∂
2
u
∂
y
2
(
x
k
,
y
m
,
t
n
)
+
O
(
τ
2
+
h
x
2
+
h
y
2
+
(
τ
h
x
)
2
+
(
τ
h
y
)
2
)
.To make (8) be consistent, it is known from (9) that the step sizes τ
,
h
x
,
h
y have to be limited such that when τ
→
0, h
x
→
0, and h
y
→
0, we have τ
/
h
x, τ
/
h
y
→
0. Letting τ
=
O
(
h
x
2
), the above requirement holds automatically.Similarly, it is easily obtained that the truncation error of the CNFD scheme isO
(
τ
2
+
h
x
2
+
h
y
2
).Remark 1.
In terms of the compatibility conditions, the constraint of time and space steps in the DFFD scheme is more demanding than that of the CNFD scheme, but it is easy for the time step to satisfy it, which is a usual method. In the following discussion of convergence, it is found that the time and space steps in the CNFD scheme have to satisfy the constraint; however, the DFFD scheme is unconditional.Remark 2.
The DFFD is a scheme in which a termμ
x
+
μ
y is added to the leading coefficient of DNFD. This extra term enhances the stability of it and greatly reduces the restriction on time step.


## 2.1. The CNFD Scheme and the DFFD Scheme

For the numerical solution the regionΩ
×
[
t
>
0
] with its boundary ∂
Ω consisting of the lines and t
=
0 being covered with a rectangular mesh, the points with coordinates (
x
,
y
,
t
)
=
(
x
k
,
y
m
,
t
n
)
=
(
L
x
0
+
k
h
x
,
L
y
0
+
m
h
y
,
n
τ
), Ω
h
=
{
(
x
k
,
y
m
)
} with k
,
m
=
0,1
,
…
,
N, and n
=
0,1
,
…. The h
x
=
(
L
x
1
-
L
x
0
)
/
N and h
y
=
(
L
y
1
-
L
y
0
)
/
N represent the space steps along x direction and y direction, while τ represents the time step. The solution of an approximating difference scheme at the same point will be denoted by u
k
,
m
n (k
,
m
=
0,1
,
…
,
N) for the purpose of analyzing stability, the numerical value actually obtained will be denoted by u
~
k
,
m
n.The CNFD scheme of (1) is
(5)
u
k
,
m
n
+
1
-
2
u
k
,
m
n
+
u
k
,
m
n
-
1
τ
2
+
ρ
u
k
,
m
n
+
1
-
u
k
,
m
n
-
1
2
τ
=
u
k
+
1
,
m
n
-
2
u
k
,
m
n
+
u
k
-
1
,
m
n
h
x
2
+
u
k
,
m
+
1
n
-
2
u
k
,
m
n
+
u
k
,
m
-
1
n
h
y
2
-
1
4
ϕ
k
,
m
(
sin
u
k
,
m
n
+
1
+
2
sin
u
k
,
m
n
+
sin
u
k
,
m
n
-
1
)
,
where ϕ
k
,
m
=
ϕ
(
x
k
,
y
m
), k
,
m
=
0,1
,
…
,
N.The DFFD scheme for (5) is obtained by replacing the term 2
u
k
,
m
n by u
k
,
m
n
+
1
+
u
k
,
m
n
-
1 in the CNFD scheme in the discretizations of ∂
u
(
t
,
x
,
y
)
/
∂
x
2 and ∂
u
(
t
,
x
,
y
)
/
∂
y
2, that is,
(6)
u
k
,
m
n
+
1
-
2
u
k
,
m
n
+
u
k
,
m
n
-
1
τ
2
+
ρ
u
k
,
m
n
+
1
-
u
k
,
m
n
-
1
2
τ
=
u
k
+
1
,
m
n
-
(
u
k
,
m
n
+
1
+
u
k
,
m
n
-
1
)
+
u
k
-
1
,
m
n
h
x
2
+
u
k
,
m
+
1
n
-
(
u
k
,
m
n
+
1
+
u
k
,
m
n
-
1
)
+
u
k
,
m
-
1
n
h
y
2
-
1
4
ϕ
k
,
m
(
sin
u
k
,
m
n
+
1
+
2
sin
u
k
,
m
n
+
sin
u
k
,
m
n
-
1
)
.
Note that
(7)
u
k
+
1
,
m
n
-
(
u
k
,
m
n
+
1
+
u
k
,
m
n
-
1
)
+
u
k
-
1
,
m
n
h
x
2
=
u
k
+
1
,
m
n
-
2
u
k
,
m
n
+
u
k
-
1
,
m
n
h
x
2
-
μ
x
2
u
k
,
m
n
+
1
-
2
u
k
,
m
n
+
u
k
,
m
n
-
1
τ
2
,
where μ
x
=
τ
/
h
x. The DFFD scheme can be expressed as follows, with μ
y
=
τ
/
h
y:
(8)
(
1
+
μ
x
2
+
μ
y
2
)
(
u
k
,
m
n
+
1
-
2
u
k
,
m
n
+
u
k
,
m
n
-
1
)
+
ρ
τ
2
(
u
k
,
m
n
+
1
-
u
k
,
m
n
-
1
)
=
μ
x
2
(
u
k
+
1
,
m
n
-
2
u
k
,
m
n
+
u
k
-
1
,
m
n
)
+
μ
y
2
(
u
k
,
m
+
1
n
-
2
u
k
,
m
n
+
u
k
,
m
-
1
n
)
-
τ
2
4
ϕ
k
,
m
(
sin
u
k
,
m
n
+
1
+
2
sin
u
k
,
m
n
+
sin
u
k
,
m
n
-
1
)
=
0
.


## 2.2. Local Truncation Error

By Tayor expansion, it can be easily obtained that for the smooth functionu
(
x
,
y
,
t
), the principal part of the local truncation error of the DFFD (8) is
(9)
(
1
+
μ
x
2
+
μ
y
2
)
×
u
(
x
k
,
y
m
,
t
n
+
1
)
-
2
u
(
x
k
,
y
m
,
t
n
)
+
u
(
x
k
,
y
m
,
t
n
-
1
)
τ
2
+
ρ
u
(
x
k
,
y
m
,
t
n
+
1
)
-
u
(
x
k
,
y
m
,
t
n
-
1
)
2
τ
-
u
(
x
k
+
1
,
y
m
,
t
n
)
-
2
u
(
x
k
,
y
m
,
t
n
)
+
u
(
x
k
-
1
,
y
m
,
t
n
)
h
x
2
-
u
(
x
k
,
y
m
+
1
,
t
n
)
-
2
u
(
x
k
,
y
m
,
t
n
)
+
u
(
x
k
,
y
m
-
1
,
t
n
)
h
y
2
=
∂
2
u
∂
t
2
(
x
k
,
y
m
,
t
n
)
+
ρ
∂
u
∂
t
(
x
k
,
y
m
,
t
n
)
-
∂
2
u
∂
x
2
(
x
k
,
y
m
,
t
n
)
-
∂
2
u
∂
y
2
(
x
k
,
y
m
,
t
n
)
+
O
(
τ
2
+
h
x
2
+
h
y
2
+
(
τ
h
x
)
2
+
(
τ
h
y
)
2
)
.To make (8) be consistent, it is known from (9) that the step sizes τ
,
h
x
,
h
y have to be limited such that when τ
→
0, h
x
→
0, and h
y
→
0, we have τ
/
h
x, τ
/
h
y
→
0. Letting τ
=
O
(
h
x
2
), the above requirement holds automatically.Similarly, it is easily obtained that the truncation error of the CNFD scheme isO
(
τ
2
+
h
x
2
+
h
y
2
).Remark 1.
In terms of the compatibility conditions, the constraint of time and space steps in the DFFD scheme is more demanding than that of the CNFD scheme, but it is easy for the time step to satisfy it, which is a usual method. In the following discussion of convergence, it is found that the time and space steps in the CNFD scheme have to satisfy the constraint; however, the DFFD scheme is unconditional.Remark 2.
The DFFD is a scheme in which a termμ
x
+
μ
y is added to the leading coefficient of DNFD. This extra term enhances the stability of it and greatly reduces the restriction on time step.


## 3. Stability Analysis of CNFD and DFFD

In this section, the stabilities of the two schemes are discussed. Following the Fourier method of analyzing stability a small error of the following form is considered:(10)
Q
k
,
m
n
=
u
k
,
m
n
-
u
~
k
,
m
n
,
with
(11)
Q
k
,
m
n
=
λ
n
e
i
(
β
k
h
x
+
γ
m
h
y
)
,
i
=
-
1
,
where λ is a complex number and β
,
γ are real. Due to von Neumann criterion for stability, the condition |
λ
|
≤
1 has to be satisfied.


### 3.1. Stability Analysis for the DFFD

Using (10)-(11) and Maclaurin’s expansion
(12)
sin
u
k
,
m
n
=
∑
j
=
0
+
∞
‍
(
-
1
)
j
(
u
k
,
m
n
)
2
j
+
1
(
2
j
+
1
)
!
,
(9) can be written as
(13)
(
1
+
μ
x
2
+
μ
y
2
+
ρ
τ
2
)
Q
k
,
m
n
+
1
+
τ
2
4
ϕ
k
,
m
S
k
,
m
n
+
1
Q
k
,
m
n
+
1
-
2
(
1
-
2
μ
x
2
sin
2
β
h
x
2
-
2
μ
y
2
sin
2
γ
h
y
2
)
Q
k
,
m
n
+
τ
2
2
ϕ
k
,
m
S
k
,
m
n
Q
k
,
m
n
+
(
1
+
μ
x
2
+
μ
y
2
-
ρ
τ
2
)
Q
k
,
m
n
-
1
+
τ
2
4
ϕ
k
,
m
S
k
,
m
n
-
1
Q
k
,
m
n
-
1
=
0
,
where
(14)
S
k
,
m
n
=
∑
j
=
0
+
∞
‍
(
-
1
)
j
(
2
j
+
1
)
!
iiiii
×
[
(
u
k
,
m
n
)
2
j
+
(
u
k
,
m
n
)
2
j
-
1
u
~
k
,
m
n
+
⋯
+
(
u
~
k
,
m
n
)
2
j
]
≈
∑
j
=
0
+
∞
‍
(
-
1
)
j
(
2
j
)
!
(
u
s
)
2
j
=
cos
⁡
u
s
.
Here we linearize the term in square brackets and u
s
=
max
⁡
k
,
m
=
0,1
,
…
,
N
-
1
u
k
,
m
0. Using (14) and Q
k
,
m
n
=
λ
n
e
i
(
β
k
h
x
+
γ
m
h
y
), we get the following stability equation:
(15)
A
λ
2
-
2
B
λ
+
C
=
0
,
where
(16)
A
=
1
+
μ
x
2
+
μ
y
2
+
τ
2
4
ϕ
k
,
m
cos
⁡
u
s
+
ρ
τ
2
,
B
=
1
-
2
μ
x
2
sin
2
β
h
x
2
-
2
μ
y
2
sin
2
γ
h
y
2
-
τ
2
4
ϕ
k
,
m
cos
⁡
u
s
,
C
=
1
+
μ
x
2
+
μ
y
2
+
τ
2
4
ϕ
k
,
m
cos
⁡
u
s
-
ρ
τ
2
.Assume thatcos
⁡
u
s
>
0, if ρ
=
0, A
>
0, C
>
0, then B
2
-
A
C
≤
0 holds for any μ
x
,
μ
y. Hence |
λ
|
≤
1 holds. Through von Neumann criterion the DFFD scheme is unconditionally stable.Assume thatρ
>
0, 1
-
ρ
τ
/
2
>
0, and A
>
0, C
>
0 and A
>
C. Note that B
2
-
A
C
≤
0 always holds. Hence the solutions,
(17)
λ
±
=
(
B
±
B
2
-
A
C
)
A
,
of (15) satisfy |
λ
±
|
≤
1, which implies that if ρ
>
0, τ
<
2
/
ρ, the scheme is stable.


### 3.2.  Stability Analysis for the CNFD

Similarly, the stability equation of CNFD (5) is
(18)
A
˘
λ
2
-
2
B
˘
λ
+
C
˘
=
0
,
where
(19)
A
˘
=
1
+
τ
2
4
ϕ
k
,
m
cos
⁡
u
s
+
ρ
τ
2
,
B
˘
=
1
-
2
μ
x
2
sin
2
β
h
x
2
-
2
μ
y
2
sin
2
γ
h
y
2
-
τ
2
4
ϕ
k
,
m
cos
⁡
u
s
,
C
˘
=
1
+
τ
2
4
ϕ
k
,
m
cos
⁡
u
s
-
ρ
τ
2
.To guarantee|
λ
|
≤
1, we need |
B
˘
|
≤
C
˘ or
(20)
-
C
˘
≤
B
˘
≤
C
˘
.The left-hand side of (20) has the form
(21)
(
1
h
x
2
sin
2
β
h
x
2
+
1
h
y
2
sin
2
γ
h
y
2
)
τ
2
+
ρ
4
τ
-
1
≤
0
,
which is true if
(22)
(
1
h
x
2
+
1
h
y
2
)
τ
2
+
ρ
4
τ
-
1
≤
0
.Letτ
~
1
,
τ
~
2 be the roots of (22) and H
=
1
/
h
x
2
+
1
/
h
y
2. Note that τ
~
1
τ
~
2
=
-
1
/
H
<
0; the roots are real and distinct. Let the positive root τ
~
2
=
(
ρ
2
+
64
H
-
ρ
)
/
8
H; we then have τ
≤
τ
~
2.The right-hand side of (20) gives
(23)
τ
2
(
ϕ
k
,
m
cos
⁡
u
s
+
4
(
1
h
x
2
sin
2
β
h
x
2
+
1
h
y
2
sin
2
γ
h
y
2
)
)
≥
τ
ρ
;
if ρ
=
0, then (23) is always satisfied, while when ρ
>
0, τ must satisfy the restriction condition
(24)
τ
≥
ρ
[
ϕ
~
+
4
(
h
x
-
2
+
h
y
-
2
)
]
-
1
.
Thus, the time step τ needs to satisfy the stability conditions (22) and (24).


## 3.1. Stability Analysis for the DFFD

Using (10)-(11) and Maclaurin’s expansion
(12)
sin
u
k
,
m
n
=
∑
j
=
0
+
∞
‍
(
-
1
)
j
(
u
k
,
m
n
)
2
j
+
1
(
2
j
+
1
)
!
,
(9) can be written as
(13)
(
1
+
μ
x
2
+
μ
y
2
+
ρ
τ
2
)
Q
k
,
m
n
+
1
+
τ
2
4
ϕ
k
,
m
S
k
,
m
n
+
1
Q
k
,
m
n
+
1
-
2
(
1
-
2
μ
x
2
sin
2
β
h
x
2
-
2
μ
y
2
sin
2
γ
h
y
2
)
Q
k
,
m
n
+
τ
2
2
ϕ
k
,
m
S
k
,
m
n
Q
k
,
m
n
+
(
1
+
μ
x
2
+
μ
y
2
-
ρ
τ
2
)
Q
k
,
m
n
-
1
+
τ
2
4
ϕ
k
,
m
S
k
,
m
n
-
1
Q
k
,
m
n
-
1
=
0
,
where
(14)
S
k
,
m
n
=
∑
j
=
0
+
∞
‍
(
-
1
)
j
(
2
j
+
1
)
!
iiiii
×
[
(
u
k
,
m
n
)
2
j
+
(
u
k
,
m
n
)
2
j
-
1
u
~
k
,
m
n
+
⋯
+
(
u
~
k
,
m
n
)
2
j
]
≈
∑
j
=
0
+
∞
‍
(
-
1
)
j
(
2
j
)
!
(
u
s
)
2
j
=
cos
⁡
u
s
.
Here we linearize the term in square brackets and u
s
=
max
⁡
k
,
m
=
0,1
,
…
,
N
-
1
u
k
,
m
0. Using (14) and Q
k
,
m
n
=
λ
n
e
i
(
β
k
h
x
+
γ
m
h
y
), we get the following stability equation:
(15)
A
λ
2
-
2
B
λ
+
C
=
0
,
where
(16)
A
=
1
+
μ
x
2
+
μ
y
2
+
τ
2
4
ϕ
k
,
m
cos
⁡
u
s
+
ρ
τ
2
,
B
=
1
-
2
μ
x
2
sin
2
β
h
x
2
-
2
μ
y
2
sin
2
γ
h
y
2
-
τ
2
4
ϕ
k
,
m
cos
⁡
u
s
,
C
=
1
+
μ
x
2
+
μ
y
2
+
τ
2
4
ϕ
k
,
m
cos
⁡
u
s
-
ρ
τ
2
.Assume thatcos
⁡
u
s
>
0, if ρ
=
0, A
>
0, C
>
0, then B
2
-
A
C
≤
0 holds for any μ
x
,
μ
y. Hence |
λ
|
≤
1 holds. Through von Neumann criterion the DFFD scheme is unconditionally stable.Assume thatρ
>
0, 1
-
ρ
τ
/
2
>
0, and A
>
0, C
>
0 and A
>
C. Note that B
2
-
A
C
≤
0 always holds. Hence the solutions,
(17)
λ
±
=
(
B
±
B
2
-
A
C
)
A
,
of (15) satisfy |
λ
±
|
≤
1, which implies that if ρ
>
0, τ
<
2
/
ρ, the scheme is stable.


## 3.2.  Stability Analysis for the CNFD

Similarly, the stability equation of CNFD (5) is
(18)
A
˘
λ
2
-
2
B
˘
λ
+
C
˘
=
0
,
where
(19)
A
˘
=
1
+
τ
2
4
ϕ
k
,
m
cos
⁡
u
s
+
ρ
τ
2
,
B
˘
=
1
-
2
μ
x
2
sin
2
β
h
x
2
-
2
μ
y
2
sin
2
γ
h
y
2
-
τ
2
4
ϕ
k
,
m
cos
⁡
u
s
,
C
˘
=
1
+
τ
2
4
ϕ
k
,
m
cos
⁡
u
s
-
ρ
τ
2
.To guarantee|
λ
|
≤
1, we need |
B
˘
|
≤
C
˘ or
(20)
-
C
˘
≤
B
˘
≤
C
˘
.The left-hand side of (20) has the form
(21)
(
1
h
x
2
sin
2
β
h
x
2
+
1
h
y
2
sin
2
γ
h
y
2
)
τ
2
+
ρ
4
τ
-
1
≤
0
,
which is true if
(22)
(
1
h
x
2
+
1
h
y
2
)
τ
2
+
ρ
4
τ
-
1
≤
0
.Letτ
~
1
,
τ
~
2 be the roots of (22) and H
=
1
/
h
x
2
+
1
/
h
y
2. Note that τ
~
1
τ
~
2
=
-
1
/
H
<
0; the roots are real and distinct. Let the positive root τ
~
2
=
(
ρ
2
+
64
H
-
ρ
)
/
8
H; we then have τ
≤
τ
~
2.The right-hand side of (20) gives
(23)
τ
2
(
ϕ
k
,
m
cos
⁡
u
s
+
4
(
1
h
x
2
sin
2
β
h
x
2
+
1
h
y
2
sin
2
γ
h
y
2
)
)
≥
τ
ρ
;
if ρ
=
0, then (23) is always satisfied, while when ρ
>
0, τ must satisfy the restriction condition
(24)
τ
≥
ρ
[
ϕ
~
+
4
(
h
x
-
2
+
h
y
-
2
)
]
-
1
.
Thus, the time step τ needs to satisfy the stability conditions (22) and (24).


## 4. Convergence and Error Estimates for DFFD and CNFD

We define a discrete inner product and its associated norm by(25)
(
u
n
,
v
n
)
=
h
x
h
y
∑
k
=
0
N
‍
∑
m
=
0
N
‍
u
k
,
m
n
v
k
,
m
,
∥
u
n
∥
2
=
(
u
n
,
u
)
,
∥
u
n
∥
∞
=
max
⁡
Ω
h
|
u
k
,
m
n
|
.
Conveniently, we note that
(26)
δ
t
u
k
,
m
n
=
u
k
,
m
n
+
1
-
u
k
,
m
n
τ
,
δ
t
¯
u
k
,
m
n
=
u
k
,
m
n
-
u
k
,
m
n
-
1
τ
,
δ
t
^
u
k
,
m
n
=
u
k
,
m
n
+
1
-
u
k
,
m
n
-
1
τ
,
δ
t
2
u
k
,
m
n
=
δ
t
u
k
,
m
n
-
δ
t
¯
u
k
,
m
n
τ
,
δ
x
u
k
,
m
=
u
k
+
1
,
m
n
-
u
k
,
m
n
h
x
,
δ
y
u
k
,
m
n
=
u
k
,
m
+
1
n
-
u
k
,
m
n
h
y
,
δ
x
2
u
k
,
m
n
=
u
k
+
1
,
m
n
-
2
u
k
,
m
n
+
u
k
-
1
,
m
n
h
x
2
,
δ
y
2
u
k
,
m
n
=
u
k
,
m
+
1
n
-
2
u
k
,
m
n
+
u
k
,
m
-
1
n
h
y
2
,
δ
x
p
u
k
,
m
n
=
δ
x
p
-
1
u
k
+
1
,
m
n
-
δ
x
p
-
1
u
k
,
m
n
h
x
(
1
≤
p
≤
+
∞
)
.The difference scheme DFFD (8) can be written in the following form:
(27)
(
1
+
μ
x
2
+
μ
y
2
)
δ
t
2
u
k
,
m
n
+
ρ
δ
t
^
u
k
,
m
n
=
δ
x
2
u
k
,
m
n
+
δ
y
2
u
k
,
m
n
-
1
4
ϕ
k
,
m
(
sin
u
k
,
m
n
+
1
+
2
sin
u
k
,
m
n
+
sin
u
k
,
m
n
-
1
)
.We have the following Lemmas. ConsiderLemma 3.
Consider(28)
(
δ
t
2
u
n
,
δ
t
^
u
n
)
=
1
2
δ
t
¯
(
∥
δ
t
u
n
∥
2
)
.Lemma 4.
Consider(29)
(
δ
x
2
u
n
,
δ
t
^
u
n
)
=
-
1
2
δ
t
^
(
∥
δ
x
u
n
∥
2
)
+
τ
2
4
δ
t
¯
(
∥
δ
x
δ
t
u
n
∥
2
)
.Lemma 5.
Lettingp
≥
p
′
+
1, for arbitrary ɛ
>
0, we have, for ∀
u
n
∈
Ω
h,
(30)
τ
2
2
∥
δ
x
p
δ
t
u
n
∥
2
≤
ɛ
(
∥
δ
x
p
u
n
∥
2
+
∥
δ
x
p
u
n
+
1
∥
2
)
+
τ
2
8
ɛ
(
4
h
x
2
)
p
′
∥
δ
x
p
-
p
′
δ
t
u
∥
2
.Proof.
Byɛ inequality,
(31)
τ
2
2
∥
δ
x
p
δ
t
u
n
∥
2
≤
ɛ
τ
2
2
∥
δ
x
p
δ
t
u
n
∥
2
+
τ
2
8
ɛ
∥
δ
x
p
δ
t
u
n
∥
2
,
the definition of the difference quotient implies
(32)
ɛ
τ
2
2
∥
δ
x
p
δ
t
u
n
∥
2
≤
ɛ
(
∥
δ
x
p
u
n
∥
2
+
∥
δ
x
p
u
n
+
1
∥
2
)
.
By using the definition of norm, we get
(33)
τ
2
8
ɛ
∥
δ
x
p
δ
t
u
n
∥
2
≤
τ
2
8
ɛ
4
h
x
2
∥
δ
x
p
-
1
δ
t
u
n
∥
2
≤
τ
2
8
ɛ
(
4
h
x
2
)
2
∥
δ
x
p
-
2
δ
t
u
n
∥
2
≤
τ
2
8
ɛ
(
4
h
x
2
)
p
′
∥
δ
x
p
-
p
′
δ
t
u
∥
2
.
Formulas (31)-(32) imply that Lemma 5 holds. In particular,
(34)
τ
2
2
(
∥
δ
x
δ
t
u
n
∥
2
+
∥
δ
y
δ
t
u
n
∥
2
)
≤
ɛ
(
{
∥
δ
y
u
n
∥
2
}
∥
δ
x
u
n
+
1
∥
2
+
∥
δ
x
u
n
∥
2
+
∥
δ
y
u
n
+
1
∥
2
+
∥
δ
y
u
n
∥
2
)
+
1
2
ɛ
(
μ
x
2
+
μ
y
2
)
∥
δ
t
u
n
∥
2
.
The proof of Lemma 5 is complete.Next we consider the convergence of the DFFD (8). Suppose e
k
,
m
n
=
u
(
x
k
,
y
m
,
t
n
)
-
u
k
,
m
n, we define the following weak form, for ∀
φ
∈
H
1
(
Ω
):
(35)
(
1
+
μ
x
2
+
μ
y
2
)
(
δ
t
2
u
n
,
φ
)
+
ρ
(
δ
t
^
u
n
,
φ
)
=
-
(
δ
x
u
n
,
δ
x
φ
)
-
(
δ
y
u
n
,
δ
y
φ
)
-
1
4
(
ϕ
(
sin
u
n
+
1
+
2
sin
u
n
+
sin
n
-
1
)
,
φ
)
.
Subtracting (8) from (27), we have
(36)
(
1
+
μ
x
2
+
μ
y
2
)
(
δ
t
2
e
n
,
φ
)
+
ρ
(
δ
t
^
e
n
,
φ
)
+
(
e
n
,
φ
)
+
(
δ
x
e
n
,
δ
x
φ
)
+
(
δ
y
e
n
,
δ
y
φ
)
=
(
e
n
,
φ
)
-
1
4
(
q
(
u
(
t
n
)
)
-
q
(
u
n
)
,
φ
)
+
(
r
n
,
φ
)
,
where
(37)
q
(
u
(
t
n
)
)
-
q
(
u
n
)
=
ϕ
[
sin
u
(
t
n
+
1
)
-
sin
u
n
+
1
+
2
(
sin
u
(
t
n
)
-
sin
u
n
)
iiiiiiiiiiiii
(
u
n
+
1
)
+
(
sin
u
(
t
n
-
1
)
-
sin
u
n
-
1
)
]
,
r
n
=
O
(
τ
2
+
h
x
2
+
h
y
2
+
(
τ
h
x
)
2
+
(
τ
h
y
)
2
)
.
In (36), letting φ
=
δ
t
^
e
n,
(38)
(
1
+
μ
x
2
+
μ
y
2
)
(
δ
t
2
e
n
,
δ
t
^
e
n
)
+
ρ
(
δ
t
^
e
n
,
δ
t
^
e
n
)
+
(
e
n
,
δ
t
^
e
n
)
(
δ
x
e
n
,
δ
t
^
δ
x
e
n
)
+
(
δ
y
e
n
,
δ
t
^
δ
y
e
n
)
=
(
e
n
,
δ
t
^
e
n
)
-
1
4
(
q
(
u
(
t
n
)
-
q
(
u
n
)
)
,
δ
t
^
e
n
)
+
(
r
n
,
δ
t
^
e
n
)
.
Applying Lemmas 3 and 4, we have
(39)
1
2
(
1
+
μ
x
2
+
μ
y
2
)
δ
t
¯
(
∥
δ
t
e
n
∥
2
)
+
ρ
∥
δ
t
^
e
n
∥
2
+
1
2
δ
t
^
(
∥
e
n
∥
2
)
-
τ
2
4
δ
t
¯
(
∥
δ
t
e
n
∥
2
)
+
1
2
δ
t
^
(
∥
δ
x
e
n
∥
2
)
-
τ
2
4
δ
t
¯
(
∥
δ
x
δ
t
e
n
∥
2
)
+
1
2
δ
t
^
(
∥
δ
y
e
n
∥
2
)
-
τ
2
4
δ
t
¯
(
∥
δ
y
δ
t
e
n
∥
2
)
=
(
r
n
,
δ
t
^
e
n
)
-
(
q
(
u
(
t
n
)
)
-
q
(
u
n
)
,
δ
t
^
e
n
)
+
(
e
n
,
δ
t
^
e
n
)
.
Summing up for n from 1 to M and multiplying 2
τ, we get
(40)
Θ
M
+
2
ρ
τ
∑
n
=
1
M
‍
∥
δ
t
^
e
n
∥
2
=
Θ
0
+
2
τ
∑
n
=
1
M
‍
[
(
r
n
,
δ
t
^
e
n
)
-
(
q
(
u
(
t
n
)
)
-
q
(
u
n
)
,
δ
t
^
e
n
)
iiiiiiiiiiiiiiiiiiiiiiiiiii
+
(
e
n
,
δ
t
^
e
n
)
]
=
Θ
0
+
2
τ
∑
n
=
1
M
‍
[
I
1
+
I
2
+
I
3
]
,
where
(41)
Θ
M
=
(
1
+
μ
x
2
+
μ
y
2
)
∥
δ
t
e
M
∥
2
+
1
2
(
∥
e
M
∥
2
+
∥
e
M
+
1
∥
2
)
-
τ
2
2
∥
δ
t
e
M
∥
2
+
1
2
(
∥
δ
x
e
M
∥
+
∥
δ
x
e
M
+
1
∥
2
)
+
1
2
(
∥
δ
y
e
M
∥
2
+
∥
δ
y
e
M
+
1
∥
2
)
-
τ
2
2
(
∥
δ
t
δ
x
e
M
∥
2
+
∥
δ
t
δ
y
e
M
∥
2
)
,
Θ
0
=
(
1
+
μ
x
2
+
μ
y
2
)
∥
δ
t
e
0
∥
2
+
1
2
(
∥
e
0
∥
2
+
∥
e
1
∥
2
)
-
τ
2
2
∥
δ
t
e
0
∥
2
+
1
2
(
∥
δ
x
e
0
∥
+
∥
δ
x
e
1
∥
2
)
+
1
2
(
∥
δ
y
e
0
∥
2
+
∥
δ
y
e
1
∥
2
)
-
τ
2
2
(
∥
δ
x
δ
t
e
0
∥
2
+
∥
δ
y
δ
t
e
0
∥
2
)
.
Note that
(42)
I
1
=
(
r
n
,
δ
t
^
e
n
)
≤
C
∥
τ
2
+
h
x
2
+
h
y
2
+
(
τ
h
x
)
2
+
(
τ
h
y
)
2
∥
∥
δ
t
^
e
n
∥
≤
C
(
τ
4
+
h
x
4
+
h
y
4
+
(
τ
h
x
)
4
+
(
τ
h
y
)
4
iiiii
(
τ
h
y
)
4
+
∥
δ
t
e
n
-
1
∥
2
+
∥
δ
t
e
n
∥
2
)
,
which follows from
(43)
|
(
ϕ
(
sin
u
(
t
n
+
1
)
-
sin
u
n
+
1
)
,
δ
t
^
e
n
)
|
≤
C
∥
u
(
t
n
+
1
)
-
u
n
+
1
∥
∥
δ
t
^
e
n
∥
≤
C
∥
e
n
+
1
∥
∥
δ
t
^
e
n
∥
≤
1
4
∥
e
n
+
1
∥
2
+
C
(
∥
δ
t
e
n
-
1
∥
2
+
∥
δ
t
e
n
∥
2
)
,
|
(
ϕ
[
2
(
sin
u
(
t
n
)
-
sin
u
n
)
(
u
n
+
1
)
iiiiiii
+
(
sin
u
(
t
n
-
1
)
-
sin
u
n
-
1
)
]
,
δ
t
^
e
n
)
|
≤
C
(
∥
e
n
∥
2
+
∥
e
n
-
1
∥
2
+
∥
δ
t
e
n
-
1
∥
2
+
∥
δ
t
e
n
∥
2
)
.
Similarly
(44)
I
2
=
(
q
(
u
(
t
n
)
)
-
q
(
u
n
)
,
δ
t
^
e
n
)
≤
1
4
∥
e
n
+
1
∥
2
+
C
(
∥
e
n
∥
2
+
∥
e
n
-
1
∥
2
+
∥
δ
t
e
n
-
1
∥
2
+
∥
δ
t
e
n
∥
2
)
,
I
3
=
(
e
n
,
δ
t
^
e
n
)
≤
C
(
∥
e
n
∥
2
+
∥
δ
t
e
n
-
1
∥
2
+
∥
δ
t
e
n
∥
2
)
.
From (42) and (44), we obtain
(45)
Θ
M
≤
Θ
0
+
2
C
τ
∑
n
=
1
M
‍
(
τ
4
+
h
x
4
+
h
y
4
+
(
τ
h
x
)
4
+
(
τ
h
y
)
4
iiiiiiiiiiiiiiiiiiiiii
(
τ
h
y
)
4
+
∥
e
n
∥
2
+
∥
δ
t
e
n
∥
2
)
.
By Lemma 5, we have
(46)
Θ
M
≥
[
1
+
μ
x
2
+
μ
y
2
-
τ
2
2
-
ɛ
4
(
μ
x
2
+
μ
y
2
)
]
∥
δ
t
e
n
∥
2
+
1
2
(
∥
e
M
+
1
∥
2
+
∥
e
M
∥
2
)
+
1
2
(
1
-
1
4
ɛ
)
(
∥
δ
x
e
M
+
1
∥
2
+
∥
δ
x
e
M
∥
2
iiiiiiiiiiiiiiiiiiiiiiii
+
∥
δ
y
e
M
+
1
∥
2
+
∥
δ
y
e
M
∥
2
)
,
(47)
Θ
0
≤
∥
δ
t
e
0
∥
2
+
1
2
(
∥
e
1
∥
2
+
∥
e
0
∥
2
)
+
1
2
(
∥
δ
x
e
0
∥
2
+
∥
δ
x
e
1
∥
2
+
∥
δ
y
e
1
∥
2
+
∥
δ
y
e
0
∥
2
)
;
in the last formula, we need to choose suitable ɛ to satisfy 1
-
1
/
4
ɛ
>
0. If the time step satisfies τ
<
2, then we will make all coefficients bigger than 0 and obtain
(48)
∥
e
M
∥
2
+
∥
δ
x
e
M
∥
2
+
∥
δ
y
e
M
∥
2
+
∥
δ
t
e
M
∥
2
≤
C
τ
∑
n
=
1
M
‍
(
τ
4
+
h
x
4
+
h
y
4
+
(
τ
2
h
x
)
4
+
(
τ
2
h
y
)
4
iiiiiiiiiiiiiiiiiiii
(
τ
2
h
y
)
4
+
∥
e
n
∥
2
+
∥
δ
x
e
n
∥
2
+
∥
δ
y
e
n
∥
2
+
∥
δ
t
e
n
∥
2
)
.By Gronwall’s inequality, we have, for allM,
(49)
∥
e
M
∥
+
∥
δ
x
e
M
∥
+
∥
δ
y
e
M
∥
+
∥
δ
t
e
M
∥
≤
C
(
τ
2
+
h
x
2
+
h
y
2
+
(
τ
h
x
)
2
+
(
τ
h
y
)
2
)
.Remark 6.
The convergence analyses for the CNFD scheme are the same as for the DFFD scheme. The difference is that there is noμ
x
2
+
μ
y
2 in the first term of CNFD scheme (5); so in (46), the convergence of the CNFD scheme needs is satisfying both 1
-
τ
2
/
2
-
(
ɛ
/
4
)
(
μ
x
2
+
μ
y
2
)
≥
σ
>
0 and 1
-
1
/
4
ɛ
>
0 for ɛ. We choose ɛ
=
1
/
2, the convergence of CDFD needs satisfying the constraint condition
(50)
τ
≤
(
1
2
+
1
8
(
1
h
x
2
+
1
h
y
2
)
)
-
1
/
2
.


## 5. The Predictor-Corrector Scheme

To avoid solving the nonlinear system arising from systems (5) and (8), the following predictor-corrector (P-C) scheme is used.


### 5.1. The Predictor of DFFD

Using an analogous scheme as in [26, 27], the predictor value u
^
k
,
m
n
+
1 was evaluated from the following three-time level explicit scheme of DFFD (8):
(51)
(
1
+
μ
x
2
+
μ
y
2
)
u
^
k
,
m
n
+
1
=
2
(
1
+
μ
x
2
+
μ
y
2
)
u
k
,
m
n
-
(
1
+
μ
x
2
+
μ
y
2
)
u
k
,
m
n
-
1
+
μ
x
2
(
u
k
+
1
,
m
n
-
2
u
k
,
m
n
+
u
k
-
1
,
m
n
)
+
μ
y
2
(
u
k
,
m
+
1
n
-
2
u
k
,
m
n
+
u
k
,
m
-
1
n
)
-
τ
2
ϕ
k
,
m
sin
u
k
,
m
n
,
for k
,
m
=
0,1
,
…
,
N.Following a similar approach for the stability analysis of the nonlinear scheme as in Section3.1, it can be proved that the characteristic equation of the predictor is given by
(52)
(
1
+
μ
x
2
+
μ
y
2
)
λ
2
-
2
(
1
-
2
(
μ
x
2
sin
2
β
k
h
x
2
+
μ
y
2
sin
2
γ
m
h
y
2
)
iiiiiiiiiiiiiiii
(
β
k
h
y
2
)
-
1
2
ϕ
k
,
m
cos
⁡
u
s
)
λ
+
1
+
μ
x
2
+
μ
y
2
=
0
.
And the scheme (51) is therefore unconditionally stable.


### 5.2. The Predictor of CNFD

The predictor of CNFD satisfies(53)
u
^
k
,
m
n
+
1
=
2
u
k
,
m
n
-
u
k
,
m
n
-
1
+
μ
x
2
(
u
k
+
1
,
m
n
-
2
u
k
,
m
n
+
u
k
-
1
,
m
n
)
+
μ
y
2
(
u
k
,
m
+
1
n
-
2
u
k
,
m
n
+
u
k
,
m
-
1
n
)
-
τ
2
ϕ
k
,
m
sin
u
k
,
m
n
.
Following the discussion in Section 3.2, the characteristic equation of (53) is
(54)
λ
2
-
2
(
1
-
2
(
μ
x
2
sin
2
β
k
h
x
2
+
μ
y
2
sin
2
γ
m
h
y
2
)
iiiiiiiiiiiiii
(
β
k
h
y
2
)
-
1
2
ϕ
k
,
m
cos
⁡
u
s
)
λ
+
1
=
0
.
Denote B
*
=
1
-
2
(
μ
x
2
sin
2
(
β
k
h
x
/
2
)
+
μ
y
2
sin
2
(
γ
m
h
y
/
2
)
)
-
(
1
/
2
)
ϕ
k
,
m
cos
⁡
u
s, the roots of (54) are
(55)
λ
±
=
B
*
±
(
B
*
)
2
-
1
.Letϕ
~
=
sup
⁡
(
x
,
y
)
∈
[
L
x
0
,
L
x
1
]
×
[
L
y
0
,
L
y
1
]
|
ϕ
(
x
,
y
)
| with |
ϕ
(
x
,
y
)
|
<
+
∞. To ensure |
λ
±
|
≤
1, we require |
B
*
|
≤
1, that is,
(56)
τ
≤
(
1
4
ϕ
~
+
1
h
x
2
+
1
h
y
2
)
-
1
/
2
.


### 5.3. The Predictor Algorithm Implementation for the DFFD Scheme

The predictor of the DFFD scheme is a three-time level explicit scheme; in order to obtain the same second order accuracy, we deal with the initial value of the DFFD scheme, that is, atn
=
0, (51) has the following form:
(57)
u
^
k
,
m
1
=
2
u
k
,
m
0
-
u
k
,
m
-
1
+
μ
x
2
(
u
k
+
1
,
m
0
-
2
u
k
,
m
0
+
u
k
-
1
,
m
0
)
+
μ
y
2
(
u
k
,
m
+
1
0
-
2
u
k
,
m
0
+
u
k
,
m
-
1
0
)
-
τ
2
ϕ
k
,
m
sin
u
k
,
m
0
(
k
,
m
=
1,2
,
…
,
N
)
;
to approximate u
-
1
(
x
,
y
) in the internal points the initial velocity is used, for this we discretize the initial velocity as
(58)
u
1
(
x
k
,
y
m
)
-
u
-
1
(
x
k
,
y
m
)
2
Δ
t
=
g
(
x
k
,
y
m
)
.
Combining (57) with (58), we have
(59)
u
^
k
,
m
1
=
u
k
,
m
0
+
2
τ
g
(
x
k
,
y
m
)
+
μ
x
2
2
(
u
k
+
1
,
m
0
-
2
u
k
,
m
0
+
u
k
-
1
,
m
0
)
+
μ
y
2
2
(
u
k
,
m
+
1
0
-
2
u
k
,
m
0
+
u
k
,
m
-
1
n
)
-
τ
2
2
ϕ
k
,
m
sin
u
k
,
m
0
,
k
,
m
=
1,2
,
…
,
N
-
1
.
Thus, (51) holds; applying the boundary conditions for ∀
n and m
=
1,2
,
…
,
N,
(60)
u
^
0
,
m
n
+
1
=
2
u
0
,
m
n
-
u
0
,
m
n
-
1
+
μ
x
2
(
u
1
,
m
n
-
2
u
0
,
m
n
+
u
-
1
,
m
n
)
+
μ
y
2
(
u
0
,
m
+
1
n
-
2
u
0
,
m
n
+
u
0
,
m
-
1
n
)
-
τ
2
ϕ
0
,
m
sin
u
0
,
m
n
,
u
^
N
,
m
n
+
1
=
2
u
N
,
m
n
-
u
0
,
m
N
-
1
+
μ
x
2
(
u
N
+
1
,
m
n
-
2
u
N
,
m
n
+
u
N
-
1
,
m
n
)
+
μ
y
2
(
u
N
,
m
+
1
n
-
2
u
N
,
m
n
+
u
N
,
m
-
1
n
)
-
τ
2
ϕ
N
,
m
sin
u
N
,
m
n
.
The second order discretization form of the boundary condition ∂
u
(
x
,
y
,
t
)
/
∂
x
=
p
(
x
,
y
,
t
) is
(61)
u
1
,
m
n
-
u
-
1
,
m
n
2
h
x
=
p
(
x
0
,
y
m
,
n
τ
)
,
u
N
+
1
,
m
n
-
u
N
-
1
,
m
n
2
h
x
=
p
(
x
N
,
y
m
,
n
τ
)
.
Substituting (61) into (60), we have
(62)
u
^
0
,
m
n
+
1
=
2
u
0
,
m
n
-
u
0
,
m
n
-
1
+
2
μ
x
2
(
u
1
,
m
n
-
h
x
p
(
x
0
,
y
m
,
n
τ
)
-
u
0
,
m
n
)
+
μ
y
2
(
u
0
,
m
+
1
n
-
2
u
0
,
m
n
+
u
0
,
m
-
1
n
)
-
τ
2
ϕ
0
,
m
sin
u
0
,
m
n
,
u
^
N
,
m
n
+
1
=
2
u
N
,
m
n
-
u
0
,
m
n
+
2
μ
x
2
(
u
N
-
1
,
m
n
+
h
x
p
(
x
N
,
y
m
,
n
τ
)
-
u
N
,
m
n
)
+
μ
y
2
(
u
N
,
m
+
1
n
-
2
u
N
,
m
n
+
u
N
,
m
-
1
n
)
-
τ
2
ϕ
N
,
m
sin
u
N
,
m
n
.Similarly, for∀
n and k
=
1,2
,
…
,
N,
(63)
u
^
k
,
0
n
+
1
=
2
u
k
,
0
n
-
u
k
,
0
n
-
1
+
μ
x
2
(
u
k
+
1,0
n
-
2
u
k
,
0
n
+
u
k
-
1,0
n
)
+
2
μ
y
2
(
u
k
,
1
n
-
h
y
q
(
x
k
,
y
0
,
n
τ
)
-
u
k
,
0
n
)
-
τ
2
ϕ
k
,
0
sin
u
k
,
0
n
,
u
^
k
,
N
n
+
1
=
2
u
k
,
N
n
-
u
k
,
N
n
+
μ
x
2
(
u
k
,
N
n
-
2
u
k
,
N
n
+
u
k
-
1
,
N
n
)
+
2
μ
y
2
(
u
k
,
N
-
1
n
+
h
y
q
(
x
k
,
y
N
,
n
τ
)
-
u
k
,
N
n
)
-
τ
2
ϕ
k
,
N
sin
u
k
,
N
n
.


### 5.4. The Correctors of DFFD and CNFD

The corrector of DFFD can be proposed as(64)
(
1
+
μ
x
2
+
μ
y
2
+
ρ
τ
2
)
u
k
,
m
n
+
1
=
-
1
4
ϕ
k
,
m
sin
u
^
k
,
m
n
+
1
+
2
u
k
,
m
n
-
1
2
ϕ
k
,
m
sin
u
k
,
m
n
-
(
1
+
μ
x
2
+
μ
y
2
-
ρ
τ
2
)
u
k
,
m
n
-
1
-
1
4
ϕ
k
,
m
sin
u
k
,
m
n
-
1
,
for k
,
m
=
0,1
,
…
,
N.The corrector of CNFD can be proposed as(65)
(
1
+
ρ
τ
2
)
u
k
,
m
n
+
1
=
-
1
4
ϕ
k
,
m
sin
u
^
k
,
m
n
+
1
+
2
u
k
,
m
n
-
1
2
ϕ
k
,
m
sin
u
k
,
m
n
-
(
1
-
ρ
τ
2
)
u
k
,
m
n
-
1
-
1
4
ϕ
k
,
m
sin
u
k
,
m
n
-
1
,
for k
,
m
=
0,1
,
…
,
N.In (64) and (65) the corrected values u
k
,
m
n
+
1 instead of the predicted u
^
k
,
m
n
+
1 were also used, and the stability analysis of the corrector is analogous to that developed in Sections 3.1 and 3.2. The initial value problem and boundary problem are solved in the same way as in Section 5.3.


## 5.1. The Predictor of DFFD

Using an analogous scheme as in [26, 27], the predictor value u
^
k
,
m
n
+
1 was evaluated from the following three-time level explicit scheme of DFFD (8):
(51)
(
1
+
μ
x
2
+
μ
y
2
)
u
^
k
,
m
n
+
1
=
2
(
1
+
μ
x
2
+
μ
y
2
)
u
k
,
m
n
-
(
1
+
μ
x
2
+
μ
y
2
)
u
k
,
m
n
-
1
+
μ
x
2
(
u
k
+
1
,
m
n
-
2
u
k
,
m
n
+
u
k
-
1
,
m
n
)
+
μ
y
2
(
u
k
,
m
+
1
n
-
2
u
k
,
m
n
+
u
k
,
m
-
1
n
)
-
τ
2
ϕ
k
,
m
sin
u
k
,
m
n
,
for k
,
m
=
0,1
,
…
,
N.Following a similar approach for the stability analysis of the nonlinear scheme as in Section3.1, it can be proved that the characteristic equation of the predictor is given by
(52)
(
1
+
μ
x
2
+
μ
y
2
)
λ
2
-
2
(
1
-
2
(
μ
x
2
sin
2
β
k
h
x
2
+
μ
y
2
sin
2
γ
m
h
y
2
)
iiiiiiiiiiiiiiii
(
β
k
h
y
2
)
-
1
2
ϕ
k
,
m
cos
⁡
u
s
)
λ
+
1
+
μ
x
2
+
μ
y
2
=
0
.
And the scheme (51) is therefore unconditionally stable.


## 5.2. The Predictor of CNFD

The predictor of CNFD satisfies(53)
u
^
k
,
m
n
+
1
=
2
u
k
,
m
n
-
u
k
,
m
n
-
1
+
μ
x
2
(
u
k
+
1
,
m
n
-
2
u
k
,
m
n
+
u
k
-
1
,
m
n
)
+
μ
y
2
(
u
k
,
m
+
1
n
-
2
u
k
,
m
n
+
u
k
,
m
-
1
n
)
-
τ
2
ϕ
k
,
m
sin
u
k
,
m
n
.
Following the discussion in Section 3.2, the characteristic equation of (53) is
(54)
λ
2
-
2
(
1
-
2
(
μ
x
2
sin
2
β
k
h
x
2
+
μ
y
2
sin
2
γ
m
h
y
2
)
iiiiiiiiiiiiii
(
β
k
h
y
2
)
-
1
2
ϕ
k
,
m
cos
⁡
u
s
)
λ
+
1
=
0
.
Denote B
*
=
1
-
2
(
μ
x
2
sin
2
(
β
k
h
x
/
2
)
+
μ
y
2
sin
2
(
γ
m
h
y
/
2
)
)
-
(
1
/
2
)
ϕ
k
,
m
cos
⁡
u
s, the roots of (54) are
(55)
λ
±
=
B
*
±
(
B
*
)
2
-
1
.Letϕ
~
=
sup
⁡
(
x
,
y
)
∈
[
L
x
0
,
L
x
1
]
×
[
L
y
0
,
L
y
1
]
|
ϕ
(
x
,
y
)
| with |
ϕ
(
x
,
y
)
|
<
+
∞. To ensure |
λ
±
|
≤
1, we require |
B
*
|
≤
1, that is,
(56)
τ
≤
(
1
4
ϕ
~
+
1
h
x
2
+
1
h
y
2
)
-
1
/
2
.


## 5.3. The Predictor Algorithm Implementation for the DFFD Scheme

The predictor of the DFFD scheme is a three-time level explicit scheme; in order to obtain the same second order accuracy, we deal with the initial value of the DFFD scheme, that is, atn
=
0, (51) has the following form:
(57)
u
^
k
,
m
1
=
2
u
k
,
m
0
-
u
k
,
m
-
1
+
μ
x
2
(
u
k
+
1
,
m
0
-
2
u
k
,
m
0
+
u
k
-
1
,
m
0
)
+
μ
y
2
(
u
k
,
m
+
1
0
-
2
u
k
,
m
0
+
u
k
,
m
-
1
0
)
-
τ
2
ϕ
k
,
m
sin
u
k
,
m
0
(
k
,
m
=
1,2
,
…
,
N
)
;
to approximate u
-
1
(
x
,
y
) in the internal points the initial velocity is used, for this we discretize the initial velocity as
(58)
u
1
(
x
k
,
y
m
)
-
u
-
1
(
x
k
,
y
m
)
2
Δ
t
=
g
(
x
k
,
y
m
)
.
Combining (57) with (58), we have
(59)
u
^
k
,
m
1
=
u
k
,
m
0
+
2
τ
g
(
x
k
,
y
m
)
+
μ
x
2
2
(
u
k
+
1
,
m
0
-
2
u
k
,
m
0
+
u
k
-
1
,
m
0
)
+
μ
y
2
2
(
u
k
,
m
+
1
0
-
2
u
k
,
m
0
+
u
k
,
m
-
1
n
)
-
τ
2
2
ϕ
k
,
m
sin
u
k
,
m
0
,
k
,
m
=
1,2
,
…
,
N
-
1
.
Thus, (51) holds; applying the boundary conditions for ∀
n and m
=
1,2
,
…
,
N,
(60)
u
^
0
,
m
n
+
1
=
2
u
0
,
m
n
-
u
0
,
m
n
-
1
+
μ
x
2
(
u
1
,
m
n
-
2
u
0
,
m
n
+
u
-
1
,
m
n
)
+
μ
y
2
(
u
0
,
m
+
1
n
-
2
u
0
,
m
n
+
u
0
,
m
-
1
n
)
-
τ
2
ϕ
0
,
m
sin
u
0
,
m
n
,
u
^
N
,
m
n
+
1
=
2
u
N
,
m
n
-
u
0
,
m
N
-
1
+
μ
x
2
(
u
N
+
1
,
m
n
-
2
u
N
,
m
n
+
u
N
-
1
,
m
n
)
+
μ
y
2
(
u
N
,
m
+
1
n
-
2
u
N
,
m
n
+
u
N
,
m
-
1
n
)
-
τ
2
ϕ
N
,
m
sin
u
N
,
m
n
.
The second order discretization form of the boundary condition ∂
u
(
x
,
y
,
t
)
/
∂
x
=
p
(
x
,
y
,
t
) is
(61)
u
1
,
m
n
-
u
-
1
,
m
n
2
h
x
=
p
(
x
0
,
y
m
,
n
τ
)
,
u
N
+
1
,
m
n
-
u
N
-
1
,
m
n
2
h
x
=
p
(
x
N
,
y
m
,
n
τ
)
.
Substituting (61) into (60), we have
(62)
u
^
0
,
m
n
+
1
=
2
u
0
,
m
n
-
u
0
,
m
n
-
1
+
2
μ
x
2
(
u
1
,
m
n
-
h
x
p
(
x
0
,
y
m
,
n
τ
)
-
u
0
,
m
n
)
+
μ
y
2
(
u
0
,
m
+
1
n
-
2
u
0
,
m
n
+
u
0
,
m
-
1
n
)
-
τ
2
ϕ
0
,
m
sin
u
0
,
m
n
,
u
^
N
,
m
n
+
1
=
2
u
N
,
m
n
-
u
0
,
m
n
+
2
μ
x
2
(
u
N
-
1
,
m
n
+
h
x
p
(
x
N
,
y
m
,
n
τ
)
-
u
N
,
m
n
)
+
μ
y
2
(
u
N
,
m
+
1
n
-
2
u
N
,
m
n
+
u
N
,
m
-
1
n
)
-
τ
2
ϕ
N
,
m
sin
u
N
,
m
n
.Similarly, for∀
n and k
=
1,2
,
…
,
N,
(63)
u
^
k
,
0
n
+
1
=
2
u
k
,
0
n
-
u
k
,
0
n
-
1
+
μ
x
2
(
u
k
+
1,0
n
-
2
u
k
,
0
n
+
u
k
-
1,0
n
)
+
2
μ
y
2
(
u
k
,
1
n
-
h
y
q
(
x
k
,
y
0
,
n
τ
)
-
u
k
,
0
n
)
-
τ
2
ϕ
k
,
0
sin
u
k
,
0
n
,
u
^
k
,
N
n
+
1
=
2
u
k
,
N
n
-
u
k
,
N
n
+
μ
x
2
(
u
k
,
N
n
-
2
u
k
,
N
n
+
u
k
-
1
,
N
n
)
+
2
μ
y
2
(
u
k
,
N
-
1
n
+
h
y
q
(
x
k
,
y
N
,
n
τ
)
-
u
k
,
N
n
)
-
τ
2
ϕ
k
,
N
sin
u
k
,
N
n
.


## 5.4. The Correctors of DFFD and CNFD

The corrector of DFFD can be proposed as(64)
(
1
+
μ
x
2
+
μ
y
2
+
ρ
τ
2
)
u
k
,
m
n
+
1
=
-
1
4
ϕ
k
,
m
sin
u
^
k
,
m
n
+
1
+
2
u
k
,
m
n
-
1
2
ϕ
k
,
m
sin
u
k
,
m
n
-
(
1
+
μ
x
2
+
μ
y
2
-
ρ
τ
2
)
u
k
,
m
n
-
1
-
1
4
ϕ
k
,
m
sin
u
k
,
m
n
-
1
,
for k
,
m
=
0,1
,
…
,
N.The corrector of CNFD can be proposed as(65)
(
1
+
ρ
τ
2
)
u
k
,
m
n
+
1
=
-
1
4
ϕ
k
,
m
sin
u
^
k
,
m
n
+
1
+
2
u
k
,
m
n
-
1
2
ϕ
k
,
m
sin
u
k
,
m
n
-
(
1
-
ρ
τ
2
)
u
k
,
m
n
-
1
-
1
4
ϕ
k
,
m
sin
u
k
,
m
n
-
1
,
for k
,
m
=
0,1
,
…
,
N.In (64) and (65) the corrected values u
k
,
m
n
+
1 instead of the predicted u
^
k
,
m
n
+
1 were also used, and the stability analysis of the corrector is analogous to that developed in Sections 3.1 and 3.2. The initial value problem and boundary problem are solved in the same way as in Section 5.3.


## 6. Numerical Results

In this section we present some numerical results of two schemes for the two-dimensional sine-Gordon equation.


### 6.1. Example 1

To observe the behavior of the numerical method, letϕ
(
x
,
y
)
=
1 in (2); it is tested on the following problem:
(66)
∂
2
u
∂
t
2
+
ρ
∂
u
∂
t
=
∂
2
u
∂
x
2
+
∂
2
u
∂
y
2
-
sin
u
,
-
7
≤
x
,
y
≤
7
,
t
>
0
,
with initial conditions
(67)
u
(
x
,
y
,
0
)
=
4
tan
-
1
(
exp
⁡
(
x
+
y
)
)
,
-
7
≤
x
,
y
≤
7
,
t
>
0
,
∂
u
(
x
,
y
,
0
)
∂
t
=
-
4
exp
⁡
(
x
+
y
)
1
+
exp
⁡
(
2
x
+
2
y
)
,
-
7
≤
x
,
y
≤
7
,
t
>
0
,
and boundary conditions
(68)
iiiiiiiiiiiiiiiiiii
∂
u
∂
x
=
4
exp
⁡
(
x
+
y
+
t
)
exp
⁡
(
2
t
)
+
exp
⁡
(
2
x
+
2
y
)
,
for
x
=
-
7
,
x
=
7
,
-
7
≤
y
≤
7
,
t
>
0
,
iiiiiiiiiiiiiiiiii
∂
u
∂
y
=
4
exp
⁡
(
x
+
y
+
t
)
exp
⁡
(
2
t
)
+
exp
⁡
(
2
x
+
2
y
)
,
for
y
=
-
7
,
y
=
7
,
-
7
≤
x
≤
7
,
t
>
0
.The theoretical solution of this problem, in which the parameterρ
=
0, is given by
(69)
u
(
x
,
y
,
t
)
=
4
tan
-
1
(
exp
⁡
(
x
+
y
-
t
)
)
.As mentioned in Section2, the truncation error of the DFFD scheme is O
(
τ
2
+
h
x
2
+
h
y
2
+
(
τ
/
h
x
)
2
+
(
τ
/
h
y
)
2
). One can easily see that if we choose the time step τ
=
h
x
2
=
h
y
2, for DFFD scheme and CNFD scheme, they are second-order convergence. Respectively by the PC-DFFD scheme and the PC-CNFD scheme we take the same space step and work out equation (66), where the initial conditions (67) and the boundary conditions (68) are employed. The error is measured by the ∥
u
E
-
u
N
∥
L
2, ∥
u
E
-
u
N
∥
L
∞ of the difference between the exact solution u
E given by (69) and the numerical solution u
N, the time is T
=
7.0. Table 1 shows that it needs a sufficiently small space and time step to keep the stability of PC-CNFD scheme while the PC-DFFD scheme is unconditionally stable, which is in consistent with the theoretical results.Table 1
Two different errors for∥
u
E
-
u
N
∥
L
2 and ∥
u
E
-
u
N
∥
L
∞ at T
=
7.
h
x
=
h
y
τ
PC-DFFD
PC-CNFD
∥
u
E
-
u
N
∥
L
2
∥
u
E
-
u
N
∥
L
∞
∥
u
E
-
u
N
∥
L
2
∥
u
E
-
u
N
∥
L
∞
1
10
1
100
0.722164
e
-
01
0.3503924
e
-
02
Unstable
Unstable
1
20
1
400
0.787769
e
-
02
0.2436561
e
-
03
Unstable
Unstable
1
40
1
1600
0.653123
e
-
04
0.1643035
e
-
05
0.249841
e
-
03
0.404723
e
-
04Take the same step lengthh
x
=
h
y
=
0.25, τ
=
0.01 as in [21]; we apply the PC-DFFD scheme to computation of the solution to (66); the absolute errors are given at times t
=
4,7
,
10,20,25. Figure 1 shows that when the space and time step are the same, compared with the absolute error at t
=
7 in [21], the accuracy of the PC-DFFD scheme is much better than the method in [21] and the absolute error becomes about 0.1 times smaller than the RBF method in [21]. Compared with [13, 15, 18–21, 26, 27], the PC-DFFD scheme is much better than the numerical algorithms presented in these articles. Figure 2 shows that the PC-DFFD scheme has better stability at t
=
25 or much longer time, where the absolute error is nearly the same as that at t
=
4. Compared with the scheme proposed in [15, 18–21, 26, 27], the PC-DFFD scheme maintains its simplicity and better stabilities. The accumulation of absolute errors can not lead to infinite increases of them; hence, the scheme can be applied in long-time numerical simulations.(a) Absolute error of the solution of RBF method [21] and (b) absolute error of the solution of the PC-DFFD scheme in t
=
7 for test problems (66)–(68).
(a)
(b)Absolute error of the solution of the PC-DFFD scheme int
=
4,10,20,25 for test problems (66)–(68).
(a)
(b)
(c)
(d)It is known (see e.g., [7–9, 14–16, 18–21, 26, 27]) that, when ρ
=
0, for the sine-Gordon equation the energy given by the following expression,
(70)
E
(
t
)
=
E
(
0
)
=
1
2
∫
‍
∫
‍
[
u
x
2
+
u
y
2
+
u
t
2
+
2
(
1
-
cos
⁡
u
)
]
d
x
d
y
,
is conserved. We also investigate this property for sine-Gordon equation; the evaluation of E
(
t
) is performed using the composite trapezoidal rule for integration.In the numerical calculations that follow, various cases involving line and ring solitons for the solution of (5) are reported. In all the following experiments, by the PC-DFFD scheme,we choose τ
=
0.001, h
x
=
h
y
=
0.1, and the boundary conditions are taken to be
(71)
∂
u
∂
x
=
0
,
∂
u
∂
y
=
0
.


### 6.2. Example 2 (Superposition of Two Line Solitons)

It has been known that superpositions of two orthogonal line solitons can be acquired forϕ
(
x
,
y
)
=
1, with initial conditions [14]
(72)
f
(
x
,
y
)
=
4
(
tan
-
1
(
exp
⁡
⁡
(
x
)
)
+
tan
-
1
(
exp
⁡
⁡
(
y
)
)
)
,
g
(
x
,
y
)
=
0
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
-
7
≤
x
,
y
≤
7
.Whenρ
=
0, the results in Figure 3 show the break up of two orthogonal line solitons which are parallel to the diagonal y
=
-
x and are moving away from each other in the direction of y
=
x, undisturbed. From t
=
0 to t
=
4, the solutions are almost not varying until while at t
=
7 a deformation has appeared, while at t
=
15 its deformation and original morphology will have a great change, but its energy E
(
t
) remains basically the same. Table 2 shows that the scheme keeps better energy conservation and reliability, compared with others presented by Dehghan and Shokri [21] and Mirzaei and Dehghan [16], Bratsos [15, 26, 27]. For a small value of ρ, the dissipative term is found to have little effect on the superposition of two line solitons. For a large value of ρ
=
1.5, Figure 4 shows that at t
=
9, even t
=
15, the function graph virtually maintains the primary condition at t
=
0, while its energy E
(
t
) decreases as ρ increases and damages its attributes of energy conservation illustrated in Table 2. However, the dissipative term is found to slow down the separation and break-up of two orthogonal line solitons as time increases. The present solutions are in good agreement with the corresponding results of [15, 16, 18–21, 26, 27].Table 2
The energyE
(
t
) of the superposition of two orthogonal line solitions.
ρ
Initial
E
(
1
)
E
(
4
)
E
(
7
)
E
(
9
)
E
(
15
)
0.0
175.5745
175.5748
175.5750
175.5737
175.5768
175.5786
0.5
175.5745
172.8013
150.4754
124.0258
108.9841
82.7631
1.5
175.5745
170.7209
112.2087
100.5008
87.3902
76.4529Figure 3
The numerical solutions at timest
=
0,1
,
2,3
,
4,7
,
9,15 for superposition of two line solitons.The numerical solutions at timest
=
9,15, with ρ
=
1.5 for superposition of two line solitons.
(a)
(b)


### 6.3. Example 3 (Circular ring Solitons)

Bogolyubskii and Makhankov [28], Bogolyubskii [29], and Christiansen and Lomdahl [12] have investigated numerically the behavior of a circular ring quasisoliton or pulsion arising from the two-dimensional sine-Gordon equation. Circular ring solitons are found for the case ϕ
(
x
,
y
)
=
1 and initial conditions [12, 14]. Consider
(73)
f
(
x
,
y
)
=
4
tan
-
1
exp
⁡
⁡
(
3
-
x
2
+
y
2
)
,
g
(
x
,
y
)
=
0
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
-
7
≤
x
,
y
≤
7
.In Figure5, the numerical solutions of circular ring solitons for the ρ
=
0 at t
=
2.8,5.6,8.4,11.2,15,18,20 are shown in terms of sin
(
u
/
2
) for three-dimensional picture. The soliton from its initial position, where it appears as two homocentric ring solitons, is shrinking until t
=
2.8 appears as a single-ring soliton. From t
=
5.6, which could be considered as the beginning of the expansion phase, a radiation appears. This expansion continues until t
=
11.2, where the soliton is almost reformed. These results are in agreement with the published ones in [15, 16, 18–21, 26, 27]. Table 3 presents the values of E
(
t
) at some selected times t, that shows the energy remains constant for a longer time as time increases. For the different ρ
=
0.05,0.5, in Figure 6, with solution changes of a smaller ρ
=
0.05 (Figure 6(a)) and a larger ρ
=
0.5 (Figure 6(b)), the results are the same as [21]. As ρ increases, the initial shrunk ring soliton was found to be changing more slowly from its initial position as time increases; the dissipative term is slowing down the evolution of the line soliton as time increases.Table 3
The energy of the circular ring solitons.
Timet
t
=
0
t
=
2.8
t
=
8.4
t
=
11.2
t
=
15
t
=
20
E
(
t
)
150.4597
150.5537
150.5604
150.5145
150.5599
150.5189Figure 5
Circular ring solitons solution att
=
0,2.8,5.6,8.4,11.2,15,18,20 with ρ
=
0.Circular ring solitons solution att
=
0,2.8,5.6,8.4,11.2,15,18,20 with (a) ρ
=
0.05 and (b) ρ
=
0.5.


### 6.4. Example 4 (Collision of Two Circular Solitons)

The collision of two expanding circular ring solitons is considered withϕ
(
x
,
y
)
=
1 and initial conditions
(74)
f
(
x
,
y
)
=
4
tan
-
1
exp
⁡
(
4
-
(
x
+
3
)
2
+
(
y
+
7
)
2
0.436
)
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiii
-
30
≤
x
≤
10
,
-
21
≤
y
≤
7
,
g
(
x
,
y
)
=
4.13
sinh
⁡
(
4
-
(
x
+
3
)
2
+
(
y
+
7
)
2
0.436
)
,
iiiiiiiiiiiiiiiiiiiiiiiiiii
-
30
≤
x
≤
10
,
-
21
≤
y
≤
7
.Numerical simulation presented in Figure7 is for sin
(
u
/
2
) at levels t = 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, and 20 with ρ
=
0, respectively. The solution which is shown includes the extension across x
=
-
10 and y
=
-
7 by symmetry properties of the problem [14, 20, 21]. Figure 7 demonstrates the collision between two expanding circular ring solitons in which two smaller ring solitons bounding an annular region emerge into a large ring soliton. The simulated solution is again precisely consistent to existing results; contour maps are given to show more clearly the movement of solitons. Though minor disturbances can be observed in middle of the numerical solution, probably due to the transactions following the symmetry features in computations, the overall simulation results match well those described in [12, 14, 16, 20, 21] with satisfaction. Simultaneously, it is found that after t
=
16, the originally formed large ring soliton is to split into two line solitons to the boundary extension.Collision of two circular solitons, att
=
0,2
,
4,6
,
8,10,12,14,16,18,20 with ρ
=
0.


### 6.5. Example 5 (Collision of Four Circular ring Solitons)

A collision of four expanding circular ring solitons is investigated, forϕ
(
x
,
y
)
=
1, and initial conditions [14]
(75)
f
(
x
,
y
)
=
4
tan
-
1
(
exp
⁡
(
4
-
(
x
+
3
)
2
+
(
y
+
3
)
2
0.436
)
)
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
-
10
≤
x
,
y
≤
10
,
g
(
x
,
y
)
=
0.436
cosh
⁡
(
exp
⁡
(
[
4
-
(
x
+
3
)
2
+
(
y
+
3
)
2
]
/
0.436
)
)
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
-
10
≤
x
,
y
≤
10
.The simulation is based on an extension acrossx
=
-
10 and y
=
-
10 due to the symmetry. The results are depicted in Figure 8 for ρ
=
0 at t = 0, 2.5, 4, 5, 6, 7.5, 8, 10, 12, and 15 in terms of sin
(
u
/
2
)
); Figure 8 demonstrates precisely the collision between four expanding circular ring solitons in which the smaller ring solitons bounding an annular region emerge into a large ring soliton. The results are in good agreement with corresponding surfaces given in [13, 14, 16, 20, 21, 26]. Similarly, as time increases, at t
=
15 the large soliton is to split and form 5 or even more solutions of varying sizes.Collision of four circular ring solitons, att
=
0,2.5,4
,
5,6
,
7.5,8
,
10,12,15 with ρ
=
0.


## 6.1. Example 1

To observe the behavior of the numerical method, letϕ
(
x
,
y
)
=
1 in (2); it is tested on the following problem:
(66)
∂
2
u
∂
t
2
+
ρ
∂
u
∂
t
=
∂
2
u
∂
x
2
+
∂
2
u
∂
y
2
-
sin
u
,
-
7
≤
x
,
y
≤
7
,
t
>
0
,
with initial conditions
(67)
u
(
x
,
y
,
0
)
=
4
tan
-
1
(
exp
⁡
(
x
+
y
)
)
,
-
7
≤
x
,
y
≤
7
,
t
>
0
,
∂
u
(
x
,
y
,
0
)
∂
t
=
-
4
exp
⁡
(
x
+
y
)
1
+
exp
⁡
(
2
x
+
2
y
)
,
-
7
≤
x
,
y
≤
7
,
t
>
0
,
and boundary conditions
(68)
iiiiiiiiiiiiiiiiiii
∂
u
∂
x
=
4
exp
⁡
(
x
+
y
+
t
)
exp
⁡
(
2
t
)
+
exp
⁡
(
2
x
+
2
y
)
,
for
x
=
-
7
,
x
=
7
,
-
7
≤
y
≤
7
,
t
>
0
,
iiiiiiiiiiiiiiiiii
∂
u
∂
y
=
4
exp
⁡
(
x
+
y
+
t
)
exp
⁡
(
2
t
)
+
exp
⁡
(
2
x
+
2
y
)
,
for
y
=
-
7
,
y
=
7
,
-
7
≤
x
≤
7
,
t
>
0
.The theoretical solution of this problem, in which the parameterρ
=
0, is given by
(69)
u
(
x
,
y
,
t
)
=
4
tan
-
1
(
exp
⁡
(
x
+
y
-
t
)
)
.As mentioned in Section2, the truncation error of the DFFD scheme is O
(
τ
2
+
h
x
2
+
h
y
2
+
(
τ
/
h
x
)
2
+
(
τ
/
h
y
)
2
). One can easily see that if we choose the time step τ
=
h
x
2
=
h
y
2, for DFFD scheme and CNFD scheme, they are second-order convergence. Respectively by the PC-DFFD scheme and the PC-CNFD scheme we take the same space step and work out equation (66), where the initial conditions (67) and the boundary conditions (68) are employed. The error is measured by the ∥
u
E
-
u
N
∥
L
2, ∥
u
E
-
u
N
∥
L
∞ of the difference between the exact solution u
E given by (69) and the numerical solution u
N, the time is T
=
7.0. Table 1 shows that it needs a sufficiently small space and time step to keep the stability of PC-CNFD scheme while the PC-DFFD scheme is unconditionally stable, which is in consistent with the theoretical results.Table 1
Two different errors for∥
u
E
-
u
N
∥
L
2 and ∥
u
E
-
u
N
∥
L
∞ at T
=
7.
h
x
=
h
y
τ
PC-DFFD
PC-CNFD
∥
u
E
-
u
N
∥
L
2
∥
u
E
-
u
N
∥
L
∞
∥
u
E
-
u
N
∥
L
2
∥
u
E
-
u
N
∥
L
∞
1
10
1
100
0.722164
e
-
01
0.3503924
e
-
02
Unstable
Unstable
1
20
1
400
0.787769
e
-
02
0.2436561
e
-
03
Unstable
Unstable
1
40
1
1600
0.653123
e
-
04
0.1643035
e
-
05
0.249841
e
-
03
0.404723
e
-
04Take the same step lengthh
x
=
h
y
=
0.25, τ
=
0.01 as in [21]; we apply the PC-DFFD scheme to computation of the solution to (66); the absolute errors are given at times t
=
4,7
,
10,20,25. Figure 1 shows that when the space and time step are the same, compared with the absolute error at t
=
7 in [21], the accuracy of the PC-DFFD scheme is much better than the method in [21] and the absolute error becomes about 0.1 times smaller than the RBF method in [21]. Compared with [13, 15, 18–21, 26, 27], the PC-DFFD scheme is much better than the numerical algorithms presented in these articles. Figure 2 shows that the PC-DFFD scheme has better stability at t
=
25 or much longer time, where the absolute error is nearly the same as that at t
=
4. Compared with the scheme proposed in [15, 18–21, 26, 27], the PC-DFFD scheme maintains its simplicity and better stabilities. The accumulation of absolute errors can not lead to infinite increases of them; hence, the scheme can be applied in long-time numerical simulations.(a) Absolute error of the solution of RBF method [21] and (b) absolute error of the solution of the PC-DFFD scheme in t
=
7 for test problems (66)–(68).
(a)
(b)Absolute error of the solution of the PC-DFFD scheme int
=
4,10,20,25 for test problems (66)–(68).
(a)
(b)
(c)
(d)It is known (see e.g., [7–9, 14–16, 18–21, 26, 27]) that, when ρ
=
0, for the sine-Gordon equation the energy given by the following expression,
(70)
E
(
t
)
=
E
(
0
)
=
1
2
∫
‍
∫
‍
[
u
x
2
+
u
y
2
+
u
t
2
+
2
(
1
-
cos
⁡
u
)
]
d
x
d
y
,
is conserved. We also investigate this property for sine-Gordon equation; the evaluation of E
(
t
) is performed using the composite trapezoidal rule for integration.In the numerical calculations that follow, various cases involving line and ring solitons for the solution of (5) are reported. In all the following experiments, by the PC-DFFD scheme,we choose τ
=
0.001, h
x
=
h
y
=
0.1, and the boundary conditions are taken to be
(71)
∂
u
∂
x
=
0
,
∂
u
∂
y
=
0
.


## 6.2. Example 2 (Superposition of Two Line Solitons)

It has been known that superpositions of two orthogonal line solitons can be acquired forϕ
(
x
,
y
)
=
1, with initial conditions [14]
(72)
f
(
x
,
y
)
=
4
(
tan
-
1
(
exp
⁡
⁡
(
x
)
)
+
tan
-
1
(
exp
⁡
⁡
(
y
)
)
)
,
g
(
x
,
y
)
=
0
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
-
7
≤
x
,
y
≤
7
.Whenρ
=
0, the results in Figure 3 show the break up of two orthogonal line solitons which are parallel to the diagonal y
=
-
x and are moving away from each other in the direction of y
=
x, undisturbed. From t
=
0 to t
=
4, the solutions are almost not varying until while at t
=
7 a deformation has appeared, while at t
=
15 its deformation and original morphology will have a great change, but its energy E
(
t
) remains basically the same. Table 2 shows that the scheme keeps better energy conservation and reliability, compared with others presented by Dehghan and Shokri [21] and Mirzaei and Dehghan [16], Bratsos [15, 26, 27]. For a small value of ρ, the dissipative term is found to have little effect on the superposition of two line solitons. For a large value of ρ
=
1.5, Figure 4 shows that at t
=
9, even t
=
15, the function graph virtually maintains the primary condition at t
=
0, while its energy E
(
t
) decreases as ρ increases and damages its attributes of energy conservation illustrated in Table 2. However, the dissipative term is found to slow down the separation and break-up of two orthogonal line solitons as time increases. The present solutions are in good agreement with the corresponding results of [15, 16, 18–21, 26, 27].Table 2
The energyE
(
t
) of the superposition of two orthogonal line solitions.
ρ
Initial
E
(
1
)
E
(
4
)
E
(
7
)
E
(
9
)
E
(
15
)
0.0
175.5745
175.5748
175.5750
175.5737
175.5768
175.5786
0.5
175.5745
172.8013
150.4754
124.0258
108.9841
82.7631
1.5
175.5745
170.7209
112.2087
100.5008
87.3902
76.4529Figure 3
The numerical solutions at timest
=
0,1
,
2,3
,
4,7
,
9,15 for superposition of two line solitons.The numerical solutions at timest
=
9,15, with ρ
=
1.5 for superposition of two line solitons.
(a)
(b)


## 6.3. Example 3 (Circular ring Solitons)

Bogolyubskii and Makhankov [28], Bogolyubskii [29], and Christiansen and Lomdahl [12] have investigated numerically the behavior of a circular ring quasisoliton or pulsion arising from the two-dimensional sine-Gordon equation. Circular ring solitons are found for the case ϕ
(
x
,
y
)
=
1 and initial conditions [12, 14]. Consider
(73)
f
(
x
,
y
)
=
4
tan
-
1
exp
⁡
⁡
(
3
-
x
2
+
y
2
)
,
g
(
x
,
y
)
=
0
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
-
7
≤
x
,
y
≤
7
.In Figure5, the numerical solutions of circular ring solitons for the ρ
=
0 at t
=
2.8,5.6,8.4,11.2,15,18,20 are shown in terms of sin
(
u
/
2
) for three-dimensional picture. The soliton from its initial position, where it appears as two homocentric ring solitons, is shrinking until t
=
2.8 appears as a single-ring soliton. From t
=
5.6, which could be considered as the beginning of the expansion phase, a radiation appears. This expansion continues until t
=
11.2, where the soliton is almost reformed. These results are in agreement with the published ones in [15, 16, 18–21, 26, 27]. Table 3 presents the values of E
(
t
) at some selected times t, that shows the energy remains constant for a longer time as time increases. For the different ρ
=
0.05,0.5, in Figure 6, with solution changes of a smaller ρ
=
0.05 (Figure 6(a)) and a larger ρ
=
0.5 (Figure 6(b)), the results are the same as [21]. As ρ increases, the initial shrunk ring soliton was found to be changing more slowly from its initial position as time increases; the dissipative term is slowing down the evolution of the line soliton as time increases.Table 3
The energy of the circular ring solitons.
Timet
t
=
0
t
=
2.8
t
=
8.4
t
=
11.2
t
=
15
t
=
20
E
(
t
)
150.4597
150.5537
150.5604
150.5145
150.5599
150.5189Figure 5
Circular ring solitons solution att
=
0,2.8,5.6,8.4,11.2,15,18,20 with ρ
=
0.Circular ring solitons solution att
=
0,2.8,5.6,8.4,11.2,15,18,20 with (a) ρ
=
0.05 and (b) ρ
=
0.5.


## 6.4. Example 4 (Collision of Two Circular Solitons)

The collision of two expanding circular ring solitons is considered withϕ
(
x
,
y
)
=
1 and initial conditions
(74)
f
(
x
,
y
)
=
4
tan
-
1
exp
⁡
(
4
-
(
x
+
3
)
2
+
(
y
+
7
)
2
0.436
)
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiii
-
30
≤
x
≤
10
,
-
21
≤
y
≤
7
,
g
(
x
,
y
)
=
4.13
sinh
⁡
(
4
-
(
x
+
3
)
2
+
(
y
+
7
)
2
0.436
)
,
iiiiiiiiiiiiiiiiiiiiiiiiiii
-
30
≤
x
≤
10
,
-
21
≤
y
≤
7
.Numerical simulation presented in Figure7 is for sin
(
u
/
2
) at levels t = 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, and 20 with ρ
=
0, respectively. The solution which is shown includes the extension across x
=
-
10 and y
=
-
7 by symmetry properties of the problem [14, 20, 21]. Figure 7 demonstrates the collision between two expanding circular ring solitons in which two smaller ring solitons bounding an annular region emerge into a large ring soliton. The simulated solution is again precisely consistent to existing results; contour maps are given to show more clearly the movement of solitons. Though minor disturbances can be observed in middle of the numerical solution, probably due to the transactions following the symmetry features in computations, the overall simulation results match well those described in [12, 14, 16, 20, 21] with satisfaction. Simultaneously, it is found that after t
=
16, the originally formed large ring soliton is to split into two line solitons to the boundary extension.Collision of two circular solitons, att
=
0,2
,
4,6
,
8,10,12,14,16,18,20 with ρ
=
0.


## 6.5. Example 5 (Collision of Four Circular ring Solitons)

A collision of four expanding circular ring solitons is investigated, forϕ
(
x
,
y
)
=
1, and initial conditions [14]
(75)
f
(
x
,
y
)
=
4
tan
-
1
(
exp
⁡
(
4
-
(
x
+
3
)
2
+
(
y
+
3
)
2
0.436
)
)
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
-
10
≤
x
,
y
≤
10
,
g
(
x
,
y
)
=
0.436
cosh
⁡
(
exp
⁡
(
[
4
-
(
x
+
3
)
2
+
(
y
+
3
)
2
]
/
0.436
)
)
,
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii
-
10
≤
x
,
y
≤
10
.The simulation is based on an extension acrossx
=
-
10 and y
=
-
10 due to the symmetry. The results are depicted in Figure 8 for ρ
=
0 at t = 0, 2.5, 4, 5, 6, 7.5, 8, 10, 12, and 15 in terms of sin
(
u
/
2
)
); Figure 8 demonstrates precisely the collision between four expanding circular ring solitons in which the smaller ring solitons bounding an annular region emerge into a large ring soliton. The results are in good agreement with corresponding surfaces given in [13, 14, 16, 20, 21, 26]. Similarly, as time increases, at t
=
15 the large soliton is to split and form 5 or even more solutions of varying sizes.Collision of four circular ring solitons, att
=
0,2.5,4
,
5,6
,
7.5,8
,
10,12,15 with ρ
=
0.


## 7. Conclusion

In this paper, two numerical methods, CNFD and DFFD, are constructed for solving the two-dimensional sin-Gordon equation. The stability, convergence, and error estimates are discussed. The DFFD scheme is unconditionally stable and convergent, while the CNFD scheme requires more critical space and time step constraints, in order to guarantee its stability and convergence. We establish two explicit PC-CNFD schemes and PC-DFFD scheme, whose stabilities are discussed. The PC-DFFD scheme is unconditionally stable. The numerical experiments indicate that the PC-DFFD scheme has better stability in comparison with the methods in [13, 15, 18–21, 26, 27] and is proper to the long-time computation.



---

*Source: 784387-2014-10-29.xml*