Next Article in Journal
Hybrid Fuzzy C-Means Clustering Algorithm, Improving Solution Quality and Reducing Computational Complexity
Previous Article in Journal
On a Matrix Formulation of the Sequence of Bi-Periodic Fibonacci Numbers
Previous Article in Special Issue
Theoretical Investigation of Fractional Estimations in Liouville–Caputo Operators of Mixed Order with Applications
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Energy-Conserving Explicit Relaxed Runge–Kutta Methods for the Fractional Nonlinear Schrödinger Equation Based on Scalar Auxiliary Variable Approach

Department of Mathematics, Northeast Forestry University, Harbin 150040, China
*
Author to whom correspondence should be addressed.
Axioms 2024, 13(9), 591; https://doi.org/10.3390/axioms13090591
Submission received: 30 July 2024 / Revised: 27 August 2024 / Accepted: 28 August 2024 / Published: 30 August 2024
(This article belongs to the Special Issue Fractional Calculus - Theory and Applications II)

Abstract

:
In this paper, we present a novel explicit structure-preserving numerical method for solving nonlinear space-fractional Schrödinger equations based on the concept of the scalar auxiliary variable approach. Firstly, we convert the equations into an equivalent system through the introduction of a scalar variable. Subsequently, a semi-discrete energy-preserving scheme is developed by employing a fourth-order fractional difference operator to discretize the equivalent system in spatial direction, and obtain the fully discrete version by using an explicit relaxed Runge–Kutta method for temporal integration. The proposed method preserves the energy conservation property of the space-fractional nonlinear Schrödinger equation and achieves high accuracy. Numerical experiments are carried out to verify the structure-preserving qualities of the proposed method.

1. Introduction

Fractional differential equations have been effectively utilized in many fields, including chemistry, physics, biology, and hydrology [1]. The classical Schrödinger equation stands as a foundation of quantum mechanics. By extending quantum mechanics based on the path integral over Lévy-like quantum-mechanical paths, the fractional Schrödinger equations can be obtained. During the past decades, the space-fractional nonlinear Schrödinger (F-NLS) equation has attracted more and more scholars’ interest. Some mathematical properties of F-NLS equations, such as the well-posedness and energy preservation, as well as many physical properties and applications, have been given in the literature, see [2,3,4] for example.
However, due to the intricacy of nonlinear terms and the intangible functions included in analytical solutions, finding an exact solution is difficult. Thus, efficiently developing accurate, reliable, and easy-to-implement numerical methods for the F-NLS equation naturally becomes an urgent topic.
In this context, we consider the F-NLS equation as follows
i u x , t t + 1 2 α 2 u x , t + ρ u x , t 2 u x , t = 0 , ( x , t ) R d × ( 0 , T ] ,
subject to the initial conditions
u x , 0 = ψ x , x R d ,
and the periodic boundary condition
u x + δ e j , t = u x , t , ( x , t ) R d × ( 0 , T ] ,
where i = 1 is the imaginary unit, u ( x , t ) is a complex-valued wave function on t 0 , T , and x R d for d = { 1 , 2 , 3 } . Equation (1) is written in multi-dimensions. The parameter ρ describes the strength of short-range (or local) nonlinear interactions: when ρ is positive, the interactions are repulsive or defocusing, whereas, if ρ is negative, the interactions are attractive or focusing [5]. e j = ( 0 , , 0 , 1 , 0 , , 0 ) , j = 1 , 2 , , n denotes an orthonormal basis of R d , where δ > 0 signifies a periodic cycle. Generally, we can select a sufficiently large domain, D, so that the solution outside D becomes negligible because of the algebraic decay of u ( x , t ) in space. This domain, D, can be modeled as a periodic region with Dirichlet boundary conditions. Therefore, we can consider an approximation to the original problem’s solution as the solution of the F-NLS equation with periodic boundaries when T is fixed [6]. In this sense, we will, without loss of generality, focus solely on the problem within the periodic region D = ( a , b ) × × ( a , b ) R d . The fractional Laplacian ( Δ ) α / 2 with 1 < α 2 represents a pseudo-differential operator, | ξ | α , in the Fourier space, which is commonly used in the context of fractional calculus to generalize the concept of the Laplace operator to non-integer orders
Δ α 2 u x = F 1 ξ α F u ,
where F denotes the Fourier transform and F 1 is its inverse. It is well known that the F-NLS (1) satisfies the specific conserved physical quantities such as mass and energy conservation [7]
M ( t ) = R d u ( x , t ) 2 d x = M ( 0 ) . E ( t ) = 1 4 R d ( Δ ) α 4 u ( x , t ) 2 + ρ u ( x , t ) 4 d x = E ( 0 ) .
In recent years, there has been an increasing focus among authors on developing efficient numerical schemes for the F-NLS equation. From a numerical standpoint, conservative algorithms exhibit superior performance compared with non-conservative algorithms because the latter are prone to nonlinear blow-up phenomena. Therefore, maintaining the properties of the original differential equations (ODEs) invariant is essential for the accuracy and reliability of numerical simulations. For the case of α = 1 , Equation (1) simplifies to the classical nonlinear Schrödinger equation. Various methods have been proposed to solve the classical NLS, such as finite difference methods [8,9], finite element methods [10,11], and time-splitting spectral methods [12]. It is widely recognized that structure-preserving algorithms offer substantial advantages for long-term integrations. These algorithms typically produce higher-quality numerical solutions with larger time steps [13,14,15,16].
For the space-fractional NLS equation, several structure-preserving algorithms have been developed in the literature. For example, Fourier pseudo-spectral methods [14], compact difference methods [17], local discontinuous Galerkin methods [18] and collocation methods [19]. However, most of the resulting schemes are costly because they are fully implicit for all nonlinearities, i.e., they require a nonlinear system that is quite complex in practical computations. To effectively address the nonlinear functions of the problem efficiently, several techniques have been derived, such as the split-step method and the linearized method [20,21]. Recently, Liu et al. proposed in [22] a class of methods combining the invariant energy quadratization (IEQ) approach with the scalar auxiliary variable (SAV), which can effectively preserve similar energy–mass laws. In [23], a new relaxation-type method for the classical nonlinear Schrödinger equation was introduced. This scheme not only preserves both the discrete mass and energy but also can efficiently avoid costly numerical treatment of the nonlinearity. Here, energy conservation refers to a discrete analogue of continuous energy, called discrete energy. These numerical methods ensure the conservation of this discrete energy, which approximates the original energy of the continuous conservative system.
Using the relaxation-type method, the nonlinear parts can be substituted by a new variable, which can linearize the nonlinearities and save computing time. Very recently, in [4], two conservative relaxation schemes were proposed for the F-NLS equation, which adopt the exponential time difference scheme and the integral factor scheme in time, respectively. In [6], an energy-preserving relaxation method developed for a class of spatial F-NLS equations with periodic boundary conditions is introduced and thoroughly analyzed. This method effectively conserves both the discrete mass and energy of the F-NLS equation and attains second-order convergence in the temporal direction. Duo and Zhang in [5] introduced a relaxation Fourier spectral method for solving the F-NLS Equation (1) and proved that the scheme demonstrates conservation in both mass and energy. However, they did not consider the unconditional stability of the relaxation scheme. Motivated by these developments, based on the SAV approach and the explicit relaxed Runge–Kutta (RRK) method, an explicit conservative scheme for the F-NLS equation is obtained, which inherits the advantages and avoids their limitations.
The remainder of this article is structured as follows. In Section 2, we transform the F-NLS equation into an equivalent system by introducing a scalar auxiliary variable. A semi-discrete conservative system is given with the weighted and shifted Lubich difference (WSLD) method employed on the reconstructed system in Section 3. In Section 4, we introduce a time discretization of the equivalent system using the invariant conservative explicit RRK method. The numerical properties of the methods are demonstrated in Section 5.

