A high-order nodal discontinuous Galerkin method for nonlinear fractional Schr\"{o}dinger type equations

TL;DR

This paper introduces a high-order nodal discontinuous Galerkin method tailored for solving nonlinear fractional Schrödinger equations, demonstrating stability, optimal convergence, and verified through numerical experiments.

Contribution

The paper develops a novel high-order discontinuous Galerkin method for nonlinear fractional Schrödinger equations, proving stability and convergence, with numerical validation.

Findings

Proves $L^{2}$ stability for the proposed method.

Achieves optimal convergence order of $O(h^{N+1})$.

Numerical experiments confirm theoretical convergence rates.

Abstract

We propose a nodal discontinuous Galerkin method for solving the nonlinear Riesz space fractional Schr\"{o}dinger equation and the strongly coupled nonlinear Riesz space fractional Schr\"{o}dinger equations. These problems have been expressed as a system of low order differential/integral equations. Moreover, we prove, for both problems, $L^{2}$ stability and optimal order of convergence $O(h^{N+1})$, where $h$ is space step size and $N$ is polynomial degree. Finally, the performed numerical experiments confirm the optimal order of convergence.

Tables

Table 1. Table 1: L 2 superscript 𝐿 2 L^{2} -Error and order of convergence for Example 7.1 with K 𝐾 K elements and polynomial order N 𝑁 N .

N	N=1          N=2          N=3
K	$L^2$-Error	order	K	$L^2$-Error	order	K	$L^2$-Error	order
64	1.57e-02	-	35	8.47e-05	-	20	1.59e-05	-
74	1.24e-02	1.63	45	3.97e-05	3.0	40	9.82e-07	4.02
84	9.2e-03	2.33	90	5.67e-06	2.81	60	2.14e-07	3.75

Table 2. Table 2: L 2 superscript 𝐿 2 L^{2} -Error and order of convergence for Example 7.2 with K 𝐾 K elements and polynomial order N 𝑁 N .

N	N=1          N=2          N=3
K	$L^2$-Error	order	K	$L^2$-Error	order	K	$L^2$-Error	order
120	1.41e-01	-	60	1.52e-04	-	40	7.02e-06	-
135	1.09e-02	2.15	80	6.54e-05	2.89	70	7.62e-07	3.97
150	8.9e-03	1.92	120	1.78e-05	3.22	90	2.6e-07	4.28

Table 3. Table 3: L 2 superscript 𝐿 2 L^{2} -Error and order of convergence for u 1 subscript 𝑢 1 u_{1} with K 𝐾 K elements and polynomial order N 𝑁 N .

N	N=1          N=2          N=3
K	$L^2$-Error	order	K	$L^2$-Error	order	K	$L^2$-Error	order
92	2.27 e-02	-	60	1.93e-04	-	50	4.1e-06	-
100	1.99e-02	1.54	90	5.60e-05	3.01	70	1.23e-06	3.58
130	1.07e-02	2.37	110	3.0e-05	3.12	100	2.98e-07	3.96

Table 4. Table 4: L 2 superscript 𝐿 2 L^{2} -Error and order of convergence for u 2 subscript 𝑢 2 u_{2} with K 𝐾 K elements and polynomial order N 𝑁 N .

N	N=1          N=2          N=3
K	$L^2$-Error	order	K	$L^2$-Error	order	K	$L^2$-Error	order
92	2.25e-02	-	60	1.7481e-04	-	50	3.87e-06	-
100	1.92e-02	1.9	90	5.03e-05	3.07	70	8.91e-07	4.37
130	1.12e-02	2.04	110	2.67e-05	3.16	100	2.4e-07	3.68

Table 5. Table 5: L 2 superscript 𝐿 2 L^{2} -Error and order of convergence for u 1 subscript 𝑢 1 u_{1} with K 𝐾 K elements and polynomial order N 𝑁 N .

N	N=1          N=2          N=3
K	$L^2$-Error	order	K	$L^2$-Error	order	K	$L^2$-Error	order
96	1.90 e-02	-	30	4.7e-04	-	40	8.68e-06	-
120	1.27e-02	2.35	60	1.47e-04	2.86	60	1.79e-06	3.89
135	9.6e-03	1.92	130	1.22e-05	3.22	80	6.03e-07	3.78

Table 6. Table 6: L 2 superscript 𝐿 2 L^{2} -Error and order of convergence for u 2 subscript 𝑢 2 u_{2} with K 𝐾 K elements and polynomial order N 𝑁 N .

N	N=1          N=2          N=3
K	$L^2$-Error	order	K	$L^2$-Error	order	K	$L^2$-Error	order
96	1.89e-02	-	40	4.18e-04	-	40	7.71e-06	-
120	1.34e-02	1.55	60	1.26e-04	2.95	60	1.47e-06	4.08
135	1.03e-02	2.22	130	1.21e-05	3.04	90	5.1e-07	3.7

Equations

$$\begin{array}{rl}& i\frac{\partial u}{\partial t}-{\lambda }_1(-\Delta {)}^{\frac{\alpha }{2}}u+{\lambda }_{2}f(|u{|}^2)u=0,\\ & u(x,0)=u_0(x),\end{array}$$

$$\begin{array}{rl}& i\frac{\partial u}{\partial t}-{\lambda }_1(-\Delta {)}^{\frac{\alpha }{2}}u+{\varpi }_{1}u+{\varpi }_{2}v+{\lambda }_{2}f(|u{|}^2,|v{|}^2)u=0,\\ & i\frac{\partial v}{\partial t}-{\lambda }_3(-\Delta {)}^{\frac{\alpha }{2}}v+{\varpi }_{2}u+{\varpi }_{1}v+{\lambda }_{4}g(|u{|}^2,|v{|}^2)v=0,\\ & u(x,0)=u_0(x),\\ & v(x,0)=v_0(x),\end{array}$$

