Numerical Simulation of Time-Fractional Nonlinear Schrödinger Equation and Time-Fractional Schrödinger–Poisson System

1. Introduction

Many natural phenomena show complex and nonlinear behavior that cannot be described accurately by classical integer-order models. Fractional calculus (FC) provides a flexible mathematical tool for studying such processes and gives a more realistic description of their behavior. For this reason, fractional models have attracted considerable interest in science and engineering during recent decades. In particular, fractional partial differential equations (FPDEs) are useful for describing systems with memory and nonlocal effects [1,2,3]. These equations have been applied to many problems in elasticity, blood flow, solid mechanics, optical fibers, signal processing, thermal radiation, fluid dynamics, medical science, and diffusion processes. Because of their nonlinear structure, exact solutions of these equations are generally difficult or impossible to obtain. Analytical solutions are available only in a limited number of special cases. Therefore, developing accurate and efficient numerical methods for NPDEs is of considerable practical importance [4]. In this paper, we focus on two important models: the nonlinear time-fractional Schrödinger equation and the time-fractional Schrödinger–Poisson system.

1.1. The Time-Fractional Nonlinear Schrödinger Equation

First, we consider the following time-fractional nonlinear Schrödinger equation:

$${\mathrm{i}}_0^{\mathrm{C}}{\mathbf{D}}_t^{\mathsf{\beta }}u(\mathbf{x},t)+\alpha \Delta u(\mathbf{x},t)+V\left(\mathbf{x}\right)u(\mathbf{x},t)+\lambda {\left|u(\mathbf{x},t)\right|}^{2}u(\mathbf{x},t)=f(\mathbf{x},t),\mathbf{x}\in \Omega ,0<t\le T,0<\mathsf{\beta }<1.$$

(1)

Equation (1) is supplemented with the boundary condition (BC)

$$u(\mathbf{x},t)=h(\mathbf{x},t),\mathbf{x}\in \partial \Omega ,0\le t\le T,$$

(2)

and the initial condition (IC)

$$u(\mathbf{x},0)=\psi \left(\mathbf{x}\right),\mathbf{x}\in \overline{\Omega },$$

(3)

where $\Omega \subset {\mathbb{R}}^2$ be a bounded computational domain, $\mathbf{x}=(x,y)\in \Omega$, and $0<t\le T$. $u:\Omega \times [0,T]\to \mathbb{C}$ is the complex-valued wave function, $\mathrm{i}=\sqrt{-1}$, $\alpha$ is the dispersion coefficient, $V\left(\mathbf{x}\right)$ is the external potential, $\lambda$ is the nonlinear interaction coefficient, and $f(\mathbf{x},t)$ denotes the source function. The Caputo fractional derivative of order $\mathsf{\beta }$ is defined as

$${}_0{}^{\mathrm{C}}\mathbf{D}_t^{\mathsf{\beta }}\psi (\mathbf{x},t)={\displaystyle \frac{1}{\Gamma (1-\mathsf{\beta })}}{\int }_0^t{\displaystyle \frac{\partial \psi (\mathbf{x},s)/\partial s}{{(t-s)}^{\mathsf{\beta }}}}ds,0<\mathsf{\beta }<1.$$

(4)

The fractional nonlinear Schrödinger equation is a nonlocal dispersive model introduced by Laskin through a generalization of Feynman’s path-integral formulation to Lévy-type trajectories [5,6]. It was later derived as the continuum limit of discrete lattice systems with long-range interactions [7]. Its nonlocal structure makes it possible to describe physical phenomena that are not captured by the classical Schrödinger equation [8,9,10,11]. Nevertheless, this nonlocality also creates considerable analytical and computational difficulties. In a different direction, time-fractional Schrödinger equations replace the classical first-order time derivative with a fractional derivative to model memory and hereditary effects.

Numerous numerical techniques have been proposed for the numerical approximation of the nonlinear time-fractional Schrödinger equation. Mohebbi et al. [12] proposed a meshless collocation method based on radial basis function (RBF). Ding and Wong [13] employed a quintic non-polynomial spline approximation to construct a high-order numerical scheme for the same class of equations. Zhang et al. [14] developed an efficient approach based on a fast approximation of the Caputo derivative for problems posed on unbounded domains. Wang et al. [15] presented two linear second-order schemes for the multidimensional nonlinear model and analyzed their stability and convergence properties. Hu et al. [16] proposed a two-grid finite element (FE) method on graded temporal meshes to handle possible initial singularities. Wang et al. [17] developed a time two-grid fitted scheme for nonsmooth solutions, reducing the computational cost while maintaining high accuracy. Huang et al. [18] introduced a meshless collocation approach for the multidimensional problem, combined with L1-type approximations of the Caputo derivative. Wang et al. [19] proposed a modified L1–Crank–Nicolson FE method and established unconditional convergence of the resulting scheme. Ma et al. [20] constructed an efficient numerical method for the two-dimensional nonlinear case and investigated its accuracy and convergence behavior. Qian et al. [21] combined an L1–Crank–Nicolson discretization on graded meshes with the virtual element method for the nonlinear fractional model. Wei et al. [22] developed an implicit fully discrete local discontinuous Galerkin method for the linear time-fractional Schrödinger equation and studied its stability and convergence. Zhang et al. [23] considered a two-dimensional formulation using temporal fractional discretization together with spectral Galerkin techniques. Serik et al. [24] proposed a high-order finite difference method combining an accurate temporal approximation with a compact spatial discretization.

1.2. The Time-Fractional Schrödinger–Poisson System

Second, we study the following nonlinear time-fractional Schrödinger–Poisson system:

$$\left\{\begin{array}{c}{\displaystyle {\mathrm{i}}_0^{\mathrm{C}}{\mathbf{D}}_t^{\mathsf{\beta }}\psi (\mathbf{x},t)+a\Delta \psi (\mathbf{x},t)+\lambda {\left|\psi (\mathbf{x},t)\right|}^2\psi (\mathbf{x},t)=\alpha \Phi (\mathbf{x},t)\psi (\mathbf{x},t)+f(\mathbf{x},t),}\hfill \\ {\displaystyle -\Delta \Phi (\mathbf{x},t)=c_D{\left|\psi (\mathbf{x},t)\right|}^2+g(\mathbf{x},t),}\hfill \end{array}\right.\mathbf{x}\in \Omega ,0<t\le T.$$

(5)

The system is supplemented with the Dirichlet BCs

$$\psi (\mathbf{x},t)=h_{\psi }(\mathbf{x},t),\Phi (\mathbf{x},t)=h_{\Phi }(\mathbf{x},t),\mathbf{x}\in \partial \Omega ,0<t\le T,$$

(6)

and the IC

$$\psi (\mathbf{x},0)={\psi }_0\left(\mathbf{x}\right),\mathbf{x}\in \overline{\Omega }.$$

(7)

where $0<\mathsf{\beta }<1,$$\psi :\Omega \times [0,T]\to \mathbb{C}$ denotes the complex-valued wave function, $\Phi :\Omega \times [0,T]\to \mathbb{R}$ represents the electrostatic potential, $\mathrm{i}=\sqrt{-1}$, and $\Delta ={\nabla }^2$ denotes the Laplace operator in d spatial dimensions, where $d=1,2,3$. The coupling constant $c_D$ may, in general, take an arbitrary value. In the present study, it is set to $c_D=2\pi$ for the two-dimensional problem and $c_D=4\pi$ for the three-dimensional problem [?]. The nonlinear coefficient $\lambda$ is assumed to take one of the values $-1$, 0, or 1. In nonlinear optics, $\lambda \ne 0$ represents the strength and type of the Kerr nonlinearity. In particular, $\lambda =1$ corresponds to an attractive local interaction, whereas $\lambda =-1$ describes a repulsive local interaction. The case $\lambda =0$ represents the absence of the local cubic nonlinear interaction. Moreover, $\mathsf{\beta }$ is the fractional order, $a>0$ is the dispersion coefficient, and $\lambda$, $\alpha$, and $c_D$ are given real parameters.

The nonlinear Schrödinger equation is an important mathematical model for wave propagation in quantum mechanics, nonlinear optics, and condensed-matter physics. Its nonlinear term describes the local self-interaction of the wave field. In many physical systems, however, the interaction is generated by a field that depends on the spatial distribution of the wave function. This leads to nonlocal models such as the Schrödinger–Poisson system, in which a nonlinear Schrödinger equation is coupled with an elliptic Poisson equation through the wave density ${\left|\psi \right|}^2$ [25,26]. Related self-consistent field models can also be obtained from mean-field approximations of many-particle quantum systems with long-range Coulomb or gravitational interactions. The Poisson equation determines a self-consistent potential generated by the wave density, and this potential, in turn, affects the evolution of the wave function. Depending on the underlying physical problem, related forms of this model are known as the Schrödinger–Newton equation or the Gross–Pitaevskii–Newton equation. Schrödinger–Poisson systems arise in the study of self-gravitating quantum particles, Bose–Einstein condensates, semiconductor devices, plasma dynamics, nonlinear optics, and scalar-field dark matter [25,27,28,29,30,31]. In the classical Schrödinger–Poisson system, the time evolution of the wave function is governed by a first-order derivative. Such a model may not fully describe processes involving memory, anomalous relaxation, or hereditary behavior. A time-fractional extension can be obtained by replacing the classical time derivative with a fractional derivative. In particular, the Caputo derivative accounts for the previous history of the solution and permits the use of initial conditions in the usual classical form. The resulting time-fractional Schrödinger–Poisson system therefore combines temporal memory, cubic nonlinearity, and a nonlocal self-consistent potential. Gu et al. [32] investigated the existence and multiplicity of solutions for a fractional Schrödinger–Poisson system involving a spatial fractional operator rather than a time-fractional derivative. Zhang et al. [33] developed a time-splitting Fourier spectral method for a space-fractional Schrödinger–Poisson model, providing an efficient spectral treatment of the spatial fractional operator.

In this paper, we employ an efficient localized meshless method for solving nonlinear time-fractional Schrödinger equations and Schrödinger–Poisson systems. The Caputo time-fractional derivative is approximated by the L1 formula, whose consistency order is $\mathcal{O}\left(\delta t^{2-\mathsf{\beta }}\right)$ for $0<\mathsf{\beta }<1$. The spatial derivatives are discretized using the localized radial basis function Hermite finite-difference (LRBF-HFD) method. The Laplacian operator is approximated over small overlapping neighborhoods, which avoids the large and ill-conditioned systems associated with global RBF methods. Finally, several numerical examples confirm the accuracy and reliability of the proposed method.

The outline of this paper is as follows. Section 2 develops the time-semi discrete schemes using the L1 formula for the nonlinear time-fractional Schrödinger equation and the Schrödinger–Poisson system. Section 3 presents the LRBF-HFD method for the spatial discretization. Section 4 reports the numerical results. Finally, Section 5 presents the main conclusions and final remarks.

2. Time Discretization

This section uses the L1 formula to approximate the time-fractional derivative. The time interval $[0,T]$ is divided into N equal subintervals, with temporal nodes $t_n=n\delta t$, where $\delta t=T/N_t$ and $n=0,1,\dots ,N_t$. For convenience, the notation $u^n\left(\mathbf{x}\right)=u(\mathbf{x},t_n)$ is used. Applying the L1 formula at $t=t_n$ leads to a semi-discrete formulation of the problem.

$$u^n\left(\mathbf{x}\right)=u(\mathbf{x},t_n),f^n\left(\mathbf{x}\right)=f(\mathbf{x},t_n),h^n\left(\mathbf{x}\right)=h(\mathbf{x},t_n).$$

We first recall the following result for the L1 approximation of the Caputo time-fractional derivative.

Lemma 1. (See [34].) Suppose that $0<\mathsf{\beta }<1$ and $g\in C^2[0,t_n]$. Then,

$$\begin{array}{cc}\hfill |& {\displaystyle \frac{1}{\Gamma (1-\mathsf{\beta })}}{\int }_0^{t_n}{\displaystyle \frac{g^{\prime }\left(s\right)}{{(t_n-s)}^{\mathsf{\beta }}}}ds\hfill \\ \hfill & -{\displaystyle \frac{\delta t^{-\mathsf{\beta }}}{\Gamma (2-\mathsf{\beta })}}\left[b_{0}g\left(t_n\right)-{\displaystyle \sum _{k=1}^{n-1}}\left(b_{n-k-1}-b_{n-k}\right)g\left(t_k\right)-b_{n-1}g\left(t_0\right)\right]|\hfill \\ \hfill & \le {\displaystyle \frac{1}{\Gamma (2-\mathsf{\beta })}}\left[{\displaystyle \frac{1-\mathsf{\beta }}{12}}+{\displaystyle \frac{2^{2-\mathsf{\beta }}}{2-\mathsf{\beta }}}-\left(1+2^{-\mathsf{\beta }}\right)\right]\underset{0\le s\le t_n}{max}\left|g^{\prime \prime }\left(s\right)\right|\delta t^{2-\mathsf{\beta }},\hfill \end{array}$$

(8)

where $b_k={(k+1)}^{1-\mathsf{\beta }}-k^{1-\mathsf{\beta }},k=0,1,2,\dots .$

2.1. Temporal Discretization of the Time-Fractional Schrödinger Equation

For $0<\mathsf{\beta }<1$, applying Lemma 1 to $\psi (\mathbf{x},t)$ at $t=t_n$, the Caputo fractional derivative can be approximated by the L1 formula as follows:

$$\begin{array}{cc}\hfill {}_0{}^{\mathrm{C}}\mathbf{D}_t^{\mathsf{\beta }}{u}^n\left(\mathbf{x}\right)& ={\displaystyle \frac{\delta t^{-\mathsf{\beta }}}{\Gamma (2-\mathsf{\beta })}}\left[b_0{u}^n\left(\mathbf{x}\right)-{\displaystyle \sum _{k=1}^{n-1}}\left(b_{n-k-1}-b_{n-k}\right)u^k\left(\mathbf{x}\right)-b_{n-1}{u}^0\left(\mathbf{x}\right)\right]+\mathcal{O}\left(\delta t^{2-\mathsf{\beta }}\right).\hfill \end{array}$$

(9)

Since $b_0=1$, introducing

$${\mu }_{\mathsf{\beta }}={\displaystyle \frac{\delta t^{-\mathsf{\beta }}}{\Gamma (2-\mathsf{\beta })}}$$

(10)

allows Eq. (9) to be written as

$$\begin{array}{cc}\hfill {}_0{}^{\mathrm{C}}\mathbf{D}_t^{\mathsf{\beta }}{u}^n\left(\mathbf{x}\right)& ={\mu }_{\mathsf{\beta }}\left[u^n\left(\mathbf{x}\right)-{\displaystyle \sum _{k=1}^{n-1}}\left(b_{n-k-1}-b_{n-k}\right)u^k\left(\mathbf{x}\right)-b_{n-1}{u}^0\left(\mathbf{x}\right)\right]+\mathcal{O}\left(\delta t^{2-\mathsf{\beta }}\right).\hfill \end{array}$$

(11)

To obtain a linear problem at each time level, the nonlinear term is treated explicitly. Assuming that u is continuously differentiable with respect to time, we have

$$u^n\left(\mathbf{x}\right)=u^{n-1}\left(\mathbf{x}\right)+\mathcal{O}\left(\delta t\right).$$

(12)

Consequently,

$${\left|u^n\left(\mathbf{x}\right)\right|}^2{u}^n\left(\mathbf{x}\right)={\left|u^{n-1}\left(\mathbf{x}\right)\right|}^2{u}^{n-1}\left(\mathbf{x}\right)+\mathcal{O}\left(\delta t\right).$$

(13)

Substituting Eqs. (11) and (13) into Eq. (1), and neglecting the truncation-error terms, gives the following time-semidiscrete formulation.

At the first time level, $n=1$, the history sum is empty. Using $u^0\left(\mathbf{x}\right)=\psi \left(\mathbf{x}\right)$, we obtain

$$\begin{array}{cc}\hfill \mathrm{i}{\mu }_{\mathsf{\beta }}{u}^1\left(\mathbf{x}\right)& +\alpha \Delta u^1\left(\mathbf{x}\right)+V\left(\mathbf{x}\right)u^1\left(\mathbf{x}\right)+\lambda {\left|\psi \left(\mathbf{x}\right)\right|}^2\psi \left(\mathbf{x}\right)=f^1\left(\mathbf{x}\right)+\mathrm{i}{\mu }_{\mathsf{\beta }}\psi \left(\mathbf{x}\right),\mathbf{x}\in \Omega .\hfill \end{array}$$

(14)

For $2\le n\le N$, the time-semidiscrete problem becomes

$$\begin{array}{cc}\hfill \mathrm{i}{\mu }_{\mathsf{\beta }}{u}^n\left(\mathbf{x}\right)& +\alpha \Delta u^n\left(\mathbf{x}\right)+V\left(\mathbf{x}\right)u^n\left(\mathbf{x}\right)+\lambda {\left|u^{n-1}\left(\mathbf{x}\right)\right|}^2{u}^{n-1}\left(\mathbf{x}\right)\hfill \\ \hfill & =f^n\left(\mathbf{x}\right)+\mathrm{i}{\mu }_{\mathsf{\beta }}{\displaystyle \sum _{k=1}^{n-1}}\left(b_{n-k-1}-b_{n-k}\right)u^k\left(\mathbf{x}\right)+\mathrm{i}{\mu }_{\mathsf{\beta }}{b}_{n-1}\psi \left(\mathbf{x}\right),\mathbf{x}\in \Omega .\hfill \end{array}$$

(15)

At every time level, Eq. (15) is supplemented with the discrete BC

$$u^n\left(\mathbf{x}\right)=h^n\left(\mathbf{x}\right),\mathbf{x}\in \partial \Omega ,n=1,2,\dots ,N,$$

(16)

and the IC

$$u^0\left(\mathbf{x}\right)=\psi \left(\mathbf{x}\right),\mathbf{x}\in \overline{\Omega }.$$

(17)

2.2. Temporal Discretization of the Time-Fractional Schrödinger–Poisson System

In what follows, we derive the time-discrete formulation of the time-fractional Schrödinger–Poisson system. For convenience, we introduce the notation

$${\psi }^n\left(\mathbf{x}\right)=\psi (\mathbf{x},t_n),{\Phi }^n\left(\mathbf{x}\right)=\Phi (\mathbf{x},t_n),$$

together with

$$f^n\left(\mathbf{x}\right)=f(\mathbf{x},t_n),g^n\left(\mathbf{x}\right)=g(\mathbf{x},t_n).$$

By applying Lemma 1 to $\psi (\mathbf{x},t)$ at $t=t_n$, we obtain

$$\begin{array}{cc}\hfill {}_0{}^{\mathrm{C}}\mathbf{D}_t^{\mathsf{\beta }}{\psi }^n\left(\mathbf{x}\right)& ={\displaystyle \frac{\delta t^{-\mathsf{\beta }}}{\Gamma (2-\mathsf{\beta })}}\left[b_0{\psi }^n\left(\mathbf{x}\right)-{\displaystyle \sum _{k=1}^{n-1}}\left(b_{n-k-1}-b_{n-k}\right){\psi }^k\left(\mathbf{x}\right)-b_{n-1}{\psi }^0\left(\mathbf{x}\right)\right]+\mathcal{O}\left(\delta t^{2-\mathsf{\beta }}\right),\hfill \end{array}$$

(18)

where $b_k={(k+1)}^{1-\mathsf{\beta }}-k^{1-\mathsf{\beta }},k=0,1,2,\dots$. Since $b_0=1$, we define ${\displaystyle {\mu }_{\mathsf{\beta }}=\frac{\delta t^{-\mathsf{\beta }}}{\Gamma (2-\mathsf{\beta })}.}$ Consequently, Eq. (18) can be written as

$$\begin{array}{cc}\hfill {}_0{}^{\mathrm{C}}\mathbf{D}_t^{\mathsf{\beta }}{\psi }^n\left(\mathbf{x}\right)& ={\mu }_{\mathsf{\beta }}\left[{\psi }^n\left(\mathbf{x}\right)-{\displaystyle \sum _{k=1}^{n-1}}\left(b_{n-k-1}-b_{n-k}\right){\psi }^k\left(\mathbf{x}\right)-b_{n-1}{\psi }^0\left(\mathbf{x}\right)\right]+\mathcal{O}\left(\delta t^{2-\mathsf{\beta }}\right).\hfill \end{array}$$

(19)

To obtain a linear problem at each time level, the nonlinear terms are treated explicitly. Assuming sufficient temporal regularity of $\psi$, we have

$${\psi }^n\left(\mathbf{x}\right)={\psi }^{n-1}\left(\mathbf{x}\right)+\mathcal{O}\left(\delta t\right).$$

Therefore,

$${\left|{\psi }^n\left(\mathbf{x}\right)\right|}^2{\psi }^n\left(\mathbf{x}\right)={\left|{\psi }^{n-1}\left(\mathbf{x}\right)\right|}^2{\psi }^{n-1}\left(\mathbf{x}\right)+\mathcal{O}\left(\delta t\right),$$

and

$${\left|{\psi }^n\left(\mathbf{x}\right)\right|}^2={\left|{\psi }^{n-1}\left(\mathbf{x}\right)\right|}^2+\mathcal{O}\left(\delta t\right).$$

At each time level, the electrostatic potential is first obtained by solving the Poisson equation

$$-\Delta {\Phi }^n\left(\mathbf{x}\right)=c_D{\left|{\psi }^{n-1}\left(\mathbf{x}\right)\right|}^2+g^n\left(\mathbf{x}\right),\mathbf{x}\in \Omega .$$

(20)

Once ${\Phi }^n$ is known, the Schrödinger equation is solved for ${\psi }^n$. At the first time level, $n=1$, the history sum is empty. Using ${\psi }^0\left(\mathbf{x}\right)={\psi }_0\left(\mathbf{x}\right)$, we obtain

$$-\Delta {\Phi }^1\left(\mathbf{x}\right)=c_D{\left|{\psi }_0\left(\mathbf{x}\right)\right|}^2+g^1\left(\mathbf{x}\right),\mathbf{x}\in \Omega ,$$

(21)

followed by

$$\begin{array}{cc}\hfill \mathrm{i}{\mu }_{\mathsf{\beta }}{\psi }^1\left(\mathbf{x}\right)& +a\Delta {\psi }^1\left(\mathbf{x}\right)+\lambda {\left|{\psi }_0\left(\mathbf{x}\right)\right|}^2{\psi }_0\left(\mathbf{x}\right)=\alpha {\Phi }^1\left(\mathbf{x}\right){\psi }^1\left(\mathbf{x}\right)+f^1\left(\mathbf{x}\right)+\mathrm{i}{\mu }_{\mathsf{\beta }}{\psi }_0\left(\mathbf{x}\right),\mathbf{x}\in \Omega .\hfill \end{array}$$

(22)

For $2\le n\le N$, the potential ${\Phi }^n$ is first determined from

$$-\Delta {\Phi }^n\left(\mathbf{x}\right)=c_D{\left|{\psi }^{n-1}\left(\mathbf{x}\right)\right|}^2+g^n\left(\mathbf{x}\right),\mathbf{x}\in \Omega .$$

(23)

The wave function at the new time level is then computed from

$$\begin{array}{cc}\hfill \mathrm{i}{\mu }_{\mathsf{\beta }}{\psi }^n\left(\mathbf{x}\right)& +a\Delta {\psi }^n\left(\mathbf{x}\right)+\lambda {\left|{\psi }^{n-1}\left(\mathbf{x}\right)\right|}^2{\psi }^{n-1}\left(\mathbf{x}\right)\hfill \\ \hfill & =\alpha {\Phi }^n\left(\mathbf{x}\right){\psi }^n\left(\mathbf{x}\right)+f^n\left(\mathbf{x}\right)+\mathrm{i}{\mu }_{\mathsf{\beta }}{\displaystyle \sum _{k=1}^{n-1}}\left(b_{n-k-1}-b_{n-k}\right){\psi }^k\left(\mathbf{x}\right)+\mathrm{i}{\mu }_{\mathsf{\beta }}{b}_{n-1}{\psi }_0\left(\mathbf{x}\right),\mathbf{x}\in \Omega .\hfill \end{array}$$

(24)

Equivalently, collecting the unknown terms on the left-hand side gives

$$\begin{array}{cc}\hfill \left(\mathrm{i}{\mu }_{\mathsf{\beta }}-\alpha {\Phi }^n\left(\mathbf{x}\right)\right){\psi }^n\left(\mathbf{x}\right)+a\Delta {\psi }^n\left(\mathbf{x}\right)& =f^n\left(\mathbf{x}\right)-\lambda {\left|{\psi }^{n-1}\left(\mathbf{x}\right)\right|}^2{\psi }^{n-1}\left(\mathbf{x}\right)+\mathrm{i}{\mu }_{\mathsf{\beta }}{\displaystyle \sum _{k=1}^{n-1}}\left(b_{n-k-1}-b_{n-k}\right){\psi }^k\left(\mathbf{x}\right)\hfill \\ \hfill & +\mathrm{i}{\mu }_{\mathsf{\beta }}{b}_{n-1}{\psi }_0\left(\mathbf{x}\right),2\le n\le N.\hfill \end{array}$$

(25)

3. Spatial Discretization by the LRBF-HFD Method

This section first provides a concise overview of the radial basis function–based Hermite interpolation method and then uses this approach to derive the LRBF–HFD formulation [35,36,37]. Let $\Omega \subset {\mathbb{R}}^2$ denote a bounded computational domain with boundary $\partial \Omega$, and let $u\left(\mathbf{x}\right)$ be defined on $\overline{\Omega }$. The domain is discretized by a set of M scattered nodes, denoted by ${\left\{{\mathbf{x}}_i\right\}}_{i=1}^M$. For each center ${\mathbf{x}}_j$, let ${\Omega }_j=\left\{{\mathbf{x}}_1^{\left(j\right)},\dots ,{\mathbf{x}}_m^{\left(j\right)}\right\}\subset \Omega$ denote a local stencil consisting of ${\mathbf{x}}_j$ and its $m-1$ nearest neighbors. The stencil size m may be fixed for all centers or allowed to vary with j. Assume that the values of u are known at the nodes in ${\Omega }_j$, whereas the values of $\mathcal{L}u$ are prescribed at a selected subset of these nodes, that is, ${\tilde{\Omega }}_j=\left\{{\tilde{\mathbf{x}}}_1^{\left(j\right)},\dots ,{\tilde{\mathbf{x}}}_l^{\left(j\right)}\right\}\subset \Omega .$ Figure 1 shows how a local stencil is formed around the central node ${\mathbf{x}}_0$ in an irregular domain. The nearby support nodes are chosen within the distance $d_{max}$, while the other interior and boundary nodes remain outside the local region.

For notational simplicity, the superscripts introduced above will be omitted hereafter. The local RBF–Hermite interpolant of $u\left(\mathbf{x}\right)$ is then written as follows:

$$\tilde{u}\left(\mathbf{x}\right)={\displaystyle \sum _{i=1}^m}{\lambda }_i\phi \left(∥\mathbf{x}-{\mathbf{x}}_i∥\right)+{\displaystyle \sum _{i=1}^l}{\vartheta }_i{\mathcal{L}}_2\phi \left(∥\mathbf{x}-{\tilde{\mathbf{x}}}_i∥\right)+{\displaystyle \sum _{i=1}^s}{\nu }_i{p}_i\left(\mathbf{x}\right).$$

(26)

The interpolation coefficients ${\left\{{\lambda }_p\right\}}_{p=1}^m$, ${\left\{{\vartheta }_{\ell }\right\}}_{\ell =1}^l$, and ${\left\{{\nu }_j\right\}}_{j=1}^s$ are determined by imposing the following conditions:

$$\begin{array}{cccc}\hfill \tilde{u}\left({\mathbf{x}}_p\right)& =u\left({\mathbf{x}}_p\right),\hfill & \hfill & p=1,2,\dots ,m,\hfill \\ \hfill \mathcal{L}\tilde{u}\left({\tilde{\mathbf{x}}}_{\ell }\right)& =\mathcal{L}u\left({\tilde{\mathbf{x}}}_{\ell }\right),\hfill & \hfill & \ell =1,2,\dots ,l,\hfill \\ \hfill {\displaystyle \sum _{p=1}^m}{\lambda }_p{p}_j\left({\mathbf{x}}_p\right)& =0,\hfill & \hfill & j=1,2,\dots ,s.\hfill \end{array}$$

The preceding relations can be assembled into the following symmetric linear system:

$$\mathbf{A}\Lambda =\mathbf{U}.$$

(27)

Here,

$$\mathbf{A}=\left[\begin{array}{ccc}\mathbf{{\rm Y}}& {\Psi }_2& \mathbf{P}\\ {\Psi }_1& {\Psi }_{12}& \mathbf{0}\\ {\mathbf{P}}^T& \mathbf{0}& \mathbf{0}\end{array}\right],\Lambda =\left[\begin{array}{c}\mathbf{\lambda }\\ \mathbf{\vartheta }\\ \mathbf{\nu }\end{array}\right],\mathbf{U}=\left[\begin{array}{c}\mathbf{u}\\ \mathcal{L}\mathbf{u}\\ \mathbf{0}\end{array}\right].$$

The entries of the block matrices are defined by

$$\begin{array}{cccc}\hfill {{\rm Y}}_{p,q}& =\phi \left(∥{\mathbf{x}}_p-{\mathbf{x}}_q∥\right),\hfill & \hfill & p,q=1,2,\dots ,m,\hfill \\ \hfill P_{p,j}& =p_j\left({\mathbf{x}}_p\right),\hfill & \hfill & p=1,2,\dots ,m,j=1,2,\dots ,s,\hfill \\ \hfill {\left({\Psi }_2\right)}_{p,\ell }& ={\mathcal{L}}_2\phi \left(∥{\mathbf{x}}_p-{\tilde{\mathbf{x}}}_{\ell }∥\right),\hfill & \hfill & p=1,2,\dots ,m,\ell =1,2,\dots ,n_h,\hfill \\ \hfill {\left({\Psi }_1\right)}_{\ell ,q}& ={\mathcal{L}}_1\phi \left(∥{\tilde{\mathbf{x}}}_{\ell }-{\mathbf{x}}_q∥\right),\hfill & \hfill & \ell =1,2,\dots ,n_h,q=1,2,\dots ,m,\hfill \\ \hfill {\left({\Psi }_{12}\right)}_{\ell ,h}& ={\mathcal{L}}_1{\mathcal{L}}_2\phi \left(∥{\tilde{\mathbf{x}}}_{\ell }-{\tilde{\mathbf{x}}}_h∥\right),\hfill & \hfill & \ell ,h=1,2,\dots ,n_h.\hfill \end{array}$$

Moreover,

$$\mathbf{\lambda }={[{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_m]}^T,\mathbf{\nu }={[{\nu }_1,{\nu }_2,\dots ,{\nu }_s]}^T,$$

$$\mathbf{u}={[u\left({\mathbf{x}}_1\right),u\left({\mathbf{x}}_2\right),\dots ,u\left({\mathbf{x}}_m\right)]}^T,\mathbf{\vartheta }={[{\vartheta }_1,{\vartheta }_2,\dots ,{\vartheta }_{n_h}]}^T,$$

and

$$\mathcal{L}\mathbf{u}={\left[\mathcal{L}u\left({\tilde{\mathbf{x}}}_1\right),\mathcal{L}u\left({\tilde{\mathbf{x}}}_2\right),\dots ,\mathcal{L}u\left({\tilde{\mathbf{x}}}_{n_h}\right)\right]}^T,$$

where the subscripts in $\mathcal{L}1$ and $\mathcal{L}2$ indicate the argument of the radial kernel with respect to which the differential operator acts. More precisely, $\mathcal{L}1$ acts on the first spatial argument, whereas $\mathcal{L}2$ acts on the second. Under the usual unisolvency assumptions, the block system in (27) admits a unique solution whenever the selected radial basis function is positive definite or conditionally positive definite and is supplemented by a polynomial space of the required degree [38,39]. To approximate the action of a linear differential operator $\mathcal{L}$ at the stencil center ${\mathbf{x}}_j$ within the LRBF–HFD approach, we first apply $\mathcal{L}$ to both sides of Equation (26) and evaluate the resulting expression at ${\mathbf{x}}_j$, yielding

$$\begin{array}{cc}\hfill \mathcal{L}\left(u\left({\mathbf{x}}_j\right)\right)=& {\displaystyle \sum _{i=1}^m}{\lambda }_i\mathcal{L}\phi \left(∥{\mathbf{x}}_j-{\mathbf{x}}_i∥\right)\hfill \\ \hfill & +{\displaystyle \sum _{i=1}^l}{\vartheta }_i\mathcal{L}{\mathcal{L}}_2\phi \left(∥{\mathbf{x}}_j-{\tilde{\mathbf{x}}}_i∥\right)+{\displaystyle \sum _{k=1}^s}{\nu }_k\mathcal{L}{p}_k\left({\mathbf{x}}_j\right)\hfill \\ \hfill =& [\mathcal{L}\phi \left(∥{\mathbf{x}}_j-{\mathbf{x}}_1∥\right),\dots ,\mathcal{L}\phi \left(∥{\mathbf{x}}_j-{\mathbf{x}}_m∥\right),\hfill \\ \hfill & \mathcal{L}{\mathcal{L}}_2\phi \left(∥{\mathbf{x}}_j-{\tilde{\mathbf{x}}}_1∥\right),\dots ,\mathcal{L}{\mathcal{L}}_2\phi \left(∥{\mathbf{x}}_j-{\tilde{\mathbf{x}}}_l∥\right),\hfill \\ \hfill & \mathcal{L}{p}_1\left({\mathbf{x}}_j\right),\dots ,\mathcal{L}{p}_s\left({\mathbf{x}}_j\right)]\Lambda \hfill \\ \hfill =& \mathbf{\psi }{\mathbf{A}}^{-1}\mathbf{U}:={\displaystyle \sum _{k=1}^m}{w}_{k}u\left({\mathbf{x}}_k\right)+{\displaystyle \sum _{k=1}^l}{\tilde{w}}_k\mathcal{L}u\left({\tilde{\mathbf{x}}}_k\right).\hfill \end{array}$$

(28)

where

$$\begin{array}{cc}\hfill \mathbf{\psi }=[& \mathcal{L}\phi \left(∥{\mathbf{x}}_j-{\mathbf{x}}_1∥\right),\dots ,\mathcal{L}\phi \left(∥{\mathbf{x}}_j-{\mathbf{x}}_m∥\right),\hfill \\ \hfill & \mathcal{L}{\mathcal{L}}_2\phi \left(∥{\mathbf{x}}_j-{\tilde{\mathbf{x}}}_1∥\right),\dots ,\mathcal{L}{\mathcal{L}}_2\phi \left(∥{\mathbf{x}}_j-{\tilde{\mathbf{x}}}_l∥\right),\hfill \\ \hfill & \mathcal{L}{p}_1\left({\mathbf{x}}_j\right),\dots ,\mathcal{L}{p}_s\left({\mathbf{x}}_j\right)].\hfill \end{array}$$

The unknown weights $w_k$, $k=1,2,\dots ,m$, and ${\tilde{w}}_k$, $k=1,2,\dots ,l$, are obtained by solving the following linear system and retaining only the components corresponding to these weights, while discarding the auxiliary coefficients ${\nu }_1,\dots ,{\nu }_s$:

$${\mathbf{A}}^T\left[\begin{array}{c}\mathbf{U}\\ \tilde{\mathbf{w}}\\ \mathbf{\nu }\end{array}\right]={\mathbf{\psi }}^T.$$

(29)

3.1. Spatial Discretization by the LRBF–HFD Method

Let $X=\left\{{\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_{N_s}\right\}\subset \overline{\Omega }$ be the set of spatial nodes, where $N_s$ denotes the total number of interior and boundary nodes. For every interior node ${\mathbf{x}}_j$, we define the local stencil ${\Omega }_j=\left\{{\mathbf{x}}_{j_1},{\mathbf{x}}_{j_2},\dots ,{\mathbf{x}}_{j_{m_j}}\right\},$ together with the Hermite derivative nodes ${\tilde{\Omega }}_j=\left\{{\tilde{\mathbf{x}}}_{j_1},{\tilde{\mathbf{x}}}_{j_2},\dots ,{\tilde{\mathbf{x}}}_{j_{l_j}}\right\}.$ By setting $\mathcal{L}=\Delta$ in the LRBF–HFD approximation, the Laplacian of a sufficiently smooth function v at ${\mathbf{x}}_j$ is approximated by

$$\Delta v\left({\mathbf{x}}_j\right)\approx {\displaystyle \sum _{r=1}^{m_j}}{w}_{j,r}v\left({\mathbf{x}}_{j_r}\right)+{\displaystyle \sum _{r=1}^{l_j}}{\tilde{w}}_{j,r}\Delta v\left({\tilde{\mathbf{x}}}_{j_r}\right).$$

(30)

The local weights are obtained from

$${\mathbf{A}}_j^T\left[\begin{array}{c}{\mathbf{w}}_j\\ {\tilde{\mathbf{w}}}_j\\ {\nu }_j\end{array}\right]={\psi }_j^T,$$

(31)

where ${\mathbf{A}}_j$ is the local augmented Hermite interpolation matrix and

$$\begin{array}{cc}\hfill {\psi }_j=[& \Delta \phi (\parallel {\mathbf{x}}_j-{\mathbf{x}}_{j_1}\parallel ),\dots ,\Delta \phi (\parallel {\mathbf{x}}_j-{\mathbf{x}}_{j_{m_j}}\parallel ),\hfill \\ \hfill & \Delta {\Delta }_2\phi (\parallel {\mathbf{x}}_j-{\tilde{\mathbf{x}}}_{j_1}\parallel ),\dots ,\Delta {\Delta }_2\phi (\parallel {\mathbf{x}}_j-{\tilde{\mathbf{x}}}_{j_{l_j}}\parallel ),\hfill \\ \hfill & \Delta p_1\left({\mathbf{x}}_j\right),\dots ,\Delta p_s\left({\mathbf{x}}_j\right)].\hfill \end{array}$$

After assembling the local approximations over all spatial nodes, we obtain

$${\mathbf{B}}_h{\mathbf{D}}^n={\mathbf{A}}_h{\mathbf{U}}^n,{\mathbf{B}}_h=\mathbf{I}-{\tilde{\mathbf{W}}}_h,{\mathbf{A}}_h={\mathbf{W}}_h,$$

(32)

where

$${\mathbf{U}}^n={\left[u^n\left({\mathbf{x}}_1\right),u^n\left({\mathbf{x}}_2\right),\dots ,u^n\left({\mathbf{x}}_{N_s}\right)\right]}^T$$

and ${\mathbf{D}}^n$ contains the corresponding approximations of the Laplacian. Therefore,

$${\mathbf{D}}^n={\mathbf{L}}_h{\mathbf{U}}^n,{\mathbf{L}}_h={\mathbf{B}}_h^{-1}{\mathbf{A}}_h.$$

(33)

In practical computations, the inverse of ${\mathbf{B}}_h$ is not formed explicitly. Instead, the matrix equation ${\mathbf{B}}_h{\mathbf{L}}_h={\mathbf{A}}_h$ is solved.

3.1.1. Space Discretization of the Nonlinear Time-Fractional Schrödinger Equation

We first define the term by

$${\mathbf{H}}_u^n={\displaystyle \sum _{k=1}^{n-1}}\left(b_{n-k-1}-b_{n-k}\right){\mathbf{U}}^k+b_{n-1}{\mathbf{U}}^0,n\ge 1.$$

(34)

For $n=1$, the summation term vanishes and hence

$${\mathbf{H}}_u^1={\mathbf{U}}^0.$$

The potential matrix is given by

$${\mathbf{V}}_h=diag\left(V\left({\mathbf{x}}_1\right),\dots ,V\left({\mathbf{x}}_{N_s}\right)\right),$$

and the nonlinear term is defined as

$$\mathcal{N}\left(\mathbf{U}\right)={\left|\mathbf{U}\right|}^2\odot \mathbf{U},$$

where ⊙ denotes componentwise multiplication. Using the relation

$${\mathbf{D}}^n={\mathbf{L}}_h{\mathbf{U}}^n,$$

the fully discrete scheme can be written as

$$\left(\mathrm{i}{\mu }_{\mathsf{\beta }}\mathbf{I}+\alpha {\mathbf{L}}_h+{\mathbf{V}}_h\right){\mathbf{U}}^n={\mathbf{F}}^n+\mathrm{i}{\mu }_{\mathsf{\beta }}{\mathbf{H}}_u^n-\lambda \mathcal{N}\left({\mathbf{U}}^{n-1}\right),n=1,2,\dots ,N,$$

(35)

where

$${\mathbf{F}}^n={\left[f^n\left({\mathbf{x}}_1\right),f^n\left({\mathbf{x}}_2\right),\dots ,f^n\left({\mathbf{x}}_{N_s}\right)\right]}^T.$$

Since ${\mathbf{B}}_h{\mathbf{L}}_h={\mathbf{A}}_h,$ multiplying Eq. (35) by ${\mathbf{B}}_h$ gives

$$\left(\mathrm{i}{\mu }_{\mathsf{\beta }}{\mathbf{B}}_h+\alpha {\mathbf{A}}_h+{\mathbf{B}}_h{\mathbf{V}}_h\right){\mathbf{U}}^n={\mathbf{B}}_h\left[{\mathbf{F}}^n+\mathrm{i}{\mu }_{\mathsf{\beta }}{\mathbf{H}}_u^n-\lambda \mathcal{N}\left({\mathbf{U}}^{n-1}\right)\right].$$

(36)

This form is convenient because ${\mathbf{L}}_h$ does not need to be computed explicitly. At the first time level, ${\mathbf{H}}_u^1={\mathbf{U}}^0$, and therefore

$$\left(\mathrm{i}{\mu }_{\mathsf{\beta }}{\mathbf{B}}_h+\alpha {\mathbf{A}}_h+{\mathbf{B}}_h{\mathbf{V}}_h\right){\mathbf{U}}^1={\mathbf{B}}_h\left[{\mathbf{F}}^1+\mathrm{i}{\mu }_{\mathsf{\beta }}{\mathbf{U}}^0-\lambda \mathcal{N}\left({\mathbf{U}}^0\right)\right].$$

Finally, the IC and BCs are imposed as ${\mathbf{U}}^0={\left[\psi \left({\mathbf{x}}_1\right),\dots ,\psi \left({\mathbf{x}}_{N_s}\right)\right]}^T,$and $u_j^n=h({\mathbf{x}}_j,t_n),{\mathbf{x}}_j\in \partial \Omega .$

3.1.2. Space Discretization of the Time-Fractional Schrödinger–Poisson System

For the coupled system, define

$${\Psi }^n={\left[{\psi }^n\left({\mathbf{x}}_1\right),\dots ,{\psi }^n\left({\mathbf{x}}_{N_s}\right)\right]}^T,{\Phi }^n={\left[{\Phi }^n\left({\mathbf{x}}_1\right),\dots ,{\Phi }^n\left({\mathbf{x}}_{N_s}\right)\right]}^T.$$

Similarly, let

$${\mathbf{F}}^n={\left[f^n\left({\mathbf{x}}_1\right),\dots ,f^n\left({\mathbf{x}}_{N_s}\right)\right]}^T,{\mathbf{G}}^n={\left[g^n\left({\mathbf{x}}_1\right),\dots ,g^n\left({\mathbf{x}}_{N_s}\right)\right]}^T.$$

The history term is defined by

$${\mathbf{H}}_{\psi }^n={\displaystyle \sum _{k=1}^{n-1}}\left(b_{n-k-1}-b_{n-k}\right){\Psi }^k+b_{n-1}{\Psi }^0.$$

(37)

We also define

$$\mathcal{R}(\Psi )={|\Psi |}^2,\mathcal{N}(\Psi )={|\Psi |}^2\odot \Psi .$$

At each time level, we first solve the Poisson equation. Using the LRBF–HFD approximation, we obtain

$$-{\mathbf{L}}_h{\Phi }^n=c_D\mathcal{R}\left({\Psi }^{n-1}\right)+{\mathbf{G}}^n.$$

(38)

Using ${\mathbf{B}}_h{\mathbf{L}}_h={\mathbf{A}}_h,$ the above equation can be written as

$$-{\mathbf{A}}_h{\Phi }^n={\mathbf{B}}_h\left[c_D\mathcal{R}\left({\Psi }^{n-1}\right)+{\mathbf{G}}^n\right].$$

(39)

After finding ${\Phi }^n$, we define

$${\mathbf{P}}^n=diag\left({\Phi }^n\left({\mathbf{x}}_1\right),\dots ,{\Phi }^n\left({\mathbf{x}}_{N_s}\right)\right).$$

The discrete Schrödinger equation is then

$$\left(\mathrm{i}{\mu }_{\mathsf{\beta }}\mathbf{I}+a{\mathbf{L}}_h-\alpha {\mathbf{P}}^n\right){\Psi }^n={\mathbf{F}}^n+\mathrm{i}{\mu }_{\mathsf{\beta }}{\mathbf{H}}_{\psi }^n-\lambda \mathcal{N}\left({\Psi }^{n-1}\right).$$

(40)

Its compact form is

$$\left(\mathrm{i}{\mu }_{\mathsf{\beta }}{\mathbf{B}}_h+a{\mathbf{A}}_h-\alpha {\mathbf{B}}_h{\mathbf{P}}^n\right){\Psi }^n={\mathbf{B}}_h\left[{\mathbf{F}}^n+\mathrm{i}{\mu }_{\mathsf{\beta }}{\mathbf{H}}_{\psi }^n-\lambda \mathcal{N}\left({\Psi }^{n-1}\right)\right].$$

(41)

For $n=1$, we have

$${\mathbf{H}}_{\psi }^1={\Psi }^0.$$

Therefore, the first time step is obtained from

$$\begin{array}{cc}\hfill -{\mathbf{A}}_h{\Phi }^1& ={\mathbf{B}}_h\left[c_D\mathcal{R}\left({\Psi }^0\right)+{\mathbf{G}}^1\right],\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill \left(\mathrm{i}{\mu }_{\mathsf{\beta }}{\mathbf{B}}_h+a{\mathbf{A}}_h-\alpha {\mathbf{B}}_h{\mathbf{P}}^1\right){\Psi }^1& ={\mathbf{B}}_h\left[{\mathbf{F}}^1+\mathrm{i}{\mu }_{\mathsf{\beta }}{\Psi }^0-\lambda \mathcal{N}\left({\Psi }^0\right)\right].\hfill \end{array}$$

(43)

The initial condition is ${\Psi }^0={\left[{\psi }_0\left({\mathbf{x}}_1\right),\dots ,{\psi }_0\left({\mathbf{x}}_{N_s}\right)\right]}^T,$ and the boundary conditions are imposed as

$${\psi }_j^n=h_{\psi }({\mathbf{x}}_j,t_n),{\Phi }_j^n=h_{\Phi }({\mathbf{x}}_j,t_n),{\mathbf{x}}_j\in \partial \Omega .$$

Thus, at each time level, we first compute ${\Phi }^n$ from Eq. (39). Then, the computed potential is used in Eq. (41) to obtain ${\Psi }^n$.

4. Numerical Results and Discussion

This section presents several numerical examples to demonstrate the accuracy and effectiveness of the LRBF-FD method for solving nonlinear time-fractional Schrödinger equations and the time-fractional Schrödinger–Poisson system. To examine the temporal convergence, the experimental convergence order ${\mathcal{C}}_t$ is calculated by

$$\begin{array}{c}\hfill {\mathcal{C}}_t={log}_2\left({\displaystyle \frac{L_{\infty }(2\delta t,h)}{L_{\infty }(\delta t,h)}}\right).\end{array}$$

Here, the maximum absolute error at the final time T is defined by $L_{\infty }=\underset{1\le j\le N-1}{max}\left|u(x_j,T)-U(x_j,T)\right|$, where $u(x_j,T)$ and $U(x_j,T)$ are the exact and numerical solutions, respectively, at the node $x_j$. For the numerical simulations, we employ the LRBF–HFD method with the polyharmonic spline radial basis function $\phi \left(r\right)=r^7$, augmented with polynomials of degree five. In Examples 4 and 5, we consider an irregular star-shaped region whose boundary is given by $r\left(\theta \right)={\displaystyle \frac{81-9cos\left(8\theta \right)}{64}},0\le \theta \le 2\pi$, and discretize the domain using Halton nodes.

Example 1. 

Let us consider the following nonlinear time-fractional Schrödinger equation [40]:

$$\begin{array}{cc}\hfill {\mathrm{i}}_0^{\mathrm{C}}{\mathbf{D}}_t^{\mathsf{\beta }}u(x,t)+{\displaystyle \frac{{\partial }^{2}u(x,t)}{\partial x^2}}+{\left|u(x,t)\right|}^{2}u(x,t)& =f(x,t),\hfill \\ \hfill 0\le x\le 1,0<t\le T,0<\mathsf{\beta }<1.\end{array}$$

The source function is defined by

$$\begin{array}{cc}\hfill f(x,t)=& -{\displaystyle \frac{2t^{2-\mathsf{\beta }}}{\Gamma (3-\mathsf{\beta })}}cos\left(2\pi x\right)+\left(-4{\pi }^2{t}^2+t^6\right)sin\left(2\pi x\right)+\mathrm{i}\left[{\displaystyle \frac{2t^{2-\mathsf{\beta }}}{\Gamma (3-\mathsf{\beta })}}sin\left(2\pi x\right)+\left(-4{\pi }^2{t}^2+t^6\right)cos\left(2\pi x\right)\right].\hfill \end{array}$$

The exact solution is $u(x,t)=t^2\left[sin\left(2\pi x\right)+\mathrm{i}cos\left(2\pi x\right)\right].$ The IC and BCs can be computed directly from the exact solution.

The LRBF–HFD method is adopted to study the numerical accuracy corresponding to various values of $\delta t$, N, and $\beta$. Table 1 gives the real, imaginary, and $L_{\infty }$ errors together with the temporal convergence rates for $\mathsf{\beta }\in \{0.4,0.7\}$ at $T=1$, showing that the rates approach first-order convergence as the time step decreases. Table 2 presents the real, imaginary, and $L_{\infty }$ errors for different numbers of spatial nodes with $\mathsf{\beta }\in \{0.25,0.85\}$ at $T=1$, showing a clear decrease in the numerical errors as the number of nodes increases. Table 3 summarizes the $L_{\infty }$, $L_2$, and RMS errors together with the CPU running times for different final times and fractional orders, showing the effects of $\mathsf{\beta }$ and T on the numerical accuracy and computational cost. Finally, Figure 2 compares the approximate and exact real and imaginary parts and their corresponding absolute errors, showing good agreement between the numerical and exact solutions.

Example 2. 

We consider the following nonlinear time-fractional Schrödinger equation [41]:

$$\begin{array}{cc}\hfill {\mathrm{i}}_0^{\mathrm{C}}{\mathbf{D}}_t^{\mathsf{\beta }}u(x,t)+{\displaystyle \frac{1}{2}}{\displaystyle \frac{{\partial }^{2}u(x,t)}{\partial x^2}}-{cos}^2\left(x\right)u(x,t)-{\left|u(x,t)\right|}^{2}u(x,t)& =0,\hfill \\ \hfill 0\le x\le 1,0<t\le T,0<\mathsf{\beta }<1.\end{array}$$

The exact solution is

$$u(x,t)=sin\left(x\right){\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \frac{{\left(-3\mathrm{i}{t}^{\mathsf{\beta }}\right)}^k}{2^k\Gamma (k\mathsf{\beta }+1)}}=sin\left(x\right)E_{\mathsf{\beta }}\left(-{\displaystyle \frac{3\mathrm{i}}{2}}{t}^{\mathsf{\beta }}\right).$$

The associated IC and BCs can be readily constructed using the given exact solution.

The LRBF–HFD method is applied to evaluate the numerical accuracy for various values of $\delta t$, N, and $\beta$. Table 4 gives the real, imaginary, and $L_{\infty }$ errors together with the temporal convergence rates for $\mathsf{\beta }\in \{0.45,0.75\}$ at $T=1$, showing that the numerical errors decrease as the time step becomes smaller. Table 5 provides similar results for $\mathsf{\beta }\in \{0.2,0.9\}$ at $T=1$, confirming the temporal convergence of the method for different fractional orders. Table 6 lists the $L_{\infty }$, $L_2$, and RMS errors along with the CPU running times for different final times and fractional orders, showing the effect of $\mathsf{\beta }$ and T on accuracy and computational cost. Finally, Figure 3 compares the approximate and exact real and imaginary parts and their absolute errors, showing good agreement between the numerical and exact solutions.

Example 3. 

Consider the following two-dimensional nonlinear time-fractional Schrödinger equation [42]:

$$\begin{array}{cc}\hfill {\mathrm{i}}_0^{\mathrm{C}}{\mathbf{D}}_t^{\mathsf{\beta }}u(x,y,t)=& -{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{{\partial }^{2}u(x,y,t)}{\partial x^2}}+{\displaystyle \frac{{\partial }^{2}u(x,y,t)}{\partial y^2}}\right)\hfill \\ \hfill & +\left[1-{sin}^2\left(x\right){sin}^2\left(y\right)\right]u(x,y,t)+{\left|u(x,y,t)\right|}^{2}u(x,y,t),\hfill \\ \hfill & (x,y)\in [0,1]\times [0,1],0<t\le T,0<\mathsf{\beta }<1.\hfill \end{array}$$

The exact solution is stated as

$$\begin{array}{cc}\hfill u(x,y,t)& =sin\left(x\right)sin\left(y\right){\displaystyle \sum _{n=0}^{\infty }}{\displaystyle \frac{{\left(-5\mathrm{i}{t}^{\mathsf{\beta }}\right)}^n}{2^n\Gamma (n\mathsf{\beta }+1)}}=sin\left(x\right)sin\left(y\right)E_{\mathsf{\beta }}\left(-{\displaystyle \frac{5\mathrm{i}}{2}}{t}^{\mathsf{\beta }}\right).\hfill \end{array}$$

The associated IC and BCs can be readily constructed using the given exact solution.

The LRBF–HFD method is implemented to examine the accuracy of the numerical solutions for different values of $\delta t$, N, and $\beta$. Table 7 gives the real, imaginary, and $L_{\infty }$ errors together with the temporal convergence rates for $\mathsf{\beta }\in \{0.35,0.95\}$, showing that the numerical errors decrease as the time step becomes smaller. Table 8 presents similar results for $\mathsf{\beta }\in \{0.4,0.6\}$, confirming the temporal convergence of the proposed method. Table 9 summarizes the $L_{\infty }$, $L_2$, and RMS errors together with the CPU times for different final times and fractional orders, showing the effects of $\mathsf{\beta }$ and T on the numerical accuracy and computational cost. Finally, Figure 4 shows the numerical real and imaginary parts, the corresponding absolute error, and the wave density ${\left|u\right|}^2$ for $\beta =0.6$ and $\delta t=1/200$, demonstrating the accuracy of the numerical solution.

Example 4. 

Next, we consider the following time-fractional nonlinear Schrödinger equation:

$$\begin{array}{cc}\hfill {\mathrm{i}}_0^{\mathrm{C}}{\mathbf{D}}_t^{\mathsf{\beta }}u(x,y,t)& -{\displaystyle \frac{{\partial }^{2}u(x,y,t)}{\partial x^2}}-{\displaystyle \frac{{\partial }^{2}u(x,y,t)}{\partial y^2}}+{\left|u(x,y,t)\right|}^{2}u(x,y,t)\hfill \\ \hfill & ={\displaystyle \frac{2\mathrm{i}{t}^{2-\mathsf{\beta }}}{\Gamma (3-\mathsf{\beta })}}sin(x+y)+t^{2}sin(x+y)\left[2+{\left|t^{2}sin(x+y)\right|}^2\right],\hfill \end{array}$$

where $0<t\le T$, $0<\mathsf{\beta }<1$. The exact solution is $u(x,y,t)=t^{2}sin(x+y).$ The IC and BCs can be obtained directly from the above exact solution.

The LRBF–HFD method is utilized to assess the accuracy of the numerical results for different values of $\Delta t$, N, and $\beta$. Table 10 gives the real, imaginary, and $L_{\infty }$ errors for $\mathsf{\beta }\in \{0.2,0.8\}$, showing approximately first-order temporal convergence as the time step decreases. Table 11 presents similar results for $\mathsf{\beta }\in \{0.35,0.75\}$, where the convergence rates also approach one for smaller time steps. Table 12 summarizes the $L_{\infty }$, $L_2$, and RMS errors together with the CPU running times for different final times and fractional orders, showing the effects of $\mathsf{\beta }$ and T on the numerical accuracy and computational cost. Figure 5 compares the numerical and exact solutions and shows the corresponding absolute error and wave density ${\left|u\right|}^2$ for $\beta =0.75$ and $\delta t=1/200$, demonstrating the accuracy of the numerical approximation. Finally, Figure 6 illustrates the numerical and exact solutions on the irregular domain, together with the associated absolute error and wave density ${\left|u\right|}^2$, for $\beta =0.5$ and $\delta t=1/200$, indicating the high accuracy of the numerical results.

Example 5. 

Finally, we consider the time-fractional Schrödinger–Poisson system

$$\begin{array}{cc}\hfill & {\mathrm{i}}_0^{\mathrm{C}}{\mathbf{D}}_t^{\mathsf{\beta }}\psi +\Delta \psi +{\left|\psi \right|}^2\psi =f,\hfill \\ \hfill & -\Delta \Phi ={\left|\psi \right|}^2+g,\hfill \end{array}$$

on $\Omega$ for $0<t\le T$. The exact solutions are

$$\begin{array}{cc}\hfill \psi (x,y,t)& =(1+t^2)sin\left(\pi x\right)sin\left(\pi y\right),\hfill \\ \hfill \Phi (x,y,t)& =t^{2}sin\left(\pi x\right)sin\left(\pi y\right).\hfill \end{array}$$

The IC and BCs are determined from the aforementioned exact solutions.

The LRBF–HFD method is employed to investigate the accuracy for different values of $\delta t$, N, and $\beta$. Table 13 gives the temporal $L_{\infty }$ errors and the corresponding convergence rates for $\psi$ and $\Phi$ with $\mathsf{\beta }\in \{0.25,0.55\}$, confirming the convergence of the proposed scheme. Table 14 presents similar results for $\mathsf{\beta }\in \{0.7,0.9\}$, showing that the errors decrease as the time step becomes smaller. Table 15 summarizes the $L_{\infty }$, $L_2$, and RMS errors for $\psi$ and $\Phi$ at different final times and fractional orders, together with the corresponding computational costs. Figure 7 shows the numerical solutions of $\psi$ and $\Phi$ and their corresponding absolute errors on the square domain for $\beta =0.5$ and $\delta t=1/200$, demonstrating the accuracy of the proposed method. Finally, Figure 8 presents the numerical solutions of $\psi$ and $\Phi$ on the irregular domain, together with their corresponding absolute errors, for $\beta =0.5$ and $\delta t=1/200$.

5. Conclusion

This paper developed an efficient localized meshless method for nonlinear time-fractional Schrödinger equations and Schrödinger–Poisson systems. The Caputo fractional derivative was approximated by the L1 formula, while the spatial derivatives were discretized using the LRBF-HFD method. The localized formulation avoided the large and ill-conditioned matrices of global RBF methods and reduced the computational cost. Several numerical examples were used to examine the accuracy and performance of the method for different fractional orders, spatial discretizations, and final times. The numerical results showed good agreement with the exact solutions and a clear decrease in the errors under mesh refinement. Therefore, the proposed method provided an accurate and efficient approach for solving the considered fractional models.