2. SAV Approach for F-NLS Equation

Compared with various typical energy-preserving schemes, the Hamiltonian system exhibits properties that ensure energy conservation and system stability [24], which are crucial for the theoretical analysis and numerical computations of the preserving relaxation-type scheme [25]. To utilize these properties, we develop the SAV approach [26] to rewrite Function (1). Similarly, we can extend such an equivalent system and method to the 2D case.
Firstly, we analyze the one-dimensional fractional F-NLS Equation (1) for x D = a , b R and t 0 , T . In this case, the fractional Laplacian on the bounded interval, D, with periodic boundary conditions can be given by the Fourier series as follows
Δ α 2 u x , t = k Z v k α u ^ k e i v k x a ,
where v k = 2 k π b a , and the Fourier coefficient, u ^ k , is given by
u ^ k = 1 b a D u x , t e i v k x a d x .
We introduce two useful lemmas provided in [27] for the subsequent theoretical analysis.
Lemma 1
([27]). For two real periodic functions, φ and ψ, we have
D φ Δ α 2 ψ d x = D ψ Δ α 2 φ d x = D Δ α 4 ψ Δ α 4 φ d x .
Lemma 2
([27]). Assuming the function G [ ϕ ] is defined as
G ϕ = D g θ , ϕ θ , Δ α 4 ϕ θ d θ ,
where g is a smooth function on D, then the variational derivative of G ϕ is expressed as follows
δ G δ ϕ θ = g ϕ + Δ α 4 g Δ α 4 ϕ .

2.1. Hamiltonian System

By setting u = p + i q , we can transform the one-dimensional system (1) into the following real-valued equations
p t = 1 2 Δ α 2 q ρ p 2 + q 2 q , q t = 1 2 Δ α 2 p + ρ p 2 + q 2 p ,
with the periodic boundary conditions
p x , t = p x + ( b a ) , t , q x , t = q x + ( b a ) , t .
Then, from Theorem 2.1 in [27], by defining the fractional Hamiltonian mass and energy equations as
M = D ( p 2 + q 2 ) d x , E = 1 4 D ( Δ ) α 4 p 2 + ( Δ ) α 4 q 2 + ρ ( p 2 + q 2 ) 2 d x ,
where D = [ a , b ] is the one-dimensional interval, the system (4) with the periodic boundary conditions obeys the fractional Hamiltonian energy and mass conservation laws, i.e.,
d d t M = 0 , d d t E = 0 .
Lemma 3
([27]). The system (4) is an infinite-dimensional canonical Hamiltonian system
p t q t = 0 1 1 0 δ E / δ p δ E / δ q ,
where the energy functional, E , is given in (5).

2.2. Reformulation of the F-NLS Equation

Now, we consider a scalar variable as follows, with the idea of the SAV approach
r = r ( p , q , t ) = f ( p , q ) , 1 + C 0 , t [ 0 , T ] ,
where f ( p , q ) = ρ 4 p 2 + q 2 2 , · , · denotes the inner product, and C 0 is a sufficiently large constant that ensures r is well-posed. Hence, we obtain the equivalent system
p t = 1 2 Δ α 2 q B 1 p , q r , q t = 1 2 Δ α 2 p + B 2 p , q r , r t = f 1 p , q , p t + f 2 p , q , q t 2 f p , q , 1 + C 0 ,
where
f 1 = ρ p 2 + q 2 p , f 2 = ρ p 2 + q 2 q , B 1 p , q = f 1 p , q f p , q , 1 + C 0 , B 2 p , q = f 2 p , q f p , q , 1 + C 0 ,
and with the initial condition
p 0 = p 0 x , 0 = R e u 0 x , q 0 = q 0 x , 0 = I m u 0 x , r 0 = f p 0 , q 0 , 1 + C 0 .
Particularly, the SAV method ensures the unconditional stability of energy and high efficiency in numerical schemes, requiring solving only decoupled equations with constant coefficients at each time step [25]. For the equivalent system (7), it exhibits both energy and mass conservation properties by defining the modified energy and mass as
E = 1 4 D Δ α 4 p 2 + Δ α 4 q 2 d x + r 2 , M = D p 2 + q 2 d x .
Theorem 1.
The equivalent system (7) satisfies the modified conservation law.
Proof. 
By taking the inner products of the first two equations in (7) with q t , p t , respectively, and adding them together, we can deduce
E = 1 4 d d t D Δ α 4 p 2 + Δ α 4 q 2 d x + f 1 p , q r , q t + f 2 p , q r , p t f p , q , 1 + C 0 = 0 .
then multiplying the third equation of (7) by 2 r , we obtain
d d t r 2 = f 1 p , q r , q t + f 2 p , q r , p t f p , q , 1 + C 0 ,
by computing with (9) and (10), the modified energy conservation law can be deduced as
d d t 1 4 D Δ α 4 p 2 + Δ α 4 q 2 d x + r 2 = 0 .
Similarly, by computing the inner product of the first two equations in (7) with p and q, the conservation of mass
d d t D p 2 + q 2 d x = 0 ,
can be proofed. □

