An attempt at obtaining stationary
solutions to nonlinear Schrödinger’s
equation via a pseudo-spectral method
A REPORT
submitted in partial fulfilment of the requirements
for the award of the dual degree of

Bachelor of Science-Master of Science
in

PHYSICS
by

CHIRAG ADWANI
(18077)
Under the guidance of

DR YANA NEC
Thompson Rivers University,
Kamloops, BC, Canada

DR AUDITYA SHARMA
Indian Institute of Science
Education and Research, Bhopal

DEPARTMENT OF PHYSICS
INDIAN INSTITUTE OF SCIENCE EDUCATION AND RESEARCH
BHOPAL
BHOPAL, MADHYA PRADESH, INDIA - 462066
April 2023

CERTIFICATE
This is to certify that Chirag Adwani, BS-MS (Physics), has worked on the
project entitled ‘An attempt at obtaining stationary solutions to nonlinear
Schrödinger’s equation via a pseudo-spectral method’ under my supervision
and guidance. The content of this report is original and has not been submitted
elsewhere for the award of any academic or professional degree.

April 2023
IISER Bhopal

Dr Yana Nec
Thompson Rivers
University, Kamloops,
BC, Canada

Committee Member

Signature

ii

Dr Auditya Sharma
Indian Institute of
Science Education and
Research, Bhopal

Date

ACADEMIC INTEGRITY AND
COPYRIGHT DISCLAIMER

I hereby declare that this project is my own work and, to the best of my knowledge, it contains no materials previously published or written by another person,
or substantial proportions of material which have been accepted for the award of
any other degree or diploma at IISER Bhopal or any other educational institution,
except where due acknowledgement is made in the document.
I certify that all copyrighted material incorporated into this document is in
compliance with the Indian Copyright Act (1957) and that I have received written
permission from the copyright owners for my use of their work, which is beyond the
scope of the law. I agree to indemnify and save harmless IISER Bhopal from any
and all claims that may be asserted or that may arise from any copyright violation.

April 2023
IISER Bhopal

Chirag Adwani

iii

ACKNOWLEDGEMENT

I wish to express immense, heartfelt gratitude towards my supervisor Dr Yana
Nec (TRU, Kamloops, BC, Canada) for her able guidance, invaluable insights and
constant support throughout the project. I am sincerely indebted to her for unending
belief in me and for investing her time in helping me become a better researcher and
an even better person. I would also like to thank my co-guide Dr Auditya (IISERB)
for all his support and insights which have aided me in the completion of this project.
I would like to extend a special thanks to Prof Jeremy Upsal from the University of
Washington for his inputs in the project.
I want to thank my dear friends (in no particular order of importance) Padu,
Kavin, Paras, Chandu, Bhavz, Mansi, Mishika, Saksham, Shivan, Hemakshi, Aarya
and Samidha for their constant love and support, even during very tough times. I
want to extend a special thanks to all my colleagues for providing me with support
whenever I have needed it and for providing me with a huge room to grow both
professionally and socially. I would like to especially thank my loving sister Suman
and my parents for their constant support in my pursuit of a career in Sciences. I
would also like to thank Dr Sheereen K Bajpai (IISERB) and the student counselling
cells at IISER Bhopal and TRU for all their help. I am also grateful to the existence
of caffeine, without which, writing this thesis would have been impossible.
Finally, I would like to thank MITACS for providing me with an opportunity to
work under the guidance of Dr Yana on the same topic, in person at the Thomson
Rivers University. My four month stay in Canada was financially aided by a MITACS
scholarship of CA$8,896 for the first three months, and by a generous scholarship of
CA$1,667 from Dr Yana Nec and TRU, for the last month, for which I am immensely
grateful.

Chirag Adwani

iv

ABSTRACT

The concept of derivatives with non-integral orders has been studied for centuries and has found numerous applications in physics and engineering. Recently,
there has been growing interest in exploring fractional versions of popular differential
equations. This study aims at finding stationary solutions to the fractional generalisation of the Nonlinear Schrödinger equation. Analytical stationary solutions to the
nonlinear Schrödinger equation (NLSE) exist and are well known in literature [1] to
have their amplitudes related to the Jacobi elliptic trigonometric functions. However, the spectrum of linearisation around the Jacobi elliptic solution is known to
be purely imaginary, making it impossible to obtain the solution with direct numerical attempts. We present a novel approach attempting to stabilise the numerical
scheme and obtaining the Jacobi elliptic solutions by replacing the NLSE with a
pair of coupled (auxiliary) PDEs inspired by the Gierer-Meinhardt (GM) model,
which is known to have stable spike-type solutions that have negative real parts in
their stability spectrum[2]. The stable nature of the solution persists even when the
Laplacian is replaced with a derivative of a non-integral order[3], which makes it a
suitable inspiration for our purposes. Results of the numerical trial of the auxiliary
PDE approach are presented. We use a pseudo-spectral scheme based on Fourier
transforms to simplify the problem from a system of partial differential equations to
ordinary differential equations and to allow for generalisation to non-integral order
derivatives. Despite extensive experimentation with various parameters, our method
fails to converge to the Jacobi elliptic solutions. Lastly, we also present a comprehensive review of the literature on fractional calculus, the Nonlinear Schrödinger
equation and the stability analysis of stationary solutions to the NLSE.

v

LIST OF SYMBOLS OR
ABBREVIATIONS

sn (x, k)

Jacobi elliptic sine function in x with elliptic modulus k

sgn(x)

The sign function. = +1 if x > 0, = −1 if x < 0 and = 0 if x = 0

Γ(x)

The gamma function

NLSE

The Nonlinear Schrödinger Equation

fNLSE

The Nonlinear Schrödinger Equation of a fractional order

GL

The Ginzburg-Landau Equation

CCGL

Complex Cubic Ginzburg-Landau Equation

AB2

Two-step explicit Adams-Basforth scheme

CN

Crank-Nicolson method

FT

Fourier Transform

DFT

Discrete Fourier Transform

FFT

The Fast Fourier Transform

vi

LIST OF FIGURES

2.1
2.2
4.1
4.2

The (analytical) trivial phase Jacobi elliptic sine solution to the NLSE
for k = 0.8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
Obtained outputs of an attempt to directly solve the NLSE. . . . . . 30
Obtained outputs of a numerical implementation of the auxiliary PDE
method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
Obtained solution for τ0 = 0.035 at different values of t . . . . . . . . 41

vii

Contents
Certificate

ii

Academic Integrity and Copyright Disclaimer

iii

Acknowledgement

iv

Abstract

v

List of Symbols or Abbreviations

vi

List of Figures

vii

1 Introduction and motivation
1.1 Dispersion relations . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 The Nonlinear Schrödinger equation . . . . . . . . . . . . . . . . . . .
1.3 Fractional calculus . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.4 The fNLSE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10
11
12
15
21

2 Theoretical background
2.1 Stationary solutions to the NLSE . . . . . . . . . . . . . . . . . . . .
2.2 Stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Numerical Schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24
24
27
29

3 Approach and Methodology
3.1 The Gierer-Meihardt model . . . . . . . . . . . . . . . . . . . . . . .
3.2 The auxiliary PDE method . . . . . . . . . . . . . . . . . . . . . . .
3.3 Implementation of Numerics . . . . . . . . . . . . . . . . . . . . . . .

33
33
35
37

4 Results

39

5 Conclusion and Future prospects

42

viii

Contents
Appendices
A Fourier Transforms
A.1 Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.2 Properties of Fourier Transforms . . . . . . . . . . . . . . . . . . . . .
A.3 Discrete Fourier Transforms . . . . . . . . . . . . . . . . . . . . . . .

45
45
46
47

B The Jacobi elliptic trigonometric functions

49

Bibliography

53

ix

Chapter 1
Introduction and motivation
There exists an idea of extending the definitions of (“usual”) derivatives of integer
order n to any arbitrary real (or even complex) order α. This extension can at
first be demanded with the same motivation that extensions to the natural number
system were demanded – purely to appease the curiosity. For instance, if we were
only aware of the integer number system, one may ask “Is there a number that can
be added to itself n times, to give 1?”, allowing us to discover the multiplicative
inverses and then eventually, the rationals. From there, one may ask “Is there a
number that can be multiplied to itself twice, to give the answer 2?”, allowing us to
expand our knowledge even beyond the rationals.
The extensions arising from mathematical curiosities have allowed us to answer
and explain a number of phenomena. The length of the hypotenuse of a right triangle
with unit base and perpendicular could only be represented with a symbol after
people accepted that numbers could exist beyond rationals. Roots to quadratic
equations like x2 + 1 = 0 could only be given when numbers beyond the reals
were accepted to exist. The argument, of course, does not apply only to number
systems. The Gamma function Γ(x) was discovered as an extension of the integer
operator known as the factorial and has also seen many applications throughout
mathematics. A number of applications of fractional derivatives also exist, as we
shall see in a later section. A common prescription to obtaining applications of
fractional derivatives involves extending the equations governing well-established
systems to their fractional counterparts. This has been a successful method in
helping us explain naturally occurring phenomena like anomalous diffusion [3, 4] and
Lévy flights. With this motivation, we shall explore the extension of the nonlinear
Schrödinger equation involving fractional derivatives.
The Nonlinear Schrödinger equation(NLSE) is a canonical equation of nonlinear
dynamics, which describes the time evolution of a wave travelling in a dispersive
10

1.1. Dispersion relations

Chapter 1. Introduction and motivation

medium. The dynamics of such systems is characterised by dispersion relations.
These relations will prove useful while motivating important questions this project
aims to answer.

1.1

Dispersion relations

At the heart Quantum Mechanics, just like any other branch of Physics, is a differential equation. The equation, termed as the Schrödinger’s equation, is a model
for one of many ways of describing quantum mechanical systems. Particularly, it
describes the time-evolution of the quantum state of a particle. It was first postulated by Erwin Schrödinger, after whom the equation has been named, in 1925.
In its usual form, the Schrödinger’s equation (1.1.1) is a linear differential equation
with its structure similar to that of the diffusion equation: consisting of a single
derivative in the time variable and double derivative in spatial variables,
ιℏ

ℏ2 ∂ 2 ψ(x, t)
∂ψ(x, t)
=−
+ V (x, t)ψ(x, t) ,
∂t
2m ∂x2

(1.1.1)

where ψ (x, t) is the wavefunction of the particle and V (x, t) is the potential that
represents the environment that the particle is in. The Schrödinger’s equation is
intended to be a representative of quantum mechanics that is analogous to Newtonian mechanics. Which is to say that it must follow the classical energy-momentum
relation
E=

