Applied Mathematics and Computation 260 (2015) 269287

Contents lists available at ScienceDirect

Applied Mathematics and Computation

journal homepage: www.elsevier.com/locate/amc

A new compact nite difference scheme for solving the complex

GinzburgLandau equation
Yun Yan a, Frederick Ira Moxley III b, Weizhong Dai a,
a
b

Mathematics & Statistics, College of Engineering & Science Louisiana Tech University, Ruston, LA 71272, USA
Hearne Institute for Theoretical Physics, Department of Physics & Astronomy Louisiana State University, Baton Rouge, LA 70803, USA

a r t i c l e

i n f o

Keywords:
Compact nite difference scheme
Complex GinzburgLandau equation
Stability

a b s t r a c t
The complex GinzburgLandau equation is often encountered in physics and engineering
applications, such as nonlinear transmission lines, solitons, and superconductivity. However,
it remains a challenge to develop simple, stable and accurate nite difference schemes for
solving the equation because of the nonlinear term. Most of the existing schemes are obtained
based on the CrankNicolson method, which is fully implicit and must be solved iteratively for
each time step. In this article, we present a fourth-order accurate iterative scheme, which leads
to a tri-diagonal linear system in 1D cases. We prove that the present scheme is unconditionally
stable. The scheme is then extended to 2D cases. Numerical errors and convergence rates of
the solutions are tested by several examples.
2015 Elsevier Inc. All rights reserved.

1. Introduction
The complex GinzburgLandau equation is often encountered in physics and engineering applications, such as nonlinear
transmission lines [15], solitons [610], and superconductivity [1117]. However, it remains a challenge to develop a simple,
stable and accurate nite difference scheme for solving this equation due to the nonlinear term, particularly in multi-dimensional
cases. There are many numerical schemes have been proposed for solving this equation [1129]. In particular, Tsertsvadze [27]
proposed a second-order CrankNicolson type of nite difference scheme for a one-dimensional (1D) GinzburgLandau equation,
and Sun and Zhu [28] proved its unconditional convergence in the l norm. Recently, Hu [29] developed several fourth-order
compact nite difference schemes by coupling the CrankNicolson method with the fourth-order compact nite difference
method [30]. The unconditional convergence in the l norm of these schemes was analyzed. However, these fourth-order
accurate nite difference schemes are implicit, particularly the implicit approximation for the nonlinear term in some schemes,
which must be solved iteratively for each time step. In this study, we present a fourth-order accurate iterative scheme which
is also obtained based on the compact nite difference method and CrankNicolson method, but leads to a tri-diagonal linear
system in 1D cases. We further prove the scheme to be unconditionally stable in the l2 norm by using the discrete energy
techniques [28,29,3138]. The obtained scheme is then extended to 2D cases. Finally, the numerical errors and convergence rates
of the solutions are tested by several examples.

Corresponding author. Tel.: +1 318 257 3301.

E-mail address: [email protected] (W. Dai).

http://dx.doi.org/10.1016/j.amc.2015.03.053
0096-3003/ 2015 Elsevier Inc. All rights reserved.

270

Y. Yan et al. / Applied Mathematics and Computation 260 (2015) 269287

2. Finite difference schemes

We rst consider the 1D complex GinzburgLandau equation with initial and periodic boundary conditions as follows:

u
2u
= (1 + i) 2 + Ru (1 + i)|u|q1 u, x (0, L), t (0, T];
t
x

(2.1a)

u(x, 0) = u0 (x),

(2.1b)

x [0, L];

u(x, t) = u(x + L, t),

t [0, T];

(2.1c)

where and are real constants, q 3 is a positive integer, i = 1, the functions u(x, t) and u0 (x) are complex valued functions,
R is a positive constant, and L is the period of u(x, t) with respect to x.
To develop a nite difference scheme for solving the above GinzburgLandau problem, we rst divide spatial interval [0, L]
L
, and the time interval [0, T] is divided into N subintervals with a time step
into M subintervals where the grid size is h = M
T
n
= N . We denote uj to be the numerical approximation of u(xj , tn ), where xj = jh, tn = n , 0 j M and 0 n N. It can be
seen from periodicity that unM = un0 , unM+1 = un1 and so on. Furthermore, we use the following rst-order and second-order nite
difference operators:

n
t uj

un+1
unj
j

n+ 12
t uj

n+ 32

uj

n+ 12

uj

n
uj+1 unj

h
x unj =

un unM

un unM

M+1
1
,
h
h

(2.2a)

0 j M 1,
(2.2b)

j = M;

n
u1 2un0 + un1
unM1 2un0 + un1

h2
h2

un 2un + un
j1
j
j+1
x2 unj =
,

2
h

un1 2unM + unM1

un 2unM + unM1

M+1

,
h2
h2

j = 0,
(2.2c)

1 j M 1,
j = M.

For the discrete space, a fourth-order compact difference scheme [30] at all grid points, xj , 0 j M, can be written as:

1
12

2 u(xj1 , t)
2 u(xj , t) 2 u(xj+1 , t)
1
+ 10
+
= 2 [u(xj1 , t) 2u(xj , t) + u(xj+1 , t)] + O(h4 ),
x2
x2
x2
h

(2.3)

which can be denoted using the second-order nite difference operator x2 as


h2 2
x
1+
12

2 u(xj , t)
= x2 u(xj , t) + O(h4 ).
x2

(2.4)

We now rewrite Eq. (2.1) as (1 + i) x2u = ut Ru + (1 + i)|u|q1 u at (xj , t), multiply it both sides by 1 +
use Eq. (2.4). This gives
2

(1 + i)x2 u(xj , t) = 1 +

h2 2

12 x

h2 2
12 x

, and then

u(xj , t)
Ru(xj , t) + (1 + i)|u(xj , t)|q1 u(xj , t) + O(h4 ).
t

(2.5)

Using the CrankNicolson technique to Eq. (2.5) with respect to time t at tn+ 1 , we obtain
2






2 u(xj , tn+1 ) + x2 u(xj , tn )
h2 2
h2 2 u(xj , tn+1 ) u(xj , tn )
x
= (1 + i) x
x u xj , tn+ 12
+R 1+
1+
12

2
12


q1 

h2 2  

x u xj , tn+ 12  u xj , tn+ 12 + O( 2 + h4 ).
(1 + i) 1 +
12

(2.6a)

Y. Yan et al. / Applied Mathematics and Computation 260 (2015) 269287

271

Note that u(xj , tn+ 1 ) must be known in order to obtain u(xj , tn+1 ) in Eq. (2.6). We may use a similar argument in obtaining
2

Eq. (2.6) from Eq. (2.5) at tn+1 and obtain


1+

h2 2

12 x

u(x , t 3 ) u(x , t 1 )
j n+
j n+
2

= (1 + i)

x2 u(xj , tn+ 32 ) + x2 u(xj , tn+ 12 )



(1 + i) 1 +

h 2

12 x



h2 2
+R 1+
x u(xj , tn+1 )
12

|u(xj , tn+1 )|q1 u(xj , tn+1 ) + O( 2 + h4 ).

(2.6b)

Thus, coupling Eq. (2.6) and then dropping the truncation error O( 2 + h4 ), we obtain an iterative fourth-order in space compact
nite difference scheme for solving Eq. (2.1) as follows:

h2 2

1+
12 x

h2 2

1+
12 x

n
t uj

= (1 + i)

2
x

un+1
+ unj
j
2

n+ 12
t uj

= (1 + i)x2

n+ 32

uj







h2 2 n+ 12
h2
n+ 1 q1 n+ 1
+R 1+
x uj (1 + i) 1 + x2 uj 2  uj 2 ,
12
12

n+ 12

+ uj
2

(2.7a)




q1

2
2
+ R 1 + h x2 un+1 (1 + i) 1 + h x2 un+1  un+1 ,
j
j
j
12
12
(2.7b)
1

