A Galerkin FEM for Riesz space-fractional CNLS

In quantum physics, fractional Schrödinger equation is of particular interest in the research of particles on stochastic fields modeled by the Lévy processes, which was derived by extending the Feynman path integral over the Brownian paths to a path integral over the trajectories of Lévy fights. In this work, a fully discrete finite element method (FEM) is developed for the Riesz space-fractional coupled nonlinear Schrödinger equations (CNLS), conjectured with a linearized Crank–Nicolson discretization. The error estimate and mass conservative property are discussed. It is showed that the proposed method is decoupled and convergent with optimal orders in L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$L^{2}$\end{document}-sense. Numerical examples are performed to support our theoretical results.


Introduction
An explosive interest has been gained to the fractional calculus and its applications during recent years. It has been known as an extension of classical calculus to non-integer orders, but compares favorably with classical calculus in modeling various applications, including image processing [1], coupled pendulums [2], capacitor microphone [3], optimal control [4], and anomalous diffusion [5,6].
In this paper, regarding the current interest in fractional PDEs as effective models, we showcase an efficient FEM for the Riesz space-fractional CNLS i∂ t u + γ ∂ 2α x u + λ |u| 2 + ρ|v| 2 u = 0, (1.1) i∂ t v + γ ∂ 2α x v + λ ρ|u| 2 + |v| 2 v = 0, t ∈ (0, T], (1.2) with i 2 = -1, real parameters γ , λ, ρ, 1/2 < α ≤ 1, and the Riesz derivative where Γ (·) denotes Euler's Gamma function. ∂ 2α x v can be analogously defined as ∂ 2α x u. Since the true solutions to Eqs. (1.1)-(1.2) subjected to suitable initial values will polynomially (α = 1) or exponentially (α = 1) decay to zero when |x| → ∞ [30], we consider the system truncated on a bounded domain Ω = (a, b) and take the initial and boundary conditions as follows: (1.4) where ∂Ω is its boundary, a 0, b 0, and u 0 (x), v 0 (x) are prescribed initial functions. In particular, it is noted that Eqs. (1.1)-(1.2) recover to the classical CNLS when α = 1, which was raised by Benney and Newell for two interacting wave trains in a weakly nonlinear conservative system [31]. The single space-fractional nonlinear Schrödinger equation (NLS) was first generated by Laskin [32,33], via extending the Feynman path integral over the Brownian paths to a new path integral over the Lévy quantum mechanical paths. Some important physical properties of self-similarity of the quantum mechanical path can be captured by this new approach. Besides, the fractional derivatives in NLS/CNLS can modify the shape of wave, while the nonlinearity and dispersion effects remain unchanged. Therefore, the spacefractional NLS/CNLS have been the topics of intense research recently. The existence and uniqueness of their global smooth solutions were investigated by [34][35][36]. The local wellposedness and ill-posedness in Sobolev spaces for power-type nonlinearities were studied in [37]. Moreover, there have been numerous works devoted to exploring the theoretical properties for those equations of different types; we refer the readers to [30,[38][39][40][41] for overall views. In the aspect of numerical algorithms for the space-fractional NLS/CNLS, Herzallah and Gepreel proposed an Adomian decomposition method to construct an approximate solution for the time-space-fractional NLS with cubic nonlinearity [42]. Atangana addressed the variable-order NLS of nonlocal type solved by the Crank-Nicolson method [43]. Amore et al. suggested a collocation approach based on little Sinc functions for the space-fractional NLS on a uniform grid [44], where good convergent behaviors are observed. In [45], the authors presented a second-order FDM for the Riesz spacefractional NLS and proved that it can conserve the discrete mass and energy. Klein et al. provided a detailed numerical study for the nonlocal NLS by Fourier spectral method [46]. Wang et al. established two mass conservative FDMs for the Riesz space-fractional CNLS [47,48], whose stability and convergence were discussed. Liu et al. solved the fractional NLS by the implicit and explicit-implicit FDMs, and the fractional CNLS by an implicit FDM [49]. Recently, a linearized compact alternating direction implicit method was constructed for the space-fractional NLS [50], by using a novel α-dependent fourth-order fractional compact difference operator, which turns out to be stable and convergent with second-order in time and fourth-order in space. The Legendre spectral and pseudospectral methods can be found in [51,52]. Zhu et al. developed an implicit FEM for spacefractional NLS [53]. Li et al. proposed a fast linearized conservative FEM for the strongly coupled space-fractional NLS [54], which uses Krylov subspace solvers with circulant preconditioners to handle the resulting Toeplitz-like systems.
Unlike the classical NLS, it is quite difficult to solve the fractional NLS analytically, let alone the space-fractional CNLS. Hence, the numerical methods are especially important. However, the existing methods are almost limited to FDM, only a few works on FEM have been done for the space-fractional NLS/CNLS. It is well known that FDMs are not suitable for complex regional problems and have a higher regularity requirement to guarantee their convergence in fractional cases, so do the spectral methods. For example, the scheme in [47] using fractional centered difference to deal with Riesz derivative requires u, v ∈ C 5 (R). As a common method, FEM compares favorably with FDM. It is widely used in the continuous medium and field problems on complex domains in engineering, and achieves higher accuracy with lower regularity requirements. The numerical simulation of the space-fractional NLS/CNLS by FEM helps to know the shape and long time behavior of the solitons governed by the fractional derivatives and their changes with different α. The main contribution of this article is to propose a fully discrete explicit-implicit FEM to solve the Riesz space-fractional CNLS on bounded domains, which is decoupled. We prove that it is conservative to mass and convergent with optimal orders in L 2 -norm. The proposed method is applied to simulate the propagation of solitons, and some meaningful physical results are found. Since the scheme is linearized, an extra iterative loop is avoided, which is beneficial to improving computing efficiency. The novelty lies in the construction of a decoupled and mass conservative FEM by decomposing u, v in Eqs. (1.1)-(1.2) into their real and imaginary parts, which can reach high accuracy with low computing burden and regularity.
The outline is organized as follows. In Sect. 2, we give some auxiliary results related to the fractional calculus, which will be used throughout the paper. Section 3 studies the weak problem along with nonlocal bilinear form. In Sect. 4, we construct the discrete scheme with FEM to discretize Eqs. (1.1)-(1.4) and prove that it is stable. In Sect. 5, the convergence is derived, and to evaluate the effectiveness of our method, numerical tests are covered in Sect. 6.

