Chebyshev pseudospectral approximation of two dimensional fractional Schrodinger equation on a convex and rectangular domain

1.   Introduction

The fractional partial differential equation was first developed as a pure mathematical theory in the 19th century ^[36]. Time- and space- nonlinear fractional Schrodinger equation(NFSE) is the fundamental equation of physics for describing nonrelativistic quantum mechanical behavior. This equation was formulated in late $1925$, and published in $1926$, by the Austrian physicist Erwin Schrodinger. In past years, the time- and space- NFSE has attracted application of various fields such as, electromagnetic waves, quantitative finance, quantum evolution of complex systems and dielectric polarization ^{[6,9,28,29,35,37]}.

Let us consider nonlinear fractional Schrodinger equation

with initial condition:

and boundary conditions:

Here, $\alpha$ and $\beta$ represent the fractional order of derivatives in time and space respectively, with values $0 < \alpha\leq 1$ and $1/2 < \beta\leq 1$. The convex domain $\Omega$ is defined as, $\Omega = (x, y)\in [x_{L}, x_{R}]\times [y_{L}, y_{R}]$, $\kappa$ and $\lambda$ are constants, $G(x, y, t)$ is a potential function, $F(x, y, t)$ is complex source term. The function $U(x, y, t)$ is assumed to be a complex wave function of time and space.

In past years, many authors have proposed various numerical methods for the solutions of time- and space- NFSE. These numerical methods are very important for understanding the physical behavior of the equations. Dong ^[11] has proposed scattering problems for two-dimensional NFSE and further, gives the mathematical expression of the Green's functions. Fan and Qi ^[14] have introduced the Galerkin finite element method for the solutions of NFSE. Some other authors ^{[15,18,22,23,25,26,42]} have used the same methods to studied numerical solutions and error analysis of two-dimensional NFSE. Zhang et al. ^[43] have proposed Galerkin-Legendre spectral schemes for the numerical solution of space NFSE. Li et al. ^[17] have discussed the numerical and stability analysis of multi-dimensional time-NFSE using the L1-Galerkin finite element method. Zhao et al. ^[44] have proposed a fourth-order compact ADI method for the numerical solutions and convergence analysis of two-dimensional space NFSE. In this method, the authors have found fourth-order of accuracy. Li et al. ^[21] have proposed a fast linearized conservative finite element method for the numerical solution of strongly coupled NFSEs and also prove that the scheme preserved both the mass and energy. Li et al. ^[24] have used Galerkin finite element method for the unconditional error analysis of NFSE. Cheng et al. ^[8] has proposed novel compact alternating direction implicit scheme for the numerical analysis of two-dimensional Riesz space fractional nonlinear reaction-diffusion equations. Bhrawy and Zaky ^[4] have studied a natural generalization of the NFSE, which is changing the variable-order NFSE to fractional quantum phenomena. Chen et al. ^[7] have introduced symplectic and multi-symplectic methods for the numerical solution of NSE. Further equation $(1.1)$ was solved by several other numerical methods, such as riccati expansion method ^[1], Legendre-Galerkin spectral method ^[19], finite element method ^[13,21], discontinuous Galerkin finite element method ^[39], momentum representation method ^[12], compact boundary value method ^[34], fractional mapping method ^[40] and finite difference method ^[20,41].

Pseudospectral methods have been developed for numerical simulation of related differential equations in many fields because of its high accuracy, especially sufficiently smooth problems. Lanczos showed the power of the Fourier series and Chebyshev polynomials in a number of problems where they had not been used before. Later in the $1970s$, Orszag introduced spectral methods again, alongside Kreiss, Oliger, and others, for the purpose of solving partial differential equations in fluid mechanics ^[38]. Many authors have used spectral method to approximate the solution of such equations ^{[2,27,30,31,32,33]}. In this paper, the authors propose a highly accurate pseudospectral method to approximate the NFSE. For the proposed method, we derive the solution of a nonlinear partial differential equation as a sum of basis functions in both space and time directions. The spectral coefficients of the sum are chosen to satisfy the solution of the nonlinear partial differential equation. The fractional derivatives matrix is considered in the modified Riemann-Liouville formula. It is a new approach to study the nonlinear partial differential equations in time and space and the method is equally important in time as it is in the space direction. The great advantage of time-space pseudospectral methods resides in their high-accuracy for a relatively small number of grid points as compared to standard time-stepping methods.

