A Conservative Difference Scheme for Solving the Coupled Fractional Schrödinger–Boussinesq System

: In this paper, a high-accuracy conservative implicit algorithm for computing the space fractional coupled Schrödinger–Boussinesq system is constructed. Meanwhile, the conservative nature, a priori boundedness, and solvability of the numerical solution are presented. Then, the proposed algorithm is proved to be second-order convergence in temporal and fourth-order spatial convergence using the discrete energy method. Finally, some numerical experiments validate the effectiveness of the conservative algorithm and confirm the accuracy of the theoretical results for different choices of the fractional-order α .


Introduction
The coupled Schrödinger-Boussinesq system (CSBS) is a basic equation in laser and plasma physics.It describes how coupled Langmuir and dust-acoustic waves propagate nonlinearly in a dusty plasma.Over the years, the theoretical results on the well-posedness and dynamic behaviors of the analytical solution to the CSBS have been widely exploited [1][2][3][4][5][6].Several methods for finding the exact solutions to the CSBS have been presented [7][8][9][10].Because the exact solutions to the CSBS often contain certain special functions, scholars have begun utilizing efficient methods to seek numerical solutions for it, such as the multi-symplectic method [11], orthogonal spline collocation method [12], radial basis function-finite difference method [13], cut-off function method [14], scalar auxiliary variable method [15], Adams prediction-correction method [16], finite-element method [17,18], and energy-preserving compact finite difference methods [19].
This paper primarily focuses on the numerical solution for the fractional coupled Schrödinger-Boussinesq system (FCSBS) where 1 < α ≤ 2, i = √ −1 and coefficients γ and ω are positive constants.The function f (v) is sufficiently smooth and real, with f (0) = 0.The complex function u(x, t) represents the electric field in Langmuir oscillations, while the real function v(x, t) characterizes low-frequency density perturbations.
Assuming v t = ϕ xx , then system FCSBS (1) transforms into the following equivalent form: where Ω = (x L , x R ) and u 0 (x), v 0 (x), ϕ 0 (x) are given smooth functions.The fractional Laplacian can be considered as the Riesz fractional derivative where x L D α x u(x, t) represents the left Riemann-Liouville fractional derivative and x D α x R u(x, t) represents the right Riemann-Liouville fractional derivative It is noteworthy that the initial-boundary value problem (3)-( 7) possesses three conservative laws, namely, the Langmuir Plasmon number the total perturbed number density and the total energy Here, F(v) > 0 is the primitive function of f (v).
Han et al. [25] investigated the local and global well-posedness in H s (R), s ≥ 1 for FCSBS, in which the nonlinear term f (v) remains undetermined.Shao and Guo [34] derived local mild solutions to the Schrödinger-damped Boussinesq system and its fractional counterpart in one dimension using the contracting mapping principle.They also presented the precise results concerning the existence and nonexistence of global mild solutions.Given the inherent nonlocality and nonlinearity of FCSBS, obtaining analytical solutions for the system is an extremely challenging task.Therefore, numerical simulation has emerged as a crucial approach for its study.Ray [35] developed a time-splitting Fourier spectral method, which has been proven to be unconditionally stable.The error norms and graphical solutions are also presented in this work.In [36], Liao et al. developed and rigorously analyzed two efficient conservative difference schemes for FCSBS.Each scheme is demonstrated to preserve two fundamental conservation laws: mass conservation and energy conservation, while converging with an accuracy of O(τ 2 + h 2 ).Compared to CSBS, numerical methods for solving FCSBS are quite scarce.Thus, the goal of this paper is to develop a new conservative scheme for solving FCSBS, while also rigorously implementing error estimates for the proposed scheme.
The contributions of this article are summarized as follows: (1) To maintain the same physical properties as the original differential equation, we develop a conservative scheme.We rigorously demonstrate that the scheme can preserve three conservative laws simultaneously, as evidenced by numerical examples.(2) The boundedness and solvability of the numerical solutions obtained from the conservative scheme are established.
(3) Importantly, the numerical solutions provided by the conservative scheme unconditionally converges to the exact solutions in the L 2 -norm with a convergence order of O(τ 2 + h 4 ).
The subsequent sections of this paper are organized as follows: In Section 2, we introduce relevant notations and auxiliary lemmas, followed by the derivation of the conservative scheme.Theoretical analyses for the proposed scheme are presented in Section 3. Numerical experiments, conducted to validate the proposed scheme, are discussed in Section 4. Finally, concluding remarks are provided in Section 5.