Auxiliary results
In the sequel, we employ Hereafter, (·, ·) will be consistently used to stand for the above inner product. In view of [18,19], for μ > 0, the left, right, and symmetric fractional semi-norms and norms are defined as , and J μ S (R), respectively. Also, let F[u] be the Fourier transform of a prescribed function u(x) ∈ L 1 (R), with variable ω, i.e., F[u] = R u(x)e -iωx dx. Then, the fractional Sobolev semi-norm and norm are defined by and H μ (R) is denoted as the closure of C ∞ (R) with regard to · H μ (R) .
. Then based on [26], there exist Observing that and via Plancherel's theorem, we immediately have In the same fashion, we obtain The proof is completed.
, and H μ (R) are equal with their equivalent semi-norms and norms.
, and H μ 0 (Ω) are equal with their equivalent semi-norms and norms.

Weak problem
Define the semi-norm and energy norm as follows: and decompose u(x), v(x) into their real and imaginary parts by then Eqs. (1.1)-(1.2) can be recast as the coupled system By using Lemma 2.1, we give the weak problem: seek u 1 , and the symmetry property.
Proof Due to Hölder's inequality, we have Considering the model only on Ω, for u, v ∈ H α 0 (Ω), we denote their extensions to R by zero beyond Ω to beũ,ṽ, respectively. Then, by the inequality we obtain where Further, for Plancherel's theorem, one gets which yields As for the coercivity, clearly, by Lemma 2.1, there holds Therefore, using Lemma 2.3 and the fractional Poincaré-Friedrichs inequalities, Λ(u, u) u 2 E is an immediate consequence, which ends its proof.

