A general conservative eighth-order compact finite difference scheme for the coupled Schrödinger-KdV equations

1.   Introduction

Nonlinear partial differential equations play an important part in various branches of science such as fluid mechanics, solid state physics, plasma physics and quantum mechanics. The coupled Schrödinger-KdV equations are put forward to model nonlinear dynamics of one-dimensional Langmuir and ion-acoustic waves in a system of coordinates moving at the ion-acoustic speed ^[18,19]. In detail, we consider the system ^[9]

where $i = \sqrt{-1}$, $m$ is a positive integer, $p, q, s, \epsilon, \alpha, \beta, \rho$ are real constants with $p\neq0$ and $\epsilon, \alpha \ge 0$. The complex-valued function $u$ and the real-valued function $v$ describe electric field of Langmuir oscillations and low-frequency density perturbation, respectively. The initial functions $\varphi$ and $\phi$ are given $l$-periodic functions. Hence, it suffices to take a single period $[0, l]$. For Eqs (1.1) and (1.2), there are four physical invariants to be considered:

The number of plasmons

The number of particles

The energy of oscillations

and the momentum

According to ^[2], these invariants may connect closely to accurate behaviors in time. Extensive numerical studies have been presented for the coupled Schödinger-KdV equations in the last decade, such as the finite element method ^[3], radial basis function (RBF) collocation method ^[4], decomposition ^[5], variational iteration ^[6], exponential time differencing three-layer implicit scheme(ETDT-P) ^[7], homotopy perturbation ^[8], and fourth-order conservative compact finite difference scheme ^[9] and so on.

In the aspect of compact difference scheme, which is well known for the narrower stencils, i.e., fewer neighboring nodes it uses, and have less truncation error comparing with typical finite difference schemes. A variety of fourth-order compact methods have been employed solving partial differential equations ^{[9,12,13,14,21,22,23,24,25,28,29]}. Furthermore, Wang ^[26] proposed a conservative eighth-order compact difference scheme for the nonlinear Schrödinger equation. In ^[27], Chen and Chen presented a conservative eighth-order compact difference scheme for the Klein-Gordon-Schrödinger equations. Motivated by ideas in ^[26,27], this article aims to construct a new general difference scheme which can deal with the conservativeness of the invariants and convergence theorem easily. In detail, there are following three advantages:

(ⅰ) The proposed scheme is compact, linearized, decoupled.

(ⅱ) The proposed scheme preserves several invariants in discrete sense.

(ⅲ) The operator form of scheme is novel and can be easily generalized from the fourth-order compact method to the eight-order method for solving other equations.

The rest of the paper is organized as follows. In Section 2, we introduce an eighth-order conservative compact finite difference scheme and apply it to solve the coupled Schödinger-KdV equations numerically. The discrete conservation properties of the proposed nonlinear scheme is analyzed and the convergence theorem of the linearized scheme is established in Section 3. Numerical experiments are presented in Section 4. Finally, a brief conclusion is given in Section 5.

2.   The eighth-order compact finite difference scheme

The domain $\Omega = \{(x, t)\, |\, 0 \leq x\leq l, \, 0 \leq t \leq T\}$ is discretized into grids described by the set $\{x_j, t_n\}$ of nodes, in which $x_j = jh, \, j = 0, 1, \ldots, J = l/h$ and $t_n = n\tau, \, n = 0, 1, \ldots, N = T/\tau$, where $h$ and $\tau$ are discretization parameters. Briefly, let $u^n_j = u(x_j, t_n)$, $v^n_j = v(x_j, t_n)$ and $\Omega_h = \{x_0, x_1, \ldots, x_J\}$. For more convenient discussion, define the following difference operators and notations:

About the approximate formulas of the first and second-order spatial derivatives at all nodes (with periodic boundary conditions) with the eighth-order accuracy, we have the following lemma. Note that we denote $u'_j = \frac{\partial u(x_j, t)}{\partial x}$ or simply denote $u'_j = (u_x)_j$ in the following lemma. Similarly, the notations $u''_j$ and $u'''_j$ are the same meaning.

Lemma 1. ^[1] For $u'$ and $u''$, we have the following approximate formulas

For the convenience to discrete and analyse the equations, we need to rewrite the relations (2.1) and (2.2) to the operator's forms.

Lemma 2. By the definition of the operators above, we have

Proof. Assume that there is an operator $\mathcal{A}_1^*u'_j = \lambda_1u_{j-1}+ \lambda_2 u_j+\lambda_1 u_{j+1}$ such that

By computation and the definition of operators above, we have $\lambda_1 = 5/42$ and $\lambda_2 = 32/42$. Hence, $\mathcal{A}_1^* = \mathcal{A}_1$ and (2.3) holds. (2.4) can be proved similarly. At last, combining (2.3) and (2.4), (2.5) follows directly.

We note that Lemma 2 shows the discrete scheme has the eighth-order accuracy if we use the operators $\mathcal{A}_1$, $\mathcal{A}_2$, $\mathcal{B}_1$ and $\mathcal{B}_2$ or their combinations to discrete the corresponding derivative values at nodes.

In the temporal discretization, we need to evaluate the function values at mid-nodes $\left(\left(n+\frac{1}{2}\right){\text{-nodes}} \right)$. The following lemma is necessary to ensure to approximate the function values at mid-nodes by values at $n$- and $(n+1)$-nodes, which can be obtained by Taylor's expansion.

Lemma 3. For any smooth function $g(t)$ and $m \in \mathbb{N}^*$, we have

where $t_{n+\frac{1}{2}} = \frac{t_n+t_{n+1}}{2}$ and

Proof. By using Taylor's expansion, we have

where $z \in \mathbb{N}^*$. Let $k < \frac{m}{2}, k \in \mathbb{N}$, from (2.11) and (2.12), we can obtain

Plugging (2.13) into (2.8), (2.7) immediately follows.

Denote the approximations of $u^n_j$ and $v^n_j$ by $U^n_j$ and $V^n_j$, respectively. Ignoring the truncation error terms in Eqs (2.3)–(2.5) and (2.7), we obtain the following implicit compact scheme with truncation error $O(\tau^2 + h^8)$ by using the Crank-Nicolson method for temporal discretization and Lemmas 2 and 3:

The schemes (2.14 and 2.15) are nonlinear and gotten by discretizing the temporal derivative with the Crank-Nicolson method, which has the second-order $O(\tau^2)$ and discretizing the special derivatives with the operators $\mathcal{B}_1$ and $\mathcal{B}_1 \mathcal{B}_2$ for (1.1) and (1.2), respectively, which has the eighth-order $O(h^8)$ by Lemma 2.

As to the linearized form of (2.14 and 2.15), we will discuss in the next section.

3.   The conservation and convergence analysis

3.1.   Notations and preliminaries

Let $H_p(\Omega_h) = \left\{u \, | \, u = \left\{u_j\right\}, \, j = 0, 1, \ldots, J \, \, {\text{and}} \, \, u_j = u_{j+J} \right\}$ denote the space of periodic real- or complex-valued grid functions defined on $\Omega_h$ with period $J$. The discrete inner product and the corresponding discrete $L^2$-norm on the grid function space $H_p(\Omega_h)$ are defined as

where $\overline{w}$ denotes the conjugate of $w$. Norm $\lvert \lvert \delta^2_x u \lvert \rvert^2 = \left < \delta^2_x u, \delta^2_x u \right >$ is well-defined with periodic condition $(u_j = u_{j \pm J})$ and the discrete $L^{\infty}$- and $H^1$-norm are defined as

The following lemmas can be easily proved:

Lemma 4. For any grid functions $u, w \in H_p(\Omega_h)$, we have

Lemma 5. For any grid functions $u \in H_p(\Omega_h)$, we have

Lemma 6. ^[20] For any grid functions $u \in H_p(\Omega_h)$, we have

Lemma 7. ^[12] For a real circulant matrix $A = C(b_0, b_1, \ldots, b_{n-1})$, all eigenvalues of $A$ are given by

where $f(x) = b_0+b_1x+b_2x^2+\ldots+b_{n-1}x^{n-1}$, and $\mu_k = { \cos}(\frac{2k\pi}{n})+i { \sin}(\frac{2k\pi}{n})$.

3.2.   Conservation and error analyses

For the compact schemes (2.14) and (2.15), we have the following conservative properties in the discrete sense. The process of proof is similar to ^[9]. Since it still has some difference, so for the convenience to read, we give the detail of proof as following:

Theorem 1. The compact schemes (2.14) and (2.15) preserve the discrete conservation laws of the numbers of plasmons and particles, i.e.,

and

where $U^n = (U^n_1, U^n_2, \ldots, U^n_J)^T$.

Proof. Setting $G^n$ is a vector with the component

then (2.14) can be written as

Computing the inner product $\left < \cdot, \cdot \right >$ on both sides of Eq (3.4) with $U^{n+\frac{1}{2}}$, $\mathcal{A}_2^{-1} G^n$, $\delta_t (\mathcal{A}_2^{-1} U^n)$, $\mathcal{A}_2^{-1} \delta_{x}^2 G^n$ and $\delta_t (\mathcal{A}_2^{-1} \delta_{x}^2 U^n)$, respectively, and applying Lemma 4, we obtain

By the Hermitian property of inner product and multiplying Eqs (3.7) and (3.9) by $i\epsilon$, we can obtain

By Lemma 4 and Eqs ($3.6$), ($3.8$), ($3.10$) and ($3.11$), it follows that

Multiplying by $p$, $\frac{20 h^2}{70}$ and $\frac{h^4}{70}$ in Eqs (3.5), (3.12) and (3.13), respectively, and eliminating the term $\left < \delta_{\overline{x}} U^{n+\frac{1}{2}}, \delta_{\overline{x}} G^n \right >$, we obtain

From the definition of $G^n$ we can see that $\left < G^n, U^{n+\frac{1}{2}} \right >$ is a real number. Hence, the imaginary part of (3.14) is zero, i.e.,

Applying Lemma 4 in (3.15), we can obtain immediately that

Computing the inner product $\left < \cdot, \cdot \right >$ on both sides of Eq (2.15) with ${\bf{1}} : = (1, 1, \ldots, 1)^T \in H_p(\Omega_h)$, we can obtain

where

using the periodic conditions, one can have the equation

i.e.,

With the periodic conditions, (3.2) immediately satisfies. The proof is finished.

Hereinafter we define

and

The compact schemes (2.14) and (2.15) are equivalent to the following matrix equations:

where ${\textbf{A}}_1$, ${\textbf{A}}_2$, ${\textbf{B}}_1$ and ${\textbf{B}}_2$ are $J \times J$ matrices corresponding to the operators $\mathcal{A}_1$, $\mathcal{A}_2$, $\mathcal{B}_1$ and $\mathcal{B}_2$, respectively,

By the properties of circulant matrices, we can see that matrices ${\textbf{A}}_1$, ${\textbf{A}}_2$, ${\textbf{B}}_1$ and ${\textbf{B}}_2$ are circulant symmetric positive definite ^[10]. Let $\mathcal{A}_1^{-1}$, $\mathcal{A}_2^{-1}$, $\mathcal{B}_1^{-1}$ and $\mathcal{B}_2^{-1}$ denote inverse operators of $\mathcal{A}_1$, $\mathcal{A}_2$, $\mathcal{B}_1$ and $\mathcal{B}_2$, respectively. The matrices corresponding to the operators $\delta_{x}^2$, $\delta_{\hat{x}}$, $\mathcal{A}_1^{-1}$, $\mathcal{A}_2^{-1}$, $\mathcal{B}_1^{-1}$ and $\mathcal{B}_2^{-1}$ are also circulant, therefore, they commute under multiplication.

The compact schemes (2.14) and (2.15) can be equivalently written as

which can be obtained by multiplying $\mathcal{B}_2^{-1}$ and $\mathcal{B}_1^{-1} \mathcal{B}_2^{-1}$ in both side of (2.14) and (2.15), respectively.

By applying Lemma 7, we can obtain the following result:

Lemma 8. ^[14] For any grid function $u \in H_p(\Omega_h)$, we have the inequalities

Define

Lemma 8 shows that $\lvert \lvert \lvert \cdot \lvert \lvert \rvert_Q$ and $\lvert \lvert \lvert \cdot \lvert \lvert \rvert_P$ are norms on $H_p(\Omega_h)$ equivalent to the discrete $L^2$-norm $\lvert \lvert \cdot \lvert \rvert$. For the proof of the following theorem, we want the relations (3.20) and (3.21).

Theorem 2. The compact schemes (2.14) and (2.15) preserve the energy of oscillations in discrete sense, i.e.,

where

Proof. Computing the inner product $\left < \cdot, \cdot \right >$ on both sides of Eq (3.20) with $\delta_tU^n$, we can obtain the following equation with the commutativity under multiplication of circulant matrices:

Then taking the real part of Eq (3.23), we obtain

Multiplying Eq (3.24) with $2\tau$ and summing from 0 to $n$, we obtain

Setting $W^n$ is a vector with the component

then Eq (3.21) can be written as

Computing the inner product $\left < \cdot, \cdot \right >$ on both sides of Eq (3.26) with $W^n$ and applying Lemma 5, we can obtain

It follows from definition (2.8) that

It is easy to see that

Multiplying Eq (3.27) with $2\tau$ and summing from 0 to $n$, we have

Since

the Eq $(3.28)$ becomes

Multiplying Eqs (3.25) and (3.29) with $\rho$ and $\frac{q}{2}$, respectively, and subtracting the results, (3.22) follows immediately.

In the following convergence analysis, we will take the symbol $C$ as a general positive constant independent of $h$ and $\tau$, not necessarily the same at different occurrences. We assume that there is a positive constant $Y^*$ such that the exact solutions $u$ and $v$ of the coupled system satisfy

Let $Y_0 = 2(Y^*+1)^2$ and define a smooth function $\Psi(r, s)\in C^{\infty}(R^2)$ as

Since schemes (2.14) and (2.15) are nonlinear, we change it into the following linearized compact scheme to reduce computational cost:

where $\hat{U}^0 = (U^{0*}+U^0)/2$, $\hat{U}^n = 3U^{n}/2-U^{n-1}/2$, and $U^{n*} = 2U^n-U^{n-1}$ for $n\ge 1$.

We can prove that the temporal and spatial convergence rates of the linearized compact schemes (3.34) and (3.35) are second- and eighth-order, respectively.

Lemma 9. Let $\left\{ y_n \right\}$ be a nonnegative real sequence, $c$ a nonnegative constant, $d$ and $\tau$ are positive constants. If

then

Theorem 3. Suppose that $u, v \in C^4(0, T;C^{11}(R))$ are the exact solutions to Eqs (1.1) and (1.2), $h^8\tau^{-1} = o(1)$, and that assumption (3.30) holds. Let $U$ and $V$ be the solutions of (3.34) and (3.35). Then there exists a constant $C = C(Y^*, T)$ such that

for $h$ and $\tau$ sufficiently small.

Proof. Let $E^n_u = u^n-U^n$ and $E^n_v = v^n-V^n$. By Eqs (1.1), (1.2), (3.34) and (3.35), and ignoring the subindex $j$, we obtain

where

By the assumption (3.30) and definition (3.31), one can see that $\Psi((v^n, v^{n+1})) = \psi((v^n, v^{n+1}))$, and hence, the truncation errors $r^n_u$ and $r^n_v$ are such that $r^n_u = O(\tau^2 + h^8)$ and $r^n_v = O(\tau^2 + h^8)$. Under the smoothness assumption of $u$ and $v$, we have

From (3.36) and (3.37), we can obtain the following equations:

where $R^n_u = \mathcal{B}_2^{-1} r^n_u$ and $R^n_v = \mathcal{B}_2^{-1} \mathcal{A}_1^{-1} r^n_u$.

We use the induction argument as in ^[15,16,17] to estimate the error bounds. To obtain our error estimate, we assume that there exists a constant $h_0 > 0$ such that, for $0 < h \le h_0$,

Since $E^0_u = E^0_v = 0$, it is easy to see that

For $n \ge 1$, by computing the inner product $\left < \cdot, \cdot \right >$ on both sides of (3.38) with $E^{n+\frac{1}{2}}_u$. we can obtain following equation by Lemma 4:

Taking the imaginary part of (3.41), we can obtain the inequality

By computing the inner product $\left < \cdot, \cdot \right >$ on both sides of (3.39) with $E^{n+\frac{1}{2}}_v$. we can obtain following equation by Lemma 4:

Since

we can have the inequality

By Lemma 5, Eq (3.43) and inequality (3.44), we have

Summing inequalities (3.42) and (3.45) side by side, and using inequality (3.44), we can have following inequality with $E^0_u = E^0_v = 0$:

By Computing the inner product $\left < \cdot, \cdot \right >$ on both sides of (3.38) with $\delta_t E^n_u$. we can obtain following equation by Lemma 4:

Taking the real part of (3.47) and summing from 0 to $k$, we can obtain

By using a method of summation by parts together with assumptions (3.30) and (3.40), we have the inequalities

By (3.48) and the above estimates, we can obtain

For any real-valued grid function $f$, an operator $\Theta$ is defined by

with $\Theta f_j = \Theta f_{j+J}$. The following results can be easily proved:

Then define a matrix ${\textbf{J}}$ corresponding to the operator $\mathcal{J}$, i.e.,

It's obvious that ${\textbf{J}}$ is invertible and ${\textbf{J}}^{-1}$ is circulant symmetric positive definite as the scale $J$ of matrix is an odd integer. By computing the inner product $\left < \cdot, \cdot \right >$ on both sides of (3.39) with $\delta_t (\mathcal{J}^{-1} \Theta E_v)^n$ and applying Lemma 4, (3.51) and (3.52), we can obtain

Since ${\textbf{J}}$, ${\textbf{A}}_1$, ${\textbf{B}}_1$ are circulant symmetric positive definite, so there exists ${\textbf{G}}$ such that ${\textbf{J}}^{-1} {\textbf{A}}^{-1}_1 {\textbf{B}}_1 = {\textbf{G}} {\textbf{G}}^T$. By (3.51) and (3.52), we can have

Summing Eq (3.54) from 0 to $k$ together with (3.55), we can obtain

By using a method of summation by parts together with assumptions (3.30) and (3.40), we have the inequalities

Noticing that $\mathcal{J}({\textbf{I}}- \frac{h^2}{4} \delta^2_x + \frac{h^4}{16} \delta^2_x \delta^2_x -\frac{h^6}{64} \delta^2_x \delta^2_x \delta^2_x) = {\textbf{I}}-\frac{h^8}{256} \delta^2_x \delta^2_x \delta^2_x \delta^2_x$, we have

By using Lemma 6 and (3.53), the above inequality can be written as

Multiplying Eq (3.56) with $2\tau$, we can obtain

Since the norms $\lvert \lvert \cdot \lvert \rvert$, $\lvert \lvert \lvert \cdot \lvert \lvert \rvert_P$, and $\lvert \lvert \lvert \cdot \lvert \lvert \rvert_Q$ are equivalent, we can have the following inequality by summing (3.46), (3.49) and (3.57):

Taking $\tau$ sufficiently small and applying Lemmas 8 and 9, we can obtain

Moreover, we need to confirm the inequality in (3.40) holds for $n = k+1$ to complete our proof. We can get the following inequalities by Lemma 6:

and similar inequalities hold for $\lvert \lvert E^{k+1}_v \lvert \rvert_{\infty}$ and $\lvert \lvert \delta_t E^{k}_u \lvert \rvert_{\infty}$. Then it's easy to see that the inequalities above hold for $n = k+1$ when $h^8 \tau^{-1} = o(1)$, i.e., $h^8 \tau^{-1} \rightarrow 0$ as $h \rightarrow 0$, and taking $h$ sufficiently small, which implies that assumption (3.40) is valid for $n = k+1$. The proof is finished.

Corollary 1. By applying Lemma 6, we can obtain the following optimal order convergence rate under the same conditions in Theorem 3:

4.   Numerical experiments

In this section, some numerical examples are presented to illustrate the conservative properties and eighth-order accuracy of the proposed compact scheme. The ultimate compact schemes (3.32)–(3.35) can be written as the following linear matrix equations:

where $E_1$, $E_2$, $F_1$ and $F_2$ are nonlinear terms. Our numerical experiments are conducted using Matlab (R2019b). The invariants $I_1, I_2, I_3$ and $I_4$ are tested by the discrete formulations:

and the errors of invariants are defined as

Moreover, the accuracy of the proposed scheme is tested by the discrete $L^2$- norm $(\lvert \lvert u-U \lvert \rvert + \lvert \lvert v-V \lvert \rvert)$ and $L^{\infty}$- norm $(\lvert \lvert u-U \lvert \rvert_{\infty} + \lvert \lvert v-V \lvert \rvert_{\infty})$.

Example 1. ^[8] We consider the following Cauchy problem:

whose exact solutions are given by $u(x, t) = {\text{exp}}(i(x+t/4))$ and $v(x, t) = 3/4$. we then compute the equations with $h = \pi/20$ and $\tau = 0.001$ in the spatial interval $[0, 2\pi]$. The errors of the numerical invariants at different time are listed in Table 1, which indicates that the proposed compact scheme preserves the conservation properties. Table 2 shows that the convergence rate of the proposed compact scheme is eighth-order in space.

Example 2. ^[3] We consider the following coupled equations:

with exact solutions

where $c$ is an arbitrary positive constant. In addition, we set the artificial boundary conditions $u(a, t) = u(b, t) = 0$ and $v(a, t) = v(b, t) = 0$ to satisfy the physical condition that $\lvert u \rvert$ and $v$ tend to zero as $\lvert x \rvert \to \infty$. Our simulations are conducted by taking $\epsilon = 1$, the traveling wave speed $c = 0.45$ and initial conditions

Table 3 lists the numerical solutions at $t = 0.001$, with $h = 0.25$, $\tau = 0.00001$ and $[a, b] = [-30, 30]$, where the scheme MECS expands $[a, b]$ to $[-150,150]$ for reducing boundary truncation error. Compared with the numerical results obtained by the fourth-order compact scheme (FCS) in ^[9] and exponential time differencing three-layer implicit scheme with Padé approximation (ETDT-P) in ^[7]. We can see that the eighth-order compact scheme (ECS) and modified eighth-order compact scheme (MECS) give better approximations. In addition, MECS gives much more accurate error estimate than ECS, which is caused by boundary truncation error. The numerical solution profiles of $\lvert U \rvert$ and $V$, as well as the contours in Figure 1 show that the waves traveling with a speed $c = 0.45$ keep the shape and hight, which are in good agreement with the exact solutions.

Example 3. ^[11] We consider the following coupled equations:

with exact solutions

where $\sigma \in (-1/6, 0)$ and $\omega \in R$. Set the artificial boundary conditions $u(a, t) = u(b, t) = 0$ and $v(a, t) = v(b, t) = 0$. Our simulations are conducted by taking $\sigma = -1/12$, $\omega = 0$, $[a, b] = [-40, 70]$, the traveling wave speed $c = (1+\sqrt{71/72})/2$ and initial conditions

The errors of the numerical invariants at different times are listed in Table 4, which indicates that the proposed compact scheme preserves the conservation properties. Table 5 shows that the convergence rate of the proposed compact scheme is eighth-order in space. The numerical solution profiles of $\lvert U \rvert$ and $V$, as well as the contours in Figure 2 show that the waves traveling with a speed $c = 0.99652$ keep the shape and hight, which are in good agreement with the exact solutions.

5.   Conclusions

In this paper, we propose an eighth-order compact finite difference scheme by constructing several circulant symmetric positive definite matrices to obtain the numerical solution of coupled Schrödinger-KdV equations. The performance of proposed compact scheme is evaluated by conservation properties and error estimate. Numerical examples demonstrate the better performance of the proposed compact scheme in accuracy compared with FCS and ETDT-P given in ^[7,9]. Since the matrices have good properties, we can discuss the possibility that the proposed compact scheme can be applied to other equations such as nonlinear Dirac equation ^[21], generalized Rosenau-RLW equation ^[22], Klein-Gordon-Schrödinger equation ^[23], coupled Gross-Pitaevskii equations ^[24] and regularized long wave equation ^[25].

Acknowledgments

This work was supported by National Natural Science Foundation of China (No. 11471092).

Conflict of interest

The authors declare no conflicts of interest.