3. Spatial Discretization Scheme

In this section, we construct a high-order energy-preserving discrete scheme in space for solving the equivalent system (7). The fractional Laplacian α / 2 has become one of the most extensively investigated nonlocal operators in recent research [7]. For a finite interval in the 1D problem, the fraction Laplacian α is equivalent to the following Riesz fractional derivative
α u x = 2 α u x x 2 α ,
where the Riesz fractional derivative in terms of order β is defined as
β u x x β = k β D x β L + D R β x u x ,
and k β = 1 / 2 cos β π / 2 . Furthermore, the left and right Riemann–Liouville fractional derivatives of order β , denoted as D x β L and D R β x , respectively, are defined as follows
D x β L u x = 1 Γ 2 β d 2 d x 2 L x x ξ 1 β u ξ d ξ , D R β x u x = 1 Γ 2 β d 2 d x 2 x R ξ x 1 β u ξ d ξ .
Here, u x defined on the interval L , R and Γ · denotes the fractional Gamma function.
Nowadays, discrete formats for second-order spatial fractional derivatives can be divided into two main types: one approach combines the central difference scheme with piecewise linear polynomial approximation, while the other involves assembling Grünwald difference operators with different weights and displacements [28]. Additionally, ref. [29] introduces the WSLD operator and its specific explanation is as follows.
According to the Lubich operator [30], we obtain L-order approximations for the α derivative α 0 or integral α < 0 with the relevant coefficients of the generating function δ α ξ , defined in [30],
δ α ξ = i = 1 L 1 i 1 ξ i α .
However, this scheme is unstable in the temporal direction. In order to address this issue, taking L = 2 in (12), for all ξ 1 , we have
3 2 2 ξ + 1 2 ξ 2 α = 2 3 α n = 0 m = 0 1 n α n 1 3 m α m ξ n + m = k = 0 q k α ξ k .
Letting k = n + m , then
q k α = 3 2 α m = 0 k 3 m g m α g k m α .
According to the following recurrence formula, defined in [31]
g 0 k = 1 , g k α = 1 α + 1 k g k 1 α , k 1 ,
the approximate solution δ α ξ can be obtained, which is evidently effective in solving fractional partial differential equations [29].
Remark 1.
The coefficient g k α = 1 k α k is for the power series expansion of the generating function 1 ξ α .
Now, let us set x i = a + i h , m i m + N x , where N x > 0 is an integer, and h = b a / N x is the spatial step size. Here, we introduce the following parameters d , d ¯ , e , e ¯ , r , r ¯ , s , s ¯ , which will be described below, then the parameter m is defined as
m = max d , d ¯ , e , e ¯ , r , r ¯ , s , s ¯ .
Assuming that μ x i , t = 0 if i = m , m + 1 , 0 and i = N x , N x + 1 , , N x + m , then the fourth-order WSLD approximation of the Riesz fractional derivative can be given by
δ x α μ x i , t = 1 h α k = 1 x 1 w i k α μ x i , t , i = 1 , 2 , , N x 1 .
Here, we introduce the coefficient
φ k θ = w d e r s w d e w d l k + d m θ + w d e r s w d e w e l k + e m θ + w d e r s w r s w r l k + r m θ + w d e r s w r s w s l k + s m θ + w d ¯ e ¯ r ¯ s ¯ w d ¯ e ¯ w d ¯ l k + d ¯ m θ + w d ¯ e ¯ r ¯ s ¯ w d ¯ e ¯ w e ¯ l k + e ¯ m θ + w d ¯ e ¯ r ¯ s ¯ w r ¯ s ¯ w r ¯ l k + r ¯ m θ + w d ¯ e ¯ r ¯ s ¯ w r ¯ s ¯ w s ¯ l k + s ¯ m θ ,
where
w d = e e d , w e = d d e , d e ,
w r = s s r , w s = s s r , r s ,
w d e = 3 r s + 2 α 3 r s d e , w r s = 3 d e + 2 α 3 d e r s , d e r s ,
w d e r s = ι z ¯ ι z ¯ ι ¯ z , w ι ¯ e ¯ r ¯ s ¯ = ι ¯ z ι ¯ z ι z ¯ , ι ¯ z ι z ¯ ,
where ι = r s p q , ι ¯ = r ¯ s ¯ d ¯ e ¯ , and d , e , r , s , d ¯ , e ¯ , s ¯ , r ¯ are integers, and
z = 6 d e r s r + s d e + 4 α r s r + s d e d + e + 9 α r s d e ,
z ¯ = 6 d ¯ e ¯ r ¯ s ¯ r ¯ + s ¯ d ¯ e + 4 α r ¯ s ¯ r ¯ + s ¯ d ¯ e ¯ d ¯ + e ¯ + 9 α r ¯ s ¯ d ¯ e ¯ .
Then, the coefficient ϖ l α = ( φ m + l α + φ m l α ) / 2 cos ( α π / 2 ) .
Theorem 2
(Fourth order approximations [32]). Let 1 < α < 2 , if the fractional derivatives D x α a μ x i , t and D x α b μ x i , t of μ x i , t together with their Fourier transform belong to L a , b , then it holds that
α u x i , t x i α = δ x α μ x i , t + o h 4 ,
for any x i , 1 i N x 1 .
Let U i t = u x i , t and U t = U 1 t , U 2 t , , U N x 1 t , the approximation (15) can be reformulated as
δ x α U t = D α U t ,
where D α = D ˜ / h α , where matrix D ˜ is a Toeplitz matrix of order N x 1 , defined as
D ˜ = ω 0 α ω 1 α ω 1 N x α ω 2 N x α ω 1 α ω 0 α ω 1 α ω 1 N x α ω 1 α ω 0 α ω N x 1 α ω 1 α ω N x 2 α ω N x 1 α ω 1 α ω 0 α .
Remark 2.
To ensure the effectiveness of the approximation (15) for space-fractional derivatives, we select the parameters p , p ¯ , q , q ¯ , r , r ¯ , s , s ¯ , and m, as specified in Lemma 2.2 in [31] or Theorem 1.12 in [29,32] so that all eigenvalues of matrix D ˜ possess positive real parts.
Now, by denoting P t = p 1 t , p 2 t , , p N x 1 t T , where p i t = p x i , t , q i t = q x i , t , and defining Q, R as the same, the spatial semi-discrete form of the F-NLS Equation (7) can be formulated as follows
P t = 1 2 D α Q B 1 P , Q R , V t = 1 2 D α Q + B 2 P , Q R , R t = 1 2 B 1 P , Q , P t + B 2 P , Q , Q t .

