Fully decoupled schemes for the coupled Schrödinger-KdV system

The coupled numerical schemes are inefficient for the time-dependent coupled Schrodinger-KdV system. In this study, some splitting schemes are proposed for the system based on the operator splitting method and coordinate increment discrete gradient method. The schemes are decoupled, so that each of the variables can be solved separately at each time level. Ample numerical experiments are carried out to demonstrate the efficiency and accuracy of our schemes.


1.
Introduction. Many phenomena such as fluid mechanics, solid-state physics, plasma physics and chemical physics, can be described by nonlinear partial differential equations. To model nonlinear dynamics of one-dimensional Langmuir and ion-acoustic waves in a system of coordinates moving at the ion-acoustic speed, the time-dependent coupled Schrödinger-KdV (CS-KdV) system was proposed [3], where i 2 = −1, is a positive constant, and the complex function E and the real function N describe electric field of Langmuir oscillations and low-frequency density perturbation, respectively. Sometimes, the following coupled Schrödinger-KdV system is used for describing the other similar phenomena. In this paper, we mainly focus on developing numerical schemes for the system (1) with initial conditions E 0 (x) = E(x, 0), N 0 (x) = N (x, 0), x ∈ [x l , x r ] and periodic boundary conditions. Our approaches can be extended to the system (2) straightforwardly. The system (1) possesses three invariants which, respectively, correspond to the number of Langmuir plasmons, the number of particles and the energy of the oscillations.
The CS-KdV system has motivated a series of studies in mathematics in the last decade. Some travelling wave solutions have been presented in [12]. Recently, numerical methods for this problem have attracted much attention, including finite element methods [5], finite difference method [2] and the other methods [27,20,6,13,22,17]. Most of them are expensive since they are fully coupled and nonlinearly implicit.
By setting the complex function E(x, t) = p(x, t) + iq(x, t), the CS-KdV system (1) can be written as an infinite Hamiltonian forṁ where u = (p, N, q) , the dot denotes ∂ t , the skew-adjoint constant linear differential operator and the functional are, respectively, Many PDEs such as the KdV equation [4,7], Schrödinger equation [11] and Maxwell's equations [21,10] can be recast into the general Hamiltonian form (4). Because symplectic or multisymplectic integrators do not preserve the nonconstant symplectic or Poisson structures like (4), one can look for energy-preserving integrators. After a suitable spatial discretization of the skew-adjoint operator and Hamiltonian in (4), one obtains the following finite dimensional Hamiltonian systeṁ where u = (p 1 , · · · , p J , N 1 , · · · , N J , q 1 , · · · , q J ) , J represents the number of grid points, D is any consistent constant skew matrix approximation to D, H∆x is any consistent approximation to H with ∆x being an increment in space. In [27], the authors used the technique introduced in [11], i.e., the time direction is integrated by the average vector field (AVF) method [23], to develop an energy-preserving scheme where τ is an increment in time, for the CS-KdV system (1). However, the scheme is fully coupled and nonlinearly implicit. Other discrete gradient (DG) methods such as the coordinate increment DG method [18] and the midpoint DG method [15] can also be employed to integrate (5), but the obtained schemes are also fully coupled and nonlinearly implicit. Although the fully implicit schemes exhibit fascinating preservation of energy of the original system, they are lack of computational efficiency. Hence, it is necessary to design efficient and robust numerical schemes. A good choice for diminishing the computational cost of the evolution equations is the operator splitting (also called fractional step or time splitting method). There have been many studies on the operator splitting methods for the Schrödinger-type equations. For example, the nonlinear Schrödinger equation can be split into linear part iu t + u xx = 0 and nonlinear part iu t + |u| 2 u = 0, respectively. The property d|u| 2 /dt = 0 implies that the nonlinear part can be integrated exactly, so that the operator splitting method for (7) is highly efficient. By making use of the splitting technique, the system (1) can also be split into a linear part and a nonlinear part. The former part is easy to be solved numerically, but the latter part cannot be integrated exactly due to dN/dt = 0. Because the coupling of variables E and N makes the nonlinear part integrated inefficiently, it is desirable to develop decoupled schemes for the coupled system (1). In general, it is not easy to make a decomposition of the coupling term directly, but one may employ a specific discretization to develop decoupled schemes for the coupled systems [14,25,8,9].
The main purpose of this study is to develop some decoupled schemes for the CS-KdV system. To this end, the Hamiltonian form (4) is decomposed into to a series of linear and nonlinear Hamiltonian subsystems which are different to those in [27]. However, the splitting method does not result in the decomposition of the variables E and N . In this work, the first-order decoupled scheme is developed by integrating the nonlinear subsystems in time with a specific DG method. And then the composing technique is employed to construct the second-order schemes which are popular in practice, This paper is organized as follows. In Section 2, a brief summary of the operator method for the evolution equation is presented. And then, some coupled schemes are presented for the CS-KdV system. In Section 3, some first-and second-order decoupled schemes are proposed by combination of the operator splitting methods and the specific DG method. Many numerical experiments are carried out to exhibit the performance of the schemes in Section 4. This paper ends with concluding remarks in Section 5.