$$-(-\Delta {)}^{\frac{\alpha }{2}}u(x,t)={\mathcal{F}}^{-1}(|\xi {|}^{\alpha }\hat{u}(\xi ,t))$$

$$\big{(}{}^{\,\,RL}_{-\infty}\mathcal{I}_{x}^{\alpha}f\big{)}(x)=\frac{1}{\Gamma(\alpha)}\int_{-\infty}^{x}\frac{f(s)ds}{(x-s)^{1-\alpha}},\quad x>-\infty,$$

$$\big{(}{}^{RL}_{\,\,\,\,x}\mathcal{I}_{\infty}^{\alpha}f\big{)}(x)=\frac{1}{\Gamma(\alpha)}\int_{x}^{\infty}\frac{f(s)ds}{(s-x)^{1-\alpha}},\quad x<\infty,$$

$$\big{(}{}^{\,\,RL}_{-\infty}\mathcal{D}_{x}^{\alpha}f\big{)}(x)=\frac{d}{dx}\big{(}{}^{\,\,RL}_{-\infty}\mathcal{I}_{x}^{1-\alpha}f\big{)}(x)=\frac{1}{\Gamma(1-\alpha)}\frac{d}{dx}\int_{-\infty}^{x}\frac{f(s)ds}{(x-s)^{\alpha}},\quad x>-\infty,$$

$$\big{(}{}^{RL}_{\,\,\,\,x}\mathcal{D}_{\infty}^{\alpha}f\big{)}(x)=\frac{-d}{dx}\big{(}^{RL}_{\,\,\,\,x}\mathcal{I}_{\infty}^{1-\alpha}f\big{)}(x)=\frac{1}{\Gamma(1-\alpha)}\bigg{(}\frac{-d}{dx}\bigg{)}\int_{x}^{\infty}\frac{f(s)ds}{(s-x)^{\alpha}},\quad x<\infty.$$

$$\big{(}{}^{\,\,RL}_{-\infty}\mathcal{D}_{x}^{\alpha}f\big{)}(x)=\bigg{(}\frac{d}{dx}\bigg{)}^{n}\big{(}{}^{\,\,RL}_{-\infty}\mathcal{I}_{x}^{n-\alpha}f\big{)}(x)=\frac{1}{\Gamma(n-\alpha)}\bigg{(}\frac{d}{dx}\bigg{)}^{n}\int_{-\infty}^{x}\frac{f(s)ds}{(x-s)^{-n+1+\alpha}},\quad x>-\infty,$$

$$\big{(}{}^{RL}_{\,\,\,\,x}\mathcal{D}_{\infty}^{\alpha}f\big{)}(x)=\bigg{(}\frac{-d}{dx}\bigg{)}^{n}\big{(}^{RL}_{\,\,\,\,x}\mathcal{I}_{\infty}^{n-\alpha}f\big{)}(x)=\frac{1}{\Gamma(n-\alpha)}\bigg{(}\frac{-d}{dx}\bigg{)}^{n}\int_{x}^{\infty}\frac{f(s)ds}{(s-x)^{-n+1+\alpha}},\quad x<\infty.$$

$$\big{(}{}^{\,\,\,\,\,\,C}_{-\infty}\mathcal{D}_{x}^{\alpha}f\big{)}(x)=\bigg{(}{}^{RL}_{-\infty}\mathcal{I}_{x}^{n-\alpha}\frac{d^{n}f}{dx^{n}}\bigg{)}(x)=\frac{1}{\Gamma(n-\alpha)}\int_{-\infty}^{x}\frac{f^{(n)}(s)ds}{(x-s)^{-n+1+\alpha}},\quad x>-\infty,$$

$$\big{(}{}^{C}_{\,x}\mathcal{D}_{\infty}^{\alpha}f\big{)}(x)=(-1)^{n}\bigg{(}{}^{RL}_{\,\,\,\,x}\mathcal{I}_{\infty}^{n-\alpha}\frac{d^{n}f}{dx^{n}}\bigg{)}(x)=\frac{1}{\Gamma(n-\alpha)}\int_{x}^{\infty}\frac{(-1)^{n}f^{(n)}(s)ds}{(s-x)^{-n+1+\alpha}},\quad x<\infty.$$

$$\frac{\partial^{\alpha}}{\partial|x|^{\alpha}}u(x,t)=-(-\Delta)^{\frac{\alpha}{2}}u(x,t)=-\frac{{}^{\,\,\,\,\,\,C}_{-\infty}\mathcal{D}_{x}^{\alpha}u(x,t)+{}^{C}_{\,x}\mathcal{D}_{\infty}^{\alpha}u(x,t)}{2\cos\big{(}\frac{\pi\alpha}{2}\big{)}}.$$

$$\Delta_{-\mu/2}u(x)=-\frac{{}^{\,\,\,\,\,\,C}_{-\infty}\mathcal{D}_{x}^{-\mu}u(x)+{}^{C}_{\,x}\mathcal{D}_{\infty}^{-\mu}u(x)}{2\cos\big{(}\frac{\pi(2-\mu)}{2}\big{)}}=\frac{{}^{\,\,\,\,\,\,C}_{-\infty}\mathcal{D}_{x}^{-\mu}u(x)+{}^{C}_{\,x}\mathcal{D}_{\infty}^{-\mu}u(x)}{2\cos\big{(}\frac{\pi\mu}{2}\big{)}}=\frac{{}^{RL}_{-\infty}\mathcal{I}_{x}^{-\mu}u(x)+{}^{RL}_{\,\,\,\,x}\mathcal{I}_{\infty}^{-\mu}u(x)}{2\cos\big{(}\frac{\pi\mu}{2}\big{)}}.$$

$$-(-\Delta)^{\frac{\alpha}{2}}u(x)=\Delta_{\frac{(\alpha-2)}{2}}\bigg{(}\frac{d^{2}u(x)}{dx^{2}}\bigg{)}.$$