where u0j = u0 (xj ), and 0 j M, 0 n N 1. To start the above iteration, uj2 , 0 j M, must be known. This can be obtained
using other methods, such as those proposed methods in [29] or the explicit pseudo-spectral methods [39] built in MATLAB. It
should be pointed out that for each time step, one needs only to solve a tri-diagonal linear system for un+1
from Eq. (2.7) once
j
n+ 12

unj and uj

n+ 32

are known, and then substitute ujn+1 into Eq. (2.7) to obtain uj

by solving another tri-diagonal linear system for

n+ 3
uj 2 .

Thus, the computation is simple and fast.

The above nite difference scheme can be extended to 2D cases. Indeed, we consider the 2D complex GinzburgLandau
equation as

u
2u 2u
= (1 + i)
+
+ Ru (1 + i)|u|q1 u, (x, y) (0, Lx ) (0, Ly ), t (0, T],
t
x2 y2

(2.8)

with initial and periodic boundary conditions, where Lx , Ly are periods of u(x, y, t) in the x and y directions, respectively. Again,
we divide [0, Lx ] into Mx subintervals, [0, Ly ] into My subintervals, and [0, T] into N subintervals with mesh sizes hx , hy and ,
respectively. We denote unj,k to be the numerical approximation of u(xj , yk , tn ), where xj = jhx , yk = khy , tn = n , 0 j Mx ,
0 k My , and 0 n N.
Instead using Eq. (2.3), we now employ a fourth-order compact nite difference scheme for Poisson equation x2u + y2u =
f (x, y) given in [40] as
2




D2y
D2y
D2x D2y
D2x D2y


D2x
D2x
+
f (xj , yk ) + O h4x + h4y ,
+
+
+
(
x
,
y
)
=
1
+
u
j
k
12
12
h2x
h2y
12h2x
12h2y

(2.9)

where the nite difference operators D2x , D2y and D2x D2y are dened as follows:

D2x uj,k

j = 0, 0 k My ,

uMx 1,k 2u0,k + u1,k ,

1 j Mx 1, 0 k My ,
= uj1,k 2uj,k + uj+1,k ,

u1,k 2uMx ,k + uMx 1,k , j = Mx , 0 k My ;

(2.10a)

D2y uj,k

0 j Mx , k = 0,

uj,My 1 2uj,0 + uj,1 ,

= uj,k1 2uj,k + uj,k+1 , 0 j Mx , 1 k My 1,

uj,1 2uj,My + uj,My 1 , 0 j Mx , k = My ;

(2.10b)

272

Y. Yan et al. / Applied Mathematics and Computation 260 (2015) 269287

D2x D2y uj,k

(u1,1 + uMx 1,1 + uMx 1,My 1 + u1,My 1 )

2(u1,0 + u0,1 + uMx 1,0 + u0,My 1 ) + 4u0,0 ,

(
u1,k+1 + uMx 1,k+1 + uMx 1,k1 + u1,k1 )

2(u1,k + u0,k+1 + uMx 1,k + u0,k1 ) + 4u0,k ,

(u1,1 + uMx 1,1 + uMx 1,My 1 + u1,My 1 )

2(u1,My + u0,1 + uMx 1,My + u0,My 1 ) + 4u0,My ,

(uj+1,1 + uj1,1 + uj1,My 1 + uj+1,My 1 )

2(uj+1,0 + uj,1 + uj1,0 + uj,My 1 ) + 4uj,0 ,

(uj+1,k+1 + uj1,k+1 + uj1,k1 + uj+1,k1 )

=

2(uj+1,k + uj,k+1 + uj1,k + uj,k1 ) + 4uj,k ,

(uj+1,1 + uj1,1 + uj1,My 1 + uj+1,My 1 )

2(uj+1,My + uj,1 + uj1,My + uj,My 1 ) + 4uj,My ,

(
u

1,1 + uMx 1,1 + uMx 1,My 1 + u1,My 1 )

2(u1,0 + uMx ,1 + uMx 1,0 + uMx ,My 1 ) + 4uMx ,0 ,

(
u
1,k+1 + uMx 1,k+1 + uMx 1,k1 + u1,k1 )

2(u + u

1,k
Mx ,k+1 + uMx 1,k + uMx ,k1 ) + 4uMx ,k ,

(u1,1 + uMx 1,1 + uMx 1,My 1 + u1,My 1 )

2(u1,My + uMx ,1 + uMx 1,My + uMx ,My 1 ) + 4uMx ,My ,

2u

j = 0, k = 0,
j = 0, 1 k My 1,
j = 0, k = My ,
1 j Mx 1, k = 0,
1 j Mx 1, 1 k My 1,

(2.10c)

1 j Mx 1, k = My ,
j = Mx , k = 0,
j = Mx , 1 k My 1,
j = Mx , k = My .

2u

We then rewrite Eq. (2.8) as (1 + i)( x2 + y2 ) = ut Ru + (1 + i)|u|q1 u at (xj , yk , t) and then apply Eq. (2.9) to it. This gives


D2y
D2x D2y
D2x D2y
D2x
(1 + i) 2 + 2 +
+
u(xj , yk , t)
hx
hy
12h2x
12h2y




D2y


u(xj , yk , t)
D2
Ru(xj , yk , t) + (1 + i)|u(xj , yk , t)|q1 u(xj , yk , t) + O h4x + h4y .
= 1+ x +
12
12
t

(2.11)

Finally, we then use the CrankNicolson method similarly and obtain the iterative fourth-order in space compact nite difference
iterative scheme for solving 2D complex GinzburgLandau equation in Eq. (2.8) as

 n+1
D2y
D2x D2y
D2x D2y uj,k + unj,k
D2x
+ 2 +
+
2
h2x
hy
12h2x
12h2y

  n+1



uj,k unj,k
 n+ 1 q1 n+ 12
D2y
D2
n+ 1
= 1+ x +
Ruj,k 2 + (1 + i) uj,k 2 
uj,k ,
12
12



1 + i

 n+ 3
n+ 1
D2y
D2x D2y
D2x D2y uj,k 2 + uj,k 2
D2x
+ 2 +
+
2
h2x
hy
12h2x
12h2y

3
1


n+
n+

q1
uj,k 2 uj,k 2
D2y
D2x


n+1
n+1
n+1

= 1+
+
Ruj,k + (1 + i) uj,k 
uj,k ,
12
12



1 + i

(2.12a)

(2.12b)

where u0j,k = u0 (xj , yk ), and 0 j Mx , 0 k My , 0 n N 1. It should be noted that Eq. (2.12) can be written as linear
systems, which avoid any nonlinear iterations. Moreover, the scheme can be solved using an ADI method. For instance, Eq. (2.12)
can be solved using a PeacemanRachford ADI scheme [41,42] as follows:










D2y
D2y 
 D2x
 D2y
  n+ 12 q1 n+ 12

D2 

n+ 12

n
u

+
1
+
i
1
+
Ru
=
1
+

1
+
i

u
u
,
1 + x 1 + i
u
j,k
j,k
 j,k 
j,k
j,k
12
2 h2x
12
2 h2y
2
12

(2.13a)


1+

D2y
12



1 + i

D2y
2 h2y


un+1
j,k








D2y
 D2x
  n+ 12 q1 n+ 12

D2x 

n+ 12

+ 1 + i
1+
Ruj,k 1 + i uj,k 
= 1+
uj,k ,
uj,k +
12
2 h2x
2
12
(2.13b)

where uj,k is an intermediate mesh function, which is equivalent to Eq. (2.12) if

are O( 2 ) and O(h2x h2y ). Similarly, we can construct an ADI scheme for Eq. (2.12).

we ignore those terms whose truncation errors

Y. Yan et al. / Applied Mathematics and Computation 260 (2015) 269287

273

3. Stability analysis
For stability, we would like to show that if solution A and solution B are different periodic solutions obtained by Eq. (2.7)
based on two different initial conditions, then the difference between these two solutions is controlled by their initial difference.
We will employ the discrete energy method [28,29,3138] to analyze the stability of the present scheme in Eq. (2.7). To this end,
we will rst introduce the vector and matrix notations, obtain some lemmas, show the numerical solution to be bounded (as
seen in Theorem 1), and nally use these properties and results to show the stability (as seen in Theorem 2).
We now introduce the vector and matrix notations as