2.
Operator splitting method and coupled scheme. Most of the existing schemes for the CS-KdV system are time-consuming since they are fully coupled and implicit. The operator splitting method is a good idea to improve the computational efficiency.
2.1. Review of the operator splitting method. We first present a brief summary of the split-step method for a general evolution equation in the form where L and N are linear and nonlinear operators, respectively. In general, the operators L and N do not commute with each other. The main idea in the splitstep method is to approximate the exact solution of original equation by solving the obtained splitting equations in a given sequential order, in which the solution of one subproblem is employed as an initial condition for the next subproblem. This may be realized by using a solution operator ϕ(τ ) that includes an appropriate combination of products of the exponential operators exp(τ L) and exp(τ N ). This produces a splitting error due to the noncommutativity of operators L and N . Splitting errors have been analyzed by a number of researchers [19,1]. For an operator equation of the form (8), the first-order solution operator is given by In the second-order version of the method, which has also received much attention in the literature, the solution operator is approximated by which is symmetric, that is, ϕ 2 (τ )ϕ 2 (−τ ) = 1. In general, L and N in the above exponential operators may be interchanged without affecting the order of the method. Approximations of higher order, that preserve the symmetry, can be constructed by a proper composition of the second-order symmetric approximation.

2.2.
Nonlinearly implicit coupled schemes. In general, one can split the Hamiltonian system (4) where They are linear and nonlinear Hamiltonian systems holding H 1 and H 2 , respectively.
In [27], the authors proposed a second-order splitting pseudo-spectral scheme (AVFS) which is similar to S-AVF-2. Let D 1 and D 2 be the first-and second-order Fourier pseudo-spectral differential matrices, whose entries are, respectively, The scheme (14b) is fully coupled and nonlinearly implicit. Other DG methods can also be employed to integrate (13), but the obtained splitting schemes are still fully coupled and nonlinearly implicit.
3. Decoupled schemes. In this section, some decoupled schemes are developed for the CS-KdV system by combination of a new operator splitting method and a specific DG method.
3.1. First-order decoupled splitting scheme. Since the operator splitting (11) is unhelpful for developing the decoupled scheme, in what follows, a new operator splitting will be considered. We split the Hamiltonian system (4) into the following three subsystemsu where One has a first-order solution operator for the CS-KdV system (1), immediately. Next, we make a spatial discretization for each of the Hamiltonian subsystems. Let D be the same as that in (12) and The spatial semi-discrete systems for the subsystems are, respectively,