$$|u|_{J_{L}^{\alpha}(\mathbb{R})}=\big{\|}{}^{RL}_{\,x_{L}}\mathcal{D}_{x}^{\alpha}u\big{\|}_{L^{2}(\mathbb{R})}.$$

$$\Vert u{\Vert }_{J_L^{\alpha }(\mathbb{R})}=(|u{|}_{J_L^{\alpha }(\mathbb{R})}^2+\Vert u{\Vert }_{L^2(\mathbb{R})}^2{)}^{\frac{1}{2}},$$

$$|u|_{J_{R}^{\alpha}(\mathbb{R})}=\big{\|}{}^{RL}_{\,\,\,\,x}\mathcal{D}_{x_{R}}^{\alpha}u\big{\|}_{L^{2}(\mathbb{R})},$$

$$\Vert u{\Vert }_{J_R^{\alpha }(\mathbb{R})}=(|u{|}_{J_R^{\alpha }(\mathbb{R})}^2+\Vert u{\Vert }_{L^2(\mathbb{R})}^2{)}^{\frac{1}{2}},$$

$$\|u\|_{J_{S}^{\alpha}(\mathbb{R})}=\big{|}\big{(}{}^{RL}_{\,x_{L}}\mathcal{D}_{x}^{\alpha}u,{}^{RL}_{\,\,\,\,x}\mathcal{D}_{x_{R}}^{\alpha}u\big{)}_{L^{2}(\mathbb{R})}\big{|}^{\frac{1}{2}},$$

$$\|u\|_{J_{S}^{\alpha}(\mathbb{R})}=\big{(}|u|_{J_{S}^{\mu}(\mathbb{R})}^{2}+\|u\|_{L^{2}(\mathbb{R})}^{2}\big{)}^{\frac{1}{2}}.$$

$$({}_{-\infty }^{RL}{\mathcal{I}}_x^{s}u,{}_x^{RL}{\mathcal{I}}_{\infty }^{s}u{)}_{\mathbb{R}}=\cos (s\pi )|u{|}_{J_L^{-s}(\mathbb{R})}^2=\cos (s\pi )|u{|}_{J_R^{-s}(\mathbb{R})}^2.$$

$$({\Delta }_{-\mu }u,u{)}_{\mathbb{R}}=|u{|}_{J_L^{-\mu }(\mathbb{R})}^2=|u{|}_{J_R^{-\mu }(\mathbb{R})}^2.$$

$$\Vert u{\Vert }_{L^2(\Omega )}\le C|u{|}_{J_{L,0}^{\alpha }(\Omega )},$$

$$\Vert u{\Vert }_{L^2(\Omega )}\le C|u{|}_{J_{R,0}^{\alpha }(\Omega )}.$$

$$\Vert {}_a^{RL}{\mathcal{I}}_x^{\mu }u{\Vert }_{L^2(\Omega )}\le K\Vert u{\Vert }_{L^2(\Omega )}.$$

$$\Vert {}_x^{RL}{\mathcal{I}}_b^{\mu }u{\Vert }_{L^2(\Omega )}\le K\Vert u{\Vert }_{L^2(\Omega )}.$$

$$\Vert {\Delta }_{-\mu }u{\Vert }_{L^2(\Omega )}\le K\Vert u{\Vert }_{L^2(\Omega )}.$$

$$\begin{array}{rl}& e={\Delta }_{(\alpha -2)/2}r,r=\frac{\partial }{\partial x}s,s=\frac{\partial }{\partial x}u,\end{array}$$

$$\begin{array}{rl}& i\frac{\partial u}{\partial t}+{\lambda }_{1}e+{\lambda }_{2}f(|u{|}^2)u=0,\\ & e={\Delta }_{(\alpha -2)/2}r,r=\frac{\partial }{\partial x}s,s=\frac{\partial }{\partial x}u.\end{array}$$

$$\begin{array}{r}u(x,t)=p(x,t)+iq(x,t),\end{array}$$

Full text

Keywords: nonlinear fractional Schrödinger equation, strongly coupled nonlinear fractional Schrödinger equations, nodal discontinuous Galerkin method, stability, error estimates.

1 Introduction

In this paper we develop a nodal discontinuous Galerkin method to solve the generalized nonlinear fractional Schrödinger equation

[TABLE]

and the strongly coupled nonlinear fractional Schrödinger equations

[TABLE]

and homogeneous boundary conditions. $f(u)$ and $g(u)$ are arbitrary (smooth) nonlinear real functions and $\lambda_{i}$, $i=1,2,3,4$ are a real constants, $\varpi_{1}$ is normalized birefringence constant and $\varpi_{2}$ is the linear coupling parameter which accounts for the effects that arise from twisting and elliptic deformation of the fiber [1]. Notice that the assumption of homogeneous boundary conditions is for simplicity only and is not essential: the method can be easily designed for nonhomogeneous boundary conditions. The fractional Laplacian $-(-\Delta)^{\frac{\alpha}{2}}$, which can be defined using Fourier analysis as [2, 3]

[TABLE]