The structure of the paper is organized as follows. In section $2$, we describe some basic definitions and notations. Discretizing and description of the methods are presented in section $3$. In section $4$, we derive the error analysis for two-dimensional time- and space-NFSE. In section $5$, we present numerical solutions and errors by the proposed scheme. In the last section, the conclusion of our work is presented.

2.   Preliminary

In this section, we present definition of modified Riemann-Liouville formula of the fractional calculus theory and Sobolev space, which will be useful throughout this paper.

Definition 1: The partial fractional derivatives of order $\alpha$ of a function $\phi_{M}(z)$, with respect to variable $z$, in the modified Riemann-Liouville formula, is defined as ^[10,16]

Here, $z$ is dummy variable which represents $x, y$ and $t$ as and when required. $\Gamma$ is the gamma function and $n = \lceil \alpha\rceil +1$ with $\lceil \alpha\rceil$ denoting the upper integer part of $\alpha$. Some basic properties of modified Riemann-Liouville fractional derivatives are as follows:

and

Construction of Chebyshev fractional differential matrix is given as

Definition 2: Sobolev space: Let $p > 0$ be an integer and $1\leq q\leq \infty$. The set of all function in $L^{q}_{w^{}}[a, b]$ such that all distribution derivatives upto order $p$ are also in $L^{q}_{w^{ }} [a, b]$ is called sobolev space, which are denoted by $W_{w^{}}^{p, q}[a, b]$ and defined by

where $w^{}$ is a weight function and norm is defined by

For $q = 2$, sobolev norm denoted by $H_{w^{}}^{p}[a, b] = W_{w^{}}^{p, 2}[a, b]$ and $L^{2}_{w^{}}[a, b]$ is a weighted space defined by

Definition 3: Kronecker products of matrices, let $P = (p_{jk})_{mn}$ and $Q = (q_{jk})_{rs}$ are two matrices, where, $m, n, r, s$ are nonnegative integers. The Kronecker product $P\otimes Q$ is a $mr\times ns$ block matrix denoted by

3.   Pseudospectral method based discretization

We seek a pseudospectral approximation $I_{M}U$, as a finite linear combination of a chosen set of Langrange basis functions in both space and time-variable with pseudospectral coefficients.

where,

here, $\phi_{i}(z)$ is a Chebyshev polynomial, $z_{j}$ is a CGL points,

normalization constant

and Chebyshev Gauss-Lobatto weight quadrature

The spectral approximation given in $(3.1)$ can be expressed as a direct product

where

Now, we can define fractional spatial derivative of equation $(3.2)$ with respect to $x$ is given by

Similarly, we can define fractional derivative of Eq. $(3.2)$ with respect to $t$ is given by

where $I_{M+1}$ denotes identity matrix of size $((M+1) \times (M+1))$.

Next, we define $U(x, y, t)$ and $F(x, y, t)$ into their real and imaginary parts

where $U_{1}(x, y, t)$, $U_{2}(x, y, t)$, $F_{1}(x, y, t)$ and $F_{2}(x, y, t)$ are real functions. Using Eq. $(3.5)$ in Eq. $(1.1)$, we obtain the following coupled system of equations

Further, we consider the following transformations which are used to transform the two-dimensional space $\Omega = (x, y)\in [x_{L}, x_{R}] \times [y_{L}, y_{R}]$ and time $[0, T]$ in to $[-1, 1]$.

We obtain the two dimensional time- and space- NFSE in new space $(x, y)\in [-1, 1]^{2}$ and time interval $t \in [-1, 1]$

with initial condition:

and boundary conditions:

Further, we consider a mapping for converting the non-homogeneous initial and boundary values to homogeneous initial and boundary values ^[5].

here corner initial and boundary values satisfy,

Define new variable $V_{k}(x, y, t)$.

the above equations, can be modified with new variable and obtain the new equations residuals,

Applying the pseudospectral method in Eqs. $(3.13)$ and $(3.14)$, we get

where $B = \left(\Psi_{[0:M]}(x)\otimes \Psi_{[0:M]}(y)\otimes \Psi_{[0:M]}(t)\right) ^{T}$ Using the CGL points in Eqs. $(3.15)$ and $(3.16)$ and obtained the algebraic equations

where, $I_{(V_{1}+\Omega_{1})}$ and $I_{(V_{2}+\Omega_{2})}$ are diagonal entry matrix

Hence,