JIAXIANG CAI, JUAN CHEN AND BIN YANG
where u = (p , N , q ) . The systems hold Hamiltonians H 1 , H 2 and H 3 , respectively, so that it is natural to integrate them with a discrete DG method. The obtained full-discrete schemes are given by where∇ is a DG of H, that is, a continuous function of (û, u) satisfying∇H(û, u) . It is clear that the choice of the AVF method in (19) is unhelpful for diminishing the computational cost. To decompose the variables p, q and N, one can consider the coordinate increment DG method 1 in (19). Combining (19) and the solution operator (16) gives a first-order scheme for the CS-KdV system. It follows from the first and third equations in (20a) that where a = 3τ 4 . It is clear that p * , q * and N * can be calculated independently. Similarly, the variables p * * , q * * and N * * can also be decomposed. For example, one can eliminate the variables q * * and N * * from the first equation of (20b) to obtain The above arguments show that the splitting scheme (20) is decoupled. Note that the matrix D 1 is cyclic. Thanks to the discrete fast Fourier transform (DFFT), we have the relationship where diag(Λ) = Fd 1 with d 1 being the first column of D 1 and F is the matrix of DFFT coefficients 2 . The relationship (24) is very helpful to improve the computational efficiency of S-CI-1. 1 The coordinate increment DG method is∇I(ȳ, y) =        I(ȳ 1 ,y 2 ,·,y d )−I(y 1 ,y 2 ,·,y d ) y 1 −y 1 I(ȳ 1 ,ȳ 2 ,·,y d )−I(ȳ 1 ,y 2 ,·,y d ) y 2 −y 2 . . .

Numerical results.
To examine the performance of our decoupled splitting schemes we now consider several different problems described below. The purpose of the numerical experiments is to verify that S-CI-1, S-CI-2â, S-CI-2a, S-CI-2b and S-CI-2b numerically: (a) exhibit the expected convergence rates in time; (b) provide accurate soliton solution; (c) preserve the discrete invariants well.
The accuracy of the solution at time t = t n is assessed by the maximal error norm e n ∞,W = max j |W (x j , t n ) − W n j | and the average error norm e n 2,W = (∆x j |W (x j , t n ) − W n j | 2 ) 1/2 , W = E, N, p, q. The conservation of invariants are examined by calculating the relative errors RE(I i )=|I n i − I 0 i |/|I 0 i |, i = 1, 2, 3, where I n 1 , I n 2 and I n 3 are computed by ∆x j |E n j | 2 , ∆x j N n j and ∆x j (3|(D 1 E n ) j | 2 + N n j |E n j | 2 + (N n j ) 3 /3 − (D 1 N n ) 2 j ), respectively. In the following simulations, we consider the Fourier pseudo-spectral discretization in space.  Table 1 lists the average and maximal solution errors of our second-order splitting schemes and the second-order implicit schemes proposed in [27,20]. As it is seen from the table, the accuracies of solution given by the present schemes are comparable to those of the schemes in [27], but more accurate than that of the scheme in [20].
Secondly, we make some numerical comparisons between our second-order schemes and the second-order implicit ones in [22,17]. The maximal errors in solution are displayed in Table 2. One can see that S-CI-2b and S-CI-2b provide more accurate solution than the others. It is interesting that the splitting schemes derived from the solution operator (27) exhibit much better performance than those originated from the solution operator (26). The results indicate that the accuracies of solution of the splitting schemes are remarkably effected by the sequential order of the  [17] 9.41e-5 2.92e-5 Table 3. The maximal solution errors for CS-KdV system (1): x ∈ [−50, 50], ∆x = 0.1, τ = 0.0001 and T = 0.1.
Finally, we intend to compare the results of S-CI-2b and S-CI-2b with the forthorder Runge-Kutta (RK-4) method [17]. The computations are done with x ∈ [−50, 50], ∆x = 0.1, τ = 0.0001 and T = 0.1. The maximal solution errors are shown in Table 3. One can see that our second-order schemes give more accurate solution than the RK-4 scheme.
Long-term simulation: Consider the simulation with x ∈ [−50, 50], ∆x = 0.1, τ = 0.1 and T = 50. Here, we only illustrate the results of S-CI-2b with Fig. 1. It is clear that the scheme simulates the problem well. Fig. 2 shows the variation of the error in solution and the relative error in invariants. One can see that, throughout simulation, the errors in solution and the changes in I 1 and I 3 are controlled well, and the I 2 invariant is captured exactly.
Convergence orders of solution and invariants: Since the Fourier pseudospectral discretization is adopted in space, the convergence order in space is infinite as the solution of original system is sufficiently smooth. In the following, only the temporal convergence orders of the schemes will be tested numerically. We carry out the simulations on the region [−60, 60] × [0, 1] with various time steps. Fig. 3 (left) shows the maximal solution error versus the time step. The slopes of lines (red: 1; blue: 2) verify that S-CI-1 and S-CI-2b are of first-and second-order accuracies, respectively.
The present splitting schemes are not I 1 -and I 3 -preserving. Fig. 3 (right) plots the changes in I 1 and I 3 versus the time step. The results indicate that, as the time step becomes small, S-CI-2b has second-order convergence in I 1 and I 3 , and S-CI-1 has first-order convergence in I 3 while second-order convergence in I 1 .
CUP time: Consider the computational region [−60, 60] × [0, 10] and J = 240. Fig. 4 displays the CPU time versus the accuracy of the solution. One can see that our second-order decoupled splitting schemes are more efficient than S-CI-1 and the second-order coupled splitting AVFS [27] and S-AVF-2 schemes. Here we do not plot the CPU time of the second-order implicit coupled AVF scheme [27] since the scheme is much more expensive than AVFS (please see the results in [27]).