where $\mathcal{F}$ is the Fourier transform. Equation (1.1) can be viewed as a generalization of the classical nonlinear Schrödinger equation. During the last decade, it has arisen as a suitable model in many application areas, such as fluid dynamics, nonlinear optics, and plasma physics [4, 5, 6]. It was first introduced by Laskin [7, 8], who derived fractional Schrödinger equation with Riesz space-fractional derivative includes a space fractional derivative of order $\alpha\,\,(1<\alpha<2)$ instead of the Laplacian in the classical Schrödinger equation, and obtained its by replacing Brownian trajectories in Feynman path integrals (corresponding to the classical Schrödinger equation) by the Lévy flights. It is generally difficult to give the explicit forms of the analytical solutions of nonlinear fractional Schrödinger equation, thus the construction of numerical methods becomes very important. In recent years, developing various numerical algorithms for solving nonlinear fractional Schrödinger equation has received much attention. For the time-fractional Schrödinger equation, Wei et al.[9] presented and analyzed an implicit fully discrete local discontinuous Galerkin (LDG) finite element method for solving the time-fractional Schrödinger equation. Hicdurmaza and Ashyralyev presented stability analysis for a first order difference scheme applied to a nonhomogeneous multidimensional time fractional Schrödinger differential equation. For the space-fractional Schrödinger equation, Wang and Huang [10] studied an energy conservative Crank-Nicolson difference scheme for nonlinear Riesz space-fractional Schrödinger equation. Yang [11] proposed a class of linearized energy-conserved finite difference schemes for nonlinear space-fractional Schrödinger equation. Galerkin finite element method for nonlinear fractional Schrödinger equations were considered [12]. Amore et.al. [13] developed the collocation method for fractional quantum mechanics.

The strongly coupled nonlinear Schrödinger system (1.2) arise in many physical fields, especially in in fluid mechanics, solid state physics and plasma waves and for two interacting nonlinear packets in a dispersive and conservative system, see, e.g.,[14, 15, 16] and reference therein. When $\alpha=2$, it represents the integer-order strongly coupled equations, and a number of conservative schemes for such case have been proposed [17, 18, 19]. When $\varpi_{1}=\varpi_{2}=0$, this system becomes the weakly coupled nonlinear fractional Schrödinger equations considered in [20, 12] and reference therein. Ran and Zhang [16] proposed a conservative difference scheme for solving the strongly coupled nonlinear fractional Schrödinger equations. A numerical study based on an implicit fully discrete LDG for the time-fractional coupled Schrödinger systems is presented [21]. To the best of our knowledge, however, the LDG method, which is an important approach to solve partial differential equations and fractional partial differential equations, has not been considered for the nonlinear Schrödinger equation and the coupled nonlinear Schrödinger equations with the Riesz space fractional derivative. Compared with finite difference methods, it has the advantage of greatly facilitates the handling of complicated geometries and elements of various shapes and types, as well as the treatment of boundary conditions.

The LDG method is a well-established method for classical conservation laws [22, 23, 24]. For application of the method to fractional problems, Mustapha and McLean [25, 26] have developed and analyzed discontinuous Galerkin methods for time fractional diffusion and wave equations. Xu and Hesthaven [27] proposed a LDG method for fractional convection-diffusion equations. They proved stability and optimal order of convergence $N+1$ for the fractional diffusion problem when polynomials of degree $N$, and an order of convergence of $N+\frac{1}{2}$ is established for the general fractional convection-diffusion problem with general monotone flux for the nonlinear term. Aboelenen and El-Hawary [28] proposed a high-order nodal discontinuous Galerkin method for a linearized fractional Cahn-Hilliard equation. They proved stability and optimal order of convergence $N+1$ for the linearized fractional Cahn-Hilliard problem. Here we propose LDG method for problems (1.1)-(1.2) with the Riesz space fractional derivative of order $\alpha$ $(1<\alpha<2)$. For $1<\alpha<2$, it is conceptually similar to a fractional derivative with an order between $1$ and $2$. We rewrite the fractional operator as a composite of first order derivatives and a fractional integral and convert the nonlinear fractional Schrödinger equation and the strongly coupled nonlinear fractional Schrödinger equations into a system of low order equations. This allows us to apply the LDG method.

The outline of this paper is as follows. In section 2, we introduce some basic definitions and recall a few central results. In section 3, we derive the discontinuous Galerkin formulation for the nonlinear fractional Schrödinger equation. In section 4, we prove a theoretical result of $L^{2}$ stability for the nonlinear case as well as an error estimate for the linear case. In section 5 we present a local discontinuous Galerkin method for the strongly coupled nonlinear fractional Schrödinger equations and give a theoretical result of $L^{2}$ stability for the nonlinear case and an error estimate for the linear case in section 6. Section 7 presents some numerical examples to illustrate the efficiency of the scheme. A few concluding remarks are offered in section 8.

2 Preliminary definitions

We introduce some preliminary definitions of fractional calculus, see, e.g.,[29] and associated functional setting for the subsequent numerical schemes and theoretical analysis.

2.1 Liouville-Caputo Fractional Calculus

The left-sided and right-sided Riemann-Liouville integrals of order $\alpha$, when $0<\alpha<1$, are defined, respectively, as

[TABLE]

and

[TABLE]

where $\Gamma$ represents the Euler Gamma function. The corresponding inverse operators, i.e., the left-sided and right-sided fractional derivatives of order $\alpha$, are then defined based on (2.1) and (2.2), as

[TABLE]

and

[TABLE]

This allows for the definition of the left and right Riemann-Liouville fractional derivatives of order $\alpha$ $(n-1<\alpha<n),\,\,n\in\mathbb{N}$ as

[TABLE]

and

[TABLE]

Furthermore, the corresponding left-sided and right-sided Caputo derivatives of order $\alpha$ $(n-1<\alpha<n)$ are obtained as

[TABLE]

and

[TABLE]

The Riesz fractional derivative is defined as

[TABLE]

If $\alpha<0$, the fractional Laplacian becomes the fractional integral operator. In this case, for any $0<\mu<1$, we define

[TABLE]

When $1<\alpha<2$, using (2.7), (2.8) and (2.10), we can rewrite the fractional Laplacian in the following form:

[TABLE]

To carry out the analysis, we introduce the appropriate fractional spaces.

Definition 2.1

(left fractional space [30]). We define the seminorm

[TABLE]

and the norm

[TABLE]