A fully discrete FEM
In this section, we propose a fully discrete FEM to approximate the solutions of Eqs. in the temporal direction, we introducē To perform FEM, a spatial mesh Ω h is given, where Π k-1 is a set of polynomials of degree at most k -1 on Ω j . To avoid iterative loop, we discretize Eqs. (1.1)-(1.4) by a linearized Crank-Nicolson scheme, where a standard extrapolation method is used to treat the nonlinear terms. An explicit-implicit FEM is derived, which reads: seek u n 1,h , u n 2,h , v n 1,h , v n 2,h ∈ V h , n = 1, 2, . . . , N , for any χ 1 1,h , subjected to the initial values As for the starting dataû where It is evident that Eqs. (4.1)-(4.10) are decoupled by marching in time alternately with u n h , v n h during the whole steps, which implies that at each temporal level, an array of linear system is only needed to be handled. This merit would be greatly beneficial to cutting the computing costs. In addition, the obtained discrete system can well preserve the mass in the discrete sense, given as follows.

Convergent analysis
In this section, we derive the convergent estimate for the linearized Crank-Nicolson FEM (4.1)-(4.10). We shall start by introducing a fractional orthogonal projection: R h : In view of [17,55], the following properties are admitted.
In what follows, we turn to the convergence of FEM (4.1)-(4.10) based on operator R h . Let u n , v n be the analytic solutions to Eqs. (1.1)-(1.4) at t = t n with their real and imaginary parts u n 1 , u n 2 , v n 1 , v n 2 respectively, i.e., u n = u n 1 + iu n 2 , v n = v n 1 + iv n 2 , and u n if ∂ 3 t f ∈ L ∞ (t n-1 , t n , L 2 ), ∂ 2 t f ∈ L ∞ (t n-1 , t n , H p ), and ∂ 2 t f ∈ L ∞ (t n-2 , t n , L 2 ). Let e n u = u nu n h = ρ n uθ n u , ρ n u = u n -R h u n , θ n u = u n h -R h u n , and u n ju n j,h = ρ n u jθ n u j ,ū Then, by the aid of mathematical induction, the convergent theorem is admitted.

Theorem 5.1 Assume that
and the same conditions hold for v 1

7)
and when α = 3/4, one has where G 0 = |u 0 1 | 2 + |u 0 2 | 2 + ρ|v 0 1 | 2 + ρ|v 0 2 | 2 , P 1 , P 2 are the truncated errors in time, which satisfy P 1 0,Ω , P 2 0,Ω O(τ ) by Taylor formula. Let ρ Setting χ 1 u,h = θ O h 2k , (5.14) As to n = 1, using the inverse inequality, one gets in the sense that û in which j,h , with j = 1, 2. In view of θ 0 u j = 0, adding these two equations together yields    Remarks 5.1 The restriction imposed on the temporal stepsize τ in Theorem 5.1 is just a need in theory. It is not necessary for us to comply with such a restriction in realistic implementation. Actually, it may be removed by the analytical technique derived in [56]; but it is beyond our scope to address this issue here.