Construction of Conservative Difference Scheme
In this section, we initially introduce relevant notations and auxiliary lemmas that will be utilized later.Subsequently, we elaborate on the establishment of a conservative difference scheme for the initial-boundary value problems (3)- (7).

Notations and Lemmas
Let h = (x R − x L )/J and τ = T/N denote the uniform step sizes in the spatial and temporal directions, respectively, where J and N are positive integers.Define For any grid function u = {u n j |u n j = u(x j , t n ), (x j , t n ) ∈ Ω h × Ω τ }, we define the following notations: Denote the space where F u(ω) = R e iωx u(x)dx is the Fourier transformation of u(x).We now introduce several lemmas crucial for constructing the conservative difference scheme.
Lemma 1 (see [22,37]).Suppose the function u ∈ L 4+α (R), 1 < α ≤ 2 and let be the fractional centered difference.Then, we have where and Z denotes the set of all integers.
Here, the A we define is consistent with the A α defined in [22,37], just with a simplified notation.
and suppose u * ∈ L 4+α (R).Then, for x ∈ (x L , x R ), we have The classical fourth-order compact approximation for standard second-order derivatives is obtained when α = 2 in Formula (11).

Derivation of the Conservative Difference Scheme
Let u n j = u(x j , t n ), ϕ n j = ϕ(x j , t n ), and v n j = v(x j , t n ).Meanwhile, U n j , Φ n j , and V n j represent the numerical approximations of u n j , ϕ n j , and v n j at the point (x j , t n ), respectively.By considering Equation (3) at both (x j , t n ) and (x j , t n+1 ), and then combining them, from Taylor expansion, we can derive Then, acting the operator A on both sides of (12), we have Using Lemma 1, we obtain Considering Equations ( 4) and ( 5) at (x j , t n ) and (x j , t n+1 ), respectively, and then combining them, from Taylor expansion, we can obtain Applying the operator B to both sides of Equations ( 15) and ( 16) and utilizing Lemma 2, we have Then, considering the discretizations ( 6) and ( 7), we have Neglecting the small terms in Equations ( 14), ( 17) and ( 18), and considering Equations ( 19) and ( 20), we can derive the following numerical scheme for (3)-( 7): Define ), , where P 1 and P 2 are square matrices with order J − 1. Denote it is easy to verify that G 1 and G 2 are symmetric positive definite matrixes.Thus, the vector forms of the fourth-order difference scheme ( 21)-( 25) and can be written as

Theoretical Analysis
Define the grid function spaces on h , the discrete inner product is defined as follows: For any two grid functions U, V ∈ V h , the discrete inner product and the associated l 2 h -norm are defined as follows: where Vj represents the complex conjugate of V j .We also define the discrete and the discrete maximum norm (l ∞ h -norm) as Provided the constant 0 ≤ σ ≤ 1, the fractional Sobolev norm ∥U∥ H σ and semi-norm |U| H σ can be defined as [39] ∥U∥ Obviously, we can obtain ∥U∥ 2 For convenience, we denote a general constant as C, which may vary across different contexts.Next, we introduce several useful auxiliary Lemmas.