*and let $J_{L}^{\alpha}(\mathbb{R})$ denote the closure of $C_{0}^{\infty}(\mathbb{R})$ with respect to $\|.\|_{J_{L}^{\alpha}(\mathbb{R})}$.

Definition 2.2

(right fractional space [30]). We define the seminorm

[TABLE]

and the norm

[TABLE]

and let $J_{R}^{\alpha}(\mathbb{R})$ denote the closure of $C_{0}^{\infty}(\mathbb{R})$ with respect to $\|.\|_{J_{R}^{\alpha}(\mathbb{R})}$.

Definition 2.3

(symmetric fractional space [30]). We define the seminorm

[TABLE]

and the norm

[TABLE]

and let $J_{S}^{\alpha}(\mathbb{R})$ denote the closure of $C_{0}^{\infty}(\mathbb{R})$ with respect to $\|.\|_{J_{S}^{\alpha}(\mathbb{R})}$.

Lemma 2.1

(see [30]). For any $0<s<1$, the fractional integral satisfies the following property:

[TABLE]

Lemma 2.2

For any $0<\mu<1$, the fractional integral satisfies the following property:

[TABLE]

Generally, we consider the problem in a bounded domain instead of $\mathbb{R}$. Hence, we restrict the definition to the domain $\Omega=[a,b]$.

Definition 2.4

Define the spaces $J_{R,0}^{\alpha}(\Omega),J_{L,0}^{\alpha}(\Omega),J_{S,0}^{\alpha}(\Omega)$ as the closures of $C_{0}^{\infty}(\Omega)$ under their respective norms.

Lemma 2.3

(fractional Poincar$\acute{e}$-Friedrichs, [30]). For $u\in J_{L,0}^{\alpha}(\Omega)$ and $\alpha\in\mathbb{R}$, we have

[TABLE]

and for $u\in J_{R,0}^{\alpha}(\Omega)$, we have

[TABLE]

Lemma 2.4

(See [31]) For any $0<\mu<1$, the fractional integration operator ${}^{RL}_{\,\,\,\,a}\mathcal{I}_{x}^{\mu}$ is bounded in $L^{2}(\Omega)$:

[TABLE]

The fractional integration operator ${}^{RL}_{\,\,\,\,x}\mathcal{I}_{b}^{\mu}$ is bounded in $L^{2}(\Omega)$:

[TABLE]

Lemma 2.5

The fractional integration operator $\Delta_{-\mu}$ is bounded in $L^{2}(\Omega)$:

[TABLE]

Proof. Combining Lemma 2.4 with (2.10), we obtain the result.

3 LDG method for nonlinear fractional

Schrödinger equation

Let us consider nonlinear fractional Schrödinger equation. To obtain a high order discontinuous Galerkin scheme for the fractional derivative, we rewrite the fractional derivative as a composite of first order derivatives and a fractional integral to recover the equation to a low order system. However, for the first order system, alternating fluxes are used. We introduce three variables $e,r,s$ and set

[TABLE]

then, the nonlinear fractional Schrödinger problem can be rewritten as

[TABLE]

For actual numerical implementation, it might be more efficient if we decompose the complex function $u(x,t)$ into its real and imaginary parts by writing

[TABLE]

where $p$, $q$ are real functions. Under the new notation, the problem (3.2) can be written as

[TABLE]

We consider problems posed on the physical domain $\Omega$ with boundary $\partial\Omega$ and assume that this domain is well approximated by the computational domain $\Omega_{h}$. We consider a nonoverlapping element $D^{k}$ such that

[TABLE]

Now we introduce the broken Sobolev space for any real number $r$

[TABLE]

We define the local inner product and $L^{2}(D^{k})$ norm

[TABLE]

as well as the global broken inner product and norm

[TABLE]

To complete the LDG scheme, we introduce the numerical flux.

The numerical traces $(p,q,s,z)$ are defined on interelement faces as the alternating fluxes [32, 24]

[TABLE]

Note that we can also choose

[TABLE]

For simplicity we discretize the computational domain $\Omega$ into $K$ non-overlapping elements, $D^{k}=[x_{k-\frac{1}{2}},x_{k+\frac{1}{2}}]$, $\Delta x_{k}=x_{k+\frac{1}{2}}-x_{k-\frac{1}{2}}$ and $k=1,...,K$. Let $p_{h},q_{h},e_{h},l_{h},r_{h},s_{h},w_{h},z_{h}\in V_{k}^{N}$ be the approximation of $p,q,e,l,r,s,w,z$ respectively, where the approximation space is defined as

[TABLE]

where $\mathbb{P}(D^{k})$ denotes the set of polynomials of degree up to $N$ defined on the element $D^{k}$. We define local discontinuous Galerkin scheme as follows: find $p_{h},q_{h},e_{h},l_{h},r_{h},s_{h},w_{h},z_{h}\in V_{k}^{N}$, such that for all test functions $\vartheta_{1},\beta_{1},\phi,\varphi,\chi,\beta_{2},\psi,\zeta\in V_{k}^{N}$,

[TABLE]

Applying integration by parts to (3.12), and replacing the fluxes at the interfaces by the corresponding numerical fluxes, we obtain

[TABLE]

4 Stability and error estimates

In the following we discuss stability and accuracy of the proposed scheme, for the nonlinear fractional Schrödinger problem.

4.1 Stability analysis

In order to carry out the analysis of the LDG scheme, we have the following results.

Theorem 4.1

*($L^{2}$ stability). The semidiscrete scheme (3.13) is stable, and $\|u_{h}(x,T)\|_{\Omega_{h}}\leq c\|u_{0}(x)\|_{\Omega_{h}}$ for any $T>0$. *