The system of nonlinear equations $(3.19)$ can be solved by using Newton Raphson method. The Jacobian matrix of this coupled nonlinear algebraic equations are as

.

4.   Theory of error estimates

In this section, we derive error estimates for two-dimensional nonlinear fractional Schrodinger equation by using time-space Chebyshev pseudospectral method.

Theorem 1: Let $U(x, y, t) \in H^{p}_{w^{}}[-1, 1]$, where $p > 0$, there exist a constant $C_{}$, such that the following estimate holds:

Proof : We approximate the solutions obtained from truncated Chebyshev expansion are defined by

where $\eta_{l}(z), \forall l\in(i, j, k)$ and $z\in(x, y, t)$ is Chebyshev polynomials, $w(z) = (\sqrt{1-z^{2}})^{-1}$ is weight function and expansion coefficients are

where $p_{l}\; \; l\in(i, j, k)$ is normalizing constant of orthogonal Chebyshev polynomials.

Next, we define the Chebyshev operator

Using Eq. $(4.3)$ in Eq. $(4.2)$, we obtain

where

using integration by parts:

again applying integration by parts, we obtain

here, $\mathcal{L}_{1} = \frac{\partial}{\partial x} \left(-\sqrt{1-x^{2}} \frac{\partial }{\partial x} \; \right)$ is Chebyshev operator,

putting the value of $\chi_{1}$ in eq(4.4), we get

Similarly integrate for space $y$ and time $t$ and obtain

In this case the order of $\mathcal{L} U$ is $6$ as $\mathcal{L}_{1}, \mathcal{L}_{2}$ and $\mathcal{L}_{3}$ are second order operator.

Repeatation of procedure for higher order concludes the results $i.e.$,

Then the Cauchy-Schwartz inequality implies

where $C$ is generic constant. Therefore, from Sobolev norm

By induction,

From Eq. $4.7$ and Eq. $4.8$, we get

Further, let us consider the discrete approximation

where $\Pi_{ijk}$ is discrete expansion coefficient and defined by

Under the assumption for sufficient smoothness functions, the aliasing error is

where $\left[., .\right]_{w^{}}$ represents the discrete weighted inner product.

We know that,

Using Eq. (4.10), Eq. (4.12) and Eq. (4.13) and we get

where the aliasing error

to simplify the expression, first interchange the summations

Since $\eta_{i}(x), \eta_{j}(y)$ and $\eta_{k}(t)$ are orthogonal, therefore, we observe that the value of Eq. $(4.15)$ is zero due to the range of summations.

Subtract the Eq. $(4.13)$ and Eq. $(3.13)$, we get

taking $L_{w^{}}^{2}[-1, 1]$ norm both side and we get

Finally putting the value of discrete inner product and we get

Combining the Eq. $(4.9)$ and Eq. $(4.17)$and we get

Take $p = 6m$ in Eq. (4.18)

Finally, putting in Eq. (4.2), we obtain

The theorem is complete.

Further, we discuss the bounds of the given Schrodinger operator. We consider the error function $E(x, y, t) = U(x, y, t)-\mathcal{I}_{M}U(x, y, t)$, where, $\mathcal{I}_{M}U(x, y, t)$ is pseudospectral approximation and $U(x, y, t)$ is exact solution. Moreover, the pseudospectral approximation satisfies the given problem, $i.e.$,

Adding both equations, we get

where $R(x, y, t) = R_{1}(x, y, t)+iR_{2}(x, y, t)$ is the residual function.

Further, we can define

where $E = E_{1}+iE_{2}$.

From theorem $1$, we get

Further, Eq. $(4.20)$ can be written in derivative form.

Using sobolev norms

Using Eq. $(4.25)$, we obtain

taking $L_{w^{}}^{2}[-1, 1]$ norm both side in Eq. $(4.22)$ and we get

Using Eq. $(4.26)$ in Eq. $(4.27)$, we get

Thus,

5.   Numerical results

In this section, we present numerical results for two-dimensional time- and space-NFSE using the Chebyshev pseudospectral method. We give three examples in this section. To demonstrate the errors in pseudospectral approximation, we consider the errors in the $L_{2}$ norms, defined by

where $\arrowvert U\arrowvert^{k}$ and $\arrowvert U \arrowvert$ are the modulus of pseudospectral approximation and modulus of analytical solution, respectively.

5.1.   Example 1

Let us consider the time- and space- NFSE on an ellipse convex region $\Omega = \left\lbrace (x, y) | \frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}} \leq 1 \right\rbrace$

with $\lambda = \kappa = 1$, $G(x, y, t) = t\left(\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}-1 \right)$.

The source term is

With $s_{1}(a) = \frac{\Gamma 5}{a^{4}\Gamma(5-2\beta)}$, $s_{2}(a) = \frac{\Gamma 4}{a^{4}\Gamma(4-2\beta)}$ and $s_{3}(a) = \frac{\Gamma 3}{a^{4}\Gamma(3-2\beta)}$, the function $h(x, a), h(y, b)$ can be given as

Initial condition

boundary conditions

and the exact solution is

where, $x_{L} = -\frac{a}{b}\sqrt{b^{2}-y^{2}}, x_{R} = \frac{a}{b}\sqrt{b^{2}-y^{2}}, y_{L} = -\frac{b}{a}\sqrt{a^{2}-x^{2}}$, $y_{R} = \frac{b}{a}\sqrt{a^{2}-x^{2}}$.

In this example, numerical solutions of the proposed method have obtained over fractional-order derivative $0 < \alpha\leq 1$ and $1/2 < \beta\leq 1$. The tabulated results of the proposed method with different fractional-order derivatives $\alpha$ and $\beta$ are presented in Table 1. The numerical results of the proposed method achieved better accuracy as the number of grid points in both space and time axis is increased. Figure 1 has shown that surface plot of the proposed method at time $T = 1.0$ with $a = 1/2, b = 1/4, \alpha = 0.8$ and $\beta = 0.85$. Contour plot also has shown the physical behavior of the equation. Moreover, Figure 2 has also shown that surface plot of the proposed method at time $T = 1.0$ with $a = 1, b = 1, \alpha = 0.8$ and $\beta = 0.85$. The numerical results are more accurate and consistent with our theoretical results.

5.2.   Example 2

In this example, let us consider the time- and space- NFSE on a rectangular region $\Omega = [0, 1]\times [0, 1]$

with $\lambda = \kappa = 1$, $G(x, y, t) = txy(1-x)(1-y)$.

The source term is

where,

Initial condition

boundary conditions

and the exact solution is

For computation purpose we have chosen space interval $\Omega = [0, 1]^{2}$ with time $T = 1$. Numerical solutions of the proposed method have been computed with different fractional-order derivative $0 < \alpha\leq 1$ and $1/2 < \beta\leq 1$. Figure 3 has shown the numerical solution of the proposed method and the contour plot has shown the physical behaviour of the proposed method. In Table 2, the tabulated results of the proposed method with different fractional-order derivatives $\alpha$ and $\beta$ are presented. The numerical results of the proposed method achieved better accuracy as the number of grid points in both space and time directions are increased. Moreover, it demonstrated that the numerical method is more efficient.

5.3.   Example 3

In this example, let us consider the time- NFSE with non-homogeneous boundary value

where $h(x, y, t) = \frac{t^{3-\alpha}(6-i(t^{6}+1)\Gamma (4-\alpha)t^{\alpha})(e^{-i(x+y)})}{\Gamma (4-\alpha)}$,

with initial condition:

and boundary conditions:

The exact solution of this problem is

In this example, error norms for two dimensional nonlinear time fractional Schrodinger equation has been calculated with different values of grids point and different fractional order derivative $0 < \alpha \leq 1$. In Table 3, it can be seen that the accuracy of the numerical results are increased along with the number of grid points and also achieved good order of accuracy. Further, proposed methods have avheived better accuracy as compared to ^[3] at $\alpha = 0.50$ with different grid points. Moreover, proposed method has obtained $14^{th}$ order of accuracy.

6.   Conclusion

In this work, we have discussed a highly accurate fully discrete time-space Chebyshev pseudospectral method for the two-dimensional time- and space-NFSE, defined on a convex and rectangular domain. The new fractional derivative matrix has been established using a modified Riemann-Liouville derivative formula at CGL points for a different order of fractional derivatives. We presented the error analysis without any dependency on time and space step restrictions of the schemes. The proposed method supports the theoretical results. To demonstrate the performance, the method has been employed on three different model problems on a convex and rectangular domain and obtained a good order of accuracy. Reported numerical results are highly accurate which shows the efficiency of the proposed method.

Acknowledgements

The first author thankfully acknowledges to the Ministry of Human Resource Development, India, for providing financial support for this research.

Conflict of interest

The authors declare no conflict of interest.