p2
+V ,
2m

(1.1.2)

where E is the total energy of the system and p is the momentum. It is clear that
the system is totally nonrelativistic. Relations such as (1.1.2) are often termed as
dispersion relations. The term dispersion relation originates from optics, where in
that context, it is the relation between the wave number k and the frequency ω.
The dispersion relation for usual optical systems, i.e., the systems that follow the
free wave equation
∂ 2 u(x, t)
1 ∂ 2 u(x, t)
=
,
c2 ∂t2
∂x2

11

(1.1.3)

1.2. The Nonlinear Schrödinger equation

Chapter 1. Introduction and motivation

can be obtained by taking the Fourier transform on both sides in all space and time
variables (refer chapter A) to get
1
F
c2

(

∂ 2 u(x, t)
∂t2

)

(
=F

x→k,t→ω

∂ 2 u(x, t)
∂x2

)
,
x→k,t→ω

1
=⇒ 2 ω 2 û = k 2 û.
c
Where û ≡ F{u}. We have used the property (A.2) of Fourier transforms to get
from the first expression above, to the second. This gives us the relation ω = ck
that we know from optics.
Similarly, if we take Fourier transforms on both sides of the Schrödinger’s equation (1.1.1) with V = 0, we get
(
ιℏF

∂ψ(x, t)
∂t

)
x→k,t→ω

ℏ2
=−
F
2m

(

∂ 2 ψ(x, t)
∂x2

)
x→k,t→ω

2 2

=⇒ ℏω ψ̂ =

ℏk
ψ̂ .
2m

Where again, we have ψ̂ ≡ F{ψ}. If we now remember the De-Broglie relations
E = ℏω and p = ℏk, we see that we have arrived at the dispersion relation (1.1.2).
In general, one may obtain the dispersion relations of a system if the basic
equations governing the system are known. However, It is quite often the case that
the dispersion relations are easier to experimentally obtain or to be realised through
physical arguments, than the basic equations themselves. And once we have a
dispersion relation governing our system, we may derive the basic equations that the
system would follow, instead of going the other way around. The whole prescription
of looking at experimental data, finding the dispersion relations and then coming
up with a basic equation of the system provides us with a phenomenological method
of doing physics. As stated earlier, we will motivate the Nonlinear Schrödinger
equation(NLSE) and even the fractional Nonlinear Schrödinger equation(fNLSE),
using dispersion relations.

1.2

The Nonlinear Schrödinger equation

The Nonlinear Schrödinger equation has been an equation of great physical and
mathematical interest for the past 60 years. The essence of NLSE can be found
in the works of Ginzburg and Landau(1950) in their study of macroscopic superconductivity. However, it was not until mid-1960s that the wider applicability of
12

1.2. The Nonlinear Schrödinger equation

Chapter 1. Introduction and motivation

the NLSE started coming into the focus. The NLSE has found applicability in a
very diverse set of fields[5], including nonlinear optics [6], plasma physics [7], fluid
dynamics, magnetic spin waves, etc. Generally speaking, it has been found
that most weakly nonlinear, dispersive, energy-conserving systems lead
to the NLSE, in an appropriate limit. Let us take a closer look at how the
NLSE appears in nonlinear optical systems.
Optical fibre technology has brought about a paradigm shift in data transmission
technology. It allows for transmission with a much larger bandwidth and over significantly longer distances as compared to the electric cables. However, as is with any
new technology, it came with its own set of limitations. Dispersion and nonlinearity
are amongst some serious limitations of optical fibres [8].
Dispersion originates from the frequency dependence of the refractive index. The
different spectral components of the signal propagate at different speeds, and as a
result, arrive at different times at their destination. This results in pulse broadening
and distortion.
Fibre nonlinearity is due to the so-called Kerr effect, where the refractive index
develops a dependence on the intensity of the beam passing through it. In the
presence of dispersion and Kerr nonlinearity, the refractive index of the material is
given by
n(ω, E) = n0 (ω) + n2 |E|2 ,

(1.2.1)

where ω and E represent the frequency and the electric field of the electromagnetic
wave, respectively; n0 (ω) is the linear, frequency-dependent part of the refractive
index; and the constant n2 is referred to as the Kerr coefficient. It is typically of
the order 10−22 m2 /W. Despite its small value, the nonlinearity can cause effects
to accumulate over large distances and result in significant pulse broadening during
propagation. Note that if we have a modulated wave, as is often the case with signal
propagation, E will represent the amplitude of the carrier wave.
The following is a sketch of the derivation. A detailed and comprehensive treatment can be found in [6] along with relevant references cited therein. We use the
motivation from the previous section to write down the dispersion relation for a
wave propagating through a medium with refractive index given by (1.2.1):
k(ω, E) =



ω
n0 (ω) + n2 |E|2 ,
c

(1.2.2)

where c denotes the speed of light in vacuum. The Taylor series expansion of k(ω, E)

13

1.2. The Nonlinear Schrödinger equation

Chapter 1. Introduction and motivation

around ω = ω0 up to second order in ω yields
ω0 n2 2
k ′′ (ω0 )
(ω − ω0 )2 +
|E| ,
k − k0 = k (ω0 )(ω − ω0 ) +
2
c
′

(1.2.3)

where the prime (′ ) denotes derivative in ω and k0 = k(ω0 ). And now, exactly in
the opposite way as in the previous section, we replace k − k0 and ω − ω0 with their
corresponding inverse Fourier operator equivalents ι∂/∂z and ι∂/∂t, respectively.
We also then let the whole expression operate on E, to obtain

ι

∂E
∂E
+ k ′ (ω0 )
∂z
∂t


−

k ′′ (ω0 ) ∂ 2 E
+ ν|E|2 E = 0,
2
∂t2

(1.2.4)

where ν ∈ R is an appropriate constant, defined by collecting all the constants in
front of the cubic term. We may identify k ′ (ω0 ) = 1/vg , where vg is the group
velocity of a given travelling wave.
To obtain the dimensionless version of (1.2.4), we first introduce retarded time
coordinate tret = t − z/vg . Assume that t∗ , z∗ and P∗ are the characteristic time,
distance and power of the system, respectively. If we choose z∗ = 1/νP∗ and t2∗ =
z∗ | − k ′′ (ω0 )| and move to the dimensionless coordinates t′ = t/t∗ , z ′ = z/z∗ and q =
√
E/ P∗ , we arrive at the NLSE
ι

∂q sgn(−k0′′ (ω0 )) ∂ 2 q
+
+ |q|2 q = 0.
2
∂z
2
∂t

(1.2.5)

There are now two cases of physical interest depending on the sign of −k ′′ (ω0 ).
When −k ′′ (ω0 ) > 0, we obtain what is termed as the defocusing NLSE, and the
resulting dispersion in this case is “normal”. While −k ′′ (ω0 ) < 0 results in the
focusing NLSE and “anomalous” dispersion. Unless specifically mentioned, we will
restrict our analysis to the defocusing NLSE, even after we change the form of the
Laplacian in the following sections. Hence, the equation of our interest is
1
ιqt = − qxx + q|q|2 .
2

(1.2.6)

There has been an extensive amount of research on the solutions of the NLSE and
their properties[1]. The NLSE is known to be a completely integrable system. In
the context of this project, we would be interested in separable solutions of the type
Ψ(x, t) = e−ιωt ϕ(x) ,
where ω is a real constant. We are specifically interested in a class of periodic

14

1.3. Fractional calculus

Chapter 1. Introduction and motivation

solutions in x, where the amplitude ϕ(x) is expressed in terms of the Jacobi elliptic
sine function sn(x, k) (refer Appendix B):
ι

2

Ψ(x, t) = k sn(x, k)e− 2 (1+k )t ,

(1.2.7)

where k ∈ [0, 1) is called the elliptic modulus. sn(x, k) is a periodic function with
the period T (k) = 4K(k), and K(k) is the elliptic integral of first kind. We will
take a deeper look at the NLSE and the Jacobi elliptic solutions in the next chapter,
where we also discuss the stability spectrum of these solutions.
Finally, we would like to point out the resemblance between the Nonlinear
Schrödinger equation (NLSE) and the complex cubic Ginzburg-Landau (CCGL) [5]
equation, which appears in the context of pattern formation processes in nonlinear
media:
ψt = gψ + (a + ιb)∇2 ψ − (d + ιc)|ψ|2 ψ,

(1.2.8)

where ψ is complex. g > 0 is known as the gain coefficient. a > 0 and b ∈ R
are coefficients controlling the diffusion and dispersion processes. c ∈ R and d > 0
account for the nonlinear dissipations. In the |a| ≪ |b| and |d| ≪ |c| limit, we
can see that the CGCL equation becomes the NLSE. So, the NLSE also makes its
appearance in the theory of pattern formation.
In summary, the above discussion has highlighted the immense physical and
mathematical significance of the Nonlinear Schrödinger equation (1.2.6). In a later
section, we will introduce a generalisation of the NLSE based on the concept of
fractional derivatives, which extends the conventional integer-order derivatives to
derivatives of arbitrary real order α. This idea falls under the realm of fractional
calculus, which is the subject matter of the next section.

1.3

Fractional calculus

The notion of fractional derivatives has existed for centuries, originally appearing in
the works of Leibniz in 1695. However, it was only after more than a century later,
in 1823, that a practical application of fractional calculus was found. Niels Abel
published his treatment of the famous tautochrone problem in 1823, at an age of
21, where he used fractional derivatives to find the solution [9]. That stood as the
sole example of an application of fractional derivatives, for almost an entire century,
that is up until the early 1920s, when Oliver Heaviside, an autodidact, found the
application of fractional derivatives in his works on the theory of transmission lines
15

1.3. Fractional calculus

Chapter 1. Introduction and motivation

[10].
Even though there was a steady progress in the subject on the side of mathematical exploration, it was not until the works of Mandelbrot on fractal geometry in the
early 1980s [11, 12] that the topic started to gain the interest of physicists too. Fractional derivatives started appearing in the contexts of fractional Brownian motion
[13] and anomalous diffusion processes [14]. Since then, fractional versions of various,
physically important differential equations have come to light and raised fundamental questions about the physical nature of the systems governed by those equations.
For instance, N. Laskin (2002) [15] presented a fractional version of the Schrödinger
equation and also proved, in the same article, that the fractional Schrödinger operator was Hermitian. Goldfain (2006) [16] demonstrated that fractional dynamics
in Minkowski spacetime mimics that of free fields in curved Riemannian spacetime.
An approach to quantise fractional free fields was also given in the same article. A
quantisation description for fractional free Klein-Gordon type fields has also been
addressed in literature [17]. Fractional calculus thus presents us with a scope to
gain novel physical insights and develop innovative methodologies and approaches.
There is a multitude of definitions of fractional derivatives in the literature [9, 10],
which are all used depending on the context they are required in. Most definitions
require that we identify derivatives and anti-derivatives as inverse operators of each
other. That is, if we define the anti-derivative (or integration) operator as
Z x
f (x0 )dx0 ,