4. Conservative Explicit SAV-RRK Method

Since we have rewritten the high-order energy function as a quadratic function by using the newly developed SAV approach, and expressed the F-NLS equation as an equivalent form of energy as a quadratic function [22], we will employ the RRK methods for time integration while using the quadratic invariance of the RRK methods to maintain the energy and mass conservation laws.
Firstly, we review the RRK methods [31,33,34] and discuss their structure-preserving characteristics. We introduce that y = ( P , Q , R ) T , y 0 = ( P 0 , Q 0 , R 0 ) T , and then the system described by (18) can be transformed as
y t = f ( y ) , t ( 0 , T ] , y ( 0 ) = y 0 ,
where, for convenience, we introduce the symbols
f = ( f 1 , f 2 , f 3 ) = 1 2 D α Q B 1 P , Q R , 1 2 D α Q + B 2 P , Q R , 1 2 B 1 P , Q , P t + B 2 P , Q , Q t .
Let M be a positive integer and 0 = t 0 < t 1 < t 2 < < t M = T with the step size τ = T / M , and let y m be the approximate value of y ( t m ) , then the s-stage explicit R K methods adopt the following form
Y m i = y m + τ j = 1 i 1 a i j f m j , i = 1 , , s , y m + 1 = y m + τ i = 1 s b i f m i ,
where f m j = f ( Y m j ) , j = 1 , , s .
Unfortunately, it is a well-established fact that only specific implicit Runge–Kutta methods can exhibit symplectic or algebraically stable properties, while no explicit Runge–Kutta methods possess these attributes. Thus, we consider the explicit RRK methods. For the one-step method in the interval [ t ^ m , t ^ m + 1 ] , m 0 , the s-stage RRK methods for (4), defined as
Y m i = y γ m + τ j = 1 i 1 a i j f m j , i = 1 , , s , y γ m + 1 = y γ m + γ m τ i = 1 s b i f m i ,
where y γ m = ( P γ m , Q γ m , R γ m ) T y ( t ^ m ) and γ m 0 , satisfy
E γ m + 1 = E γ m , i f i = 1 s b i f m i 0 , γ m = 1 , i f i = 1 s b i f m i = 0 ,
with E γ m defined as
E γ m = 1 4 D α 2 P γ m 2 + 1 4 D α 2 Q γ m 2 + R γ m 2 .
Like the RK methods, the RRK methods can also be expressed by using the Butcher tableau as follows
c 1 c 2 c s 0 0 0 a 21 0 0 a s 1 a s 2 0
γ m b 1 γ m b 2 γ m b s
The RRK methods are explicit, from which we can compute the relaxation parameter γ m explicitly, and this is the advantage of using explicit methods.
Then, from the definition in (22) and scheme in (20), we can obtain that the discrete energy satisfies
E γ m + 1 = 1 4 D α 2 P γ m + 1 2 + 1 4 D α 2 Q γ m + 1 2 + R γ m + 1 2 = 1 4 P γ m + γ m τ i = 1 s b i f m i 1 , D α P γ m + γ m τ i = 1 s b i f m i 1 1 4 Q γ m + γ m τ i = 1 s b i f m i 2 , D α Q γ m + γ m τ i = 1 s b i f m i 2 + R γ m + γ m τ i = 1 s b i f m i 3 2
with continuous calculation, we obtain
E γ m + 1 = 1 4 D α 2 P γ m 2 + 1 4 D α 2 Q γ m 2 + R γ m 2 1 4 [ P γ m , D α γ m τ i = 1 s b i f m i 1 + Q γ m , D α γ m τ i = 1 s b i f m i 2 + γ m τ i = 1 s b i f m i 1 , D α P γ m + γ m τ i = 1 s b i f m i 2 , D α Q γ m + D α i = 1 s b i f m i 1 , P γ m + D α i = 1 s b i f m i 2 , Q γ m 8 R γ m γ m τ i = 1 s b i f m i 3 ] + γ m τ i = 1 s b i f m i 3 2 + 1 4 D α 2 γ m τ i = 1 s b i f m i 1 2 + 1 4 D α 2 γ m τ i = 1 s b i f m i 2 2 = E γ m 1 4 γ m τ [ P γ m , D α i = 1 s b i f m i 1 + Q γ m , D α i = 1 s b i f m i 2 + i = 1 s b i f m i 1 , D α P γ m + i = 1 s b i f m i 2 , D α Q γ m 8 R γ m i = 1 s b i f m i 3 ] + 1 4 γ m 2 τ 2 4 i = 1 s b i f m i 3 2 + D α 2 i = 1 s b i f m i 1 2 + D α 2 i = 1 s b i f m i 2 2 .
Then, in order to ensure E γ m = E γ m + 1 , we have that γ m satisfies
γ m = P γ m , D α i = 1 s b i f m i 1 + Q γ m , D α i = 1 s b i f m i 2 + D α i = 1 s b i f m i 1 , P γ m + D α i = 1 s b i f m i 2 , Q γ m 8 R γ m i = 1 s b i f m i τ 4 i = 1 s b i f m i 3 2 + D α 2 i = 1 s b i f m i 1 2 + D α 2 i = 1 s b i f m i 2 ) 2
This implies that the value of γ m can be explicitly calculated by the above formula. Since lim τ 0 γ m = 1 , the explicit SAV-RRK methods (20) are valid. Since the relaxation coefficient, γ m , varies with each step, it can be estimated similarly at different steps. As τ 0 , the relaxation coefficient approaches 1, indicating that the relaxation methods can be viewed as a small perturbation of the standard methods. For simplicity, we assume that the relaxation coefficient, γ m , is computed exactly below. In [22], the implications of the estimates of the relaxation coefficient for the accuracy of the explicit RRK methods were investigated, and similar conclusions can be drawn for the present method.
What is more, it can be demonstrated that the explicit SAV-RRK methods have relative energy-conservation properties.
Theorem 3
(Energy conservation). The fully discrete scheme obtained by using fourth-order WSLD approximation and SAV-RRK methods whose order is at least two satisfies the relative energy conservation, i.e., E γ 0 = E γ m for m = 0 , 1 , , M , where E γ m is defined in (22).
Proof. 
When i = 1 s b i f m i = 0, the second equation in (20) implies y γ m + 1 = y γ m , satisfying this theorem naturally. In the case where i = 1 s b i f m i 0 , we also deduce it from the corresponding condition (21). □
Remark 3.
Actually, the energy conservation property of the WSLD-SAV-RRK method described in Theorem 3 should be relatively energy conservative, which indicates that the proposed schemes can preserve the original energy with high precision at the discrete level. The precision of energy conservation is correlated with the convergence order of the RRK method chosen in the present method and the round-off error generated by computing the relaxation coefficient, γ m .
For the convergence of the WSLD-SAV-RRK methods, due to the properties of the WSLD difference operator, it can be inferred that the WSLD-SAV-RRK methods can reach the fourth-order accuracy in the spatial direction, and the convergence order in the temporal direction is the same as that of RRK methods applied to ordinary differential equations. For the stability of the SAV-RRK method, we give the following theorem, which can be simply obtained by the properties of the RRK method without proof.
Theorem 4.
(Stability) For any SAV-WSLD-RRK methods, it is satisfies that p m C , q m C for any m = 0 , 1 , , M , where C is a positive constant independent of τ and h.