Lemma 3 ([23]
).For any two grid functions U, V ∈ V h , a linear operator Λ α exists such that
Theorem 1.The scheme (26)-( 30) is conservative in the sense Proof.Computing the inner product of ( 26) with (U n+1 + U n ) and taking the imaginary part, we have By virtue of the first identity (32) of Lemma 4 and direct computation, we can deduce Then, substituting the above equalities into (34), one can obtain Making the inner product of ( 27) with Then, according to the first identity of Lemma 3, we obtain that By computing the inner product of ( 26) with 2τδ t U n and taking the real part, we have Calculating directly, we have According to (33) of Lemma (32), we obtain Thus, and using the relation By making the inner product of (28) with 2τδ t V n , we have Taking into account the relations Summing this equation and ( 35) and (36), we arrive at the formula Noting that

A Priori Bound
Lemma 5 ([23]).For any nonsingular matrix G and U ∈ V, two positive integers C 0 and C 1 exist such that where C 0 = min{λ j }, C 1 = max{λ j }, λ j is the singular value of the nonsingular matrix G.

Solvability
In this section, we discuss the solvability of the finite difference scheme ( 26)- (30).
, and z = (a T , b T , c T ) T ; then, z is a 3(J − 1)-dimensional vector or a point of 3(J − 1)-dimensional Euclidean space R 3(J−1) .Now, we use the Schauder fixed point to prove the existence of the solutions for the finite difference scheme ( 26)- (30).For this purpose, we construct a mapping T λ : R 3(J−1) −→ R 3(J−1) of the 3(J − 1)-dimensional Euclidean space into itself, with a parameter λ ∈ (0, 1) Obviously, the mapping T λ (z) defined here is continuous and there is a fixed point satisfying T 0 (z 0 ) = z 0 .Now, we prove the boundedness of all the possible solutions to the mapping.Making the inner product of ( 48) with (a + U n ) and taking the real part, by virtue of Lemma 4, we obtain and thus a is uniformly bounded.Now, we prove the boundedness of b and c.
Computing the inner product of ( 49) and (50) with γ(b + V n ) and (c + Φ n ), respectively, we obtain The addition of (52) to (53) yields Using Young's inequality Lemma 6 with p = q = 2 and ϵ = 1, we have Applying Taylor's Theorem, Young's inequality Lemma 6 with p = q = 2 and ϵ = 1, and Theorem 2, we obtain By substituting ( 55)-( 57) into (54), we have It follows from Theorem 2 that we have If τ is sufficiently small, we have According to the definition of the norm ∥ • ∥ l p h and (51), we can obtain the following estimates: Obviously, ( 51) and (60) imply that ∥a∥, ∥b∥, and ∥c∥ are uniformly bounded.Thus, ∥z∥ is uniformly bounded.It follows from the Schauder fixed-point theorem [43] that the conclusion of Theorem 3 holds.This completes the proof.