4.2.
Numerical results for the CS-KdV system (2). The system can also be rewritten as the Hamiltonian form (4) with u = (p, N, q) and possesses the same invariants I 1 and I 2 in (3) while the energy invariant I 3 = H. The methodology introduced in the previous section can be applied to this system analogously. The obtained splitting schemes are named as those in the previous section. For simplicity, we only test the performance of S-CI-2b. The CS-KdV system (2) admits an analytical solution where ξ = x + 2γt, θ = γx + (γ + 3γ 2 ± δ)t/3, δ = √ 2c 1 c 3 + 4c 2 2 with c 1 , c 3 and c 2 > 0 being arbitrary constants. In the following numerical tests, c 1 = c 3 = 0 and c 2 = 0.01 are used.
Next, we carry out the simulation on the region [−30, 30] × [0, 1] with parameters τ = 0.01, ∆x = 0.25 and γ = 0.1 or γ = 1. Figs. 5 and 6 indicate that the numerical solutions dovetail with the exact ones well for both cases γ = 0.1 and γ = 1. The accuracies of the solution and the preservation of the invariants for the case γ = 1 are illustrated with Fig. 7 (left).
The real and imaginary parts of the solution E are oscillatory for the large γ. The larger the γ is, the higher the frequency of oscillation is. In general, one should use a scheme which has a sufficient resolution in space to simulate the high oscillation case, however, it is quite time-consuming. Thanks to the decoupled feature of the scheme and the DFFT technique, our scheme deals with the highly oscillatory cases comfortably. We choose the initial value with γ = 10 and x ∈ [−60, 60], and the computational parameters τ = 0.001, ∆x = 0.1 and T = 1. Fig. 8 shows that the numerical solutions coincide with the exact ones. The errors in solution and the relative changes in invariants are plotted in Fig. 7 (right). The results indicate that our scheme is robust for the highly oscillatory case.
5. Concluding remarks. The CS-KdV system is a very important nonlinear dynamics model, which has attracted much attention. The system is difficult to solve numerically since the variables E and N are coupled. To our best knowledge, almost the schemes in the literature are coupled and implicit. These features lead that they are expensive in many simulations, especially in the simulation of highly oscillatory problem. In this study, we firstly develop a decoupled splitting scheme for the CS-KdV system by combination of the operator splitting method and coordinate increment DG method. And then, some second-order decoupled splitting schemes are proposed based on the composing method. The obtained schemes are  easy to use and highly efficient since each of the variables can be solved independently at each time level. Numerical results show the excellent performance of our schemes.