5. Numerical Experiments and Discussions

In this context, we will primarily demonstrate the algorithmic process, convergence, and energy conservation. We utilize the WSLD method for spatial discretization and combine it with the high-order approximation scheme of SAV-RRK. This approach is used to solve a certain class of one-dimensional and two-dimensional F-NLS equations through the following numerical examples.

5.1. One-Dimensional F-NLS Equation Problem

Firstly, we consider the 1D fractional problem
i u x , t t + 1 2 α 2 u x , t + ρ u x , t 2 u x , t = 0 , x D ,
with the initial value
u x , 0 = sech x · exp 2 i x ,
with the parameters ρ = 2 . If we take α = 2 and set the space interval D = ( , ) , the exact solution of Equation (24) is
u x , t = sech x 4 t · exp i 2 x 3 t .
Because the initial value given in the scheme in (25) tends to decay to zero exponentially as x moves away from the origin, the wave function can be identified as negligible when x outside [ 8 π , 8 π ] . Consequently, let u ( a , t ) = u ( b , t ) = 0 for a 0 and b 0 .
In Figure 1, the surfaces of the numerical solution | u ( x , t ) | computed by SAV-WSLD-RRK methods with τ = 1 / 2 9 and h = π / 2 10 are depicted for different values of order α . We observe from Figure 1 that the values of order α affect the shape of the soliton. What is more, in Figure 2, we plot the numerical solutions | u ( x , t ) | for Equation (24) for time t = 1 , 3 , 5 for α = 1.5 , 1.8 with τ = 1 / 2 9 and h = π / 2 10 .
It can be seen that there are two turning points. From Figure 1 and Figure 2, it can be observed that the shape of the soliton changes more rapidly as α decreases. As α approaches 2, the numerical solutions of the F-NLS Equation (24) converge to the exact solutions of the Integer-order Schrödinger equation.
To estimate the error and verify the convergence of the new methods, we define the discrete L 2 -norm error of the numerical solution for the 1D problem as follows
E r r o r h τ = j = 0 N x 1 h U j M u ( x j , T ) 2 ,
where U j M | 0 j N x 1 is the numerical solution with time grid τ = T / M and space grid h = ( b a ) / N x at time T, and u ( x j , T ) is the analytical solution for the F-NLS equations. However, most of the exact solutions of fractional differential equations are hard to obtain. For this reason, as in [35,36,37], the errors in the temporal and spatial directions can be computed by
E r r o r τ = j = 0 N x 1 h U j M U j 2 M 2 , E r r o r h = j = 0 N x 1 h U j M U 2 j M 2 ,
with sufficiently small h and τ , respectively.
Then, the experimental order of convergence (EOC) can be calculated by
E O C τ = log 2 ( E r r o r τ 1 / E r r o r τ 2 ) log 2 ( τ 1 / τ 2 ) , E O C h = log 2 ( E r r o r h 1 / E r r o r h 2 ) log 2 ( h 1 / h 2 ) ,
if h and τ are sufficiently small, respectively.
In Table 1, we verify the convergence orders in the temporal and spatial directions for the present method. Table 1 shows the errors at T = 1 and the EOCs in temporal and spatial directions of Equation (24) with α = 1.5 and 2. We find that the present method, when using the second-order RRK method, exhibits an approximate second-order EOC in the temporal direction and is close to fourth-order in the spatial direction.
Next, we verify the energy conservation property. To this end, we define the discrete energy error, E r r o r E ( t j ) = E γ j E γ 0 / E γ 0 , where E γ j is defined in (22).
Figure 3 displays the discrete energy errors for Equation (24) with different fractional orders, α , for time t [ 0 , 6 ] . It can be seen from Figure 3 that the proposed scheme preserves the energy at a discrete level, which coincides with Theorem 3 and Remark 3, and the energy error reaches the order of 10 3 , which has a good conservation performance.

