Simulating DNLS models
Mario Mulansky
University of Potsdam, Department of Physics and Astronomy

arXiv:1304.1608v1 [physics.comp-ph] 5 Apr 2013

Louisiana State University, Center for Computation and Technology (mmulansky@cct.lsu.edu)

Abstract
We present different techniques to numerically solve the equations of motion for the
widely studied Discrete Nonlinear Schrödinger equation (DNLS). Being a Hamiltonian
system, the DNLS requires symplectic routines for an efficient numerical treatment. Here,
we introduce different such schemes in detail and compare their performance and accuracy
by extensive numerical simulations.
Keywords: Symplectic Integration; Performance; DNLS; Hamiltonian System

1. Introduction
The Discrete Nonlinear Schrödinger Equation (DNLS) has been known to physicists,
biologists, chemists and mathematicians for more than 30 years now. It describes a simple
model of coupled unharmonic oscillators defined by the complex valued phase/amplitude
ψn ∈ C at lattice site n. The dynamics of the oscillators is governed by the DNLS that
writes in one spatial dimension:
i

∂
ψn = Vn ψn + ψn+1 + ψn−1 + β|ψn |2 ψn ,
∂t

(1)

where i is the imaginary unit, Vn ∈ R is the local potential and β ∈ R the nonlinear
strength. These equations of motion can be derived from the following Hamilton function:
H=

X

∗
Vn |ψn |2 + ψn+1
ψn + ψn+1 ψn∗ +

n

β
|ψn |4 .
2

(2)

The Hamiltonian character immediately gives an integral of motion which is usually
called energy and denoted as E := P
H = const. However, the DNLS possesses another
conserved quantity, the norm N := n |ψn |2 , usually set to unity N = 1. The existence
of two conserved quantities makes this equation particularly challenging for numerical
approaches, as will be explained later.
This model first appeared as a description of polarons by Holstein [21]. Later, the
DNLS was used by Davydov in his studies of protein dynamics [13, 49] as well as in
the context of local modes in small molecules [48]. In recent years, however, two new
applications of the DNLS gained increasing attention: Bose-Einstein-Condensates (BEC)
and coupled optical wave guides.
Preprint submitted to Journal of Computational Physics

April 8, 2013

In the context of Bose-Einstein condensates, the DNLS describes the mean field approximation of a weakly interacting Bose gas at zero temperature as shown e.g. in [35].
The predictions of this description have been, to some extent, experimentally verified
for condensates in optical traps and a periodic potential, e.g. [2, 8]. A very interesting
phenomenon that can be studied in BECs is Anderson localization, which denotes the
“absence of diffusion” [3] in disordered systems. Several experimental realizations of this
effect have been reported, e.g. [47, 46, 27, 4, 45]. The corresponding DNLS accompanied
by a random local potential is often called the Discrete Anderson Nonlinear Schrödinger
Equation (DANSE) and has been subject of heavy numerical studies in the past years.
Moreover, the DNLS gives a very accurate description of the propagation of light in
optical waveguides with a nonlinear material [25, 11]. As for BECs, optical waveguides
were also used to study Anderson localization by imposing a disordered potential. Especially the fact that the nonlinear strength can be controlled quite precisely in terms of
the properties of the nonlinear material makes this system very appealing for studying
Anderson localization and nonlinearity experimentally. This has been done recently by
Schwartz et al. in [47]. A more extensive review on the history and applications of the
DNLS can be found in [14].
The increasing experimental accessibility of systems described by the DNLS has led to
numerous numerical investigations of such models. Special attention has been given to the
case where Anderson localization, induced by disorder, is accompanied by nonlinearity. In
terms of the DNLS this can be modeled by choosing the local potential Vn to be a random
variable, typically as independent and identically distributed (iid) from some interval of
size W : Vn ∈ [−W/2, W/2]. A large number of numerical investigations studied the slow
destruction of Anderson localization by the nonlinear interactions between the localized
eigenmodes [43, 16, 39, 28, 24, 52]. In all these works, subdiffusive spreading of initially
localized modes is reported up to computationally accessible times, recently confirmed
by an experimental study [32]. Additionally, similar spreading laws have been observed
in quasi-periodic nonlinear lattices [29] and the nonlinear Stark ladder [17, 26]. Recent
results on the scaling properties of chaos [42] and on the scaled spreading in a reduced
model [38] raised some questions on the asymptotic validity of the observed subdiffusive
spreading. Several attempts to find a theoretical description of the spreading have been
proposed, mainly based on an effective noise theory [16], but a full understanding of
the interplay between disorder and nonlinearity is still lacking [15]. Lately, the theory
of chaotic diffusion has been developed [37] and successfully applied to a different class
of systems with disorder and nonlinearity [40, 41], but can possibly be extended to the
DNLS/DANSE models as well. Here, however, we are not going to further address
the very interesting problem of asymptotic spreading, but rather focus on the different
algorithmic techniques that are used to obtain those numerical results.
This article is organized as follows: After this introduction, we will review the general
properties of symplectic integrators and the idea of operator splitting to construct such
methods in section 2. In section 3 we introduce the different numerical approaches for
simulating the DNLS model and in section 4 we compare their performance. Finally, we
end with our conclusions in section 5.

2

2. Symplectic Schemes
Solving the DNLS model here means to find an approximate solution ψ(t) of the
equations of motion (1) from an initial condition ψ(t = 0). Being a Hamiltonian system,
the DNLS has a symplectic flow map and exhibits conservation of energy E = const.
This requires the usage of symplectic routines that preserve the symplectic nature of
the system [31, 19]. Formally, such a solution can be written in terms of the Liouville
operator etLH , defined via Poisson brackets [18] and acting on the initial condition ψ(t) =
etLH ψ(0). If the action of this operator is known explicitly, the system is called integrable
and the solution ψ(t) can be written analytically. In most cases, however, this solution
can only be approximated by numerical methods.
A very common technique to find a symplectic scheme for a non-integrable system is
the so-called operator splitting [34]. For this, the Hamilton function has to be separable,
which means it can be written as H = A + B, where the action of the operators etLa and
etLB are known explicitly. Then, the integration of the solution from time t to t + ∆t for
some small time step ∆t can be approximated by:
e∆tLH = e∆t(LA +LB ) = e∆tLA e∆tLB + O(∆t2 ).