a Ix (f ) :=
a

then we must identify this as the inverse of derivative operator a Dx ,
−1
a Dx ≡ a Ix .

In the above notation the superscript on the right denotes the order of the operator,
the subscript to the left denotes the lower limit of the integral in the definition of the
operator and the subscript on the right denotes the independent variable that the
operator is defined with. Now, if we look at higher order anti-derivative operators,
n
a Ix for n ∈ N, we see that it assumes the form of a repeated integral
n
a Ix (f ) =

Z x

Z xn−1
dxn−2 · · ·

dxn−1
a

Z x2

a

Z x1
dx1

a

f (x0 ) dx0 ,
a

where we assume that the function f is of a sufficiently regular behaviour for all
required integer integrals and derivatives to exist. The above repeated integral can
be written in terms of a single integral using the Cauchy’s formula for repeated
16

1.3. Fractional calculus

Chapter 1. Introduction and motivation

integrals.
Theorem 1.1 (Cauchy’s formula for repeated integrals). For a continuous, real
function f , the nth repeated integral
n
a Ix (f ) =

dxn−2 · · ·

dxn−1

f (x0 ) dx0 ,

dx1
a

a

a

Z x1

Z x2

Z xn−1

Z x

a

equals the single integral
n
a Ix (f ) =

1
(n − 1)!

Z x

(x − t)n−1 f (t) dt.

a

Proof. The proof follows from mathematical induction.
For n = 1, the formula reads
1
a Ix (f ) =

1
(0)!

Z x

Z x

0

(x − t) f (t) dt =
a

Z x
f (t)dt =

a

f (x0 ) dx0 ,
a

which satisfies the base case without much effort, if we incorporate the extension of
the factorial to have 0! = 1. Next, assume that the hypothesis holds for k ∈ N. We
also note that

 Z
Z x
1
d 1 x
k
(x − t) f (t) dt =
(x − t)k−1 f (t) dt.
(1.3.1)
dx k! a
(k − 1)! a
using the Leibniz Integral rule1 . Now,
k+1
(f ) =
a Ix

Z x

Z xk

Z x2

Z xk−1

Z x1

dx1
f (x0 )dx0
dxk−1
dxk−2 · · ·
a
a
a
Z x
Z xk
1
=
dxk
(xk − t)k−1 f (t) dt (from induction hypothesis)
(k
−
1)!
 Z a xk

Za x
d
1
k
=
dxk
(xk − t) f (t)dt
(from equation (1.3.1))
dxk k! a
a
Z
1 x
=
(x − t)k f (t) dt,
k! a
(1.3.2)
dxk

a

a

which completes the proof.
1

The Leibniz Integral rule for differentiating an input function with an integral sign:
!
Z b(x)
Z b(x)
d
∂
g(x, t) dt = g(x, b(x))b′ (x) − g(x, a(x))a′ (x) +
g(x, t) dt.
dx
∂x
a(x)
a(x)

17

1.3. Fractional calculus

Chapter 1. Introduction and motivation

We can then write an expression for a higher order anti-derivative as
−n
n
a Dx (f ) = a Ix (f ) =

1
(n − 1)!

Z x

(x − t)n−1 f (t) dt,

(1.3.3)

a

where n ∈ N. The above expression provides us with an approach of defining integrals
of more general orders. If we allow the formula to hold for any α ∈ R and replace the
factorial with its analytic continuation, Γ(α), we can define the Riemann-Liouville
fractional Integral:
Definition 1.2 (The Riemann-Liouville fractional integral). For α ∈ R, the RiemannLiouville fractional integral can be defined as
−α
α
a Dx (f ) = a Ix (f ) =

1
Γ(α)

Z x

(x − t)α−1 f (t) dt,

a

where the constant a is the lower bound of the integration domain and needs to be
specified a priori.
Now, we again emphasise the fact that the integration operators and the derivative operators conceptually are the inverses of each other2 , allowing us to write
Dα = Dm D−(m−α) ,
dm
= m I m−α
dx

m∈N
(1.3.4)

where we dropped the subscripts for convenience. We shall make the choice m = ⌈α⌉,
which denotes the smallest integer ≥ α. The above may be interpreted as a definition
of a fractional derivative of order α, wherein we are required to first compute a
Riemann-Liouville fractional integral of the order (m − α), and then take the usual
derivative of order m. It must be noted that when α is an integer, ⌈α⌉ will be the
same as α and we will just be left with computing the usual derivative (i.e. we
want to identify the integral operator of order 0, a Ix0 , with the identity operator so
that Dm ≡ dm /dxm , for all m ∈ N). Thus, we arrive at the Riemann-Liouville
definition of a fractional derivative:
Definition 1.3 (Riemann-Liouville fractional derivative). For α ∈ R and m = ⌈α⌉,
the Riemann-Liouville fractional derivative can be defined as
Z x
dm
dm
1
m−α
α
(f ) = m
(x − t)m−α−1 f (t) dt,
a Ix
a Dx (f ) =
m
dx
dx Γ(m − α) a
2

Note that this is not exactly the case! A beautiful discussion on this is given in [10]. However,
it is beyond the scope of this thesis.

18

1.3. Fractional calculus

Chapter 1. Introduction and motivation

where again the constant a that defines the domain, will have to be specified a priori.
One may also have given a definition of fractional derivatives by inverting the
order in which the derivative and integral operators appear in equation (1.3.4). This
gives us the Caputo fractional derivative.
Definition 1.4 (Caputo fractional derivative). For α ∈ R+ and m = ⌈α⌉, the
Caputo fractional derivative can be defined as
m
C α
m−α d
D
(f
)
=
I
(f ) =
a
a
x
x
m

dx

1
Γ(m − α)

Z x

(x − t)m−α−1 f (m) (t) dt,

a

where once again, the constant a will have to be specified a priori. We introduce the
α
left superscript in C
a Dx to distinguish the Caputo derivative (1.4) from the RiemannLiouville derivative (1.3). One of the eminent advantages of the Caputo derivative
is that the derivative of a constant function is 0, which is not the case for the RL
derivative and which may pose problems in contexts where this feature important,
such as equilibria of dynamical systems.
It can be noted that the above definitions of fractional derivatives suggest that
the fractional operators are nonlocal. Which means that the fractional derivative
of a function f at a point x in its domain, is dependent upon the shape of the
function elsewhere in the domain too (this is also the reason why we have been
carrying the left subscript in our notation). This is quite visibly the case for the
Riemann-Liouville and Caputo fractional derivatives as they are both dependent on
the computation of an integral involving the function itself, over the whole domain.
So, if the function (or its derivative in the case of Caputo) were to be modified
continuously elsewhere, it would affect the value of derivative of non-integral order
α, at the point x. However, for integer values of α, the derivative completely loses
its property of nonlocality and such a change will not make a difference in the value
of derivative at the point x. There are many such interesting properties of fractional
derivatives beyond the scope of this discussion.
Both Riemann-Liouville and Caputo derivatives provide us with plausible definitions of the fractional derivatives. However, they are not the most practical when it
comes to numerical implementations. One will not only have to worry about using
a good quadrature rule to carry out the integration process, but also a good scheme
to compute the derivative. We will also have to worry about the cases for which α
is an integer, separately.
It may be noted that the above definitions were developed based on generalisations of the properties of the usual, integer order derivatives. Specifically, both arise
19

1.3. Fractional calculus

Chapter 1. Introduction and motivation

from the composition property of derivatives exploited in (1.3.4). One may look for
other such properties of usual derivatives to generalise and thence define fractional
derivatives. To this end, we recall the Fourier transform property of derivatives
Property 1.1 (Fourier transform property of derivatives). For a derivative operator
dm /dxm of integer order m in the variable x, the following relation holds true for a
Fourier-transformable function
 m


d
m
F
f
(x)
=
(ιξ)
F
f (x) x−→ξ ,
(1.3.5)
dxm
x−→ξ
where F denotes the Fourier transform.
It is desirable to define a similar property for non-integral orders. This is where
the Riesz fractional derivative comes to the rescue [18].
Definition 1.5 (Riesz fractional derivative). For α ∈ R+ , we define the Riesz
fractional derivative as the Riemann-Liouville derivative with a = −∞,
α
−∞ Dx (f ) =

dm
1
dm
m−α
I
(f
)
=
−∞
x
m
m
dx
dx Γ(m − α)

Z x

(x − t)m−α−1 f (t) dt.

(1.3.6)

−∞

It can be proven that the Riesz fractional derivative as defined above, possesses
the following Fourier transform property:
F



α
−∞ Dx (f )


:= (ιξ)α F f (x) ,

However, using this as the Fourier transform for derivatives of any order α ̸= 2, α > 0
will bring in another issue as the derivative operator itself does not possess the
reflection symmetry x 7→ −x. This symmetry is essential to keep the expression
of the derivative real, under a Fourier transform. The Riemann-Liouville and the
Caputo derivative also do not posses this symmetry either, whilst all highest order
derivatives of popular equations (like ∇2 , ∇4 , etc.) always do seem to possess it.
To obtain a definition of a fractional derivative that has this symmetry, we will first
make a variable shift in 1.3.6
Z ∞
dm
1
dm
m−α
α
(f ) = m
(t − x)m−α−1 f (t) dt,
(1.3.7)
xI
∞ Dx (f ) =
dxm ∞
dx Γ(m − α) x
the Fourier transform of which, would look like
F



α
∞ Dx (f )


= (−ιξ)α F f (x) .

20

1.4. The fNLSE

Chapter 1. Introduction and motivation

This allows us to define fractional Laplacian, which possesses reflection symmetry
(note that we omit α = 2m + 1, m ∈ Z as that would correspond to the case which
should typically not possess the reflection symmetry that we are trying to achieve.):
Definition 1.6 (The fractional Laplacian). For α ∈ R+ , α ̸= 2m + 1, m ∈ Z, we
define the fractional Laplacian ( −(−∇2 )α/2 ) as
2 α/2

−(−∇ )


(f ) :=

−

α
α
−∞ Dx + ∞ Dx

2 cos(απ/2)


(f )

(1.3.8)

Now, if we keep in mind the fact that
|ξ|α =