Proof. Set $(\vartheta_{1},\beta_{1},\phi,\varphi,\chi,\beta_{2},\psi,\zeta)=(p_{h},-r_{h}+e_{h},p_{h},-z_{h},q_{h},l_{h}-w_{h},-q_{h},s_{h})$ in (3.13), and consider the integration by parts formula \big{(}u,\frac{\partial r}{\partial x}\big{)}_{D^{k}}+\big{(}r,\frac{\partial u}{\partial x}\big{)}_{D^{k}}=[ur]_{{}_{x_{k-\frac{1}{2}}}}^{x_{k+\frac{1}{2}}}, we get

[TABLE]

with entropy fluxes

[TABLE]

Employing Young’s inequality and Lemma 2.5, we obtain

[TABLE]

Recalling Lemma 2.3, provided $c_{i},\,\,i=1,2,3,4$ are sufficiently small such that $c_{i}\leq 1$, we obtain that

[TABLE]

we notice that, with the definition (3.9) of the numerical fluxes and with simple algebraic manipulations and summing over all elements (4.4), we easily obtain

[TABLE]

This implies that

[TABLE]

Hence

[TABLE]

Employing Gronwall’s inequality, we obtain $\|u_{h}(x,T)\|_{\Omega_{h}}\leq c\|u_{0}(x)\|_{\Omega_{h}}$.$\quad\Box$

4.2 Error estimates

We consider the linear fractional Schrödinger equation

[TABLE]

It is easy to verify that the exact solution of the above (4.8) satisfies

[TABLE]

Subtracting (4.9), from the linear fractional Schrödinger equation (3.13), we have the following error equation

[TABLE]

For the error estimate, we define special projections, $\mathcal{P}^{-}$ and $\mathcal{P}^{+}$ into $V_{h}^{k}$. For all the elements, $D^{k}$, $k=1,2,...,K$ are defined to satisfy

[TABLE]

Denoting

[TABLE]

For the special projections mentioned above, we have, by the standard approximation theory [33], that

[TABLE]

where here and below $C$ is a positive constant (which may have a different value in each occurrence) depending solely on u and its derivatives but not of $h$.

Lemma 4.1

[TABLE]

where

[TABLE]

Proof. From the Galerkin orthogonality (4.10), we get

[TABLE]

We take the test functions

[TABLE]

we obtain

[TABLE]

Summing over $k$, simplify by integration by parts and (3.9). This completes the proof. $\quad\Box$

Theorem 4.2

Let $u$ be the exact solution of the problem (4.8), and let $u_{h}$ be the numerical solution of the semi-discrete LDG scheme (3.13). Then for small enough $h$, we have the following error estimates:

[TABLE]

where the constant $C$ is dependent upon $T$ and some norms of the solutions.

Proof. Integrating both sides of the above identity Lemma 4.1 with respect to $t$ over $(0,T)$, we get

[TABLE]

Next we estimate the term $\int_{0}^{T}Q_{i}dt$, $\,\,i=1,...,4$. So we employ Young’s inequality (4.15) and the approximation results (4.13), we obtain

[TABLE]

Using the definition of the numerical traces, (3.9), and the definitions of the projections $\mathcal{P}^{+},\mathcal{P}^{-}$ (4.11), we get

[TABLE]

So

[TABLE]

From the approximation results (4.13) and Young’s inequality, we obtain

[TABLE]

Combining (6.19), (4.23) and (4.20), we obtain

[TABLE]

Recalling Lemmas 2.3, we obtain

[TABLE]

provided $c_{i},\,\,i=1,2,3,4$ are sufficiently small such that $c_{i}\leq 1$, we obtain

[TABLE]

Employing Gronwall’s lemma, we can get (4.19). $\quad\Box$

5 LDG method for strongly nonlinear coupled

fractional Schrödinger equations

In this section, we present and analyze the LDG method for the strongly coupled nonlinear fractional Schrödinger equations

[TABLE]

To define the local discontinuous Galerkin method, we rewrite (5.1) as a first-order system:

[TABLE]

We decompose the complex functions $u(x,t)$ and $v(x,t)$ into their real and imaginary parts. Setting $u_{1}(x,t)=p(x,t)+iq(x,t)$ and $u_{2}(x,t)=\upsilon(x,t)+i\theta(x,t)$ in system (5.1), we can obtain the following coupled system

[TABLE]

We define local discontinuous Galerkin scheme as follows: find $p_{h},q_{h},e_{1},r_{h},s_{h},l_{1},w_{h},z_{h}$,

$\upsilon_{h},\theta_{h},e_{2},\rho_{h},\varpi_{h},l_{2},\xi_{h}$,$\varrho_{h}\in V_{k}^{N}$, such that for all test functions $\vartheta_{1},\beta_{1},\phi,\varphi,\chi,\beta_{2},\psi$, $\zeta,\gamma,\beta_{3},\delta,\varsigma,o,\beta_{4},\omega,\kappa\in V_{k}^{N}$,

[TABLE]

[TABLE]

Applying integration by parts to (5.4), and replacing the fluxes at the interfaces by the corresponding numerical fluxes, we obtain

[TABLE]

The numerical traces $(p,q,s,z,\upsilon,\theta,\varpi,\varrho)$ are defined on interelement faces as the alternating fluxes

[TABLE]

6 Stability and error estimates

In the following we discuss stability and accuracy of the proposed scheme, for the nonlinear fractional coupled Schrödinger problem.

6.1 Stability analysis

In order to carry out the analysis of the LDG scheme,

Theorem 6.3

*($L^{2}$ stability). The semidiscrete scheme (5.5) is stable, and

$\|u_{h}(x,T)\|_{\Omega_{h}}+\|v_{h}(x,T)\|_{\Omega_{h}}\leq c(\|u_{0}(x)\|_{\Omega_{h}}+\|v_{0}(x)\|_{\Omega_{h}})$ for any $T>0$. *