(3)

The splitting scheme above is the most simple and only accurate in first order of ∆t. In
a more general way, the approximation can be written as a product of j operators:
e∆tLH =

j
Y

eai ∆tLA ebi ∆tLB + O(∆tp+1 ).

(4)

i=1

The operators eai ∆tLA and ebi ∆tLB represent the exact integrations according to Hamiltonians A and B over times ai ∆t and bi ∆t respectively. The parameters ai , bi are chosen
such that the resulting product is an exact representation of e∆tLH of order ∆tp . Such a
splitting scheme effectively computes a trajectory not of the original system with Hamiltonian H = A + B, but of a new Hamiltonian K = A + B + O(∆tp ). Thus, the symplectic property of the dynamics is preserved. Moreover, the computed trajectory, obtained
from many subsequent steps (4), conserves a new energy Ẽ := K = E + O(∆tp ), which
means that along this trajectory the original energy is not exactly conserved, but only
fluctuates around its exact value with a magnitude of O(∆tp ). With a non-symplectic integration scheme this is not the case, as there the numerical error of the energy increases
at every step by an amount O(∆tp ), which accumulates along the trajectory. Hence,
especially for long integration times the symplectic routines are essential as they allow
reasonable energy conservation even for large time steps ∆t which makes them much
more efficient than non-symplectic methods. Thus, symplectic integrators have become
standard for Hamiltonian systems and several splitting routines of different orders have
been developed in recent years, see e.g. [53, 7, 33, 10, 9, 30].
In the following, we will use two second order (p = 2) symmetric splittings: the
SBAB1 and SBAB2 scheme [30]. Symmetric splittings are especially interesting as they
naturally lead to self-adjoint algorithms. The time evolution according to those splittings

3

for a separated Hamilton function H = A + B are defined as follows:
SBAB1 : e∆tLH = e∆t/2LB e∆tLA e∆t/2LB + O(∆t3 )

(5)

2

K = A + B + O(∆t · |B|)
∆tLH

SBAB2 : e

= eb1 ∆tLB ea1 ∆tLA eb2 ∆tLB ea1 ∆tLA eb1 ∆tLB + O(∆t3 )
2

(6)

2

K = A + B + O(∆t · |B| ),
with a1 = 1/2 , b1 = 1/6 and b2 = 2/3 in the last line and K denoting the formal
Hamiltonian integrated exactly by the given scheme. The symmetric composition ensures
that these schemes are identical to their inverse with negative timestep: SBAB1/2 (∆t) =
(SBAB1/2 (−∆t))−1 , a property denoted as self-adjoint.
3. Numerical Methods
Here, we will present different algorithms to numerically calculate approximate trajectories of the DNLS equation. At first we focus on methods based on the two-part
operator splitting described above, all of which have a different approach to solve the
non-local term in (1). The first method uses Fourier transform for the non-local part and
has been applied to the DNLS in [51, 16, 28, 24], it will be abbreviated as “FT” here.
The second algorithm employs an implicit Crank-Nicolson scheme and has been used
in [36, 39] to integrate the DNLS, it will be referred to as “CN” throughout this text.
The third method is the so-called “PQ” scheme introduced in [5], where the non-local
part is again separated into two integrable operators. Finally, a new multi-symplectic
scheme is introduced based on the ideas of Bridges and Reich [6] that was already applied
to a DNLS without local potential [54] and will be called Euler-Box “EB” scheme here.
3.1. Fourier Method
The obvious way to treat the coupling term in the DNLS is to use a spectral method.
This involves forward and backward Fourier transforms at every time step and is thus
called “FT” in this text. It has already been applied to the DNLS as described in [51].
For this scheme, the Hamilton function (2) is split into two parts H = AF T + BF T with:
X
∗
AF T =
ψn−1 ψn∗ + ψn ψn−1
n

BF T =

X
n

Vn |ψn |2 +

β
|ψn |4 .
2

(7)

Here, A contains only the coupling and B consists of the linear and nonlinear local
potential. This implies that the action of the time evolution operator e∆tLAF T becomes
local when applied to the Fourier transform of the state ψ̃q . Specifically, e∆tLAF T acts
as follow [51]:

N
X



ψ̃
=
ψn e2πiq(n−1)/N

q



n=1

e∆tLAF T : ψ̃q0 = e−2i∆t cos(2π(q−1)/N ) ψ̃q
(8)



N


1 X 0 −2πin(q−1)/N

 ψn0 =
ψ̃ e
,

N n=1 q
4

where the first and the last operation are forward- and backward Fourier transforms of
the state ψn and the second step is the time evolution in Fourier space. The action of
e∆tLBF T can be written explicitly as it is fully local in real space:
2

e∆tLBF T : ψn0 = e−i∆t(Vn +β|ψn | ) ψn .

(9)