(ιξ)α + (−ιξ)α
,
2 cos(απ/2)

the fractional Laplacian will posses the following Fourier transform property
F



2 α/2

− (−∇ )


(ιξ)α + (−ιξ)α 
(f ) = −
F f (x) = −|ξ|α F f (x) .
2 cos(απ/2)

Hence, we now have a reflection symmetric fractional derivative. This also gives
us a fairly straightforward way to implement fractional derivatives numerically. In
contrast to Riemann-Liouville, Caputo and Riesz derivatives, we would only need to
perform a Fourier transform, and then invert it after multiplying by some factors, to
implement the fractional Laplacian numerically. This process can be achieved with
a far greater efficiency, owing to the existence of Fast Fourier Transforms (FFTs).
The ability to implement fractional derivatives numerically, allows us to investigate a wide range of previously unexplored problems. The specific problem that we
shall be interested in, is the fractional generalisation of the Nonlinear Schrödinger
equation, as discussed at the end of the previous section. We will take a closer look
at the motivations to investigate this problem, in the next section.

1.4

The fractional Nonlinear Schrödinger equation (fNLSE)

In the previous section, we saw that fractional derivatives have historically been
used to provide a generalisation for equations governing various branches of physics
and engineering, allowing us to explain a wide variety of novel physical phenomena
[5]. The fractional generalisation has been particularly successful in explaining large
classes of diffusion-related phenomena. For instance, in the case of free diffusion

21

1.4. The fNLSE

Chapter 1. Introduction and motivation

processes, whenever one finds a deviation from Fick’s law of diffusion3 , fractional
calculus might be employed [4]. The aesthetic similarity between the diffusion equation and the Schrödinger equation has been a motivation for enquiries on fractional
quantum mechanics and fractional linear Schrödinger equations of various types
[15, 19, 20]. We have seen in section 1.2 that similarities exist between the Complex
Cubic Ginzburg-Landau (CCGL) equation and the Nonlinear Schrödinger equation,
and because there are fractional generalisations of the CCGL [4, 21], there is motivation to look at the fractional generalisation of the Nonlinear Schrödinger equation
as well.
The Nonlinear Schrödinger equation was introduced in section 1.2, wherein a
sketch of a derivation of the NLSE was also given. We can follow a similar process to
derive the fractional Nonlinear Schrödinger equation as well, by essentially replacing
the dispersion relation (1.2.3) by the following (as done in [21]):
k = ωk ′ + |ω|α

k ′′
+ ν|E|2 ,
2

(1.4.1)

where α ∈ R. The next step in the derivation was to replace the powers of ω with
their inverse Fourier operator equivalents. We can recall that the equivalent of |ω|α
would be the fractional Laplacian (definition 1.6), multiplied by a factor of −1. If
we then keep following the rest of the derivation exactly as in the above section, we
will arrive at the fractional Nonlinear Schrödinger equation:
ι

1
∂Ψ (x, t)
2
= − Dαx Ψ (x, t) + Ψ (x, t) Ψ(x, t) ,
∂t
2

(1.4.2)

where Ψ(x, t) is a complex valued function and α ∈ R+ , α ̸= 2m + 1, m ∈ Z is the
order of the fractional Laplacian (Dαx = −(−∇2 )α/2 ). The α = 2 case corresponds
to the usual nonlinear Schrödinger equation (1.2.6). Thus, we are now at a place
where we can establish the objectives of this project.
Objective 1. An objective of this project is to make an attempt at numerically
obtaining stationary solutions of the nonlinear Schrödinger’s equation of a fractional
order (fNLSE) (1.4.2).
To make an approach towards achieving objective 1, we may first attempt to
3

The usual Fickian diffusion processes will have the Mean Square Displacement (MSD) linearly
proportional to t:
⟨r2 (t)⟩ ∝ t.
However, various types of ‘anomalous’ diffusion processes have emerged over the years that have
⟨r2 (t)⟩ ∝ tγ , where γ can be ≶ 1. Fractional calculus has been highly successful in explaining such
processes[4].

22

1.4. The fNLSE

Chapter 1. Introduction and motivation

obtain the Jacobi elliptic solutions (1.2.7) of the usual NLSE (α = 2), that were
introduced in section (1.2). Once we have a working scheme to obtain the Jacobi
elliptic solutions, we may proceed to generalise the algorithm to work for fractional
order derivatives. However, to achieve our intended generalisation, we will try to
obtain the Jacobi elliptics by using the Fourier transform property (1.3.5) of derivatives in the numerical implementation, as it is the closest to the definition of Riesz
fractional derivatives. Doing so will also aid us in another way; as after we transform
the Laplacian, we would only be left to solve an ordinary Initial Value Problem in
t, as compared to a partial differential equation in two variables that we started off
with. Thus, we have
Objective 2. To attain objective 1, we will first aim to obtain the Jacobi elliptic solutions (1.2.7) numerically using a scheme that utilises the Fourier transform
property 1.3.5 of derivatives.
In this chapter, we presented an introduction to the NLSE, its stationary solutions, and the various places that it arises in literature. An elaborate discussion
on fractional derivatives was also given, wherein we saw various definitions of the
fractional derivatives, before eventually settling for the one that is most useful to us:
the Riesz fractional derivative. We remarked the recent upsurge in the interest in
fractional derivatives, due to various discoveries of its usefulness being made across
many branches of Physics and Engineering. Finally, we established the objectives
1 and 2, of the project. In the subsequent chapters, we will make progress towards
achieving those objectives by first studying the Jacobi elliptic stationary solutions of
NLSE, in greater depth and then developing schemes to obtain them numerically.

23

Chapter 2
Theoretical background
In chapter 1, the (defocusing) NLSE was introduced as an equation governing the
amplitude of a wave in dispersive media. It was established that the equation has a
class of stationary solutions related to the Jacobi elliptic functions. In this chapter,
will look at the Jacobi elliptic solutions in a greater detail: how they arise and
whether they can be obtained numerically or not. While answering the later, we will
see what happens to perturbations around stationary solutions of the NLSE. This
chapter will also present the first numerical trial done towards a primary objective
of the project (2), which is to obtain the Jacobi elliptic numerically.

2.1

Stationary solutions to the NLSE

The (defocusing) NLSE is given by
1
ιΨt = − Ψxx + Ψ|Ψ|2 ,
2

(2.1.1)

where Ψ(x, t) is intensity of the carrier wave propagating in a dispersive medium.
We aim to look for separable solutions of the type
Ψ(x, t) = e−ιωt ϕ(x),

(2.1.2)

where ω ∈ R is a constant. Plugging this into (2.1.1), we obtain an equation for
ϕ(x)
1
ωϕ(x) = − ϕxx + ϕ(x)|ϕ(x)|2 .
2

24

(2.1.3)

2.1. Stationary solutions to the NLSE

Chapter 2. Theoretical background

Decomposing the solution further, we write ϕ(x) as an amplitude and a phase
ϕ(x) = R(x)eιθ(x) ,
where both R(x) and θ(x) are now real-valued functions. Plugging this into (2.1.3)
and separating the real and imaginary parts, we get the differential equations for
R(x) and θ(x)
1
(2.1.4a)
Rx θx + Rθxx = 0,
2
Rxx = 2R3 − 2ωR + Rθx2 .

(2.1.4b)

Equation (2.1.4a) can be solved right away by noting that it implies
θxx
−2Rx
=
θx
R
=⇒ ln(θx ) = −2ln(R) + c1 ,
Z x
c
dy,
=⇒ θ =
2
0 R (y)

(2.1.5)

where c1 and c are some integration constants. Plugging this in (2.1.4b)
Rxx = 2R3 − 2ωR +

c2
.
R3

(2.1.6)

Equation (2.1.6) can be solved using standard techniques of elliptic integration [1],
and it can be shown that R(x) is given by
R2 (x) = k 2 sn2 (x, k) + b,

(2.1.7)

where sn(x, k) is the Jacobi elliptic sine function, k ∈ [0, 1) is the elliptic modulus
and b is another real constant. A discussion on Jacobi elliptic functions can be
found in Appendix B. The Jacobi sine function sn(x, k) is a periodic function with
the period T (k) = 4K(k), where K(k) is the elliptic integral of the first kind
Z π/2
K(k) =
0

dy
p
.
1 − k 2 sin2 y

For (2.1.7), we also have
1
3
ω = (1 + k 2 ) + b,
2
2

c2 = b(b + 1)(b + k 2 ).

25

(2.1.8)

2.1. Stationary solutions to the NLSE

Chapter 2. Theoretical background

Noting that R(x), θ(x) and all of the associated constants must be real, we get b ≥ 0
(to keep c in R). If we choose to look at the solution with b = c = 0, ω = 12 (1 + k 2 ),
we will get a class of solutions with θ(x) = 0, i.e., the trivial phase solutions:
ι

2

Ψ(x, t) = k sn(x, k)e− 2 (1+k )t ,

(2.1.9)

up to symmetry transformations permitted by the original equation (2.1.1). Fig-

Figure 2.1: The Jacobi elliptic sine solution to the NLSE for k = 0.8.
ure 2.1 provides an example of a trivial phase solution (2.1.9) for k = 0.8. Only the
amplitude |ψ(x)| has been plotted in the figure as a function of x, alongside a scaled
sine function. We can see that the Jacobi sine function bears a striking resemblance
to the regular sine function, except for a flatter top.
To this end, we recall that an objective of the project is to solve the NLSE (2.1.1)
numerically and obtain the Jacobi elliptic solution (2.1.9) as a stationary solution of
the system (i.e. we see where the solution is tending towards, in the large t limit).
Furthermore, we place emphasis on employing a scheme that involves performing a
Fourier transform on the Laplacian. This is done with the intention of facilitating
the generalisation to fractional derivatives that we eventually plan to pursue. Doing
this also helps us in transforming a partial differential equation in two variables
(say x and t, just like in our case), into an ordinary differential equation in just one
variable (in our case, t), thus greatly reducing the complexity of the problem. We
can then use standard numerical integration techniques for ODEs to progress in t
and then eventually invert the Fourier transform to reach to the solution. However,
it turns out that there is a fundamental issue in trying to directly use numerical
integration methods on (2.1.1), owing to the stability of its solutions.
Numerical implementations inherently introduce perturbations to the system being analysed, that are outside of our control. These perturbations can be attributed
26

2.2. Stability analysis

Chapter 2. Theoretical background