T

un = un0 , un1 , . . . , unM1 ,

T
1
n+ 1
n+ 1
n+ 1
un+ 2 = u0 2 , u1 2 , . . . , uM12
,


|un |q1 un = |un0 |q1 un0 , . . . , unM1 q1 unM1



 n+ 12 q1 n+ 12
u

u

,
T




 n+ 12 q1 n+ 12
 n+ 12 q1 n+ 12
= u0 
u0 , . . . , uM1 
uM1 ,

where (...)T is the transpose of the vector (...), and

10 1
1 10

..
1
..
.
H=
.
12
.
.
.
1

..

..

10

..

10

MM

It should be pointed out that because of the periodic mesh function, un0 = unM , we do not include unM in the vectors.
Thus, the present scheme in Eq. (2.7) can be written into vector form as

t un = (1 + i)H1 x2

un+1 + un
2


t un+ 2 = (1 + i)H1 x2
1




1
1 q1
1

+ Run+ 2 (1 + i) un+ 2  un+ 2 ,
1

un+ 2 + un+ 2
2

(3.1a)


+ Run+1 (1 + i)(|un+1 |q1 un+1 ).

(3.1b)

It can be seen that if H = I, an identity matrix, then Eq. (3.1) will reduce to a second-order CrankNicolson scheme.
We dene the inner product of two vectors and norms as

u , v
=
n

M1


unj vnj