Numerical experiments
In this part, numerical experiments are performed to support our theoretical results as well as gauge the actual performance of FEM (4.1)-(4.10). As stated earlier, the exact solutions to Eqs. (1.1)-(1.2) will polynomially or exponentially decay to zero when |x| tends to infinity, a large Ω will be truncated so that the wave functions are negligible outside of Ω and thereby (1.4) can successfully be enforced. In the tests below, we always set V h to be a piecewise linear space.
Here, γ = 1, λ = 2 are taken. Then, when α = 1, the soliton solution is Let Ω = (-20, 20); we simulate the propagation of a single soliton with fixed τ = 0.01, h = 0.1 and different α. Fig. 1 displays the evolutions of the amplitudes of single solitons for α = 0.6, 0.8, 0.9, and 1, respectively, where the solitons of different width and height are observed. In particular, when α = 1, it collapses to the classical one. As depicted in [47], the order α (α = 1) can create two turning points in the solitary wave solution; we see the same phenomena appear here from Fig. 2. Table 1 reports the convergent results at t = 1 when α = 1 by selecting τ = 0.1 * h. Table 2 lists the values of conserved quantity Q n u at Figure 2 The numerical solutions at different t for α = 0.9; red stars: t = 1, blue circles: t = 2, green triangles: t = 3  different T for α = 0.55, 0.85, and 1. It is drawn that the mass is well preserved throughout the process and our scheme is available for long-time simulation.
Example 6.2 Consider the following space-fractional CNLS: subjected to the initial values As α = 1 and ρ = 1, the model is known as the Manakov system, which is integrable, i.e., a collision without any reflection, transmission, trapping, and creation of new solitary waves would be produced, when two solitons intersect each other. Let γ = 1, x 0 = 5, and c = 3; we consider the double solitons collision and the impact may be caused by different fractional orders on the wave propagation within Ω = (-20, 20) under the parameters τ = h = 0.1. In Fig. 3, we display the interaction of two solitons and the evolution of the amplitude of each soliton when α = 0.7, 0.8, 0.9, and 1. In Fig. 4, we present the state solutions at t = 1 and t = 4 for α = 0.75. It is observed that the width and height of the solitons have been greatly changed when α = 1. More precisely, the smaller α we take, the bigger change will be created. The time point when the solitons intersect is also varying with α. Besides, the shapes of the solitons may not be retained after they collide with each other when α = 1. Table 3 reports the evolution of the quantity Q n u for α = 0.65, 0.75, and 1, respectively, by which we easily see that our scheme well conserves the initial mass, where the quantity Q n v is omitted since Q n v = Q n u .
At last, we compare the computing efficiency of the explicit-implicit FEM and the Crank-Nicolson FEM with Newton's iteration (termed "Newton Method"). Table 4 gives the detailed time costs for those two methods, from which we conclude that the explicitimplicit FEM is more efficient than the Newton Method.
Remarks 6.1 In the tests, Newton's procedure was terminated by reaching a solution with tolerant error 1.0E-12; all the algorithms were run by using Matlab R2012a on a Lenovo PC with Intel Pentium G2030 3.00 GHz CPU and 4 GB RAM.

Conclusion
In this research, a fully discrete mass conservative FEM is proposed to solve the Riesz space-fractional CNLS. Its convergence is discussed and the convergent orders are showed   to be optimal in L 2 -sense. Since the scheme is linearized and decoupled, an extra Newton iterative loop is naturally avoided, which can significantly improve the computing efficiency in practice. The codes are also utilized to simulate the wave propagation of solitons, and a detailed numerical investigation on the double solitons collision is done with different fractional orders. The test results indicate the soliton shapes are greatly changed with α, while the conserved quantities are well preserved. It is noteworthy that the problem discussed here can possibly be solved by other numerical methods like the lattice fractional equation methods in [28,29], operational matrix methods in [57,58], or the fourth-order nonstandard compact and sixth-order weighted essentially non-oscillatory (WENO) finite difference methods in [59,60], which can reduce computing effort by converting the original problem to a system of differential/algebraic equations, or achieving high-order convergence without spurious oscillations. In future works, we will use the Riesz fractional differences and operational matrix techniques to deal with Riesz derivative or apply the ideas of the high-order nonstandard compact and WENO schemes in the above citations to derive more efficient numerical methods for Eqs. (1.1)-(1.2). Although some difficulties will be encountered, it is meaningful to extend these works to the space-fractional NLS/CNLS. The extensions of these methods to the CNLS with both time-and spacefractional derivatives will also be considered.