to several factors [22], including the discrete nature of the implementation, limitations in the floating point representation system, errors caused by truncation and
rounding off, and inaccuracies that arise when reading or writing data to the disk.
Hence, the qualitative nature of stationary solutions becomes important when we
wish to obtain them numerically. The stability of a solution is determined, as we will
see in the next section, by its linearisation spectrum. It is essential that the solution
we are trying to obtain be stable, meaning that small perturbations near them decay
to zero. In order to determine whether a solution is stable or not, certain conditions
must be met, which we will explore further. We will also discuss the aforementioned
issue associated with obtaining stationary solutions of the NLSE.

2.2

Stability analysis

The notion of stability of stationary solutions falls under the realm of qualitative
stability analysis of differential equations. The idea involves linearising the system
around a stationary solution, and then analysing how the perturbations evolve over
time, while being restricted to only the first order terms. To see how this process
works out for our problem, we first transform (2.1.1) to incorporate solutions of the
type
Ψ(x, t) = e−ιωt ψ(x, t).
Using this, we get that ψ(x, t) follows
1
ιψt = −ωψ − ψxx + ψ|ψ|2 .
2

(2.2.1)

Stationary solutions will, by definition, have ψt ≡ 0. Remember that in the previous
section, we had shown that such solutions exist for the NLSE. Now, we introduce
perturbations to these solutions, in the following manner
ψ(x, t) = eιθ(x) (R(x) + ϵu(x, t) + ιϵv(x, t)) + O(ϵ2 ),

(2.2.2)

where u(x, t) and v(x, t) are real-valued functions, and |ϵ| ≪ 1. R(x) and θ(x) are the
same as in (2.1.2). Note that because there are no restrictions on u or v as functions
of space and time, the above form represents the most general perturbation. Any
changes in the phase can be factored inside. Substituting (2.2.2) in (2.2.1) and

27

2.2. Stability analysis

Chapter 2. Theoretical background

dropping any terms of order ϵ2 or higher, we get
1
eιθ (ιϵut − ϵvt ) = −ωeιθ (R + ϵu + ιϵv) − [(Reιθ )xx + ϵ(ueιθ )xx + ιϵ(veιθ )xx ]
2
+ eιθ (R3 + 3ϵR2 u + ιϵR2 v). (2.2.3)
From here, we collect all the order zero terms in ϵ and realise that they give us back
(2.1.3), vanishing thereby. Next, we separately collect the real and imaginary parts
of O(ϵ) terms to get
∂
∂t

u
v

!
=L

!
u
=J
v

L+ S
−S L−

!

!
u
,
v

(2.2.4a)

where L+ , L− and S are linear operators given by


1 2
c2
2
L+ ( · ) = − ∂x ( · ) + 3R (x) − ω +
( · ),
2
2R4 (x)


c2
1 2
2
( · ),
L− ( · ) = − ∂x ( · ) + R (x) − ω +
2
2R4 (x)


c
(·)
S( · ) =
∂x
,
R(x)
R(x)

(2.2.4b)
(2.2.4c)
(2.2.4d)

and
J=

!
0 1
.
−1 0

(2.2.4e)

In equation (2.2.4a), we can separate the variables and consider solutions of the
form
!
!
u(x, t)
U
(x)
= eλt
,
(2.2.5)
v(x, t)
V (x)
essentially reducing the system to an eigenvalue problem
U
L
V

!

U
=λ
V

!
.

(2.2.6)

By looking at equations (2.2.5) and (2.2.6), we can conclude that the spectrum of
L (denoted σ (L)) is what determines the stability of the solutions. Negative real
parts in an eigenvalue would imply stability in that eigenspace, and vice versa. And
in case the real part happens to be zero, i.e. the eigenvalue is purely imaginary, the
28

2.3. Numerical Schemes

Chapter 2. Theoretical background

perturbations would neither grow nor die out to zero.
The stability spectrum of Jacobi elliptic solution has been extensively studied in
the 2011 article by Bottman et al [1] and references cited therein. All the numerical
and analytical enquiries suggest that the Jacobi elliptic solution has a purely imaginary spectrum. Consequently, although small perturbations around the solution
remain small forever, they never approach zero. An illustration of this fact is given
in fig. 2.2. Even for large times, the solution does not appear to converge to the
Jacobi elliptic (orange), denoted by ψJE , as small perturbations seem to persist.
The obtained numerical solution (blue), denoted as ψnum , deviates away from ψJE
and then returns closer to it periodically, as can be seen from the first four panels
of the figure. The graph of relative error in the last panel confirms this behaviour
as the error oscillates with an approximately constant amplitude over time. This
output establishes that we cannot obtain the Jacobi elliptic solution by tackling the
problem head on, just as we expected from theory.
To obtain the solution here, we performed a Fourier transform on the Laplacian
in (2.1.1) as we had emphasised earlier;
n
o
1 2
dψ̂
2
= − q ψ̂ + F ψ ψ ,
ι
dt
2

(2.2.7)

where ψ̂ := F{ψ} and q is the transform variable. We then use a combination of two
different numerical integration algorithms to progress in time. The nonlinear terms
were integrated by the explicit two-step Adams-Bashforth scheme, while the linear
terms were progressed using the implicit Crank-Nicolson method. A description of
these methods can be found in the next section (2.3). Finally, we invert the Fourier
transform to obtain the solution.
The computation was carried out in Wolfram Mathematica (v12.0.0.0) on a Mac
OS X x86 (64-bit) platform. Number of points in the spatial domain, npts, was set
to 101. While the time variable ran from tmin = 0 to tmax = 10, in steps of dt =
0.0001. The plots were obtained using Mathematica’s default Plot[] and ListPlot[]
commands.

2.3

Numerical Schemes

In the previous section, we saw how Fourier transforming the Laplacian turns the
partial differential equation (2.1.1) into a first order ordinary differential equation
(2.2.7). It is still, however, not analytically solvable, just like most differential equations known to humanity. Nonetheless, for most practical applications, numerical
29

2.3. Numerical Schemes

Chapter 2. Theoretical background

(a) t = 0.0 units

(b) t = 1.6 units

(c) t = 8.4 units

(d) t = 9.9 units

(e) Relative error vs t

Figure 2.2: Time progression of the numerical solution ψnum obtained from a direct implementation of the NLSE, plotted alongside the Jacobi elliptic solution
ψJE = 0.8 sn(x, 0.8), which is also taken as the initial conditions for the computation
of ψnum . The spatial domain ranges from x = 0 to x = 2K(0.8), where K(0.8) is
the complete elliptic integral of the first kind with elliptic modulus k = 0.8. Panels (a)-(d) shows the time progression of ψnum . Panel (e) shows the relative error
(ψnum − ψJE ) /ψJE at the mid point of the domain, x = K(0.8).
approximations to the solutions are sufficient. Hence, over the years, algorithms to
obtain those numerical solutions have been developed. All of the algorithms focus on
30

2.3. Numerical Schemes

Chapter 2. Theoretical background

approximating the shape of the solution at discrete steps, given an initial condition
and an integrand.
The algorithms for finding numerical solutions to first order ordinary differential equations can be broadly categorised into two types: linear multistep methods
and predictor-corrector methods methods. The predictor-corrector methods, like
the famous Runge-Kutta methods, calculate intermediate steps between the current
and next step in order to approximate the next step. However, these methods do
not make use of any information from any of the previous steps. Linear multistep
methods, on the other hand, use one or more previously computed steps to approximate the next step. Although both algorithms can give accurate solutions up to an
arbitrary order, linear multistep methods require us to evaluate the integrand only
once per step, whereas the predictor-corrector methods access the integrand for as
many times as the number of intermediate steps. The choice of method to be used,
depends on the specific requirements of the problem at hand.
For our specific problem, the integrand is the right-hand side of (2.2.7). The
expression involves computation of a Fourier transform. Computing the Fourier
transform at multiple intermediate steps, including the repeated evaluation of the
integrand, would entail prohibitive computational cost. Hence, we will rely exclusively on linear multistep methods. There are two methods that must be pointed
out, which would be utilised, in combination, for all numerical implementations in
this project:
• Second order Adams-Bashforth (AB2) method,
• Trapezoidal rule or the Crank-Nicolson (CN) method.
Let us first give a precise definition of multistep methods [22]:
Definition 2.1 (m-step multistep methods). Given the initial value problem (IVP)
y ′ = f (t, y),

a ≤ t ≤ b,

y(a) = α,

(2.3.1)

an m-step multistep method approximates the solution yi+1 at ti+1 , as
yi+1 = am−1 yi + am−2 yi−1 + · · · + a0 yi−m+1

(2.3.2)

+ h [bm f (ti+1 , yi+1 ) + bm−1 f (ti , yi ) + · · · + b0 f (ti−m+1 , yi−m+1 )],
for i = m − 1, m, . . . , N − 1. Here, a0 , a1 , . . . , am−1 , b0 , b1 , . . . , bm are constants
and h is the grid spacing (b − a)/N .

31

2.3. Numerical Schemes

Chapter 2. Theoretical background

If bm = 0, the method is termed explicit as the next step is explicitly given by
only the terms that come before it. And similarly, if bm ̸= 0, the method is termed
implicit as the equation for yi+1 is implicit in it. Generally, implicit methods have
better accuracy than their explicit counterparts, since they are stable.
The aforementioned AB2 and CN methods are different types of multistep methods. Their corresponding definitions are as follows. The second order AdamsBashforth (AB2) is an explicit, two-step method given by
Definition 2.2 (Two-step Adams-Bashforth method).
h
yn+1 = yn + [3 f (ti , yi ) − f (ti−1 , yi−1 )],
2
where i = 1, 2, . . . , N − 1. The local truncation error for this algorithm is O(h2 ).
The trapezoidal rule is an implicit, two-step method given by
Definition 2.3 (Trapezoidal rule).
h
yn+1 = yn + [f (ti+1 , yi+1 ) + f (ti , yi )],
2
where i = 1, 2, . . . , N − 1. The local truncation error for this algorithm is O(h3 ).
This rule is equivalent to the implicit Crank-Nicolson (CN) scheme, when applied
to our initial value problem (2.3.1). Hence, we may use the names exchangeably.
To conclude, it has been established that stationary solutions to the NLSE exist
in the form of Jacobi elliptic solutions (2.1.9) and that their stability spectrum is
purely imaginary. As a consequence of the later, we saw that it is impossible to
obtain them directly by using numerical integration techniques. So to acquire the
Jacobi elliptic solution from the NLSE numerically, we seek an alternative system
that differs from NLSE in concept, but remains sufficiently similar to enable a conversion to NLSE. In the following chapter, we will explore a potential candidate for
the same.

32