5.2. Two-Dimensional F-NLS Equation Problems

Firstly, we consider the following two-dimensional F-NLS equation problem [38]
i u t + 1 2 x α 2 u + y α 2 u + ρ u 2 u = 0 , x , y D = 20 , 20 2 , t 0 , T u x , y , t = 0 , x , y D , t 0 , T , u x , y , 0 = 2 π exp x 2 + y 2 , x , y D D .
The numerical solutions | u | of problem (26) with α = 1.3 , 1.7 , 1.9 and ρ = 2 at time t = 0.3 and t = 0.8 obtained by the present numerical method are shown in Figure 4. It can be seen from Figure 4 that the wave peak lies at the origin of the coordinate plane, and that the amplitude of the wave will diminish with the extension of time. Moreover, from Figure 4, it also can be seen that the wave disperses faster for α = 1.7 , 1.9 , than that for α = 1.3 , i.e., the larger α is, the greater is the dispersion effect on the shape of the wave.
Secondly, we consider the following two-dimensional problem [39]
i u t + 1 2 x α 2 u + y α 2 u + ρ u 2 u = 0 , x , y D = 0 , 2 π 2 , t 0 , T , u x , y , 0 = 1 + sin x 2 + sin y , x , y D D ,
with periodic boundary conditions on D .
The numerical solutions | u | of problem (27) with α = 1.3 , 1.4 , 1.5 , 1.8 , 1.9 and ρ = 2 at time t = 0.108 obtained by the present numerical method are shown in Figure 5. It can be seen from Figure 5 that there is a singularity near the boundary occurring at t = 0.108 . Moreover, from Figure 5, it also can be seen that the singularity becomes weaker as α becomes smaller.

6. Conclusions

This paper presented an energy-preserving method for space F-NLS equations based on the SAV approach and RRK methods, in which the spatial direction is discretized by a fourth-order approximation. The energy conservation of the novel method is proven. Finally, the effectiveness and the conservation of the numerical method are verified by several numerical experiments in both one- and two-dimensional cases.

Author Contributions

Conceptualization, Y.L.; methodology, Y.L. and Y.Z.; software, J.Z. and Y.C.; validation, Y.Z. and J.Z.; formal analysis, Y.C. and Y.Z; investigation, Y.L. and Y.C.; writing—original draft preparation, Y.Z.; writing—review and editing, J.Z. and Y.Z.; funding acquisition, Y.Z. and Y.L.; visualization, Y.C. All authors have read and agreed to the published version of the manuscript.

Funding

This work is supported by the College Students Innovations Special Project funded by Northeast Forestry University of China (No. DC-2023166) and the Natural Science Foundation of Heilongjiang Province of China (Grant Number: LH2022A002).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All the data were computed using our algorithm.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:
F-NLSFractional nonlinear Schrödinger
SAVScalar auxiliary variable
WSLDWeighted and shifted Lubich difference
RRKRelaxed Runge–Kutta