Proof. Set $(\vartheta_{1},\beta_{1},\phi,\varphi,\chi,\beta_{2},\psi,\zeta,\gamma,\beta_{3},\delta,\varsigma,\beta_{4},o,\omega,\kappa)=(p_{h},T_{h}-r_{h},p_{h},-z_{h},q_{h},H_{h}-w_{h},-q_{h},s_{h},\upsilon_{h},L_{h}-\rho_{h},\upsilon_{h},-\varrho_{h},\theta_{h},E_{h}-\xi_{h},-\theta_{h},\varpi_{h})$ in (3.13), and consider the integration by parts formula \big{(}u,\frac{\partial r}{\partial x}\big{)}_{D^{k}}+\big{(}r,\frac{\partial u}{\partial x}\big{)}_{D^{k}}=[ur]_{{}_{x_{k-\frac{1}{2}}}}^{x_{k+\frac{1}{2}}}, we get

[TABLE]

Summing over all elements (6.1), employing Young’s inequality and using the definition of the numerical traces, (5.6), we obtain

[TABLE]

Recalling Lemma 2.3 and provided $c_{i},\,\,i=1,2,...,8$ are sufficiently small such that $c_{i}\leq 1$, we obtain that

[TABLE]

Hence

[TABLE]

Employing Gronwall’s inequality, we obtain

[TABLE]

6.2 Error estimates

We consider the linear fractional coupled Schrödinger system

[TABLE]

It is easy to verify that the error equations of the above (6.6) satisfies

[TABLE]

Theorem 6.4

Let $u$ and $v$ be the exact solutions of the linear coupled fractional Schrödinger equations (6.6), and let $u_{h}$ and $v_{h}$ be the numerical solutions of the semi-discrete LDG scheme (5.5). Then for small enough $h$, we have the following error estimates:

[TABLE]

where the constant $C$ is dependent upon $T$ and some norms of the solutions.

Proof. We donate

[TABLE]

From the Galerkin orthogonality (6.7), we get

[TABLE]

We take the test functions

[TABLE]

we obtain

[TABLE]

Summing over $k$, simplify by integration by parts and (5.6), we get

[TABLE]

Now, we estimate $T_{i}$ term by term.

[TABLE]

Employing Young’s inequality, we obtain

[TABLE]

and

[TABLE]

Employing Young’s inequality and Lemma 2.5, we obtain

[TABLE]

and

[TABLE]

and

[TABLE]

Using the definition of the numerical traces, (5.6), and the definitions of the projections $\mathcal{P}^{+},\mathcal{P}^{-}$ (4.11), we get

[TABLE]

Combining (6.15), (6.17), (6.20) and (6.13), we obtain

[TABLE]

Recalling Lemma 2.3, we get

[TABLE]

provided $c_{i},\,\,i=1,2,...,8$ are sufficiently small such that $c_{i}\leq 1$, we obtain

[TABLE]

An integration in $t$ plus the standard approximation theory then gives the desired error estimates.

7 Numerical examples

In this section we will present several numerical examples to illustrate the previous theoretical results. Before that, we adopt the nodal discontinuous Galerkin methods for the full spatial discretization using a high-order nodal basis set of orthonormal Lagrange-Legendre polynomials of arbitrary order in space on each element of computational domain as a more suitable and computationally stable approach As shown by Aboelenen and El-Hawary [28]. We use the high-order Runge-Kutta time discretizations [34], when the polynomials are of degree $N$, a higher-order accurate Runge-Kutta (RK) method must be used in order to guarantee that the scheme is stable. In this paper we use a fourth-order non-Total variation diminishing (TVD) Runge-Kutta scheme [35]. Numerical experiments demonstrate its numerical stability

[TABLE]

where $\mathbf{u}_{h}$ is the vector of unknowns, we can use the standard fourth-order four stage explicit RK method (ERK)

[TABLE]

to advance from $\mathbf{u}_{h}^{n}$ to $\mathbf{u}_{h}^{n+1}$, separated by the time step, $\Delta t$. In our examples, the condition $\Delta t\leq C\Delta x^{\alpha}_{min}\,\,\,(0<C<1)$ is used to ensure stability.

Example 7.1

As the first example, we consider the linear fractional Schrödinger equation

[TABLE]

with the initial condition $u_{0}(x)=x^{6}$ and the corresponding forcing term $g(x,t)$ is of the form

[TABLE]

to obtain an exact solution $u(x,t)=e^{-it}x^{6}$ with $\nu=1.2,\,\lambda=\frac{\Gamma(8-\nu)}{2\Gamma(8)}$. The errors and order of convergence are listed in Table 1, confirming optimal $O(h^{N+1})$ order of convergence across.

Example 7.2

Consider the following nonlinear fractional Schrödinger equation

[TABLE]

with the initial condition $u_{0}(x)=x^{7}$ and the corresponding forcing term $g(x,t)$ is of the form

[TABLE]

The exact solution $u(x,t)=e^{-it}x^{7}$ with $\nu=1.1,\,\lambda=\frac{\Gamma(8-\nu)}{\Gamma(8)}$. The errors and order of convergence are listed in Table 2, confirming optimal $O(h^{N+1})$ order of convergence across.

Example 7.3

We consider the nonlinear fractional Schrödinger equation

[TABLE]

with the initial condition $u_{0}(x)=(x^{2}-1)^{6}$ and the corresponding forcing term $g(x,t)$ is of the form

[TABLE]

to obtain an exact solution $u(x,t)=e^{-it}(x^{2}-1)^{6}$ with $\nu=1.5,\,\lambda=\frac{0.2\Gamma(13-\nu)}{\Gamma(13)}$. We consider cases with $N=2,3$ and $K=20,30,40,50$. The numerical orders of convergence are shown in Figure 1, showing an $O(h^{N+1})$ convergence rate for all orders.

Example 7.4

We consider the nonlinear fractional Schrödinger equation (1.1) with initial condition,

[TABLE]