Chapter 3
Approach and Methodology
In this chapter we shall present a candidate for stabilising the Jacobi elliptic solution.
As we saw in the last chapter, we cannot hope to solve the NLSE directly and get
the Jacobi elliptics. Hence, we instead replace the NLSE with a system of partial
differential equations inspired by the well known Gierer-Meinhardt system. We shall
also present the details of the numerical implementation for solving the system on
a computer.

3.1

The Gierer-Meinhardt model and spike-type
solutions

In 1972 Geirer and Meinhardt [23] introduced a model for reaction-diffusion system
involving activators and inhibitors. The system has found applicability in biological
pattern formation [24], cell differentiation [25], and pattern formation on the shells
of sea animals [26]. The model is a system of coupled partial differential equations.
In their dimensionless form the 1D equations at the limit of slow activator and fast
inhibitor are [3, 27]
∂t a = ϵ2 ∂x2 a − a +

ap
,
hq

t > 0,

1 am
, −1 < x < 1, t > 0,
ϵ hs
∂x h(x = ±1, t), a(x, 0) = a0 (x), h(x, 0) = h0 (x).

τ0 ∂t h = D ∂x2 h − µh +
∂x a(x = ±1, t),

−1 < x < 1,

(3.1.1a)
(3.1.1b)
(3.1.1c)

Here a(x, t), h(x, t) denote the activator and the inhibitor concentrations, respectively. Also ϵ, D, µ > 0 and τ ≥ 0 are the activator diffusivity, inhibitor diffusivity,
inhibitor decay rate and the reaction time constant, respectively. The exponents

33

3.1. The Gierer-Meihardt model

Chapter 3. Approach and Methodology

(p, q, m, s) satisfy [3, 23, 27]
p > 1,

q > 0,

m > 0,

s ≥ 0,

qm
− (s + 1) > 0.
p−1

(3.1.2)

The system has been known to have a stable (spike-type) stationary solution for
0 < ϵ ≪ 1 and 0 < D < D1 (ϵ) ∼ ϵ2 e2/ϵ [2, 3, 27], for over two decades now,
and is well studied in the literature. Moreover, the eigenvalues of the the system
obtained by the linearisation of (3.1.1) around a stationary solution are proven to
have negative real parts [2]. If we recall from the previous chapter, that is a property
we desired the NLSE system to have, in order to obtain them numerically.
There also exists a generalisation of the above system to fractional derivatives,
as studied in [3]:
∂t a = ϵγ − (−∇2 )γ/2 a − a +

ap
,
hq

−1 < x < 1,

t > 0,

1 am
, −1 < x < 1, t > 0,
ϵ hs
∂x h(x = ±1, t), a(x, 0) = a0 (x), h(x, 0) = h0 (x),

τ0 ∂t h = D ∂x2 h − µh +
∂x a(x = ±1, t),

a(x, t) = a(x + 2k, t),

h(x, t) = h(x + 2k, t),

x ∈ R, t > 0, k ∈ Z.

(3.1.3a)
(3.1.3b)
(3.1.3c)
(3.1.3d)

In [3], the range 1 ≤ γ ≤ 2 of the order of fractional derivatives was explored.
It was found that the spike type stationary solution still existed, although with
slightly different properties. For 1 ≤ γ < 2, spike type stationary solution produced
algebraically decaying tails, compared to an exponential decay in the γ = 2 case.
The 1 ≤ γ < 2 solutions were also comparatively less stable, but were found to be
stable nonetheless.
It must be noted that in (3.1.3d), we take values of x from all of R to make
sure that the fractional derivative can be evaulated. Fractional derivatives, as can
be recalled from section (1.3), are nonlocal operators and require functions to be
defined over the interval of their definition, which is R for the above case. The
functions a(x, t) and h(x, t) are made periodic to agree with the domain of interest.
The Neumann boundary conditions (3.1.3c) makes sure that the function is smooth
throughout.
For eqns. (3.1.3), one can find an equation that the spike-type solution follows
by constructing an asymptotic expansion of the solution near the centre of the spike
(the spike inner region). The analysis of the region eventually leads to the ground
state equation −(−∇2 )γ/2 u − u + up = 0, which assumes a form similar to that of
the NLSE. Now, if one tries to solve this equation directly using the same method
34

3.2. The auxiliary PDE method

Chapter 3. Approach and Methodology

as in section 2.3, it is seen that the solution diverges very quickly. To stabilise the
solution, a new system of augmented PDEs is introduced [3], which provides us with
the conceptual basis for the process which we can use to stabilise our system.
The Gierer-Meinhardt system, hence, would be an excellent candidate to take
an inspiration from, in order to construct a system that stabilises the Jacobi elliptic
solutions. This can only prove useful though, if we are able to connect it to the
NLSE. In the following section, a system of auxiliary coupled partial differential
equations will be introduced, that will be of the similar form as system (3.1.1).
This will be done in the hope to obtain a stable stationary solution from them and
connect it to the Jacobi elliptic solution of the NLSE.

3.2

The auxiliary PDE method

Recall that the Nonlinear Schrödinger equation is given by
1
ιΨt = − Ψxx + Ψ|Ψ|2 .
2

(3.2.1)

As in the section (2.2), we transform the NLSE to incorporate solutions of the type
Ψ(x, t) = e−ιωt ψ(x, t).
Hence, ψ(x, t) satisfies
1
ιψt = −ωψ − ψxx + ψ|ψ|2 ,
2
and is known to have stationary solutions related to the Jacobi elliptic functions
(2.1.9), which we are trying to obtain numerically. We expect ω = (1 + k 2 )/2, as in
the Jacobi elliptic solution, cf. equation (2.1.8). To achieve this, we introduce the
following system of auxiliary PDEs in two functions, ψ and ϕ:
2

ι

ψψ
∂ψ
1 ∂2
+ ωψ ψ +
ψ
−
2 = 0,
∂t
2 ∂x2
ϕ

(3.2.2a)

∂ϕ
D ∂2
2
ιτ0
+ ωϕ ϕ +
ϕ − ψ ψ = 0.
(3.2.2b)
2
∂t
2 ∂x
Next, following an identical construction for the Gierer-Meinhardt system [27], we
assume an asymptotic expansion of ψ and ϕ in 1/D, and take D ≫ 1 to limit the
analysis to linear terms. The expansions are as follows:
ψ ∼ ψ0 +

ψ1
1
+ O( 2 ),
D
D
35

ϕ ∼ ϕ0 +

ϕ1
1
+ O( 2 ).
D
D

(3.2.3)

3.2. The auxiliary PDE method

Chapter 3. Approach and Methodology

The following table contains the O(1) and O(D) terms for both eqs. (3.2.2a) and (3.2.2b):

Equation (3.2.2a)
Equation (3.2.2b)

O(D)

O(1)

-

ι∂t ψ0 + ωψ ψ0 + 21 ∂x2 ψ0 − ψ0 |ψ0 |2 /|ϕ0 |2 = 0
ιτ0 ∂t ϕ0 + ωϕ ϕ0 + 21 ∂x2 ϕ1 − ψ0 |ψ0 |2 = 0

∂x2 ϕ0 = 0

Now, The O(D) term of eq. (3.2.2b) along with the no flux boundary condition
implies that ϕ0 is only ϕ0 (t), i.e., it is a constant function in x. O(1) term of equation
(3.2.2b) at t → ∞ implies
1 ∂2
ϕ1 = ψ0 |ψ0 |2 − ωϕ ϕ0 ,
2 ∂x2
on which we now impose the solvability condition, which yields
1
ϕ0 =
2K(k) ωϕ

Z 2K(k)

ψ0 |ψ0 |2 dx,

(3.2.4)

0



where 0,2K(k) is the interval of
 integration, as per the definition of the inner
R
a
1
product ⟨u, v⟩ = (b−a)
u · v dx . K(k) is the complete elliptic integral of the
b
first kind, with elliptic modulus k. It is also the half-period of the Jacobi sine
sn(x, k).
Now, O(1) term of equation (3.2.2a) implies
ωψ ψ0 = −

ψ0 |ψ0 |2
1 ∂2
ψ
+
.
0
2 ∂x2
|ϕ0 |2

Since ϕ0 is just a constant at steady state, we can simplify the expression by substituting ψ0 = cη and make the equation independent of ϕ0 by choosing c = ϕ0 .
This yields that η then must follow
ωψ η = −

1 ∂2
2
η+η η ,
2
2 ∂x

(3.2.5)

which is exactly the NLSE. Hence, we now have a method to connect the NLSE with
the auxiliary system (3.2.2). Once we obtain a solution to the auxiliary system, we
just need to normalise it with c = |ϕ0 | to get a solution to the NLSE, which we
expect to be the Jacobi elliptic.
A numerical trial was attempted for the above problem, details of which will
be discussed in the next section. It must be noted that the parameter values were
chosen so that they closely resembled the ones from (3.1.1). We allowed ourselves to
explore a wide range of values of τ0 . Results of the numerical trial will be presented

36

3.3. Implementation of Numerics

Chapter 3. Approach and Methodology

in the next chapter. Before that, we discuss our numerical implementation of the
problem.

3.3

Implementation of Numerics

The numerical schemes that were used to solve system (3.2.2), were once again a
combination of two-step Adams-Bashforth and Crank-Nicholson algorithms. This
is the same setup we followed in the direct numerical trial of the NLSE (figure 2.2)
given in the previous chapter. We recall the equations governing the above schemes
are as follows.
The two-step Adams-Bashforth scheme for an ODE of the shape y ′ = f (y, t) is given
by
yn+1 = yn +



h
′
3yn′ − yn−1
,
2

(3.3.1)

and the Crank-Nicholson method, given by
yn+1 = yn +



h ′
yn+1 + yn′ .
2

(3.3.2)

The subscript n denotes the evaluation of the term at nth time step. We shall stick
with this notation hereafter.
We start our approach, once again, by performing a Fourier transform on the
Laplacian and turn our system into an ordinary differential equation in t. The
Fourier transformed versions of equations (3.2.2a) and (3.2.2b) are:
∂ ψ̂
1
ι
= −ωψ ψ̂ − q 2 ψ̂ + F
∂t
2

(

ψψ
|ϕ|2

2

)
,

(3.3.3a)

n
o
∂ ϕ̂
D 2
2
ιτ0
= −ωϕ ϕ̂ − q ϕ̂ + F ψ ψ ,
(3.3.3b)
∂t
2
where hats denote Fourier transformed quantities and q is the transform variable.
We choose to progress the nonlinear terms using the two-step Adams-Bashforth algorithm (3.3.1) and the remaining using Crank-Nicholson algorithm (3.3.2). Hence,
using (3.3.3) as the integrand for (3.3.1) and (3.3.2), we get the following as the
time-step progressions:

37

3.3. Implementation of Numerics

Chapter 3. Approach and Methodology

"
1
h
1
ιψ̂n+1 = ιψ̂n +
− ωψ ψ̂n+1 − ωψ ψ̂n − q 2 ψ̂n+1 − q 2 ψ̂n
2
2
2
(
)
(
)
2
2
ψψ
ψψ
+ 3F
−F
|ϕ|2
|ϕ|2
n

#
,

n−1

"
D
D
h
− ωϕ ϕ̂n+1 − ωϕ ϕ̂n − q 2 ϕ̂n+1 − q 2 ϕ̂n
ιτ0 ϕ̂n+1 = ιτ0 ϕ̂n +
2
2
2
#

o
o
n
n
2
2
−F ψ ψ
+ 3F ψ ψ

.

n−1

n

Collecting like terms and simplifying, we get

ψ̂n+1 =

"

1
ι + h4 q 2 + h2 ωψ


h 2 h
ι − q − ωψ ψ̂n +
4
2
(
)
)
(
2
2
ψψ
ψψ
h
−F
3F
2
|ϕ|2
|ϕ|2




n

ϕ̂n+1 =

"

1
ωϕ
ιτ0 + h4 q 2 + hD
2




,

n−1


hD 2 h
q − ωϕ ϕ̂n +
ι−
4
2
(
)
(
)
h
2
2
3F ψ ψ
−F ψ ψ
2
n

!#

!#
.

n−1

At the end, we shall invert the Fourier transform and normalise the solution according to the prescription given in section 3.2 and thereby recover a solution to the
NLSE. As can be noted from the algorithms, we require values from the nth and the
(n − 1)th step for the computation of the (n + 1)th step. This is because we chose to
use the two-step Adams-Bashforth method. However, this feature will imply that
if the initial conditions are at n = 1, then we do not have enough previous steps
to compute the step at n = 2. This problem is circumvented by progressing the
nonlinear terms with the Euler algorithm at n = 2. This completes all the details
that will be needed to implement a numerical trial to solve 3.2.2. We shall present
the results of this numerical trial in the next chapter.

38

Chapter 4
Results
In the previous chapter we set up the details of the numerical trial for finding
solutions to the auxiliary PDE system 3.2.2. We first converted the system of PDEs
in x and t into ODEs in t, by performing a Fourier transform on the variable x.
Then, we made use of a combination of AB2 and CN algorithms , resulting in an
algorithm to progress in t. In this chapter, we shall present the solutions obtained
from this setup for various parameter values.
In anticipation of the Jacobi elliptic solutions, we set up the initial conditions of
the algorithm as the following
ψn=1 = ν k sn(x, k),
Z 2K(k)
1
ψ|ψ|2 dx .
ϕn=1 =
2K(k) ωϕ 0

(4.1)
(4.2)

Here, ν ≈ 1 but ̸= 1 as we do not wish to start exactly at the solution itself, but
slightly away from it, in order to be able to observe a clear convergence towards the
desired solution. We also choose k = 0.8, D = 1000, ν = 0.5
and ωϕ = ωψ = 12 (1+k 2 ),
0.8
as in (2.1.8). The number of grid points in the spatial domain, npts, was set to 210
while the time variable ran from tmin = 0 to tmax = 10, in steps of dt = 0.01. All
of the computation was carried out in Wolfram Mathematica (v12.0.0.0) on a Mac
OS X x86 (64-bit) platform. The plots were obtained using Mathematica’s default
Plot[] and ListPlot[] commands.
Figure 4.1 shows the output of the algorithm for different values of τ0 . Even
complex values of τ0 were tried out. However, all of the obtained solutions are
convergent to the trivial instead of the Jacobi elliptic, marked by the orange curves
which implies that the given prescription has failed to stabilise the Jacobi elliptic
solution and a different candidate is needed for the same. Nonetheless, certain values
of τ0 produced interesting evolution scenarios. One such case is listed in fig. 4.2
39

Chapter 4. Results

(a) τ0 = 0.1

(b) τ0 = 0.02

(c) τ0 = −ι

(d) τ0 = −0.01ι

(e) τ0 = −1 − 0.01ι

(f) τ0 = 0.01 − 0.03ι

Figure 4.1: Obtained outputs for a numerical implementation of the auxiliary PDE
method. Panels (a) - (f) show the evolution for different values of τ0 . In each
panel, ψnum (x, t) is shown in blue at different values of t, each indicated above their
corresponding curves. The Jacobi elliptic solution ψJE (x) = 0.8 sn(x, 0.8) is also
displayed alongside in orange, for comparison. The spatial domain is 0 ≤ x ≤
2K(0.8).
(τ0 = 0.035). As can be clearly seen, the graph of ψnum (x, t) seems to converge to 0
by t = 9.3, before attaining large values again only to be seen going down to 0 by
t = 99.9.

40

Chapter 4. Results

(a) t = 0

(b) t = 0.82

(c) t = 0.93

(d) t = 1.04

(e) t = 2.44

(f) t = 9.99

Figure 4.2: Panels (a)-(f) display the obtained solution (blue) ψnum (x, t) for τ0 =
0.035 at different values of t. The Jacobi elliptic solution ψJE (x) = 0.8 sn(x, 0.8) is
also displayed alongside in orange, for comparison. The spatial domain is 0 ≤ x ≤
2K(0.8).
Hence, it is apparent from the results given in figs. 4.1 and 4.2 that achieving
stabilisation via this scheme is unlikely. The NLSE must be sufficiently different
from the GM model to elude this attempt.

41

Chapter 5
Conclusion and Future prospects
Although objectives 1 and 2 of the computational aspect of the project were not
achieved in the allocated duration of this project, a wide range of topics in the field
of dynamics were explored in review of the literature required for the project.
We came across multiple famous partial differential equation models like the
nonlinear Schrödinger equation, the complex Ginzburg-Landau equation and the
Gierer-Meinhardt system, all bearing tremendous usefulness in modelling of natural
phenomena. An extensive discussion on the nonlinear Schrödinger equation was
presented. The Jacobi elliptic solutions were introduced as stationary solutions of
the NLSE and their (purely imaginary) linear stability spectrum was discussed.
We also presented an enquiry into the intriguing notion of fractional order derivatives, wherein we saw four out of a multitude of its plausible definitions that appear
in literature. We also noted its immense success in various practical applications
and the recent upsurge in the interest in exploring fractional versions of popular
differential equations. This is where our interest in exploring the solutions to the
nonlinear Schrödinger equation of a fractional order was conceived. To do that, we
first wished to obtain the Jacobi elliptic solutions of the NLSE, numerically. However, the spectral nature of the Jacobi solutions prevented us from doing so directly.
Hence, a novel scheme, inspired by the Gierer-Meinhardt system (3.1.1), was developed in hopes of stabilising the Jacobi elliptic solution. The scheme involved
development of an auxiliary system of PDEs (3.2.2), which resembled the GM system. The auxiliary PDE system was then solved numerically for various values of
different parameters.
A comprehensive investigation of an adequate range of parameters ν and τ0 was
made and all of them resulted in a convergence to the trivial solution. So, we
are forced to conclude that the mentioned approach might require a much deeper
numerical or, even better, an analytical insight to glean any details about the exis42

Chapter 5. Conclusion and Future prospects
tence and/or the location of the basins of attraction of the stationary solutions of
the auxiliary PDE system (3.2.2). Such a basin, of course, may not even exist as
(3.2.2) merely bears resemblance to the well established Gierer-Meinhardt system
(3.1.1) on visual grounds. A better auxiliary system may be found if one is able to
circumvent the issue created by the imaginary factor ι on the LHS of the NLSE,
while attempting to connect the NLSE and the Gierer-Meinhardt system (3.1.1) in
a fashion similar to the approach in section (3.2). One can also search for better
candidates for processes that can stabilise the Jacobi elliptic solution. Once a working algorithm for obtaining the Jacobi elliptic from the NLSE is found, one can
generalise the process to obtain solutions to the fractional Nonlinear Schrödinger
equation (fNLSE).

43

Appendices

44

Appendix A
Fourier Transforms
Fourier transform (FT) is one of many available types of integral transforms in
mathematics. It converts a function of time (space) into a form that describes the
frequencies (wave numbers) present in the original function (i.e. the spectrum). FTs
have found applicability in various branches of mathematics, physics and engineering. In the context of this thesis, FTs have primarily been used as a tool to convert
difficult-to-solve partial differential equations into ordinary differential equations,
which are far easier to deal with. We shall take a look at the definition and some
properties of FTs along with a brief introduction to Discrete Fourier Transforms.

A.1

Definitions

There are multiple conventions that are followed to define FTs, depending on the
context. We shall stick to the convention that has the formulas of FT and the
inverse Fourier Transform symmetric (as it makes the definitions easy to remember
and follow).
Definition A.1 (Fourier Transform). The Fourier Transform fˆ(ξ) of a function
f (x); x, ξ ∈ R, is given by
1
fˆ(ξ) = √
2π

Z +∞

f (x)e−ιξx dx

(A.1.1)

−∞

provided that f is Lebesgue-measurable on the real line [28], i.e.,

R∞
−∞

|f (x)|dx < ∞.

The inverse Fourier Transform can also be defined in a similar fashion
Definition A.2 (Inverse Fourier Transform). The inverse Fourier Transform

45

A.2. Properties of Fourier Transforms

Appendix A. Fourier Transforms

f (x) of a function fˆ(ξ); x, ξ ∈ R, is given by
1
f (x) = √
2π

Z +∞

fˆ(ξ)eιξx dξ

(A.1.2)

−∞

The pair of functions f and fˆ is often termed as a Fourier transform pair. There
are a few commonly used short-hand notations to denote FTs

fˆ(ξ) = F f (x) x→ξ


or simply fˆ(ξ) = F f (x) ,

and similarly for Inverse FT

f (x) = F −1 fˆ(ξ) ξ→x


or f (x) = F −1 fˆ(ξ) .

Note that Fourier Transforms may be defined in any number of space (and\or time)
dimensions. In the case where both space and time variables are present,
it is customary to define the Fourier transform in a manner that treats
the time variables distinctively from all spatial variables. For instance, here
is how one may define a FT in 3 + 1 dimensions
Definition A.3 (Fourier Transform in 3 + 1 dimensions). The Fourier Transform
fˆ(k, ω) of a function f (x, t); x, k ∈ R3 ; t, ω ∈ R, is given by
fˆ(k, ω) =

1
(2π)2

ZZ

f (x, t)e−ι(k·x−ωt) d3 xdt

(A.1.3)

R4

and similarly
Definition A.4 (Inverse Fourier Transform in 3 + 1 dimensions). The Inverse
Fourier Transform f (x, t) of a function fˆ(k, ω); x, k ∈ R3 ; t, ω ∈ R, is given by
1
f (x, t) =
(2π)2

ZZ

fˆ(k, ω)eι(k·x−ωt) d3 kdω.

(A.1.4)

R4

It may be noted that the time variable in the above definitions, is being treated
in exactly the opposite fashion as the spatial variables.

A.2

Properties of Fourier Transforms

There are a few properties of FTs that are of particular interest for the context of
this thesis, they are (in no particular order of importance)
46

A.3. Discrete Fourier Transforms

Appendix A. Fourier Transforms

Property A.1 (Linearity). If a and b are any two given complex scalars, then



F ag(x) + bh(x) = aF g(x) + bF h(x) .
Property A.2 (Differentiation). Suppose f (x) is continuously differentiable and
both f and its derivative f ′ are integrable, then
F


d
f (x) = (ιξ)F f (x) ,
dx

and in general

 dn

n
f
(x)
=
(ιξ)
F
f (x) .
dn x
where n ∈ N. In fact, the Fourier transform properties of the Riesz fractional
derivatives enables us to generalise this property even to derivatives of non-integral
orders:


F R Dαx f (x) = −|ξ|α F f (x) ,
F

where R Dαx denotes the Riesz fractional derivative of order α ∈ R, in the variable x.
Property A.3 (Convolution theorem). If f (x) and g(x) are integrable functions,
then



F f (x) · F g(x) = F (f ∗ g)(x) ,
R∞
where ∗ denotes convolution (f ∗ g)(x) = −∞ f (y)g(x − y)dy.

A.3

Discrete Fourier Transforms

For FTs to have any practical viability, there must exist a way of performing the
transform (and its inverse), numerically. And since numerics can only be performed
on discrete sets, we need a discrete transform that converts a list of numbers into
another, such that the desired spectral information about the former list is captured
in the latter. This can be achieved by defining
Definition A.5 (Discrete Fourier Transform). The Discrete Fourier Transform
(DFT) vs of a list ur of length n is defined to be
n

1 X
vs = √
ur e2πι(r−1)(s−1)/n .
n r=1
and

47

(A.3.1)

A.3. Discrete Fourier Transforms

Appendix A. Fourier Transforms

Definition A.6 (Discrete Inverse Fourier Transform). The inverse Discrete Fourier
Transform (inverse DFT) ur of a list vs of length n is defined to be
n

1 X −2πι(r−1)(s−1)/n
vs e
.
ur = √
n r=1

(A.3.2)

It must be noted that just like FTs and inverse FTs, many definitions of DFT
and inverse DFT exist (again often involving different conventions). The definition
A.5 here, is the same as in the documentation of Wolfram Mathematica [29]. This
definition has been chosen here because Wolfram Mathematica was used for all
computation aspects of this project.

48

Appendix B
The Jacobi elliptic trigonometric
functions
The Jacobi elliptic functions [30] are a set of special functions which appear as
solutions to differential equations of the type
d2 y
= A + By + Cy 3 ,
dx2
where all of A, B, and C are real constants. The three basic Jacobi elliptic trig
functions are sn(x, k), cn(x, k) and dn(x, k), where k is known as the elliptic modulus.
These may also be defined using the elliptic integral of the first kind,
Z ϕ
x = F (ϕ, k) =
0

dt
p
.
1 − k 2 sin2 t

Using this, we define
ϕ = F −1 (x, k) =: am(x, k),
where am(x, k) is the Jacobi amplitude. Now, sn(x, k), cn(x, k) and dn(x, k) may be
defined using am(x, k) in the following fashion
sn(x, k) := sin ϕ = sin(am(x, k)),
and similarly
cn(x, k) := cos ϕ = cos(am(x, k)).

49

Appendix B. The Jacobi elliptic trigonometric functions
Then we have
p
1 − k 2 sn2 (x, k)
q
= 1 − k 2 sin2 (am(x, k)) .

dn(x, k) :=

From the above definitions, it must be clear that if one defines the complete elliptic
integral of the first kind as
Z π/2
K(k) =
0

dt
p
,
1 − k 2 sin2 t

then we have am(K, k) = π/2. Hence, the Jacobi elliptic functions sn(x, k), cn(x, k) and
dn(x, k) will all be periodic in K(k) w.r.t. their first argument as
sn(x + 2mK, k) = (−1)m sn(x, k),
cn(x + 2mK, k) = (−1)m cn(x, k)
and
dn(x + 2mK, k) = dn(x, k),
where m is a positive integer.

50

Bibliography
1. Bottman, N., Deconinck, B. & Nivala, M. Elliptic solutions of the defocusing
NLS equation are stable. Journal of Physics A: Mathematical and Theoretical
44, 285201. doi:10.1088/1751-8113/44/28/285201 (June 2011).
2. Iron, D., Ward, M. & Wei, J. The stability of spike solutions to the onedimensional Gierer–Meinhardt model. Physica D: Nonlinear Phenomena 150,
25–62. doi:10.1016/S0167-2789(00)00206-2 (Mar. 2001).
3. Nec, Y. Spike-Type Solutions to One Dimensional Gierer–Meinhardt Model
with Lévy Flights. Studies in Applied Mathematics 129, 272–299. doi:10.1111/
j.1467-9590.2012.00550.x (2012).
4. Nec, Y. Pattern formation in systems with anomalous diffusion PhD thesis
(Technion – Israel Institute of Technology, 2009).
5. Malomed, B. Encyclopedia of Nonlinear Science (Scott, Alwyn (ed.), New York:
Routledge, pp. 639–643, 2005).
6. Hasegawa, A. & Kodama, Y. Solitons in optical communications Akira Hasegawa
and Yuji Kodama English, xiv, 320 p. : (Clarendon Press ; Oxford University
Press Oxford : Oxford ; New York, 1995).
7. Zakharov, V. E. Collapse of Langmuir Waves. Soviet Journal of Experimental
and Theoretical Physics 35, 908 (Jan. 1972).
8. Ablowitz, M. J., Prinari, B. & Trubatch, A. D. Discrete and Continuous Nonlinear Schrödinger Systems doi:10.1017/CBO9780511546709 (Cambridge University Press, 2003).
9. Herrmann, R. Fractional Calculus – An Introduction for Physicists (3rd revised
and extended Edition), World Scientific Publishing, Singapore, September 2018,
636 pp, ISBN: 978-981-3274-57-0 636 (Sept. 2018).
10. Oldham, K. B. & Spanier, J. The fractional calculus; theory and applications
of differentiation and integration to arbitrary order by Keith B. Oldham and
Jerome Spanier English, xiii, 234 p. (Academic Press New York, 1974).
51

Bibliography

Bibliography

11. Mandelbrot, B. B. The fractal geometry of nature (Freeman, San Francisco,
CA, 1982).
12. Mandelbrot, B. B. Self-affine fractals and fractal dimension. Phys. Scr. 32, 257
(1985).
13. Metzler, R. & Klafter, J. The random walk’s guide to anomalous diffusion: a
fractional dynamics approach. Physics Reports 339, 1–77. doi:10.1016/S03701573(00)00070-3 (2000).
14. Eliazar, I. & Klafter, J. Anomalous is ubiquitous. Annals of Physics 326, 2517–
2531. doi:10.1016/j.aop.2011.07.006 (2011).
15. Laskin, N. Fractional Schrödinger equation. Phys. Rev. E 66, 056108. doi:10.
1103/PhysRevE.66.056108 (5 Nov. 2002).
16. Goldfain, E. Complexity in quantum field theory and physics beyond the standard model. Chaos, Solitons and Fractals 28. Sir Hermann Bondi 1919–2005,
913–922. doi:10.1016/j.chaos.2005.09.012 (2006).
17. Lim, S. Fractional derivative quantum fields at positive temperature. Physica
A: Statistical Mechanics and its Applications 363, 269–281. doi:10.1016/j.
physa.2005.08.005 (2006).
18. Bayın, S. Ş. Definition of the Riesz derivative and its application to space
fractional quantum mechanics. Journal of Mathematical Physics 57, 123501.
doi:10.1063/1.4968819 (Dec. 2016).
19. Naber, M. Time fractional Schrödinger equation. Journal of Mathematical
Physics 45, 3339–3352. doi:10.1063/1.1769611 (Aug. 2004).
20. Laskin, N. ”Fractional quantum mechanics and Lévy path integrals”. Physics
Letters A 268, 298–305. doi:10.1016/s0375-9601(00)00201-2 (Apr. 2000).
21. Some applications of fractional equations. Communications in Nonlinear Science and Numerical Simulation 8. Chaotic transport and compexity in classical and quantum dynamics, 273–281. doi:10.1016/S1007-5704(03)00049-2
(2003).
22. Burden, R., Faires, J. & Burden, A. Numerical Analysis 10th Edition (Cengage
Learning, 2015).
23. Gierer, A. & Meinhardt, H. A theory of biological pattern formation. Kybernetik
12, 30–39. doi:10.1007/BF00289234 (1972).
24. Meinhardt, H. Models of Biological Pattern Formation (Academic Press, London, 1982) (May 1982).
52

Bibliography

Bibliography

25. Huang, X., Dong, Y. & Zhao, J. HetR homodimer is a DNA-binding protein
required for heterocyst differentiation, and the DNA-binding activity is inhibited by PatS. Proceedings of the National Academy of Sciences 101, 4848–4853.
doi:10.1073/pnas.0400429101 (2004).
26. Meinhardt, H. & Klingler, M. A model for pattern formation on the shells
of molluscs. Journal of Theoretical Biology 126, 63–89. doi:10.1016/S00225193(87)80101-7 (1987).
27. WARD, M. J. & WEI, J. Hopf bifurcation of spike solutions for the shadow
Gierer–Meinhardt model. European Journal of Applied Mathematics 14, 677–
711. doi:10.1017/S0956792503005278 (2003).
28. Boas, M. L. Mathematical methods in the physical sciences; 3rd ed. (Wiley,
Hoboken, NJ, 2006).
29. Wolfram language function: “Fourier”, Wolfram Language and System Documentation Center. Wolfram Research. https://reference.wolfram.com/
language/ref/Fourier.html.
30. Weisstein, E. W. “Jacobi Elliptic Functions” From MathWorld–A Wolfram
Web Resource. https://mathworld.wolfram.com/JacobiEllipticFunctions.
html.

53