References

  1. Baleanu, D.; Karaca, Y.; Vázquez, L.; Macías-Díaz, J.E. Advanced fractional calculus, differential equations and neural networks: Analysis, modeling and numerical computations. Phys. Scr. 2023, 98, 110201. [Google Scholar] [CrossRef]
  2. Guo, X.; Xu, M. Some physical applications of fractional Schrödinger equation. J. Math. Phys. 2006, 47, 082104. [Google Scholar] [CrossRef]
  3. Muhammad, J.; Younas, U.; Nasreen, N.; Khan, A.; Abdeljawad, T. Multicomponent nonlinear fractional Schrödinger equation: On the study of optical wave propagation in the fiber optics. Partial. Differ. Equ. Appl. Math. 2024, 11, 100805. [Google Scholar] [CrossRef]
  4. Xu, Z.; Fu, Y. Two novel conservative exponential relaxation methods for the space-fractional nonlinear Schrödinger equation. Comput. Math. Appl. 2023, 142, 97–106. [Google Scholar] [CrossRef]
  5. Duo, S.; Zhang, Y. Mass-conservative Fourier spectral methods for solving the fractional nonlinear Schrödinger equation. Comput. Math. Appl. 2016, 71, 2257–2271. [Google Scholar] [CrossRef]
  6. Cheng, B.; Wang, D.; Yang, W. Energy preserving relaxation method for space-fractional nonlinear Schrödinger equation. Appl. Numer. Math. 2020, 152, 480–498. [Google Scholar] [CrossRef]
  7. Lischke, A.; Pang, G.; Gulian, M.; Song, F.; Glusa, C.; Zheng, X.; Mao, Z.; Cai, W.; Meerschaert, M.M.; Ainsworth, M.; et al. What is the fractional Laplacian? A comparative review with new results. J. Comput. Phys. 2020, 404, 109009. [Google Scholar] [CrossRef]
  8. Simos, T.; Williams, P. A finite-difference method for the numerical solution of the Schrödinger equation. J. Comput. Appl. Math. 1997, 79, 189–205. [Google Scholar] [CrossRef]
  9. Bao, W.; Cai, Y. Uniform error estimates of finite difference methods for the nonlinear Schrödinger equation with wave operator. SIAM J. Numer. Anal. 2012, 50, 492–521. [Google Scholar] [CrossRef]
  10. Karakashian, O.; Makridakis, C. A space-time finite element method for the nonlinear Schrödinger equation: The continuous Galerkin method. SIAM J. Numer. Anal. 1999, 36, 1779–1807. [Google Scholar] [CrossRef]
  11. Zhang, H.; Shi, D.; Li, Q. Nonconforming finite element method for a generalized nonlinear Schrödinger equation. Appl. Math. Comput. 2020, 377, 125141. [Google Scholar] [CrossRef]
  12. Thalhammer, M. Convergence analysis of high-order time-splitting pseudospectral methods for nonlinear Schrödinger equations. SIAM J. Numer. Anal. 2012, 50, 3231–3258. [Google Scholar] [CrossRef]
  13. Macías-Díaz, J.E. A structure-preserving method for a class of nonlinear dissipative wave equations with Riesz space-fractional derivatives. J. Comput. Phys. 2017, 351, 40–58. [Google Scholar] [CrossRef]
  14. Wang, P.; Huang, C. Structure-preserving numerical methods for the fractional Schrödinger equation. Appl. Numer. Math. 2018, 129, 137–158. [Google Scholar] [CrossRef]
  15. Macías-Díaz, J.E. An explicit dissipation-preserving method for Riesz space-fractional nonlinear wave equations in multiple dimensions. Commun. Nonlinear Sci. Numer. Simul. 2018, 59, 67–87. [Google Scholar] [CrossRef]
  16. Bai, G.; Hu, J.; Li, B. High-Order mass-and energy-conserving methods for the nonlinear Schrödinger equation. SIAM J. Sci. Comput. 2024, 46, A1026–A1046. [Google Scholar] [CrossRef]
  17. Hu, D.; Gong, Y.; Wang, Y. On convergence of a structure preserving difference scheme for two-dimensional space-fractional nonlinear Schrödinger equation and its fast implementation. Comput. Math. Appl. 2021, 98, 10–23. [Google Scholar] [CrossRef]
  18. Castillo, P.; Gómez, S. On the conservation of fractional nonlinear Schrödinger equation’s invariants by the local discontinuous Galerkin method. J. Sci. Comput. 2018, 77, 1444–1467. [Google Scholar] [CrossRef]
  19. Bhrawy, A.H.; Zaky, M. An improved collocation method for multi-dimensional space–time variable-order fractional Schrödinger equations. Appl. Numer. Math. 2017, 111, 197–218. [Google Scholar] [CrossRef]
  20. Zhai, S.; Wang, D.; Weng, Z.; Zhao, X. Error analysis and numerical simulations of Strang splitting method for space fractional nonlinear Schrödinger equation. J. Sci. Comput. 2019, 81, 965–989. [Google Scholar] [CrossRef]
  21. Wang, D.; Xiao, A.; Yang, W. A linearly implicit conservative difference scheme for the space fractional coupled nonlinear Schrödinger equations. J. Comput. Phys. 2014, 272, 644–655. [Google Scholar] [CrossRef]
  22. Liu, Y.; Ran, M. Arbitrarily high-order explicit energy-conserving methods for the generalized nonlinear fractional Schrödinger wave equations. Math. Comput. Simul. 2024, 216, 126–144. [Google Scholar] [CrossRef]
  23. Besse, C. A relaxation scheme for the nonlinear Schrödinger equation. SIAM J. Numer. Anal. 2004, 42, 934–952. [Google Scholar] [CrossRef]
  24. McLachlan, R. Featured review: Geometric numerical integration: Structure-preserving algorithms for ordinary differential equations. SIAM Rev. 2003, 45, 817–821. [Google Scholar]
  25. Leimkuhler, B.; Reich, S. Simulating Hamiltonian Dynamics; Number 14; Cambridge University Press: Cambridge, UK, 2004. [Google Scholar]
  26. Li, D.; Li, X.; Zhang, Z. Linearly implicit and high-order energy-preserving relaxation schemes for highly oscillatory Hamiltonian systems. J. Comput. Phys. 2023, 477, 111925. [Google Scholar] [CrossRef]
  27. Fu, Y.; Hu, D.; Wang, Y. High-order structure-preserving algorithms for the multi-dimensional fractional nonlinear Schrödinger equation based on the SAV approach. Math. Comput. Simul. 2021, 185, 238–255. [Google Scholar] [CrossRef]
  28. Sousa, E.; Li, C. A weighted finite difference method for the fractional diffusion equation based on the Riemann–Liouville derivative. Appl. Numer. Math. 2015, 90, 22–37. [Google Scholar] [CrossRef]
  29. Chen, M.; Deng, W. Fourth order accurate scheme for the space fractional diffusion equations. SIAM J. Numer. Anal. 2014, 52, 1418–1438. [Google Scholar] [CrossRef]
  30. Lubich, C. Discretized fractional calculus. SIAM J. Math. Anal. 1986, 17, 704–719. [Google Scholar] [CrossRef]
  31. Zhao, J.; Li, Y.; Xu, Y. An explicit fourth-order energy-preserving scheme for Riesz space fractional nonlinear wave equations. Appl. Math. Comput. 2019, 351, 124–138. [Google Scholar] [CrossRef]
  32. Li, Y.; Shan, W.; Zhang, Y. High-Order Dissipation-Preserving Methods for Nonlinear Fractional Generalized Wave Equations. Fractal Fract. 2022, 6, 264. [Google Scholar] [CrossRef]
  33. Ketcheson, D.I. Relaxation Runge–Kutta methods: Conservation and stability for inner-product norms. SIAM J. Numer. Anal. 2019, 57, 2850–2870. [Google Scholar] [CrossRef]
  34. Ranocha, H.; Sayyari, M.; Dalcin, L.; Parsani, M.; Ketcheson, D.I. Relaxation Runge–Kutta methods: Fully discrete explicit entropy-stable schemes for the compressible Euler and Navier–Stokes equations. SIAM J. Sci. Comput. 2020, 42, A612–A638. [Google Scholar] [CrossRef]
  35. Zhao, X.; Sun, Z.Z.; Hao, Z.P. A fourth-order compact ADI scheme for two-dimensional nonlinear space fractional Schrodinger equation. SIAM J. Sci. Comput. 2014, 36, A2865–A2886. [Google Scholar] [CrossRef]
  36. Weng, Z.; Zhuang, Q. Numerical approximation of the conservative Allen–Cahn equation by operator splitting method. Math. Methods Appl. Sci. 2017, 40, 4462–4480. [Google Scholar] [CrossRef]
  37. Willoughby, M.R. High-Order Time-Adaptive Numerical Methods for the Allen-Cahn and Cahn-Hilliard Equations. Ph.D. Thesis, University of British Columbia, Vancouver, BC, Canada, 2011. [Google Scholar]
  38. Wang, P.; Huang, C. Split-step alternating direction implicit difference scheme for the fractional Schrödinger equation in two dimensions. Comput. Math. Appl. 2016, 71, 1114–1128. [Google Scholar] [CrossRef]
  39. Liang, X.; Khaliq, A.Q. An efficient Fourier spectral exponential time differencing method for the space-fractional nonlinear Schrödinger equations. Comput. Math. Appl. 2018, 75, 4438–4457. [Google Scholar] [CrossRef]