= (u

n T

vn ,

 u  u  2
n

= h

j=0

u n 

 
 
= max unj  ,
0jM1

12

1
 n 2
= un , hun
2 ,
uj 

M1


(3.2a)

j=0

un p = h

1p

 n p
,
uj 

M1


p 1,

(3.2b)

j=0

where vnj stands for the conjugate of complex vnj . To analyze the stability of the present scheme in Eq. (3.1), we need the following
six lemmas.
Lemma 1 ([36,37]). For any periodic mesh functions un , vn , it holds that

x2 un , vn = x un , x vn
.

(3.3)

Lemma 2. For any periodic mesh function un , it holds that

 


Re H1 x2 un , hun = H1 x2 un , hun = Qx un 2 ,

Re H1 x2

!
un+1 + un
, t hun = 21 (Qx un+1 2 Qx un 2 ),
2

where Q = Chol(H1 ), the Cholesky factorization, and Re denotes the real part of complex value.

(3.4a)

(3.4b)

274

Y. Yan et al. / Applied Mathematics and Computation 260 (2015) 269287

Proof. The proof is similar to that in [31]. Since H is a real, symmetric and positive denite matrix, H1 is also a real, symmetric
and positive denite matrix. By the Cholesky factorization, there exists a real upper-triangular matrix Q such that H1 = QT Q.
Hence, we obtain from Lemma 1 that

 

H1 x2 un , hun = x2 H1 un , hun

= x H1 un , hx un

= H1 x un , hx un

= QT Qx un , hx un

= Qx un , hQx un

= Qx un 2 .

(3.5)

Here, we have used two facts that (1)

x u = H
n

=
=

 n n
T 
T
u1 , u2 , . . . , unM un0 , un1 , . . . , unM1

h

T

T
H1 un1 , un2 , . . . , unM H1 un0 , un1 , . . . , unM1
h

x H1 un ,

(3.6)

and similarly, H1 x2 un
obtain from Eq. (3.5)
Since

= x2 H1 un ; and (2) Qun , Qun

= (Qun )T Qun = (un )T QT Q un = (QT Qun )T un =
that Re H1 x2 un , hun
= H1 x2 un , hun
= Qx un 2 and hence Eq. (3.4) holds.

QT Qun , un
. Thus, we

!
un+1 + un
, t hun
2

1  2 1 n+1
=
x H (u + un ), h(un+1 un )
2
1
=
x H1 (un+1 + un ), hx (un+1 un )

2
1
=
H1 x (un+1 + un ), hx (un+1 un )

2
1
=
QT Qx (un+1 + un ), hx (un+1 un )

2
1
=
Qx (un+1 + un ), hQx (un+1 un )

2
1
=
Qx un+1 + Qx un , hQx un+1 hQx un

2
1
1
=
[ Qx un+1 , hQx un+1
Qx un , hQx un
]
[ Qx un , hQx un+1
Qx un+1 , hQx un
]
2
2
1
1
=
[Qx un+1 2 Qx un 2 ]
[ Qx un , hQx un+1
Qx un+1 , hQx un
]
2
2

H1 x2

(3.7)

and the fact that Re( u, v

v, u
) = 0 for any vectors u and v, we obtain

Re H1 x2

un+1 + un
, t hun
2

1
1
[Qx un+1 2 Qx un 2 ]
[Re( Qx un , hQx un+1
Qx un+1 , hQx un
)]
2
2
1
=
(Qx un+1 2 Qx un 2 ),
2

(3.8)

and hence Eq. (3.4) holds. 

Lemma 3. The following inequalities hold


Re

z1
z2

z
z2

| 1| =

|z1 |
, Re(z1 z2 ) |z1 z2 | = |z1 ||z2 |;
|z2 |

||z1 | |z2 || |z1 z2 | |z1 | + |z2 |;

(3.9a)

(3.9b)

Y. Yan et al. / Applied Mathematics and Computation 260 (2015) 269287

275

||z1 |n z1 |z2 |n z2 | |z1 | + |z2 | n |z1 z2 |, n 0;

| un , hvn
| un  vn 
" n
"
" u + vn "2
"
"
"
"
2

u n 2
2

u n 2
2

vn 2

vn 2
2

(3.9c)

(3.9d)

(3.9e)

where z1 and z2 are complex numbers.

Proof. We prove only for Eq. (3.9) and omit the detailed proofs for others since they are straightforward. It can be seen that
Eq. (3.9) holds if |z1 | = |z2 |. For the case |z1 | > |z2 | , we obtain

 


 
 |z1 |n z1 |z2 |n z2   |z1 |n (z1 z2 ) + |z1 |n |z2 |n z2 

=


 

z1 z2
z1 z2


 

n
n
n
 |z1 | (z1 z2 )   |z1 | |z2 | z2 




 + 

z1 z2
z1 z2


n
n
|z1 |n |z1 z2 | |z1 | |z2 |  |z2 |
+
=
|z1 z2 |
|z1 z2 |

n
n
z

z
|
|
|
|
|z2 |
1
2
|z1 |n +
|z1 | |z2 |
n1


n
k
nk1
= |z1 | + |z2 |
| z1 | | z2 |
=

 n


k=0

|z1 |k |z2 |nk

k=0

(|z1 | + |z2 |)n ,

and hence Eq. (3.9) holds. Similarly, we can prove that Eq. (3.9) holds for the case |z1 | < |z2 |.
Lemma 4 ([42,43]). For any
n 1

un p (u 
with =

1
2

un

(3.10)


and integer p 2, it holds that

x un  + un )

(3.11)

1p , where is a constant independent of p and h.

Lemma 5. For any x 0, y 0 and integer n 1, it holds that

(x + y)n 2n1 (xn + yn ).

(3.12)

Proof. We use the mathematical induction to prove it. It is obviously

that Eq. (3.12) is true for n = 1. Assume that Eq. (3.12)

holds for n up to k. Then, we have (x + y)k (x + y) 2k1 xk + yk (x + y) . Moreover, it can be seen that

2k1 (xk + yk )(x + y) 2k (xk+1 + yk+1 )

= 2k1 (xk+1 + yk+1 + xk y + xyk ) 2k (xk+1 + yk+1 )
= 2k1 [(xk y + xyk ) (xk+1 + yk+1 )]
= 2k1 [xk (y x) + yk (x y)]
= 2k1 [(y x)(xk yk )]

k1

k1
l kl1
=2
(y x)(x y) x y

= 2k1 (x y)2

l=0
k1


xl ykl1

l=0

0,

(3.13)

implying that (x + y)k+1 2k1 (xk + yk )(x + y) 2k (xk+1 + yk+1 ). By the induction, we conclude that Eq. (3.12) holds for any
integer n 1. 

276

Y. Yan et al. / Applied Mathematics and Computation 260 (2015) 269287

Lemma 6. For H1 = QT Q, where Q is a real upper-triangular matrix, then there exists two positive constants Ca and Cb such that

"
"
Ca un  "Qun " Cb un 

(3.14)

Proof. Note that H1 is a real, symmetric and positive denite matrix. By Schurs Lemma [44], there exists a unitary matrix U
T ) such that U1 H1 U = D , where D is a real diagonal matrix and the diagonal elements are the eigenvalues of H1 .
(i.e., U1 = U
Thus, we have H1 = UDU1 and

Qun 2 = h Qun , Qun

= h QT Qun , un

= h H1 un , un

= h UDU1 un , un

= h DU1 un , U1 un

= h Dwn , wn
= h

M1


dj |wj |2 ,

(3.15)

j=0

implying that dmin wn 2 Qun 2 dmax wn 2 , where wn = U1 un , dj is the eigenvalue of H1 . Since U is a unitary matrix,
we have wn 2 = Uun 2 = un 2 , and hence dmin un 2 Qun 2 dmax un 2 . By the denition of matrix H, we obtain that
1
(10 + 2 cos jM ), j = 0, . . . , M 1, indicating that 23 j 1. Since the eigenvalues of H and H1
the eigenvalues of H are j = 12
are reciprocal, we obtain that dmin 1 and dmax 32 . Hence, we may choose Ca = 1 and Cb = 32 . 
Theorem 1. For the scheme in Eq. (3.1), its solution satises these priori estimates as follows:

"

"

1
un  C , ""un+ 2 "" C ;

(3.16a)

"

"

1
x un  C , ""x un+ 2 "" C ;

"

(3.16b)

"

1
un  C , ""un+ 2 "" C .

(3.16c)

where 0 n T , C , C , C are constants which are independent on both h and .

Proof. We use the mathematical induction method to prove it. For n = 0, from the initial condition we should have u0  C ,
1

x u0  C , and u0  C . Since u 2 must be obtained using another method, we may choose a stable numerical method to
1
1
1
obtain the solution at t 1 . Thus, we should have u 2  C , x u 2  C , u 2  C . Assume that Eq. (3.16) holds for n up to
2

m 1 T 1. We would like to show Eq. (3.16) to be true for n = m.

1
m
m1
. This gives
To estimate um  and um+ 2 , we rst take an inner product of Eq. (3.1) at n = m 1 with h u +u
2

t um1 , h

 m
!

!
!
um + um1
um + um1
u + um1
um + um1
1
= (1 + i) H1 x2
,h
+ R um 2 , h
2
2
2
2
m
m1 !
u +u
1
1
.
(1 + i) |um 2 |q1 um 2 , h
2

Taking the real part of the above equation and using Lemmas 2 and 3, one may obtain

!
um + um1
2


!
!
 m
um + um1
um + um1
u + um1
1
= Re (1 + i) H1 x2
,h
+ RRe um 2 , h
2
2
2



m
m1 !
 m 12 q1 m 12 u + u
u
,h
+ Re (1 + i) u

2



"
 m
!
2
m1 "
 m 1 um + um1 
"
"
u +u

"

2 ,h
Re (1 + i) "
+
R
u

Q

" x
"

2
2
 
!

 
um + um1 
1 q1
1
um 2 , h
+ |(1 + i)|  um 2 

2
"
""
" #
"
" m
"

"
m1
" m 1 q1 m 1 " " um + um1 "
"
u +u
" m 12 " "
2
"
"
"
"
"
"
2
2
u
R "u
""

""
" + 1 + "u
"
2
2

1
(um 2 um1 2 ) = Re
2

t um1 , h

(3.17)

Y. Yan et al. / Applied Mathematics and Computation 260 (2015) 269287


2 12
M1
"
" " um + um1 "
  m 1 q1 m 1 
" #
" m 12 " "
" + 1 + 2 h
= R "u
uj 2 
uj 2 
""
"
"


2
j=0

12
M1
"
" " um + um1 "
  m 1 2q
" #
" m 12 " "
2
2

"
"

= R "u
h
""
uj

"+ 1+
2
j=0

277

" m
"
" u + um1 "
"
"
"
"
2

" m
"
" u + um1 "
"
"
"
"
2

2q1 q
" m
"
M1
"
" " um + um1 "
  m 1 2q
" u + um1 "
" #
" m 12 " "
2
u

"
" + 1 + 2 h
"
= R "u
""
 j

"
"
"
"
2
2
j=0

"
" " m + um1 "
"
"q "
"
m1 "
" #
" m
1 ""u
"
m 12 " " u + u
" + 1 + 2 "
"
R "um 2 " "
"u
"
"
"
"
2q "
2
2


" m
" m
"2  $
"2 
"
"
"
"
m1
"
u +u
u + um1 "
1 + 2 " m 1 "2q "
R " m 1 "2 "
"
"
"
"
2
2
" +"
" +"
"u
"u
" +
"
2q
2
2
2
2


$

"
"
R "
1
1 + 2 "
" m 12 "2 1
" m 12 "2q 1 m 2 1 m1 2
2

"u
"u
" + um  + um1 2 +
" + u  + u
2q
2
2
2
2
2
2


"
"
"
2q
2
c1 "
1 "
"
" m 12 "
(3.18)
" + "um 2 " + um 2 + um1 2 ,
"u
2q
2

$
m
m1
m
m1
where c1 = max{R, 1 + 2 }. Here, we have used the fact that Re((1 + i)Qx ( u +u
)2 ) = Qx ( u +u2 )2 0. By
2
1

Lemma 4 with p = 2q and Lemma 5 as well as the assumption um 2  C ,x um 2  C , we have

"
"
"
"1 1 "
"1 1
" 2q
"
1 " 2 + 2q "
1 " 2 2q
1"
"
"
" m 12 "2q
+ "um 2 "
" "um 2 "
"x um 2 "
"u
2q
"
"
"
"
"
"
1 "q+1 "
1 "q1
1 "2q
"
"
2q 22q1 "um 2 "
+ "um 2 "
"x um 2 "

2q 22q1 Cq+1 Cq1 + C2q c2 .

(3.19)
1

Substituting Eq. (3.19) into Eq. (3.18) and then using the assumption um1  C , um 2  C , we obtain

"

"

"

"

1 2q
1 2
(um 2 um1 2 ) c1 ""um 2 "" + ""um 2 "" + um 2 + um1 2

2q



c1 c2 + C2 + C2 + um 2
= c3 + c1 um 2 ,

(3.20)

where c3 = c1 c2 + 2c1 C2 . Thus, we obtain

u m  2

1
[c3 + um1 2 ]
1 c1

2[c3 + C2 ],

(3.21)

if is small enough such that 1 c1

1
2

and 1, implying that um  is bounded independently on h and . Once um  is

1

bounded, we can obtain from Eq. (3.1) that um+ 2  is also bounded using a similar argument.
m+ 12

To estimate x um  and x u
gives
2

2 , we rst take the inner product of Eq. (3.1) at n = m 1 with ht um1 /(1 + i). This

!
um + um1
, ht um1
2
!


 (1 + i) 
R
1
 m 12 q1 m 12
um 2 , ht um1
+
u
, ht um1 .
u

(1 + i)
(1 + i)

1
um1 , ht um1
= H1 x2
(1 + i) t

(3.22)

278

Y. Yan et al. / Applied Mathematics and Computation 260 (2015) 269287

Taking the real part of the above equation and then using Lemmas 2 and 4, we obtain

(1 + 2 )
=

t um1 2

(1 + 2 )

t um1 , ht um1

1
t um1 , ht um1

(1 + i)

 m
!




u + um1
R
1
= Re H1 x2
um 2 , ht um1
, ht um1 + Re
2
(1 + i)

!
(1 + i)  m 12 q1 m 12
u
, ht um1
+ Re
u

(1 + i)

 
 
!
q1
 


  m 1

1
R
m 12
m1 
  u 2 , ht um1  +  (1 + i)   um 12 
u
,
h

(Qx um 2 Qx um1 2 ) + 
t
 (1 + i)  

2
(1 + i) 
$
"
"
"
"

2
q1
R
1+ "
1
1 "
" m 12 "
m1
"um 12 
um 2 "


(Qx um 2 Qx um1 2 ) +
"u
" t um1  +
"
" t u
2
1 + 2
1 + 2
$

"
"
" 1
1
R
1 + 2 "
" m 12 "2q 2
" m 12 "

t um1 
(Qx um 2 Qx um1 2 ) +
"u
"u
" t um1  +
"
2q
2
1 + 2
1 + 2
"
1
R2 "
1
1 2
 t um1 2

(Qx um 2 Qx um1 2 ) + ""um 2 "" + 

2
2
2 1 + 2


"
1 + 2 "
1
" m 12 "2q
 t um1 2 .
+
(3.23)
" + 
"u
2q
2
2 1 + 2
= Re

Thus, moving the term 21 (Qx um 2 Qx um1 2 ) to the left-hand-side in Eq. (3.23) and dropping the term (1+1 2 ) t um1 2 ,
we simplify Eq. (3.23) to

"

"2

"

"2q

(Qx um 2 Qx um1 2 ) R2 ""um 2 "" + (1 + 2 ) ""um 2 ""

1

2q

R C + (1 + )c2 ,
2

implying that Qx

um 2

Qx

um1 2

+ R2 C

+ (1 + 2 )c

(3.24)
if 1. By Lemma 6, we obtain

1
Qx um 2
Ca
1

(Qx um1 2 + R2 C + (1 + 2 )c2 )
Ca

1 2
C x um1 2 + R2 C + (1 + 2 )c2

Ca b

1 2 2
C C + R2 C + (1 + 2 )c2 c5 .

Ca b

x um 2

(3.25)

"
"
1 "2
"
Using a similar argument, we may obtain from Eq. (3.1) that "x um+ 2 " c5 .
Finally, we have by Lemmas 4 and 5 that

um 2 um  2 x um  2 + um 


2 2

um 2 + x um 2
2

&2

2 2 (um x um  + um 2 )

 2

C + C2
m 2
2
2
+ u  2
+ C c6 .
2

(3.26)

Using a similar argument, we also obtain um+ 2  c6 , Thus, by the mathematical induction, Eq. (3.16) holds for any n, and
hence we have completed the proof. 
Theorem 2. Assume that ua and ub are two different periodic solutions obtained by Eq. (3.1) based on two different initial conditions.
1

Letting en = (ua )n (ub )n and en+ 2 = (ua )n+ 2 (ub )n+ 2 , then en and en+ 2 satisfy

"
"
" 1 "2

" n+ 32 "2
" "
"e
" + en+1 2 C "e 2 " + e0 2 ,

for any n and small , where C is a constant, implying that the scheme is unconditionally stable.

(3.27)

Y. Yan et al. / Applied Mathematics and Computation 260 (2015) 269287

279

Proof. It can be seen from Eq. (3.1) that for any n > 0, en and en+ 2 satisfy

t en = (1 + i)H1 x2

en+1 + en
2


t en+ 2 = (1 + i)H1 x2
1






1 q1
1
1 q1
1
1

(ua )n+ 2 (ub )n+ 2  (ub )n+ 2 ,
+ Ren+ 2 (1 + i) (ua )n+ 2 

en+ 2 + en+ 2
2



+ Ren+1 (1 + i) |(ua n+1 )|q1 (ua )n+1 |(ub )n+1 |q1 (ub )n+1 ,

n+1 +en

where e0 = (ua )0 (ub )0 and e 2 = (ua ) 2 (ub ) 2 . Taking the inner products of Eq. (3.28) with h e
tively, we obtain

t en , h

'

t e

n+ 3
2

and h e

en+ 2 + en+ 2
,h
2

'

en+ 2 + en+ 2
2

en+ 2 + en+ 2
,h
2

'

n+ 1
2

, respec-

(3.29a)

en+ 2 + en+ 2
= (1 + i)
+ R e ,h

2
(
'
3
1
q1
q1


en+ 2 + en+ 2


(1 + i) (ua )n+1 
.
(ua )n+1 (ub )n+1  (ub )n+1 , h
2
H1 x2

(3.28b)

+e
2

!

!
!
 n+1
en+1 + en
en+1 + en
en+1 + en
e
+ en
1
= (1 + i) H1 x2
,h
+ R en+ 2 , h
2
2
2
2
!




n+1
q1
q1
e
+ en
1 
1
1
1

,
(ua )n+ 2 (ub )n+ 2  (ub )n+ 2 , h
(1 + i) (ua )n+ 2 
2
3

n+ 12

(3.28a)

n+1

(3.29b)

Since (ua )n+ 2  C , (ub )n+ 2  C by Theorem 1, we may simplify the jth component of the nonlinear vector in
Eq. (3.29) based on Eq. (3.9) as




q1
q1
 

 

n+ 12 
n+ 12
n+ 12 
n+ 12
(
u
)
(
u
)

(
u
)
(
u
)
 b


  a

a
b

j
j
 
 q1 


1
1 



n+ 1 
n+ 1 
(ua )n+ 2 (ub )n+ 2 
(ua )j 2  + (ub )j 2 


j
j



q1  n+ 1 

 e 2 
2C

j





1
= c7  en+ 2  ,
j

(3.30)

q1

where c7 = 2C
. Similarly, we have





q1
q1
 



 
n+1 
n+1 
n+1
n+1 
(ua )
(ub ) 
(ub )
 (ua ) 
 c7  en+1  .

j
j
j

(3.31)

Taking the real part of Eq. (3.29), using Lemmas 2 and 3, and then using Eq. (3.30), we obtain

1
(en+1 2 en 2 )
2
"
 n+1
"2
#
"
" " en+1 + en "
"
" " en+1 + en "
"
"
"
e
+ en "
" n+ 12 " "
" n+ 12 " "
2
"
"
"
"
"
Q

+
R
1
+

+
c
"e
""
"e
""
7
" x
"
"
"
"
2
2
2
"
" n+1
"
 #
"
+ en "
1 ""e
"
"
c7 1 + 2 + R "en+ 2 " "
"
"
2

$
"
"2 
"
"
n+1
+ en "
c7 1 + 2 + R " n+ 1 "2 "
e
"
"
2

" +"
"e
"
2
2
$

"


c7 1 + 2 + R "
" n+ 12 "2 1

en+1 2 + en 2
" +
"e
2
2
$


"
"
2
c7 1 + + R " n+ 1 "2
n+1 2
n 2
2

+

e

+
e
.


"
"e
2

(3.32)

280

Y. Yan et al. / Applied Mathematics and Computation 260 (2015) 269287

which can be simplied to

"
"
"2
"
1 2
(en+1 2 en 2 ) c8 ""en+ 2 "" + ""en+1 "" + en 2 ,

(3.33)

$

where c8 =
1 + 2 c7 + R . Using a similar argument for Eq. (3.29), we obtain

"

"

"
"
"
"
"
"
3 "2
1 "2
"
" n+ 32 "2 " n+ 12 "2
"
"e
" "e
" c8 "en+ 2 " + "en+ 2 " + en+1 2 .

(3.34)

Summing Eqs. (3.33) and (3.34) together gives

"

"
"
"2 "
1 "2
"
" n+ 32 "2 "
2
"e
" + "en+1 " "en+ 2 " en 

"

"
"
"
"2
"
3 "2
1 "2
"
"
c8 "en+ 2 " + 2 "en+1 " + 2 "en+ 2 " + en 2
"

"
"
"2 "
"
3 "2
1 "2
"
"
c9 "en+ 2 " + "en+1 " + "en+ 2 " + en 2 ,
1

(3.35)

where c9 = 2c8 , implying that

(En+1 En ) c9 (En+1 + En ),

(3.36)

where En = en+ 2 2 + en 2 . Thus, we obtain from Eq. (3.36) that

1 + c9
En
1 c9



1 + c9 n+1 0

E
1 c9

En+1

e2 (n+1)c9 E0
c10 E0 ,

(3.37)

and (n + 1) T if is suciently small such that 1 c9

where c10 =
proof of Theorem 2 is completed. 
e2Tc9

1
2.

Hence, Eq. (3.27) has been obtained and the

Furthermore, the stability analysis for the 2D scheme is similar to the above analysis, but much more complicated. We omit
the detailed derivations here because of the limitation of pages.

4. Numerical examples
To test the accuracy of our numerical schemes, we rst considered a 1D complex GinzburgLandau equation with initial and
periodic boundary conditions as follows:

2 2
u
2u
, t (0, 1];
= (1 + i) 2 + Ru (1 + i)|u|2 u, x 0,
t
x
R
u(x, 0) =

2R i( 2R x)
e 2 ,
2


x 0,

2 2
;



2 2
2R i( (+ )R t)
2
e
u(0, t) = u
,
,t =
2
R

2R

(4.1a)

(4.1b)

(4.1c)

( + )R

where the exact solution is u (x, t) = 22R ei( 2 x 2 t).

We employed our present scheme in Eq. (2.7) to solve the above problem, in which the Thomas algorithm [45] was used for
solving the obtained tri-diagonal linear systems. The maximum of l -norm errors of the numerical solutions, as compared with
the analytical solution, were computed for 0 t 1 based on the formula

, h = max

0n 1





max unj u(xj , tn ) .

0jM

(4.2)

To obtain the convergence rate with respect to the spatial variable, we may assume that e ( , h) = O( p + hq ). Thus, e (2q/p , 2h)


= O[(2q/p )p + (2h)q ] = 2q O p + hq . Consequently, e (2q/p , 2h)/e ( , h) = 2q and hence q = log2 [e (2q/p , 2h)/e ( , h)]

Y. Yan et al. / Applied Mathematics and Computation 260 (2015) 269287

281

Table 1
1
Maximum error and convergence rate in the rst example (R = 1, = = 1, and u 2 was obtained
using Eq. (4.3)).

(h, )

e (h, )

Rate

(22 /5, 1/5)

(22 /10, 1/20)
(22 /20, 1/80)
(2 2 /40, 1/320)

0.00485726115701
2.946157986007180e-004
1.820597872425982e-005
1.139391147327777e-006

4.0432
4.0164
3.9981

Table 2
1
Maximum error and convergence rate in the rst example (R = 2, = = 1, and u 2 was obtained
using Eq. (4.3)).

(h, )

e (h, )

Rate

(2 /5, 1/5)
(2 /10, 1/20)
(2 /20, 1/80)
(2 /40, 1/320)

0.01610696193415
0.00106759559593
6.710783491111773e-005
4.230900881952949e-006

3.9152
3.9917
3.9874

Table 3
1
Maximum error and convergence rate in the rst example (R = 1, = = 1, and u 2 was obtained
using the analytical solution).

(h, )

e (h, )

Rate

(22 /5, 1/5)

(22 /10, 1/20)
(22 /20, 1/80)
(2 2 /40, 1/320)

0.00489551131950
2.943305960467438e-004
1.819984324406890e-005
1.139197836597871e-006

4.0560
4.0154
3.9980

Table 4
1
Maximum error and convergence rate in the rst example (R = 2, = = 1, and u 2 was obtained
using the analytical solution).

(h, )

e (h, )

Rate

(2 /5, 1/5)
(2 /10, 1/20)
(2 /20, 1/80)
(2 /40, 1/320)

0.01550766547724
0.00105786656594
6.696133491033216e-005
4.228439981965586e-006

3.8738
3.9817
3.9851

is the convergence rate with respect to the spatial variable. For our present scheme, we expect to have p = 2 and
q = 4.
Since the present scheme needs those values at time levels n = 0 and n = 12 , the values at the time level n = 12 must be
obtained using other methods. Here, we used the fully implicit method developed in [29] as

h2 2

1+
12 x

u 12 u0
j
j

/2

= (1 + i)

2
x

uj2 + u0j
2


u 21 + u0
2
j
+ R 1 + h x2 j

12
2
1
2

h2 2 |uj
(1 + i)(1 +
)
12 x

|q1 + |u0j |q1 uj2 + u0j

2

(4.3)

Note that the analytical solution is known in this example, we also used the exact values at the time level n = 12 to check if our
scheme is fourth-order. In our computation, we chose two different values of R. Results are shown in Tables 14. It can be seen
from these tables that the convergence rate of the scheme is approximately fourth-order with respect to the spatial variable,
which coincides with the theoretical analysis. Fig. 1 shows the solution proles obtained based on three different meshes as
compared with the analytical solutions. Results do not show much of a difference between the analytical solution and the
numerical solution.

282

Y. Yan et al. / Applied Mathematics and Computation 260 (2015) 269287

1
--- analytical solution
present scheme, 10 grid points
o present scheme, 20 grid points
+ present scheme, 40 grid points

0.8
0.6
0.4

Re(u)

0.2
0
-0.2
-0.4
-0.6
-0.8
-1

x
(a) R =1

--- analytical solution

present scheme, 10 grid points
o present scheme, 20 grid points
+ present scheme, 40 grid points

1
0.8
0.6
0.4

Re(u)

0.2
0
-0.2
-0.4
-0.6
-0.8
-1
0

x
(b) R=2

Fig. 1. Proles of the real part of the numerical solution at t = 1 obtained using three different meshes and corresponding = 1/20, 1/80, 1/320, respectively,
as compared with the analytical solution in the rst example.

Table 5
1
Maximum error and convergence rate in the second example (R = 1, = = 1, and u 2 was obtained using the analytical solution).




hx , hy ,
e hx , hy ,
Rate



2 2 /5, 2 2 /5, 1/5



2 2 /10, 2 2 /10, 1/20



2 2 /20, 2 2 /20, 1/80



2 2 /40, 2 2 /40, 1/320

0.00110542709763

6.406423411185432e-005

4.1089

3.996877803767565e-006

4.0026

2.482725014645710e-007

4.0089

We then considered a 2D complex GinzburgLandau equation with initial and periodic boundary conditions as

2 2
u
2u 2u
= (1 + i)
+
+ Ru (1 + i)|u|2 u, x, y 0,
, t (0, 1];
t
x2 y2
R
u(x, y, 0) =

2R i( R x+ R y)
2
e 2
,
2


x, y 0,

2 2
;

(4.4a)

(4.4b)

Y. Yan et al. / Applied Mathematics and Computation 260 (2015) 269287

283

Table 6
1
Maximum error and convergence rate in the second example (R = 2, = = 1, and u 2 was obtained using the analytical solution).




hx , hy ,
e hx , hy ,
Rate


2 /5, 2 /5, 1/5


2 /10, 2 /10, 1/20


2 /20, 2 /20, 1/80


2 /40, 2 /40, 1/320

0.01239497829524
8.097105045953942e-004
4.978160430930223e-005
3.117868510497469e-006

3.9362
4.0237
3.9970

t=3

t=1

0.14
0.08
0.12

1
0.07

0.8

0.8

0.1

0.06
0.6
0.05
0.4

|u|

|u|

0.6

0.08
0.4

0.04

0.06

0.2

0.2
0.03
0
20

0.02
10

0.04

0
20
10

20
10

20

-20

-10

-10
-20

-10
-20

0.02

10

-10

0.01

(a)

-20

(b)

t=5

t=7
0.8
0.4

1
0.8

0.25

0.4

|u|

|u|

0.6

0.8

0.3

0.6

0.7

0.35

0.6

0.5

0.4

0.4

0.2

0.3

0.2
0.2
0.15
0
20

0
20

0.1
10
10

0.05

-10

0.2
10

20

20

-20

-10

-10
-20

(c)

0.1

10

0
-10
-20

-20

(d)

Fig. 2. Solution distributions at various times at (a) t = 1, (b) t = 3, (c) t = 5, and (d) t = 7 obtained based on the mesh 40 40 and = 0.1 in the third example.

u(x, 0, t) = u(x, Ly , t) =

2R i( R x (+ )R t)
2
e 2
,
2

u(0, y, t) = u(Lx , y, t) =

2R i( R y (+ )R t)
2
,
e 2
2

t [0, 1]

(4.4c)

t [0, 1],

(4.4d)

( + )R

where the analytical solution is u(x, y, t) = 22R ei( 2 x+ 2 y 2 t). For simplicity, we used the exact solution for u 2 in the
computation. Results are shown in Tables 5 and 6. Again, as expected, it can be seen from these two tables that the convergence
rate of the scheme is approximately fourth-order.
R

284

Y. Yan et al. / Applied Mathematics and Computation 260 (2015) 269287

1

1
-- pseudo-spectral method
o present scheme

0.8

0.8

0.7

0.7

0.6

0.6

0.5
0.4

0.5
0.4

t=7
t=5
t=3
t=1

0.3

t=7
t=5
t=3
t=1

0.3

0.2

0.2

0.1

0.1

0
-20

-15

-10

-5

10

15

0
-20

20

t=7
t=5
t=3
t=1

-15

t=7
t=5
t=3
t=1

-10

-5

(a)

(b)

10

15

20

1
-- pseudo-spectral method
o present scheme

0.9

-- present scheme

0.9

0.8

0.8

0.7

0.7

0.6

0.6

|u|

|u|

-- pseudo-spectral method
o present scheme

0.9

|u|

|u|

0.9

0.5
0.4

0.5
0.4

t=7
t=5
t=3
t=1

0.3

t=7
t=5
t=3
t=1

0.3

0.2

0.2

0.1

0.1

0
-20

-15

-10

-5

10

15

20

0
-20

t=7
t=5
t=3
t=1

-15

t=7
t=5
t=3
t=1

-10

-5

(c)

10

15

20

(d)

Fig. 3. Solution proles along the y-axis at various times obtained based on four meshes of (a) 20 20, (b) 30 30, (c) 40 40, (d) 50 50, and = 0.1 in the
third example.

Finally, we considered a more complex example with no periodic boundary condition, which we could not nd the exact
solution, as follows:

u
2u 2u
= (1 + i)
+
+ u (1 + i)|u|2 u, x, y (, +) , t > 0;
t
x2 y2

(4.5a)

2
u(x, y, 0) =

(4.5b)

(x + iy)e(x +y ), x, y (, +).
2

2
2
The example is related to the non-equilibrium condensate. For the coecient of the diffusion term x2u + y2u is 1 + i, we expect
to have some diffusions in the solution. In our computation, the domain was taken to be 20 x, y 20, where the boundary
condition was set to be zero. The number of grid points in both x and y was chosen to be 40 40 with the time step = 0.1.
1

For this case, u 2 was computed using the fourth-order accurate and explicit pseudo-spectral method [39] built in the software
MATLAB. Fig. 2 shows the simulation of the solution at various times in 0 t 7. As expected, we see from Fig. 2 that the vortex
grows and diffuses toward the boundary.
We then chose four different meshes of 20 20, 30 30, 40 40, and 50 50 with the time step = 0.1 and compared
with the fourth-order accurate and explicit pseudo-spectral method, as shown in Figs. 3 and 4. It can be seen from these two
gures that the solutions obtained based on these two methods are not signicantly different for the meshes of 20 20, 30 30,
and 40 40, as shown in Figs. 3(a)(c) and 4(a)(c). However, for the mesh of 50 50, the pseudo-spectral method produces a
divergent solution in which the maximum value of u(x, y) = 4.670916294752263e180 + i8.621058722691569e180 at t = 0.8 and
u(x, y) = NaN at t = 0.9. On the other hand, our method still produces a stable solution, as shown in Figs. 3(d) and 4(d). This
indicates that our method has a better stability condition.

Y. Yan et al. / Applied Mathematics and Computation 260 (2015) 269287

1

1
-- pseudo-spectral method
o present scheme

0.8

0.8

0.7

0.7

0.6

0.6

0.5
0.4

0.5
0.4

t=7
t=5
t=3
t=1

t=7
t=5
t=3
t=1

0.3

0.2

0.2

0.1

0.1

0
-20

-- pseudo-spectral method
o present scheme

0.9

|u|

|u|

0.9

0.3

-15

-10

-5

10

15

0
-20

20

t=7
t=5
t=3
t=1

-15

t=7
t=5
t=3
t=1

-10

-5

(a)

10

15

20

(b)
-- pseudo-spectral method
o present scheme

-- present scheme

0.9

0.8

0.8

0.7

0.7

0.6

0.6

|u|

|u|

0.9

0.5
0.4

0.5
0.4

t=7
t=5
t=3
t=1

t=7
t=5
t=3
t=1

0.3

0.2

0.2

0.1

0.1

0
-20

0.3

285

-15

-10

-5

10

15

20

0
-20

t=7
t=5
t=3
t=1

-15

t=7
t=5
t=3
t=1

-10

-5

(c)

(d)

10

15

20

Fig. 4. Solution proles along the x-axis at various times obtained based on four meshes of (a) 20 20, (b) 30 30, (c) 40 40, (d) 50 50, and = 0.1 in the
third example.

5. Conclusion
In this study, we have developed a new, simple, and accurate nite difference iterative scheme for solving the 1D and 2D
complex GinzburgLandau equations with initial and periodic boundary conditions, respectively. Coupled with the Crank
Nicolson nite difference technique and the fourth-order compact nite difference method for spatial variables, the new
scheme is proved to be unconditionally stable and provides fourth-order accurate numerical solutions with respect to the
spatial variables. Numerical errors and convergence rates of the solutions have been tested by several examples. Results
show that the maximal l -norm errors are small as expected, and the convergence rates of the numerical solutions are
fourth-order with respect to the spatial variables. Further research will be focused on the applications of the new method
to practical physics and engineering problems, such as the phenomenological GinzburgLandau complex superconductivity
model [22]:

286

Y. Yan et al. / Applied Mathematics and Computation 260 (2015) 269287

+A
+ ik +
t
k

A
+ + curl curlA+R
t

+ ||2 = 0, in  (0, T );

i
+ A
k

(5.1a)

= 0, in  (0, T );

(5.1b)

with the initial and boundary conditions

i
+ A n= 0, curlA =H, on  (0, T );
k


x,0 = 0 (x) ,

 
A x,0 = A0 (x) ,

in ;

(5.1c)

(5.1d)

where Eq. (5.1) may be solved using the present method. Here, is a complex valued function and is referred to as the
2
order parameter so that || gives the relative density of the superconducting electron pairs, and the normal and the pure
2
2
superconducting states are characterized accordingly by || = 0 and || = 1. stands for the complex conjugate of . A is a
real vector potential for the total magnetic eld and is a real scalar function called electric potential. H is the applied magnetic
eld that points out of the (x1, x2 )-plane. , k are positive constants which are related to the known physical quantities.
Acknowledgments
The research was supported by a grant from NASA EPSCoR & LaSPACE, Louisiana (NASA/LEQSF (2013)-RAP-2) and was partially
supported by the National Science Foundation (NSF EPS 1003897). The authors appreciate referees for their valuable suggestions.
References
[1] M. Lin, W. Duan, Wave packet propagating in an electrical transmission line, Chaos Solitons Fractals 24 (2005) 191196.
[2] F. Ndzana, A. Mohamadou, T. Kofane, Modulated waves and chaotic-like behaviours in the discrete electrical transmission line, J. Phys. D: Appl. Phys. 40
(2007) 32543262.
[3] E. Kengne, R. Vailancourt, Propagation of solitary waves on lossy nonlinear transmission lines, Int. J. Modern Phys. B 1 (2009) 118.
[4] E. Kengne, C. Tadmon, T. Nguyen-Ba, R. Vaillancourt, On the dissipative complex Ginzburg-Landau equation governing the propagation of solitary pulses in
dissipative nonlinear transmission lines, Chin. J. Phys. 47 (2009) 8090.
[5] E. Kengne, R. Vaillancourt, 2D Ginzburg-Landau system of complex modulation for coupled nonlinear transmission lines, J. Infrared Millimeter Waves 30
(2009) 679699.
[6] J. Duan, P. Holmes, Fronts, domain walls and pulses in a generalized Ginzburg-Landau equations, Proc. Edinb. Math. Soc. 38 (1995) 7797.
[7] J.M. Soto-Crespo, N. Akhmediev, Exploding solition and front solutions of the complex cubic-quintic Ginzburg-Landau equation, Math. Comput. Simul. 69
(2005) 526536.
[8] E.N. Tsoy, A. Ankiewicz, N. Akhmediev, Dynamical models for dissipative localized waves of the complex Ginzburg-Landau equation, Phys. Rev. E 73 (2006)
036621.
[9] N.K. Efremidis, D.N. Christodoulides, K. Hizanidis, Two-dimensional discrete Ginzburg-Landau solitions, Phys. Rev. A 76 (2007) 043839.
[10] W.H. Renninger, A. Chong, F.W. Wise, Dissipative solitons in normal-dispersion ber lasers, Phys. Rev. A 77 (2008) 023814.
[11] Q. Du, M.D. Gunzburger, J.S. Peterson, Analysis and approximation of the Ginzburg-Landau model of superconductivity, SIAM Rev. 34 (1992) 5481.
[12] Q. Du, M.D. Gunzburger, J. Deang, Finite element approximation of a periodic Ginzburg-Landau model for type-II superconductors, Numer. Math. 64 (1993)
85114.
[13] Q. Du, Discrete gauge invariant approximations of time dependent Ginzburg-Landau model of superconductivity, Math. Comput. 67 (1998) 965986.
[14] Q. Du, R.A. Nicolaides, X. Wu, Analysis and convergence of covolume approximation of the Ginzburg-Landau model of superconductivity, SIAM J. Numer.
Anal. 35 (1998) 10491072.
[15] J. Deang, Q. Du, M.D. Gunzburger, Modeling and computation of random thermal uctuations and material defects in the Ginzburg-Landau model for
superconductivity, J. Comput. Phys. 181 (2002) 4567.
[16] Q. Du, L. Ju, Numerical simulations of the quantized vortices on a thin superconducting hollow sphere, J. Comput. Phys. 201 (2004) 511530.
[17] Q. Du, L. Ju, Approximations of a Ginzburg-Landau model for superconducting hollow spheres based on spherical centroidal voronoi tessellations, Math.
Comput. 74 (2005) 12571280.
[18] A. Doelman, Finite-dimensional models of the Ginzburg-Landau equation, Nonlinearity 4 (1991) 231250.
[19] Q. Wang, Z.D. Wang, Simulating the time-dependent dx2 y2 Ginzburg-Landau equations using the nite-element method, Phys. Rev. B 54 (1996) R15645.
[20] G.J. Lord, Attractors and inertial manifolds for nite-difference approximations of the complex Ginzburg-Landau equation, SIAM J. Numer. Anal. 34 (1997)
14831512.
[21] P. Takac, A. Jungel, A nonstiff Euler discretization of the complex Ginzburg-Landau equation in one space dimension, SIAM J. Numer. Anal. 34 (1997)
292328.
[22] Z. Chen, S. Dai, Adaptive Galerkin methods with error control for a dynamical Ginzburg-Landau model in superconductivity, SIAM J. Numer. Anal. 38 (2001)
19611985.
[23] J. Willers, E.H. Twizell, A nite-difference solution of the Ginzburg-Landau equation, J. Differ. Eq. Appl. 9 (2003) 10591068.
[24] Y. Zhang, W. Bao, Q. Du, Numerical simulation of vortex dynamics in Ginzburg-Landau-Schrdinger equation, Eur. J. Appl. Math. 18 (2007) 607630.
[25] P. Degond, S. Jin, M. Tang, On the time splitting spectral method of the complex Ginzburg-Landau equation in large time and space scale, SIAM J. Sci. Comput.
30 (2008) 24662487.
[26] L. Zhang, Long-time behavior of nite difference approximations for the two-dimensional complex Ginzburg-Landau equation, Numer. Funct. Anal. Optimization 31 (2010) 11901211.
[27] G.Z. Tsertsvadze, On the convergence of difference schemes for the Kuramoto-Tsuzuki equation and for systems of reaction diffusion type, J. Comput. Math.
Math. Phys. 31 (1991) 698707.
[28] Z.Z. Sun, Q. Zhu, On Tsertsvadzes difference scheme for the Kuramoto-Tsuzuki equation, J. Comput. Appl. Math. 98 (1998) 289304.
[29] X. Hu, Numerical methods on compact difference schemes for 1D nonlinear Kuramoto-Tsuzuki equation, Numer. Methods Partial Differ. Eq., in press.
[30] S.K. Lele, Compact nite difference schemes with spectral-like resolution, J. Computat. Phys. 103 (1992) 1642.
[31] T. Wang, B. Guo, Unconditional convergence of two conservative compact difference schemes for non-linear Schrdinger in one dimension, Sci. Sin. Math.
41 (2011) 207233.

Y. Yan et al. / Applied Mathematics and Computation 260 (2015) 269287

287

[32] T. Wang, B. Guo, Q. Xu, Fourth-order compact and energy conservative difference schemes for the nonlinear Schrdinger equation in two dimensions, J.
Comput. Phys. 243 (2013) 382399.
[33] T. Wang, Optimal point-wise error estimate of a compact difference scheme for the coupled Gross-Pitaevskii equations in one dimension, J. Sci. Comput. 59
(2014) 158186.
[34] T. Wang, Optimal point-wise error estimate of a compact difference scheme for the coupled nonlinear Schrdinger equations, J. Comput. Math. 32 (2014)
5874.
[35] T. Wang, Optimal point-wise error estimate of a compact difference scheme for the Klein-Gordon-Schrdinger equation, J. Math. Anal. Appl. 412 (2014)
155167.
[36] M. Lees, Alternating direction and semi-explicit difference methods for parabolic partial differential equations, Numer. Math. 3 (1961) 398412.
[37] B. Guo, On discrete energy method (i), Chin. Univ. J. Numer. Math. 4 (1984) 327344.
[38] W. Dai, An unconditionally stable three-level explicit difference scheme for the Schrdinger equation with a variable coecient, SIAM J. Numer. Anal. 29
(1992) 174181.
[39] J. Yang, Nonlinear Waves in Integrable and Non-integrable Systems, SIAM, 2010.
[40] R.K. Mohanty, V. Gopal, High accuracy arithmetic average type discretization for the solution of two-space dimensional nonlinear wave equations, Int. J.
Model. Simul. Sci. Comput. 3 (2012) 1150005.
[41] D.W. Peaceman, H.H. Rachford Jr., The numerical solution of parabolic and elliptic differential equations, J. Soc. Ind. Appl. Math. 3 (1955) 2841.
[42] Z.Z. Sun, Numerical Methods of Partial Differential Equations, 2nd edition, Science Press, Beijing, 2012 (in Chinese).
[43] Y. Zhou, Application of Discrete Functional Analysis to the Finite Difference Methods, International Academic Publishers, New York, 1990.
[44] K.E. Atkinson, An Introduction to Numerical Analysis, 2nd edition, John Wiley & Sons, New York, 1989.
[45] K.W. Morton, D.F. Mayers, Numerical Solution of Partial Differential Equations, Cambridge University Press, London, 1994.