with parameters $\lambda_{1}=\lambda_{2}=1$ and $\,x\in[-20,20]$. We consider cases with $N=2$ and $K=80$ and solve the equation for several different values of $\alpha$. The numerical solution $u_{h}(x,t)$ for $\alpha=1.1,\,1.4,\,1.8,\,2.0$ is shown in Figure 2. We observe that the order $\alpha$ will affect the shape of the soliton case. When $\alpha$ becomes smaller, the shape of the soliton will change more quickly. This property of the fractional Schrödinger equation can be used in physics to modify the shape of wave without change of the nonlinearity and dispersion effects. The numerical solutions of the fractional equation are convergent to the solutions of the classical non-fractional equation when $\alpha$ tends to $2$.

Example 7.5

Consider the linear coupled fractional Schrödinger equations

[TABLE]

and the corresponding forcing terms $g_{1}(x,t)$ and $g_{2}(x,t)$ are of the form

[TABLE]

The exact solutions $u_{1}(x,t)=e^{-it}x^{7}$ and $u_{2}(x,t)=e^{-it}x^{7}$ with $\nu=1.1,\,\lambda_{1}=\frac{\Gamma(8-\nu)}{\Gamma(8)}$, $\lambda_{2}=\frac{\Gamma(8-\nu)}{\Gamma(8)}$. The errors and order of convergence are listed in Tables 3 and 4, confirming optimal $O(h^{N+1})$ order of convergence across.

Example 7.6

We consider the nonlinear coupled fractional Schrödinger equations

[TABLE]

and the corresponding forcing terms $g_{1}(x,t)$ and $g_{2}(x,t)$ are of the form

[TABLE]

to obtain an exact solutions $u_{1}(x,t)=e^{-it}x^{7}$ and $u_{2}(x,t)=e^{-it}x^{7}$ with $\nu=1.2,\,\lambda_{1}=\frac{\Gamma(8-\nu)}{2\Gamma(8)}$, $\lambda_{2}=\frac{\Gamma(8-\nu)}{2\Gamma(8)}$. The errors and order of convergence are listed in Tables 5 and 6, confirming optimal $O(h^{N+1})$ order of convergence across.

Example 7.7

Consider the following nonlinear coupled fractional Schrödinger equations

[TABLE]

and the corresponding forcing terms $g_{1}(x,t)$ and $g_{2}(x,t)$ are of the form

[TABLE]

The exact solutions $u_{1}(x,t)=e^{-it}(x^{2}-1)^{6}$ and $u_{2}(x,t)=e^{-it}(x^{2}-1)^{6}$ with $\nu=1.3,\,\lambda_{1}=\frac{\Gamma(13-\nu)}{2\Gamma(13)}$, $\lambda_{2}=\frac{\Gamma(13-\nu)}{2\Gamma(13)}$. We consider cases with $N=2,3$ and $\log_{10}(h)$. The numerical orders of convergence are shown in Figure 3, showing an $O(h^{N+1})$ convergence rate for all orders.

Example 7.8

We consider the following weakly coupled problem

[TABLE]

subject to the initial conditions

[TABLE]

when $\beta=1$ and $\alpha=2$, the problem collapses to the Manakov equation, and the solitary waves collide elastically see Figure 4. The exact solutions are given by

[TABLE]

where $r_{1}=1$, $r_{2}=1$, $V_{0}=0.4$, $D=10$ and $x\in[-40,40]$. The Figures 5 and 6 present the numerical solutions for different values of order $\alpha$ and $\beta$. From these figures it is obvious that the collision of solitons are inelastic. In particular, the colliding particles stick together after interaction when $\alpha=1.8$, which means that there may occur a completely inelastic collision see Figure 6.

Example 7.9

Finally, we consider the strongly coupled system as follows

[TABLE]

subject to the initial conditions

[TABLE]

where $r_{1}=r_{2}=1$, $V_{0}=0.4$, $D=10$ and $x\in[-40,40]$.

Elastic collisions: The collision of the solitary waves is elastic [36] when $\varpi_{1}=1$, $\alpha=2$ see Figure 7. We observe that the two waves emerge without any changes in their shapes and velocities after collision. Taking $\varpi_{1}=1$, we compute the numerical solutions for different values of $\alpha$, which are depicted in Figures 8 and 9. From these figures, for any $1<\alpha\leq 2$, the collision is always elastic. When $\alpha$ tends to $2$, the shape of the solitons will change more slightly and the waveforms become closer to the classical case with $\alpha=2$.

Inelastic collision: The collision is inelastic [36] when $\varpi_{1}=0.0175$ and $\alpha=2$ see Figure 10. It is clear that the shapes and directions of two waves have changed after interaction. The observation is in accordance with the known result.

The Figures 11 and 12 present the numerical solutions for different values of order $\alpha$ for fixed $\varpi_{1}=0.0175$. From these figures it is obvious that the collision is always inelastic.

8 Conclusions

In this work, we developed and analyzed a nodal discontinuous Galerkin method for solving the nonlinear fractional Schrödinger equation and the strongly coupled nonlinear fractional Schrödinger equations, and have proven the stability of these methods. They are discretized using high-order nodal basis set of orthonormal Lagrange-Legendre polynomials as a more suitable and computationally stable approach. Numerical experiments confirm that the optimal order of convergence is recovered. As a last two examples, the weakly coupled nonlinear fractional Schrödinger equations with initial conditions are solved for different values of $\alpha$ and results show that the collision of solitons are inelastic when $\alpha\neq 2$ and the results of the strongly nonlinear fractional Schrödinger equations are the shape of the soliton will change slightly as $\alpha$ increase, with the classical case $\varpi_{1}=1$ and $\alpha=2$ as the limit. When $\varpi_{1}=1$ and $\alpha\neq 2$, the collision is always elastic and the collision is inelastic when $\varpi_{1}=0.0175$ and $1<\alpha\leq 2$.