Convergence
In this subsection, we will first introduce two important Lemmas and then prove the convergence of the conservative scheme.
Lemma 11 ([43]).Suppose that g(x) Lemma 12 ([43]).Suppose that the discrete time sequence {w n |n = 0, 1, • • • , N; Nτ = T} satisfies the recurrence formula where A, B, and where τ is sufficiently small, such that From the construction of the scheme, we obtain the following result.
We present the error estimation of the conservative difference scheme ( 26)- (30) in the following theorem.Readers can refer to [23,24,32,33] for the same style of analysis method regarding the proof of convergence.
Theorem 5. Suppose that u(x, t), v(x, t), and ϕ(x, t) are the sufficiently smooth solutions to the problem (3)- (7).Then, the solutions U n j , V n j , and Φ n j of the scheme (26)-( 30) converges to the solutions u n j , v n j , and ϕ n j of the problem (3)-( 7) with order O(τ 2 + h 4 ).
When α = 2, f (v) = θv 2 , 3γ ̸ = θ, and 4γb 1 ̸ = d 1 , the FCSBS (3)-( 7) has the exact solitary wave solutions Here, and M, δ are free parameters.In the following simulations, the solutions at t = 0 are taken as the initial conditions, and the parameters are chosen as follows: Furthermore, we ensure that the computational domain x ∈ [x l , x r ] is sufficiently large to minimize errors introduced by the boundary conditions relative to the entire spatial domain.Firstly, we measure the accuracy of the proposed scheme by computing errors and convergence rates through Since exact solutions for α ∈ (1, 2) are not available, we obtain reference "exact" solutions U and V using the derived scheme with h = 1/80 and τ = 1/160.Additionally, we set −x l = x r = 40 and T = 2. Tables 1 and 2 demonstrate that the numerical scheme ( 26)- (30) achieves second-order accuracy in time and fourth-order accuracy in space for numerical solutions U n and V n with α = 2. Table 3 presents the l 2 h -norm errors and convergence orders for different α ∈ (1, 2).The table indicates that the scheme is second-order accurate in time.Next, we maintain a fixed time step of τ = 1/160 to assess spatial errors and convergence orders for varying α ∈ (1, 2), with the results summarized in Table 4.It is evident that the scheme achieves fourth-order accuracy in space.These convergence findings align with the theoretical expectations.Secondly, we compute the discrete conservation laws.For this test, we set h = 1/2, τ = 1/100, −x l = x r = 200, and T = 30, leaving α as the only free parameter.Tables 5-7 present the Langmuir plasmon number I n U , the total perturbed number density I n V , and the total energy I n E , respectively, where It is observed that schemes (26)- (30) effectively preserve these quantities, making them suitable for long-term simulation.Notably, the Langmuir plasmon number and the total perturbed number density remain unaffected by α, while the total energy varies selectively with α.To enhance conservation accuracy, a smaller iteration tolerance can be applied, albeit at the expense of increased computational cost.Next, we simulate the solitary wave solution with x ∈ [−100, 100], t ∈ [0, 30], h = 1/10, and τ = 1/2.Figures 1-3 depict the waveforms for |U| and V at varying α values.These figures demonstrate that variations in the order α directly impact the shapes of the solitons.As α decreases, small oscillations are observed near the solitary wave of |U| and V, with the amplitudes of these oscillations in V being larger than those in |U|.These characteristics mirror the numerical simulations of the space fractional Schrödinger system, which are utilized in physics to alter waveforms without altering nonlinearity and dispersion effects.Finally, we present a numerical comparison between our scheme ( 26)-( 30) and Scheme I from [36].Table 8 lists the errors and computational times for both schemes with different values of α.Clearly, our scheme provides a more accurate solution than Scheme I in [36], with only a slightly higher computational time cost.

Conclusions
In this paper, based on the compact difference approach, we derived a novel conservative numerical algorithm for solving the space fractional coupled Schrödinger-Boussinesq system.The conservative property, boundedness, solvability, and convergence of the numerical solution are evidenced.Finally, numerical experiments for different fractional-order α illustrated that the derived algorithm can guarantee conservation and convergent with order O(τ 2 + h 4 ).The results presented in this paper are applicable for numerical solutions to the classical coupled Schrödinger-Boussinesq system.In further work, we will try to discuss more complex and higher-dimensional fractional partial differential equations.

Figure 1 .
Figure 1.The wave forms of the numerical solution for |U| and V with α = 2.

Figure 2 .
Figure 2. The wave forms of the numerical solution for |U| and V with α = 1.5.

Figure 3 .
Figure 3.The wave forms of the numerical solution for |U| and V with α = 1.2.

Foundation
of the Department of Education of Hunan Province (2022JJ40021), and the Educational Department of Hunan Province of China (21B0722).

Table 5 .
The values of I n U at different times with τ = 1/100 and h = 1/2.

Table 6 .
The values of I n V at different times with τ = 1/100 and h = 1/2.

Table 7 .
The values of I n E at different times with τ = 1/100 and h = 1/2.