Figure 1. Surfaces of the numerical solutions | u | of Equation (24) for different α with τ = 1 / 2 9 and h = π / 2 10 . (a) Numerical solution for α = 1.1 ; (b) numerical solution for α = 1.2 ; (c) numerical solution for α = 1.3 ; (d) numerical solution for α = 1.5 ; (e) numerical solution for α = 1.7 ; (f) numerical solution for α = 2 .
Figure 1. Surfaces of the numerical solutions | u | of Equation (24) for different α with τ = 1 / 2 9 and h = π / 2 10 . (a) Numerical solution for α = 1.1 ; (b) numerical solution for α = 1.2 ; (c) numerical solution for α = 1.3 ; (d) numerical solution for α = 1.5 ; (e) numerical solution for α = 1.7 ; (f) numerical solution for α = 2 .
Axioms 13 00591 g001
Figure 2. Numerical solutions for Equation (24) for time t = 1 , 3 , 5 for different α with τ = 1 / 2 9 and h = π / 2 10 . (a) Numerical solution for time t = 1 , 3 , 5 for α = 1.5 ; (b) numerical solution for time t = 1 , 3 , 5 for α = 1.8 .
Figure 2. Numerical solutions for Equation (24) for time t = 1 , 3 , 5 for different α with τ = 1 / 2 9 and h = π / 2 10 . (a) Numerical solution for time t = 1 , 3 , 5 for α = 1.5 ; (b) numerical solution for time t = 1 , 3 , 5 for α = 1.8 .
Axioms 13 00591 g002
Figure 3. Energy errors for Equation (24) for time t [ 0 , 6 ] for different α with τ = 1 / 2 9 and h = π / 2 10 . (a) Energy error for α = 1.3 ; (b) energy error for α = 1.5 ; (c) energy error for α = 2 .
Figure 3. Energy errors for Equation (24) for time t [ 0 , 6 ] for different α with τ = 1 / 2 9 and h = π / 2 10 . (a) Energy error for α = 1.3 ; (b) energy error for α = 1.5 ; (c) energy error for α = 2 .
Axioms 13 00591 g003
Figure 4. Numerical solutions | u | of Equation (26) for different α at t = 0.3 and t = 0.8 using the present method with τ = 1 / 100 and h = 1 / 2 5 .
Figure 4. Numerical solutions | u | of Equation (26) for different α at t = 0.3 and t = 0.8 using the present method with τ = 1 / 100 and h = 1 / 2 5 .
Axioms 13 00591 g004
Figure 5. Numerical solutions | u | of (27) for different α at t = 0.108 using the present method with τ = 1 / 100 and h = 1 / 2 5 .
Figure 5. Numerical solutions | u | of (27) for different α at t = 0.108 using the present method with τ = 1 / 100 and h = 1 / 2 5 .
Axioms 13 00591 g005
Table 1. The errors in the temporal and spatial directions and convergence orders of Equation (24) using the proposed method.
Table 1. The errors in the temporal and spatial directions and convergence orders of Equation (24) using the proposed method.
α τ h E r r o r h E O C h h τ E r r o r τ E O C τ
1.5 1 / 2 10 π / 2 6 1.1991 × 10 1 π / 2 9 1 / 2 5 5.2266 × 10 2
π / 2 7 7.9256 × 10 3 3.9192 1 / 2 6 1.2603 × 10 2 2.0521
π / 2 8 8.5764 × 10 4 3.2081 1 / 2 7 5.8514 × 10 3 1.1069
π / 2 9 5.3900 × 10 5 3.9920 1 / 2 8 1.8240 × 10 3 1.6817
2 1 / 10 3 π / 2 6 4.0759 × 10 2 π / 2 10 1 / 2 5 2.2878 × 10 2
π / 2 7 3.2603 × 10 3 3.6440 1 / 2 6 7.3454 × 10 3 1.6390
π / 2 8 1.6732 × 10 4 4.2843 1 / 2 7 2.9819 × 10 3 1.3006
π / 2 9 1.0642 × 10 5 3.9748 1 / 2 8 1.0803 × 10 3 1.4648
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Zhao, Y.; Li, Y.; Zhu, J.; Cao, Y. Energy-Conserving Explicit Relaxed Runge–Kutta Methods for the Fractional Nonlinear Schrödinger Equation Based on Scalar Auxiliary Variable Approach. Axioms 2024, 13, 591. https://doi.org/10.3390/axioms13090591

AMA Style

Zhao Y, Li Y, Zhu J, Cao Y. Energy-Conserving Explicit Relaxed Runge–Kutta Methods for the Fractional Nonlinear Schrödinger Equation Based on Scalar Auxiliary Variable Approach. Axioms. 2024; 13(9):591. https://doi.org/10.3390/axioms13090591

Chicago/Turabian Style

Zhao, Yizhuo, Yu Li, Jiaxin Zhu, and Yang Cao. 2024. "Energy-Conserving Explicit Relaxed Runge–Kutta Methods for the Fractional Nonlinear Schrödinger Equation Based on Scalar Auxiliary Variable Approach" Axioms 13, no. 9: 591. https://doi.org/10.3390/axioms13090591

APA Style

Zhao, Y., Li, Y., Zhu, J., & Cao, Y. (2024). Energy-Conserving Explicit Relaxed Runge–Kutta Methods for the Fractional Nonlinear Schrödinger Equation Based on Scalar Auxiliary Variable Approach. Axioms, 13(9), 591. https://doi.org/10.3390/axioms13090591

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop