Quintic non-polynomial spline for time-fractional nonlinear Schrödinger equation

In this paper, we shall solve a time-fractional nonlinear Schrödinger equation by using the quintic non-polynomial spline and the L1 formula. The unconditional stability, unique solvability and convergence of our numerical scheme are proved by the Fourier method. It is shown that our method is sixth order accurate in the spatial dimension and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$(2-\gamma )$\end{document}(2−γ)th order accurate in the temporal dimension, where γ is the fractional order. The efficiency of the proposed numerical scheme is further illustrated by numerical experiments, meanwhile the simulation results indicate better performance over previous work in the literature.


Introduction
In the past few decades, fractional differential equations have gained much importance due to their usefulness in modeling phenomena in various areas such as physics, engineering, finance, biology and chemistry [12,33,46,51]. To cite some recent developments: in 2019 Jajarmi and Baleanu [26] studied a general form of fractional optimal control problems involving fractional derivative with singular or non-singular kernel; Jothimani et al. [30] discussed an exact controllability of nondensely defined nonlinear fractional integrodifferential equations with the Hille-Yosida operator; Valliammal et al. [61] studied the existence of mild solutions of fractional-order neutral differential system with statedependent delay in Banach space. In 2020 Jajarmi et al. [28] investigated a fractional version of SIRS model for the HRSV disease; Baleanu et al. [2] proposed a new fractional model for the human liver involving the Caputo-Fabrizio fractional derivative; Baleanu et al. [3] studied the fractional features of a harmonic oscillator with position-dependent mass; Sajjadi et al. [54] discussed the chaos control and synchronization of a hyperchaotic model in both the frameworks of classical and of fractional calculus; Jajarmi and Baleanu [27] proposed a new iterative method to generate the approximate solution of nonlinear fractional boundary value problems in the form of uniformly convergent series; Shiri et al. [56] employed discretized collocation methods for a class of tempered fractional differential equations with terminal value problems; Tuan et al. [60] tackled the problem of finding the solution of a multi-dimensional time-fractional reaction-diffusion equation with nonlinear source from the final value data; Li et al. [44] proposed a new approximation for the generalized Caputo fractional derivative based on WSGL formula and solved a generalized fractional sub-diffusion problem; Gao et al. [19] studied the epidemic predictability for the novel coronavirus (2019-nCoV) pandemic by analyzing a time-fractional model and finding its solution by a q-homotopy analysis transform method; Gao et al. [20] investigated the infection system of the novel coronavirus (2019-nCoV) with a nonlocal operator defined in the Caputo sense; Gao et al. [21] tackled the fractional Phi-four equation by using a q-homotopy analysis transform method numerically; Sabir et al. [53] presented a novel meta-heuristic computing solver for solving the singular three-point second-order boundary value problems using artificial neural networks.
The subject of the present work, the Schrödinger equation, was first proposed by the Austrian physicist Erwin Schrödinger in 1926 [55]. It is a fundamental equation in quantum physics that describes the evolution of the position-space wave function of a particle. In fact, the nonlinear Schrödinger equations describe a wide class of physical phenomena such as models of protein dynamics, self-focusing in laser pulses and nonlinear fiber optics [1,11,17,58]. The Schrödinger equation has also been generalized to fractional differential equations. In 2000, Laskin [34] generalized the non-fractional Schrödinger equation to a space-fractional Schrödinger equation by using the Feynman path integrals over the Lévy trajectories and replacing the quantum Riesz derivative with the Laplace operator. Later, in 2004 Naber [50] proposed a different generalization by changing the first order timederivative to a Caputo fractional derivative-this time-fractional Schrödinger equation has been used to describe fractional quantum mechanical behavior. In 2010, Muslih et al. [49] obtained a fractional Schrödinger equation by using a fractional variational principle and a fractional Klein-Gordon equation. In 2017, Gómez-Aguilar and Baleanu [23] presented an alternative model of fractional Schrödinger equation involving Caputo-Fabrizio fractional operator.
Many researchers pay attention to the numerical treatment of fractional Schrödinger equations. For the space-fractional Schrödinger equation: a linear implicit conservative difference scheme of order O(τ 2 + h 2 ) has been proposed in [62] for the case of a coupled nonlinear Schrödinger equation, where τ is the temporal step size and h is the spatial step size; Zhao et al. [65] have used a compact operator to approximate the Riesz derivative, and the proposed linearized difference scheme for a two-dimensional nonlinear spacefractional Schrödinger equation can achieve O(τ 2 +h 4 ), whereh = max{h 1 , h 2 }, h 1 and h 2 are the spatial step sizes in the x and y dimensions, respectively; Wang and Huang [64] have presented a conservative linearized difference scheme, which can achieve the order of O(τ 2 + h 2 ); a collocation method has been applied to a multi-dimensional spacetime variable-order fractional Schrödinger equation in [5]; a fourth-order implicit timediscretization scheme based on the exponential time differencing approach together with a fourth-order compact scheme in space have been proposed in [31], the method is of order O(τ 4 + h 4 ); Li et al. [37] have used a fast linearized conservative finite element method to solve the coupled type equation; Wang and Xiao [63] have proposed an efficient conservative scheme for the fractional Klein-Gordon-Schrödinger equation with central difference and Crank-Nicolson scheme, their method can achieve O(τ 2 + h 2 ). It is also noted that Hashemi and Akgül [24] have utilized Nucci's reduction method and the simplest equation method to extract analytical solutions specially of soliton kinds of nonlinear Schrödinger equations in both time and space fractional terms.
For the time-fractional Schrödinger equation: Khan et al. [32] have applied the homotopy analysis method; Mohebbi et al. [47] have used the Kansa approach to approximate the spatial derivative and L1 discretization to approximate the Caputo time-fractional derivative; a Krylov projection method has been developed in [22]; a Jacobi spectral collocation method has been applied to a multi-dimensional time-fractional Schrödinger equation in [4]; a quadratic B-spline Galerkin method combined with L1 discretization scheme has been proposed in [16]; a linearized L1-Galerkin finite element method has been used in [35] for a multi-dimensional nonlinear time-fractional Schrödinger equation; a cubic non-polynomial spline method combined with L1 discretization has been proposed in [36] and the stability has been shown by the Fourier method, the convergence order is not proved but is observed from numerical experiments to be O(τ 2-γ + h 4 ).
Motivated by the above research, in this paper we consider the following time-fractional nonlinear Schrödinger equation: is the Caputo fractional derivative of order γ ∈ (0, 1) defined by [12,51] We shall employ a quintic non-polynomial spline together with L1 discretization to solve (1.1). The stability, unique solvability and convergence of our numerical scheme are then proved by the Fourier method-we note that this method of proof is rare for numerical methods of time/space-fractional Schrödinger equation, especially in establishing the convergence order; on the other hand, the energy method has been commonly used to show the convergence of numerical methods for space-fractional Schrödinger equation [31,[62][63][64][65]. By the Fourier method, it is shown that our method is of order O(τ 2-γ + h 6 )-this improves the spatial convergence achieved by other methods for timefractional Schrödinger equation. Further, on the choices of our tools, we have observed in several different problems that a non-polynomial spline usually exhibits a better approximation than a polynomial spline because of its parameter [13, 15, 25, 38-43, 52, 57]; while L1 discretization is a stable and widely used approximation for the Caputo fractional derivative [15,18,29,38,48]. The organization of this paper is as follows. We derive the numerical scheme in Sect. 2. The stability, unique solvability and convergence are established by the Fourier method in Sects. 3, 4 and 5 respectively. In Sect. 6, we present three examples to verify the efficiency of our numerical scheme and to compare with other methods in the literature. Finally, a conclusion is drawn in Sect. 7.

Derivation of the numerical scheme
In this section, we shall develop a numerical scheme for problem (1.1) by using quintic non-polynomial spline and L1 discretization. The details of quintic non-polynomial spline will be presented first. LetP be uniform meshes of the spatial interval [0, L] with step size h = L M and the temporal interval [0, T] with step size τ = T N , respectively. For any given function y(x, t), we denote y(x j , t n ) by y n j , and y(x j , t) by y j for fixed t. Let u(x, t) denote the exact solution of (1.1) and U n j denote the numerical approximation of u n j . We shall set U n j to be the value of the quintic non-polynomial spline at (x j , t n ). We define the quintic non-polynomial spline as follows.

Definition 2.1 ([57]
) Let t = t n , 1 ≤ n ≤ N be fixed. For a given meshP, we say P n (x) is the quintic non-polynomial spline with parameter k (> 0) if P n (x) ∈ C (4) [0, L], P n (x) has the form span{1, x, x 2 , x 3 , sin(kx), cos(kx)} and its restriction P j, From the above definition, we can express P j,n (x) on [x j , x j+1 ], 0 ≤ j ≤ M -1 as Denote ω = kh. Using (2.2), a direct computation gives a n j = U n j -  Using the continuity of the first and third derivatives of the spline at x = x j+1 , i.e., P (1) j,n (x j+1 ) = P (1) j+1,n (x j+1 ) and P (3) j,n (x j+1 ) = P (3) j+1,n (x j+1 ), we obtain the following relations for Note that the consistency relation for (2.5)(b) will lead to 2(α + β) = 1. Manipulating (2.5)(a) and (2.5)(b), we can easily get Using the quintic non-polynomial spline to approximate the exact solution u(x, t) of (1.1), the spline relation (2.7) leads to where ϒ n j is the local truncation error in the spatial dimension. The next lemma gives a result on this error.

Lemma 2.1 For any fixed t
then the local truncation error ϒ n j associated with the spline relation (2.9) satisfies Proof We carry out the Taylor expansion at x = x j in (2.9), this gives To achieve the highest order of O(h 6 ), we set Z 1 = Z 2 = Z 3 = 0, which together with the consistency relation 2(α + β) = 1, gives (2.10) and (2.11).
A similar result to Lemma 2.1 has also been obtained in [38].
Remark 2.1 In order to compute the numerical solution U n j , 1 ≤ j ≤ M -1, we need another two equations besides (2.7) or (2.9). We consider the following equations which incorporate the boundary conditions in (1.1): where ϒ n 1 and ϒ n M-1 are local truncation errors in the spatial dimension, and the constant a i and b i have to be computed such that By carrying out Taylor expansion at x = x 2 and x = x M-2 in (2.12)(a) and (2.12)(b), respectively, we get the following which satisfies (2.13): Remark 2.2 If we let the spline parameter α = 1 10 in (2.9), then (2.9)| j=2 and (2.9)| j=M-2 are simply the same as (2.12)(a) and (2.12)(b), respectively. Therefore, we should have α = 1 10 .
To simplify the notations of the spline relations (2.9) and (2.12), we introduce the following definition.

Definition 2.2
For y = (y 1 , . . . , y M-1 ), we define the operators ∧ and ∧ 1 by and  .12) can be presented as The next lemma gives the L1 discretization for the Caputo fractional derivative.
where τ is the temporal step size, We are now ready to derive the numerical scheme for (1.1). To begin, we discretize (1.1) at (x j , t n ) to get Using the L1 discretization (2.17) in (2.20), we obtain Next, applying the operator ∧ to (2.21) yields Noting (2.16), it follows that Upon rearranging the above relation, we have After omitting the error and replacing the exact solution u with the numerical solution U, we get the numerical scheme (2.24) Remark 2.4 It is obvious that our scheme (2.23) is a nonlinear scheme. To linearize (2.23), following [45] we shall use the iterative algorithm where U n(r) j is the rth iterate of U n j , and (2.26) Note that the scheme (2.25) is linear in U n(r+1) . Hence, instead of (2.23)-(2.24), in practice we shall employ the linearized iterative scheme (2.25) with (2.27) For practical purposes, we would certainly need a stopping criterion to get U n(r+1) j to a desired accuracy. An example of such a stopping criterion is In fact, we shall use the above stopping criterion with the small constant as 1 × 10 -6 for our numerical simulations in Sect. 6.