Having found two explicit representations of the two parts of the time evolution one
can now construct integrators by consecutively applying e∆tLAF T and e∆tLBF T . The
simplest form is given by the first order approximation e∆tLH = e∆tLAF T e∆tLBF T +
O(∆t2 ). As described above, symmetric second order schemes can be obtained by using
the SBAB1 or SBAB2 splitting (5), (6). A numerical study on the behavior of the energy
error for the simple first order and the two second order splitting schemes of the FT
method is shown in Figure 1. We compared the three splitting methods by integrating the
DNLS system (1) with a random potential chosen independent and identically distributed
Vn ∈ [−2, 2], a nonlinear strength of β = 1 and N = 128 lattice sites – the DNLS with a
random potential is usually referred to as DANSE model. As initial condition
we chose
P
|ψn |2 = 1. We
a Gaussian in the lattice center with width σ = 10, ensuring N =
performed the numerical time evolution using each of the above schemes up to the time
T = 100 and repeated for decreasing step sizes ∆t = 1.0 . . . 10−6 . As quantification for
the accuracy we calculate the mean square energy error ∆E and norm ∆N along the
trajectory:
s
p
p
1 X
2
(Em − E0 )2 ,
∆N = h(Nm − N0 )2 i
(10)
∆E = h(Em − E0 ) i =
M m
where the index m denotes the time and E0 is the exact energy and N0 = 1 the exact
norm, fixed to unity in these simulations. Figure 1 shows the behavior of the energy error
for the different splittings for decreasing stepsize. One clearly observes the second order
behavior of the SBAB1 and SBAB2 , compared to the first order trivial splitting (3).
Additionally, the SBAB2 splitting shows a significantly smaller energy error than the
simpler ABC and SBAB1 schemes. For stepsizes smaller than ∆t ∼ 10−4 the finite double
precision starts limiting the accuracy of the computation, hence smaller stepsizes ∆t <
10−4 are not advisable for the SBAB1/2 FT schemes. Note, that the FT based splitting
schemes exhibit numerically exact norm conservation so no analysis of ∆N is required.
3.2. Crank-Nicolson (CN)
As the Schrödinger equation is of such great importance, many numerical methods
have been developed to find trajectories in this system. One of the most important approaches is the “Crank-Nicolson” scheme [12], a second order symplectic scheme to treat
the usual (linear) Schrödinger equation, also used for solving other partial differential
equations like the wave equation [44]. To use this scheme, we will split the Hamilton
function (2) into a linear and a nonlinear part and then use the ideas of Crank and
Nicolson to solve the linear part. Hence, we write H = ACN + BCN with:
X
∗
ACN =
Vn |ψn |2 + ψn−1 ψn∗ + ψn ψn−1
n

BCN =

Xβ
n

2

(11)
|ψn |4 .
5

∆E

10-1
10-3
10-5
10-7
10-9
10-11
10-13
10-1510-6

FT
FT SBAB1
FT SBAB2

10-5

10-4

10-3

10-2

∆t

∼ ∆t2

10-1

100

Figure 1: Behavior of the energy error for the FT scheme and different splitting methods.
The simulation was done for a DNLS with N = 128 lattice sites, β = 1 and a random
potential Vn ∈ [−2, 2]. ∆E is the averaged energy error along the trajectory (10) up to
a total integration time T = 100.
The Liouville operators for the two individual terms are then defined by their actions
on ψ as:
ψ̇n = (LACN ψ)n = −i(ψn+1 + ψn−1 + Vn ψn )
2

ψ̇n = (LBCN ψ)n = −iβ|ψn | ψn .

(12)
(13)

The time evolution operator e∆tLBCN is nonlinear, but can then be written explicitly due
to its local character:
2
e∆tLBCN : ψn0 = e−i∆tβ|ψn | ψn ,
(14)
where ψ 0 denotes the time evolution of one time step ∆t according to the Hamilton
function B only.
Unfortunately, e∆tLACN can only be written explicitly for a constant or periodic potential. As we want to treat the generic case here, we will use the Crank-Nicolson scheme,
a second order, norm preserving approximative method. Applying it to the linear part
LA (12) (we skip the index CN her for simplicity), the exponential is approximated by
1+LA ∆t/2
, which leads to an implicit equation for the evolved
Caley’s formula: e∆tLA ≈ 1−L
A ∆t/2
0
state ψ :
e∆tLA : (1 − LA ∆t/2)ψ 0 = (1 + LA ∆t/2)ψ.
(15)
From (12) can be seen that LA is a band matrix with −iVn on the diagonal and −i
on the upper/lower sub-diagonal, this means (15) is a set of linear equations which
can be solved by a Gauß-algorithm. Moreover, due to the band structure of LA the
computational requirement is only linear in the system size.
Using the symplectic integrators for the two parts ACN and BCN one can now construct a symplectic scheme as described in section 2. Namely, we have implemented the
usual Crank-Nicolson splitting scheme of order p = 1: CN after (3), the second order
6

10-1
10-3

∆E

10-5
10-7
10-9
CN
CN SBAB1
CN SBAB2

10-11
10-13
10-1510-6

10-5

10-4

10-3
∆t

10-2

∼ ∆t2

10-1

100

Figure 2: Behavior of the energy error for the CN scheme and different splitting methods.
The simulation setup is as described in Figure 1. The dashed line shows exemplarily the
second order behavior and is located at the same position as in Figure 1 for comparison.
CNSBAB1 scheme following (5) and the CNSBAB2 scheme from (6). Remember that all
these schemes are symplectic and norm preserving. We also note that for weak nonlinearities the last method CNSBAB2 represents a considerable improvement over CNSBAB1 ,
as the error gets significantly smaller when |BCN |  1. In a typical spreading setup as
chosen in the numerical examples of this article, the total error in the energy reduces by
about one order of magnitude, as seen in Figure 2.
The results are presented in Figure 2 in terms of the energy error ∆E, as the norm
again is preserved exactly. The plots illustrate the respective p = 1, 2 behavior of the
methods, as well as the improvement of the CNSBAB2 over the more simple CNSBAB1 .
For both methods we found an optimal stepsize ∆t ≈ 10−4 that leads to minimal energy
errors of ∆E ≈ 10−12 . For stepsizes below these values, the finite double precision again
leads to a linear increase of the error. Furthermore, the error of the CNSBAB2 is clearly
smaller than for the FTSBAB2 shown in Figure 1, as easily seen from comparing with the
dashed black lines that are plotted at the same position in both graphs.
3.3. PQ-Splitting
While in the methods above Fourier transforms or the Crank-Nicolson scheme were
used as an integrator of the linear coupling term in H, Bodyfelt et al. presented a new
idea to treat the nonlocal term – the PQ scheme [5]. They proposed to apply another
operator splitting of the linear part such that the final evolution operator consists of
three separate mappings.
To derive this method, we first have to divide the state into real and imaginary part:
ψn = an + ibn . The Hamilton function (2) is then:
H=

X

Vn (a2n + b2n ) + 2(an an+1 + bn bn+1 ) +

n

7

β 2
(a + b2n )2 .
2 n

(16)