Stability analysis of the numerical scheme
In this section, we shall analyze the stability of the scheme ( has also been used in [14,15]. With the linearization, we can rewrite (2.23)-(2.24) as (3.1) Consider the perturbed system of (3.1) with perturbation in the initial values Let U n j be the solution of (3.1) and letÛ n j be the solution of the perturbed system (3.2). To apply the Fourier method, we expand ρ n to a piecewise constant function ρ n (x), where Then ρ n (x) can be expanded as a Fourier series, To carry out the Fourier stability analysis, it is sufficient to consider an individual harmonic of the form [6] ρ n j (m) = d n (m)e i2π mjh/L . Proof For notational simplicity, let d n ≡ d n (m) for a fixed wave number m. We shall use mathematical induction to complete the proof. First, we consider n = 1. Upon substituting (3.8) into (3.3), we obtain After a series of computations, (3.10) leads to and (3.14) It is obvious that | | ≤ 1, therefore from (3.11) we have Next, assume that |d | ≤ |d 0 |, 1 ≤ ≤ n -1. We shall show that |d n | ≤ |d 0 |. For any n ≥ 2, we substitute (3.8) into (3.3) to get Using Lemma 2.3, |d | ≤ |d 0 | for 1 ≤ ≤ n -1, and | | ≤ 1, it follows from (3.17) that  Hence, we have completed the proof of (3.9).

Theorem 3.1 (Stability) The numerical scheme (2.25)-(2.27) or equivalently (3.1) is unconditionally stable with respect to the initial data.
Proof From the definition of the L2 norm (3.4) and (3.5), we find which shows that the numerical scheme (3.1) is robust to perturbation of initial data.

Solvability of the numerical scheme
In this section, we shall investigate the solvability of the numerical scheme (2.25)-(2.27) or equivalently (3.1). It is clear that the corresponding homogeneous system of (3.1) is (4.1) Similar to the proof of Theorem 3.1, we can verify that the solution of (4.1) satisfies U n 2 ≤ U 0 2 , 1≤ n ≤ N.
Since U 0 = 0, it follows that U n = 0 for 1 ≤ n ≤ N . Hence, (4.1) has only the trivial solution and we obtain the following theorem.

Convergence of the numerical scheme
In this section, we shall establish the convergence of the numerical scheme (2.25)-(2.27) or equivalently (3.1). Recall that u(x, t) is the exact solution of (1.1) and U n j is the numerical approximation of u n j obtained from (3.1). Let the error E n j at (x j , t n ) be E n j = u n j -U n j , 0 ≤ j ≤ M, 0 ≤ n ≤ N .
Clearly, from (3.1) we obtain the error equation where T n j is the local truncation error, In view of (2.16), (2.17) and (2.22), we see that T n j = O(h 6 + τ 2-γ ) and there is a constant C 1 such that It follows that Next, similar to (3.5)-(3.7), we define the piecewise constant functions E n (x) and T n (x) by and Then E n (x) and T n (x) have the Fourier series expansions Further, noting (5.2), we see that ∞ m=-∞ |δ n (m)| 2 converges and there is a positive constant C 2 ≥ 1 [7] such that δ n (m) ≤ C 2 δ 1 (m) , 1≤ n ≤ N. where θ = 2πm/L. Denote n ≡ n (m) and δ n ≡ δ n (m) for a fixed wave number m. Upon substituting (5.9) into (5.1), we get Now, we shall present two lemmas which are essential in the proof of the convergence result. The proof is completed.
Remark 5.1 As shown in Theorem 5.1, the numerical scheme (2.25)-(2.27) achieves sixth order convergence in the spatial dimension. This improves the work of [36] where the spatial convergence order is observed to be four through numerical experiment. Furthermore, in [36] the authors have only proved the stability of their scheme, while we have proven the stability, unique solvability and convergence of our method by the Fourier method. We remark that the Fourier method is rarely used in the analysis of numerical methods of time/space-fractional Schrödinger equation, especially in establishing the convergence order, so we have successfully illustrated the analytical technique of the Fourier method in this work. It is noted that the energy method has been commonly used to show the convergence of numerical methods for the space-fractional Schrödinger equation [31,[62][63][64][65].

Numerical examples
In this section, we shall present three numerical examples to verify the efficiency of the scheme (2.25)-(2.27) and to compare with other methods in the literature. Note that our theoretical convergence result is in L2 norm, nonetheless it would also be interesting to see the convergence in maximum norm. In fact, for a fixed pair of (h, τ ), we shall compute the following accuracy indicators: maximum L2 norm error E 2 (h, τ ), maximum modulus error E ∞ (h, τ ), maximum real part absolute error E Re ∞ (h, τ ) and maximum imaginary part absolute error E Im ∞ (h, τ ), defined by Im u n j -U n j . (6. 3) The convergence orders in temporal and spatial dimensions can be computed by where 0 < γ < 1 and The exact solution is u(x, t) = t 2 [sin(2πx) + i cos(2πx)].
Let α = 1/20 be fixed. We shall apply the scheme (2.25)-(2.27) to compute the errors and convergence orders, and comparisons are made with other methods in the literature. There follows a brief description of the numerical simulation in Tables 1-4

and Figs. 1-3:
(i) In Table 1, we fix h = 1/1000 and let τ vary. Applying (2.25)-(2.27), we compute E 2 (h, τ ), E ∞ (h, τ ) and the respective temporal convergence orders. The numerical results indicate that our method is of order (2γ ) in the temporal dimension, thus verifying the theoretical temporal convergence order in Theorem 5.1. (ii) In Table 2, we fix τ = 1/5000 and let h vary. We present E 2 (h, τ ), E ∞ (h, τ ) and the spatial convergence orders of our scheme (2.25)-(2.27) and those of the cubic non-polynomial spline (CNS) method [36]. From Table 2, we see that our method can achieve at least O(h 6 ), while the CNS method is O(h 4 ). Moreover, our method obtains smaller errors in all the cases. The observation also confirms the theoretical spatial convergence order in Theorem 5.1. (iii) Let τ = 1/512 be fixed. In Table 3, we compare E Re ∞ (h, τ ) and E Im ∞ (h, τ ) of our scheme (2.25)-(2.27) with those of the meshless collocation (MC) method [47] and the CNS method [36]. The numerical results indicate that our method gives the smallest errors in all the cases. (iv) In Table 4 [16] and the CNS method [36]. Once again, the numerical results show that our scheme outperforms these methods.
Let α = 1/20 in the implementation of the scheme (2.25)-(2.27). In Table 5, we fix τ = 1/5000 and present the spatial convergence orders of our scheme and those of the CNS method [36]. In Table 6, we fix τ = 1/512 and compare the maximum real/imaginary part absolute errors of our numerical scheme with those of the CNS method [36]. Furthermore, in Fig. 4 we plot the absolute modulus error |u n j -U n j | obtained from our scheme for γ = 0.3 and (h, τ ) = (1/20, 1/200).
From the numerical simulation and plot, once again we demonstrate that the theoretical spatial convergence order of our scheme is at least six (Theorem 5.1) and our scheme performs better than the CNS method.    In the next example, we shall investigate the effect of α on the actual maximum L2 norm error E 2 (h, τ ). + -4π 2 t γ +2t 2 + t γ +2t 2 3 sin(2πx) + i cos(2πx) . (6.11) The exact solution is u(x, t) = (t γ +2t 2 )[sin(2πx) + i cos(2πx)]. First, we let α = 1/20 and τ = 1/5000 and present the spatial convergence orders of our scheme (2.25)-(2.27) and those of the CNS method [36] (Table 7). To visualize the efficiency of our method, in Figs. 5-7 we plot the real/imaginary parts of the numerical and exact solutions and the absolute modulus error for γ = 0.5 and (h, τ ) = (1/16, 1/200). We observe that the theoretical spatial convergence of our method is indeed at least O(h 6 ) (Theorem 5.1) and our method gives a good approximation of the exact solution.

Conclusion
In this paper, we derive a numerical scheme to solve the time-fractional nonlinear Schrödinger equation (1.1) of fractional order γ ∈ (0, 1). Our tools include the quintic non-polynomial spline and L1 discretization. The unconditional stability, unique solvability and convergence are proved by the Fourier method. It is shown that our method can achieve sixth order convergence in space and (2γ )th order convergence in time.
Three numerical examples are presented to verify the theoretical results and to compare with other methods in the literature.