This can then be split into the three integrable parts:
B=

X

Q=

X

Vn (a2n + b2n ) +

n

β 2
(a + b2n )2
2 n

2an an+1

(17)

n

P =

X

2bn bn+1 .

n

As above, the time evolution operator can now be approximated:
PQ : e∆tH = e∆tLB e∆tLQ e∆tLP + O(∆t2 ),

(18)

where with defining αn = 2Vn + β(a2n + b2n ) the individual steps are given as:
(
a0 = an cos(αn ∆t) + bn sin(αn ∆t)
e∆tLB : 0n
bn = bn cos(αn ∆t) − an sin(αn ∆t)
(
a0 = an
e∆tLQ : 0n
bn = bn − 2∆t(an−1 + an+1 )
(
a0 = an + 2∆t(bn−1 + bn+1 )
e∆tLP : 0n
bn = bn ,

(19)

(20)

(21)

This concatenation of three solvable evolution operators defines a first order symplectic scheme which is an integrator of the formal Hamiltonian K = B + P + Q + O(∆t).
However, as pointed out in section 2, the error of the splitting can be improved by using
better approximations than (18). The simplest symmetric second order splitting is given
by:
PQABC : e∆tLH = e

∆t
2 LB

e

∆t
2 LP

e∆tLQ e

∆t
2 LP

e

∆t
2 LB

+ O(∆t3 )

(22)

2

K = B + P + Q + O(∆t · |B||P ||Q|).
Surprisingly, this obvious three-operator concatenation has been discussed only very
recently by Skokos et al. [50] where it also was applied to the DNLS. It was called
“ABC” scheme there, so we denote it as PQABC here. Besides that, one can also use the
two-operator splitting presented above and apply that two times. This idea was used
in [5] where the Hamiltonian H is split in two steps as follows: first H = A + B followed
by A = P + Q. Applying the Leapfrog method SBAB1 to the two splittings we get the
following scheme:
PQSBAB1 : e∆tLH = e

∆t
4 LB

e

∆t
2 LP

e

∆t
4 LB

e∆tLQ e

∆t
4 LB

e

∆t
2 LP

e

∆t
4 LB

+ O(∆t3 )

(23)

2

K = B + P + Q + O(∆t · |B||P ||Q|).
Applying the SBAB2 method individually to the two splittings also results in a second
8

∆N

10-1
10-3
10-5
10-7
10-9
10-11
10-13
10-1510-6

PQ
PQ ABC
PQ SBAB1
PQ SBAB2
EB

10-5

10-4

10-3

10-2

∆t

∼ ∆t2

10-1

100

Figure 3: Behavior of the norm error for different splitting methods for the PQ scheme
and the Euler-Box method. The simulation setup is as described in Figure 1. The dashed
line shows exemplarily the second order behavior and is located at the same position as
in Figures 1 and 2 for comparison.
order symmetric scheme given by the following 13 successive simple mappings:
2

2

PQSBAB2 : e∆tH = ed1 ∆tLB ed1 c2 ∆tLP ec2 ∆tLQ ed2 c2 ∆tLP ec2 ∆tLQ
2

× ed1 c2 ∆tLP ed2 ∆tLB ed1 c2 ∆tLP ec2 ∆tLQ
d2 c2 ∆tLP c22 ∆tLQ d1 c2 ∆tLP d1 ∆tLB

×e

e

e

e

(24)
3

+ O(∆t )

K = B + P + Q + O(∆t2 · |B|2 |P |2 |Q|)
with d1 = 1/6, d2 = 2/3 and c2 = 1/2 [30]. This is again symplectic and has order p = 2.
For small B this would lead to reduction of the error. However, for the splitting (17), the
norms of the operators B , P , Q are all ∼ 1. Hence, it is not obvious that this scheme
leads to a better error behavior, but a numerical investigation shows indeed that the
energy error ∆E decreases by more than one order of magnitude when using the SBAB2
method compared to SBAB1 . This is shown in Figure 3, where we again integrated the
DNLS (1) as above with N = 128 sites, an idd. random potential Vn = −2 . . . 2 and
β = 1 up to the time T = 100 using decreasing time steps ∆t = 0.5 . . . 10−6 . As initial
conditions we again used a Gaussian at the center with width σ = 10. Because the
norm is not exactly conserved by the PQ scheme, we will here use ∆N as quantification
for the accuracy of the method for the sake of variety. One finds that the PQABC and
the PQSBAB1 produce very similar errors. The Euler-Box scheme, presented in the next
chapter, exhibits a slightly smaller error, but by far the best splitting is the PQSBAB2 .
Comparing Figures 2 and 3 one finds that the CN based schemes still produce the smallest
errors. However, the PQ schemes are computationally much simpler which allows smaller
time steps ∆t at the same CPU time compared to the CN methods. Hence, a detailed
performance study is required to identify which method is more efficient. This will be
done in section 4.
9

3.4. Multi-symplectic Euler-Box Scheme
The above schemes relied on the splitting of the Hamilton function into integrable
parts. In this section, we will present a different approach that is based on a multisymplectic formulation of the equations. Multi-symplectic schemes were introduced by
Bridges and Reich [6] for Hamiltonian PDEs and they are known to locally conserve
symplectivity. As the DNLS system can be thought of as a discretized version of a
continuous nonlinear Schrödinger equation, the application of the multi-symplectic theory
is reasonable [23]. Actually, this has been done already by Wang et al. [54] where a
multi-symplectic scheme for the DNLS without local potential (Vn = 0) was developed.
Here, we will adopt this scheme closely following the calculations in [54] for the general
DNLS (1). As result, we will obtain a locally implicit, self-adjoint multi-symplectic
scheme of second order.
To derive this scheme we have to start at the continuous wave equation for the complex
valued, space- and time-dependent wave function ψ(x, t):
iψt = ψxx + 2ψ + V (x)ψ + β|ψ|2 ψ.

(25)

Note, that from a straight forward spatial discretization with grid size ∆x = 1 one
immediately obtains (1). By setting ψ(x, t) = a + ib and v = ax , w = bx this can be
formulated as:
M zt + Kzx = ∇z S(z),
(26)
with

1
β
2+V 2
(a + b2 ) + (v 2 + w2 ) + (a2 + b2 )2
2
2
4
and M and K being anti-symmetric matrices and z containing the state:
 




a
0 0 −1 0
0 −1 0 0
b
0 0 0 −1
1 0 0 0

,

z=
K=
M =
v  .
1 0 0
0 0 0 0 ,
0
w
0 1 0
0
0 0 0 0
S=

(27)

(28)

This formulation emphasizes the multi-symplectic nature of the system. By introducing
a spatial and timely discretization xn , n = 1, 2, . . . and tm , m = 1, 2, . . . with grid size ∆x
and time step ∆t the discretized state is defined znm := z(xn , tm ). The multi-symplectic
Euler-Box scheme for this discrete state variable is then given as:
M+ δ + znm + M− δ − znm + K+ δ+ znm + K− δ− znm = ∇z S(znm ),

(29)

with δ +− , δ+− are the forward/backward difference operators acting on time and space,
m
e.g. δ + znm = (znm+1 − znm )/∆t and δ− znm = (znm − zn−1
)/∆x. M+ , M− and K+ , K− are
splittings of the symplectic structure matrices M = M+ + M− and K = K+ + K− , where
T
T
the conservation of symplectivity demands M+
= −M− and K+
= −K− . This matrix
splitting is not unique. However, in [54] it was found that all possible splittings lead to
one of two fundamental schemes. The first one is given by taking M+ and K+ as upper
triangular matrices:




0 −1 0 0
0 0 −1 0
0 0 0 0
0 0 0 −1
,

.
M+ = 
K
=
(30)
+
0 0 0 0
0 0 0
0
0 0 0 0
0 0 0
0
10

Submitting this into (29) we obtain, after substituting vnm and wnm , the following scheme:
∆t m
m
m
(a
− 2am
n + an−1 ) − ∆t(Vn + 2)an
∆x2 n+1
2
m 2 m
− ∆tβ((am
n ) + (bn ) )an
∆t m+1
m+1
= am
(b
− 2bm+1
+ bm+1
n +
n
n−1 ) + ∆t(Vn + 2)an
∆x2 n+1
+ ∆tβ((am+1
)2 + (bm+1
)2 )bm+1
.
n
n
n

bm+1
= bm
n
n −

am+1
n

(31)

This is a first order, multi-symplectic scheme [54, 22]. The scheme is locally implicit
as the second equation for am+1
is given implicitly and requires to solve a quadratic
n
equation. However, no implicit dependence on neighboring lattice sites is included and
hence the scheme does not involve solving a set of nonlinear equations which makes it
locally implicit. From here on, we will set ∆x = 1 which makes this method solving the
DNLS (1).
The second fundamental splitting is obtained by exchanging M+ ↔ M− while leaving
K+ and K− as is. This gives another scheme similar to the one above but with a locally
. It can be shown that the two schemes are
instead of am+1
implicit equation for bm+1
n
n
adjoint to each other with respect to time reversal [54]. The adjoint method of a one-step
method z m+1 = Φ∆t z m is defined as the inverse original map with reversed time −∆t,
that is Φ∗∆t := Φ−1
−∆t and thus can be obtained by solving the following equations for
z m+1 : z m = Φ−∆t (z m+1 ). Given the two adjoint schemes above by composition theory
a new second-order, self-adjoint method can be created by constructing the mapping
Φ∗∆t/2 Φ∆t/2 [31], which is given by the following equations:
m+ 21

an

m+ 21

bn

∆t m
m
m
(b
− 2bm
n + bn−1 ) − ∆t(Vn + 2)bn
2 n+1
2
m 2 m
− ∆tβ((am
n ) + (bn ) )bn
∆t m+ 12
m+ 1
m+ 1
m+ 1
= bm
(an+1 − 2an 2 + an−12 ) − ∆t(Vn + 2)an 2
n −
2
= am
n −

m+ 12 2

− ∆tβ((an

m+ 21 2

) + (bn

m+ 12

) )an

(32)

m+ 12

bm+1
= 2bn
n
am+1
n

− bm
n
∆t m+1
m+ 12
m+1
= an
−
(b
− 2bm+1
+ bm+1
n
n−1 ) − ∆t(Vn + 2)bn
2 n+1
− ∆tβ((am+1
)2 + (bm+1
)2 )bm+1
.
n
n
n

This scheme consists of two explicit (first and third line) and two locally implicit steps
(second and fourth line) and defines the multi-symplectic Euler-Box (EB) scheme. The
behavior of the energy error for this scheme is plotted in Figure 3 (black stars). As expected, it is second order accurate, and the error lies between the PQSBAB1 and PQSBAB2
schemes. However, the Euler-Box method requires to evaluate two root functions for each
lattice site at each step, which is computationally quite demanding. So again, only a detailed performance study can reveal if this method presents some advantage over the
other schemes introduced above. It should be mentioned that this scheme is still a particularly interesting approach for solving the DNLS, as it can be easily generalized to two11

10-4

10-1
10-3
10-5
10-7
10-9
10-11
10-13
10-1510-3

∆E

∆E

101
10-1
10-3
10-5
10-7
10-9
10-11
10-13

EB4
CN4

10-3

10-2
∆t

10-1

∼ ∆t4

100

(a) Fourth order Yoshida composition.

EB6
CN6

10-2

∆t

10-1

∼ ∆t6

100

(b) Sixth order Yoshida composition.

Figure 4: Energy error of the fourth (a) and sixth (b) order composition methods based
on the SBAB2 Crank-Nicolson (CN4/CN6) and the Euler-Box (EB4/EB6) scheme. The
numerical setup is the same as for the results above.
or even higher-dimensional lattices. A clear advantage over the other methods described
above, for which a generalization to higher dimensions is not straight forward.
3.5. High-Order Schemes
All the methods above were accurate to at most second order. Of course, one often
requires higher order methods, especially when a good accuracy is wanted. Fortunately,
there exists a generic way of constructing higher order methods from symmetric, second
order schemes such as the SBAB1, SBAB2 or the Euler-Box algorithm [31]. Many such
schemes have been developed, for example by Yoshida [53], but here we restrict ourselves
to presenting two higher-order methods, one of order four and one of order six. Assume
we have a symmetric second order scheme defined by the mapping ψ 0 = Φ(∆t)ψ with
the stepsize ∆t as parameter. In this context, a scheme is called symmetric if it is selfadjoint Φ(∆t) = Φ∗ (∆t), where the adjoint scheme is the inverse mapping with negative
stepsize: Φ∗ (∆t) = Φ−1 (−∆t). For the SBAB1/2 schemes it can be easily seen that they
are indeed symmetric, but this also holds for the Euler-Box scheme, by construction, as
explained in section 3.4 and in [54].
Following [53], we construct a fourth order scheme Φ4 (∆t) as follows:
Φ4 (∆t) = Φ(α1 ∆t)Φ(α0 ∆t)Φ(α1 ∆t),

(33)

with α0 = −21/3 /(2 − 21/3 ) and α1 = 1/(2 − 21/3 ) [53]. Similarly, a sixth order scheme
can be constructed:
Φ6 (∆t) = Φ(µ3 ∆t)Φ(µ2 ∆t)Φ(µ1 ∆t)Φ(µ0 ∆t)Φ(µ1 ∆t)Φ(µ2 ∆t)Φ(µ3 ∆t).

(34)

The parameters µ0...3 are not defined uniquely in this case. Here we will use the set
called “solution A” in [53]:
µ1 = −1.17767998417887

µ2 = 0.235573213359357

µ3 = 0.784513610477560

µ0 = 1 − µ1 − µ2 − µ3 .
12

CN SBAB2
FT SBAB2
PQ SBAB2
EB

102
101

103

CN
CN4
CN6
RK

CPU time (minutes)

CPU time (minutes)

103

102
101

100

10-110-8 10-7 10-6 10-5 10-4 10-3 10-2 10-1

0
1010
-1110-1010-9 10-8 10-7 10-6 10-5 10-4 10-3 10-2 10-1

(a) Performance of second order methods.

(b) Performance of higher order methods.

∆E

∆E

Figure 5: Performance of the different second order schemes (a) and higher order schemes
based on CNSBAB2 (b). These graph show the CPU time to integrate a system with
N = 1024 lattice sites up to T = 105 with a given averaged energy error ∆E on an Intel
Core i7, 2.93GHz machine. In (b), additionally the results for a non-symplectic 8-th
order Runge-Kutta scheme (RK) are plotted.
Note, that these methods of constructing higher order schemes based on the PQABC
scheme have recently been discussed by Skokos et al. in [50], where they also used the
DNLS with random potential as benchmarking example. In Figure 4 we examplarily show
the energy error for the schemes Φ4 (a) and Φ6 (b) based on the symplectic CNSBAB2
and the Euler Box schemes. The O(∆t4 ) and O(∆t6 ) behavior is clearly visible, again
with a smaller error for the CN based schemes.
4. Performance
After having introduced and discussed several numerical methods to solve the Discrete Nonlinear Schrödinger equation, we will now study the efficiency of the different
approaches. We will first consider only the second order schemes, as the most efficient of
these can then also be used to construct the most efficient higher order scheme. For all
of the methods described above the computational effort increases linearly with system
size1 , hence no study on different system sizes is required to compare the different methods. Note especially that although the Crank-Nicolson (CN) scheme involves solving a
linear system of N equations, which formally is of complexity O(N 2 ), one can make use
of the band structure of the coupling term LACN with only three non-zero diagonals and
thus implement this algorithm with complexity O(N ) as well.
We analyze the performance of the different methods by computing a trajectory
for a lattice with N = 1024 sites and a random potential Vn ∈ [−2, 2] and nonlinear
strength β = 1. This is the typical setup for studying spreading in the DANSE model [43,
16, 39]. We again start the time evolution from an initial Gaussian with width σ =
1 The FT method has an N log N behavior, so it might become more disadvantageous for very large
systems

13

10 and total norm N = 1. We compute a trajectory up to a time T = 105 with
decreasing stepsizes ∆t. As before, we use the mean squared energy error ∆E as given
in (10) as quantification of the accuracy, but here the energy is only computed every 10
timesteps, so the major computational effort lies in the time evolution. The results for
the different second order schemes are shown in Figure 5a, where the CPU time is plotted
in dependence of the mean energy error ∆E. The FTSBAB2 and the EB exhibit a similar
performance, while the PQSBAB2 scheme introduced by Skokos et al. [51] represents a
clear improvement. However, one gets the best performance when using the CNSBAB2
scheme. It requires about the same computational effort as the FTSBAB2 scheme, but
its different splitting gives a much smaller energy error, as already seen from comparing
Figures 1 and 2. We also checked the performance of other splittings (SBAB1 and ABC),
but the SBAB2 versions were always superior.
In Figure 5b we compare the performance of higher order schemes based on the most
efficient second order scheme CNSBAB2 . Namely, we show results for the Yoshida 4 and
Yoshida 6 compositions introduced in (33) and (34). For comparison, we also plot the
second order results for CNSBAB2 . Interestingly, the higher order schemes are only more
efficient when a high accuracy is required. If an accuracy ∆E > 10−5 is sufficient, the
second order scheme CNSBAB2 provides the best performance. If one requires smaller
errors, a higher order scheme should be chosen. In Figure 5b also the performance of a
non-symplectic 8-th order Runge-Kutta scheme with stepsize control [20] is presented,
as implemented in the Boost.Odeint library [1]. We found that for the total integration
time chosen here, T = 105 , this non-symplectic scheme is competitive and even the most
efficient for very high accuracy. However, one has to keep in mind that this scheme is
non-symplectic, hence the energy error ∆E will increase linearly in time. Thus, for very
long time scales T > 105 , the symplectic schemes will outperform the controlled RungeKutta scheme. For the DANSE model that was used here examplarily, the integration
times go up to T = 109 [36], where the higher order symplectic schemes will surely beat
the RK method.
5. Conclusions
We introduced and described several numerical schemes to compute approximate
trajectories for the Discrete Nonlinear Schrödinger equation (1). The DNLS is a very
important model with many applications, hence identifying the most efficient numerical
scheme is of great interest. For our numerical performance test we relied on the DANSE
setup that is used extensively in the past to study the interplay between disorder and
nonlinearity. Of the mostly used second order splitting schemes, the CNSBAB2 method
was found to show the best performance in terms of the least CPU time for a given energy
error ∆E. Hence, we conclude that this scheme is superior to the others and should be the
first choice for numerically treating DNLS models when a moderate precision is required.
A new and particularly interesting approach is the Euler-Box scheme as it can also be
used for higher dimensional systems. But the implicit treatment of the nonlinearity in
this case strongly restricts the possibility to treat different local nonlinear potentials (e.g.
higher powers). The FT, CN and PQ splitting schemes, on the other hand, can be used
for any local nonlinear potential in the DNLS and are thus more flexible in this regard.
Note, that all presented schemes here are suitable to additionally integrate the set of
14

linearized equations to obtain Lyapunov exponents. For the CN scheme this has been
already done in [36].
Finally, we showed that when a high accuracy is required, higher order schemes
provide better performance than the standard second order schemes. At an integration
time of T ≈ 105 , even non-symplectic schemes are competitive. But this changes with
increasing T . We note that performance results are always to be taken with care as they
might vary greatly for different compilers, operating systems or hardware. However, we
believe that our results are a helpful guide for which simulation code to choose when
dealing with DNLS models.
Acknowledgements
I thank A. Pikovsky and B. Mulansky for fruitful discussions. I also thank the CCT at
Louisiana State University for hospitality and financial support as well as the DAAD for
financial support under the project number 50015188. Finally, financial support under
the Project HPC-EUROPA2 (Project number 228398), with the support of the European
Community - under the FP7 “Research Infrastructure” Programme is acknowledged.
References
[1] K. Ahnert and M. Mulansky, Odeint – Solving ordinary differential equations in C++, AIP
Conference Proceedings, 1389 (2011), pp. 1586–1589, www.odeint.com.
[2] B. Anderson and M. A. Kasevich, Macroscopic quantum interference from atomic tunnel arrays,
Science, 282 (1998), pp. 1686–1689.
[3] P. W. Anderson, Absence of diffusion in certain random lattices, Physical Review Letters, 109
(1958), pp. 1492–1505.
[4] J. Billy, V. Josse, Z. Zuo, A. Bernard, B. Hambrecht, P. Lugan, D. Clément, L. SanchezPalencia, P. Bouyer, and A. Aspect, Direct observation of Anderson localization of matter
waves in a controlled disorder, Nature, 453 (2008), pp. 891–894.
[5] J. Bodyfelt, T. Laptyeva, C. Skokos, D. Krimer, and S. Flach, Nonlinear waves in disordered
chains: Probing the limits of chaos and spreading, Physical Review E, 84 (2011), p. 016205.
[6] T. J. Bridges and S. Reich, Multi-symplectic integrators: numerical schemes for Hamiltonian
PDEs that conserve symplecticity, Physics Letters A, 284 (2001), pp. 184–193.
[7] J. Candy and W. Rozmus, A symplectic integration algorithm for separable Hamiltonian functions,
Journal of Computational Physics, 92 (1991), pp. 230–256.
[8] F. Cataliotti, S. Burger, C. Fort, P. Maddaloni, F. Minardi, A. Trombettoni, A. Smerzi,
and M. Inguscio, Josephson junction arrays with Bose-Einstein condensates, Science, 293 (2001),
pp. 843–846.
[9] P. Channell and C. Scovel, Symplectic integration of Hamiltonian systems, Nonlinearity, 3
(1999), p. 231.
[10] S. A. Chin, Symplectic integrators from composite operator factorizations, Physics Letters A, 226
(1997), pp. 344–348.
[11] D. Christodoulides and R. Joseph, Discrete self-focusing in nonlinear arrays of coupled waveguides, Optics Letters, 13 (1988), pp. 794–796.
[12] J. Crank and P. Nicolson, A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type, Advances in Computational Mathematics, 6 (1996),
pp. 207–226.
[13] A. Davydov, The theory of contraction of proteins under their excitation, Journal of Theoretical
Tiology, 38 (1973), pp. 559–569.
[14] J. C. Eilbeck and M. Johansson, The discrete nonlinear schrödinger equation-20 years on, in
Conference on Localization and Energy Transfer in Nonlinear Systems, 2003, p. 44.
[15] S. Fishman, Y. Krivolapov, and A. Soffer, The nonlinear Schrödinger equation with a random
potential: results and puzzles, Nonlinearity, 25 (2012), p. R53.

15

[16] S. Flach, D. Krimer, and C. Skokos, Universal spreading of wave packets in disordered nonlinear
systems, Physical Review Retters, 102 (2009), p. 24101.
[17] I. Garcia-Mata and D. L. Shepelyansky, Nonlinear delocalization on disordered Stark ladder,
European Physical Journal B, 71 (2009), pp. 121–124.
[18] H. Goldstein, Classical Mechanics, 3/E, Pearson Education, 2002.
[19] E. Hairer, C. Lubich, and G. Wanner, Geometric Numerical Integration: Structure-Preserving
Algorithms for Ordinary Differential Equations, Springer Series in Computational Mathematics,
Springer London, Limited, 2006.
[20] E. Hairer, S. Nørsett, and G. Wanner, Solving Ordinary Differential Equations I: Nonstiff
Problems, Springer, 1993.
[21] T. Holstein, Studies of polaron motion: Part II. The small polaron, Annals of Physics, 8 (1959),
pp. 343–389.
[22] J. Hong, H. Liu, and G. Sun, The multi-symplecticity of partitioned Runge-Kutta methods for
Hamiltonian PDEs, Mathematics of Computation, 75 (2006), p. 167.
[23] A. Islas, D. Karpeev, and C. Schober, Geometric integrators for the nonlinear schrödinger
equation, Journal of Computational Physics, 173 (2001), pp. 116–148.
[24] M. Ivanchenko, T. Laptyeva, and S. Flach, Anderson localization or nonlinear waves: A matter
of probability, Physical Review Letters, 107 (2011), p. 240602.
[25] S. Jensen, The non-linear coherent coupler: A new optical logic device, in Digest of the Topical
Meeting on Integrated and Guided-Wave Optics, Incline Village, Nev., p. MB4-1, 1980.
[26] D. O. Krimer, R. Khomeriki, and S. Flach, Delocalization and spreading in a nonlinear Stark
ladder, Physical Review E, 80 (2009), p. 036201.
[27] Y. Lahini, A. Avidan, F. Pozzi, M. Sorel, R. Morandotti, D. N. Christodoulides, and
Y. Silberberg, Anderson localization and nonlinearity in one-dimensional disordered photonic
lattices, Physical Review Letters, 100 (2008), p. 013906.
[28] T. Laptyeva, J. Bodyfelt, D. Krimer, C. Skokos, and S. Flach, The crossover from strong
to weak chaos for nonlinear waves in disordered systems, EPL (Europhysics Letters), 91 (2010),
p. 30001.
[29] M. Larcher, T. V. Laptyeva, J. D. Bodyfelt, F. Dalfovo, M. Modugno, and S. Flach,
Subdiffusion of nonlinear waves in quasiperiodic potentials, arXiv preprint arXiv:1206.0833, (2012).
[30] J. Laskar and P. Robutel, High order symplectic integrators for perturbed Hamiltonian systems,
Celestial Mechanics and Dynamical Astronomy, 80 (2001), pp. 39–62.
[31] B. Leimkuhler and S. Reich, Simulating hamiltonian dynamics, vol. 14, Cambridge University
Press, 2005.
[32] E. Lucioni, B. Deissler, L. Tanzi, G. Roati, M. Zaccanti, M. Modugno, M. Larcher, F. Dalfovo, M. Inguscio, and G. Modugno, Observation of subdiffusion in a disordered interacting
system, Physical Review Letters, 106 (2011), p. 230403.
[33] R. I. McLachlan, On the numerical integration of ordinary differential equations by symmetric
composition methods, SIAM Journal on Scientific Computing, 16 (1995), pp. 151–168.
[34] R. I. McLachlan and G. R. W. Quispel, Splitting methods, Acta Numerica, 11 (2002), pp. 341–
434.
[35] R. Mishmash and L. Carr, Ultracold atoms in 1d optical lattices: mean field, quantum field,
computation, and soliton formation, Mathematics and Computers in Simulation, 80 (2009), pp. 732–
740.
[36] M. Mulansky, Localization properties of nonlinear disordered lattices. Diploma Thesis, University
of Potsdam, 2009.
[37] M. Mulansky, Chaotic Diffusion in Nonlinear Hamiltonian Systems, PhD thesis, University of
Potsdam, 2012.
[38] M. Mulansky, K. Ahnert, and A. Pikovsky, Scaling of energy spreading in strongly nonlinear
disordered lattices, Physical Review E, 83 (2011), p. 026205.
[39] M. Mulansky and A. Pikovsky, Spreading in disordered lattices with different nonlinearities,
Europhysics Letters, 90 (2010), p. 10015.
[40] M. Mulansky and A. Pikovsky, Scaling of energy spreading in strongly nonlinear lattices, arXiv
preprint arXiv:1205.3592, (2012).
[41] M. Mulansky and A. Pikovsky, Scaling properties of energy spreading in nonlinear Hamiltonian
two-dimensional lattices, Physical Review E, 86 (2012), p. 056214.
[42] A. Pikovsky and S. Fishman, Scaling properties of weak chaos in nonlinear disordered lattices,
Physical Review E, 83 (2011), p. 025201.
[43] A. Pikovsky and D. Shepelyansky, Destruction of Anderson localization by a weak nonlinearity,

16

Physical Review Letters, 100 (2008), p. 94101.
[44] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in
C++, Cambridge University Press, 2002.
[45] G. Roati, C. DErrico, L. Fallani, M. Fattori, C. Fort, M. Zaccanti, G. Modugno, M. Modugno, and M. Inguscio, Anderson localization of a non-interacting Bose–Einstein condensate,
Nature, 453 (2008), pp. 895–898.
[46] L. Sanchez-Palencia, D. Clément, P. Lugan, P. Bouyer, G. V. Shlyapnikov, and A. Aspect,
Anderson localization of expanding Bose-Einstein condensates in random potentials, Physical Review Letters, 98 (2007), p. 210401.
[47] T. Schwartz, G. Bartal, S. Fishman, and M. Segev, Transport and Anderson localization in
disordered two-dimensional photonic lattices, Nature, (2007), pp. 52–55.
[48] A. Scott and J. Eilbeck, On the ch stretch overtones of benzene, Chemical Physics Letters, 132
(1986), pp. 23–28.
[49] A. C. Scott, Dynamics of Davydov solitons, Physical Review A, 26 (1982), p. 578.
[50] C. Skokos, E. Gerlach, J. Bodyfelt, G. Papamikos, and S. Eggl, High order three part split
symplectic integrators: Application to the disordered discrete nonlinear Schrödinger equation, arXiv
preprint arXiv:1302.1788 [physics.comp-ph], (2013).
[51] C. Skokos, D. Krimer, S. Komineas, and S. Flach, Delocalization of wave packets in disordered
nonlinear chains, Physical Review E, 79 (2009), p. 056211.
[52] H. Veksler, Y. Krivolapov, and S. Fishman, Spreading for the generalized nonlinear Schrödinger
equation with disorder, Physical Review E, 80 (2009), p. 037201.
[53] H. Yoshida, Construction of higher order symplectic integrators, Physics Letters A, 150 (1990),
pp. 262–268.
[54] W. Yu-Shun, L. Qing-Hong, and S. Yong-Zhong, Two new simple multi-symplectic schemes for
the nonlinear Schrödinger equation, Chinese Physics Letters, 25 (2008), pp. 1538–1540.

17

