SAV Galerkin-Legendre spectral method for the nonlinear Schrödinger-Possion equations

1.   Introduction

In this paper, we present a fully-discrete structure-preserving scheme for the following nonlinear Schrödinger-Possion equations

where ${\rm{i}} = \sqrt{-1}$, $\Omega = (-1, 1)$. $\psi(x)$ is the real-valued potential function that represents the external field. $\beta\in R$ is a coupling constant that represents the relative strength of the Poisson potential, $\beta\ge0$ holds in the case of attracting forces and $\beta < 0$ holds in the case of repulsive forces. The Schrödinger-Possion system is a local single particle approximation of time-dependent Hartree-Fock system. This nonlinear system (1.1) has important applications in many quantum systems^[1,2], where the complex-valued function $u(x, t)$ denotes the single particle wave function, and $v(x, t)$ represents the Poisson potential relative to the boundary condition.

In recent years, much attention has been paid to developing efficient numerical schemes for solving the Schrödinger-Possion equations (1.1). For instance, Soler and Ringhofer ^[3] proposed a modified Crank-Nicolson scheme, which is mass- and energy-conserving in the discrete level. Dong ^[4] improved the numerical methods for the general 3D case. And the computational cost was significantly reduced. Zhang and Dong ^[5] applied the backward Euler and time-splitting sine pseudo-spectral methods to study the ground state and dynamics of 3D system in different setups. Auzinger et al. ^[6] used operator splitting methods combined with finite element spatial discretizations to solve the problem. Cheng et al. ^[7] proposed a fast spectral element method which can reach exponential accuracy for the Schrödinger-Poisson system. Lubich ^[8] gave an error analysis of Strang-type splitting integrators for the Schrödinger-Poisson equations. Furthermore, there exists many other numerical schemes for Schrödinger-type equations. Li et al. ^[9] developed efficient numerical schemes for the coupled fractional Klein-Gordon-Schrödinger equation by combining the Crank-Nicolson/leap-frog difference methods for the temporal discretization and the Galerkin finite element methods for the spatial discretization. Li et al. ^[10] constructed the conservative linearized Galerkin finite element methods (FEMs) for the nonlinear Klein-Gordon-Schrödinger equations. Li et al. ^[11] proposed a fully discrete scheme for the nonlinear fractional Schrödinger equation with wave operator by combining the Crank-Nicolson method in temporal direction with the Galerkin finite element method in the spatial direction. Antoine et al. ^[12] applied the finite difference time domain methods and time-splitting spectral method to solve the nonlinear Schrödinger/Gross-Pitaevskii equations. More details about the related papers can be found in ^{[13,14,15,16,17,18,19]}. Up to now, most schemes are proved to be mass- and energy-conserving. The convergence results are missing in most references.

In this paper, we aim to develop structure-preserving schemes as well as their error analysis for the Schrödinger-Poisson system. We firstly introduce a scalar auxiliary variable (SAV) and rewrite the equations as a new family of partial differential equations. Then, we apply the Legendre-Galerkin spectral method in the spatial discretization and the Crank-Nicolson method in the temporal discretization. We show the fully-discrete scheme is mass- and energy-conserved. Moreover, we obtain the boundedness of the numerical solutions based on the Gagliardo-Nirenberg and some Sobolev inequalities.

In terms of the bounded numerical solutions, we present a rigourous error analysis of the fully discrete numerical scheme. The convergence results indicate that the fully-discrete scheme is of order $2$ in the temporal direction and decreases exponentially in the spatial direction.

Compared with other numerical schemes for the nonlinear Schrödinger-Possion equations, the proposed scheme is proved to be mass- and energy-conserved without any time steps restrictions. Furthermore, by applying spectral method in space, the proposed scheme can reach exponential convergence in space. We also give the unconditional convergence in the $L^\infty$-norm in the final section.

It is remarkable that the key to developing structure-preserving method is the so called scalar auxiliary variable (SAV) approach. The approach was firstly proposed by Shen et al. in ^[20,21] and has a successful application in developing energy stable or energy-conserving schemes for time-dependent partial differential equations, such as gradient flows ^[22,23,24], wave equations ^[25,26,27], Schrödinger equations ^[28,29] and so on ^[30]. Up to now, most numerical results are obtained in the real spaces and much attention is paid on energy-conserving properties of the scheme. In this paper, the study is extended to a coupled system in the complex space and the unconditional converge results are investigated.

The rest of the paper is organized as follows. In Section 2, a fully discrete scheme for the generalized nonlinear Schrödinger-Possion equation is introduced. In Section 3, the convergence of the proposed scheme is discussed. In Section 4, some numerical experiments are provided to demonstrate the theoretical results. Finally, some concluding remarks are presented in Section 5.

Throughout the paper we use $C$ and $C_i(i\in N)$ to denote positive constants, which could have different values in different places.

2.   Fully discrete scheme

In this section, we present the fully discrete method based on SAV approach.

2.1.   Preparation

Through the paper, we use the following notations:

where $0\le l\le k, \, k\in N$ and $\bar{v}$ denotes the conjugate of $v$. Firstly, we assume that the exact solutions of Schrödinger-Possion equations (1.1) satisfy

where $l_1 = 0, \, 1, \, 2, \, 1_2 = 0, 1, \, \sigma\ge2$ and $K^* > 0$ is a positive constant.

We also let $\phi_k(x) = L_k(x)-L_{k+2}(x)$, where $\{L_k(x)\}_{k = 0}^N$ represents the Legendre polynomials. Define

Denote the polynomial space by

For all $u\in L^2(\Omega)$, we define $P_L:L^2(\Omega)\mapsto X_N$ as:

Let $\nu = P_L u$ and $\nu = \Delta P_L u$ respectively, we can obtain the following stability properties:

Thus $\Vert P_L u\Vert_{H^1}\le\Vert u\Vert_{H^1}, \, \forall u\in H_0^1(\Omega)$. Meanwhile, we define $\Pi_N:H_0^1(\Omega)\mapsto X_N$ as:

Then we have the following lemma.

Lemma 1. (^[31])For all $u\in H_0^1(\Omega)\cap H^m(\Omega)(m\ge 1)$, there holds

where $C > 0$ is a constant independent of $N$.

We collect some lemmas here. They play an important role in the proof of the main results.

Lemma 2. (^[32])Suppose that $v(x)\in C^1[a, b]$, and $v(a) = v(b) = 0$, then

Lemma 3. (^[33])Let m be a nonnegative integer, and assume $g\in H^m(\Omega)$ and $\partial\Omega\in C^{m+2}$. Suppose that $\nu\in H_0^1(\Omega)$ is the unique solution of the boundary problem

Then $\nu\in H^{m+2}(\Omega)$ and the following estimate holds

where $C$ depending on $m, \, \Omega$ and $\alpha, \, \beta$.

Lemma 4. (^[32])Suppose that $u(t)\in C^3[0, T]$, then we have

where C is a constant independent of $t_{n+1}-t_n$.

Lemma 5. (^[34])Suppose that the discrete mesh function $\{\varpi^n|n = 1, \, 2, \, \cdots, \, M;M\tau = T\}$ is nonnegative and satisfies recurrence formula

where A, B and ${C_n(n = 1, \, 2, \, \cdots, \, M)}$ are nonnegative constants. Then

where $\tau$ is small such that $(A+B)\tau\le\frac{M-1}{2M}(M > 1)$.

Lemma 6. (^[35])The Legendre polynomials $\{L_n(x)\}$ satisfy the following three-term recurrence relation

and the orthogonality relation

where $k, \, j\ge 1$.

2.2.   The SAV-CN scheme

In system (1.1), a scalar auxiliary variable can be introduced:

where

Thus (1.1) can be rewritten as

where $b(u) = (f(|u|^{2})u+\psi(x)u)/{\sqrt{\mathcal{H}(t)+C_0}}$ and $\Re(\cdot)$ means to take the real part of a complex number.

Then, we present the fully-discrete scheme. Let $t_n = n\tau, n = 0, 1, \cdots, N_t$, where $\tau = \frac{T}{N_t}$ and $N_t$ is a positive integer. For a sequence of functions $\{\varphi^i\}_{i = 0}^{N_t}$, we denote

The SAV-CN scheme is to find $u_N^{n}, v_N^n\in X_N$ and $w_N^{n}\in R^1$ such that, for $n = 0, \, 1, \, 2, \, \cdots, \, N_t-1$:

where $b({u}_N^{n+\frac{1}{2}}) = (f(|u_N^{n+\frac{1}{2}}|^2){u}_N^{n+\frac{1}{2}}+\psi(x){u}_N^{n+\frac{1}{2}})/\sqrt{H(t_{n+\frac{1}{2}})+C_0}$.

Firstly, we have the following continuous and discrete mass and energy conservation laws, respectively.

Theorem 1. Suppose that $u, \, v$, and $w$ are the solutions of system $(2.2)–(2.4)$, then the continuous mass and energy conservation laws can be obtained as follows:

Besides, suppose that $u_N^n, \, v_N^n$ and $w_N^n$ are the solutions of system $(2.5)–(2.7)$, we have

Proof. Case (a): Proof of (I) and (II).

Taking the inner product of Eq (2.2) with $\bar u$, we can obtain that

Meanwhile, multiplying Eq (2.3) with $-\beta\bar v$, and integrating it over $\Omega$, we can get

Summing up Eqs (2.8) and (2.9), we can get

Taking the imaginary part, we can obtain

Besides, multiplying Eq (2.2) with $\bar u_t$, and integrating it over $\Omega$, we can get:

Taking the real part, we can obtain

which is equivalent to

Since $(v, |u|^2) = \Vert v_x\Vert^2$, we can obtain the mass and energy conservation laws

Case (b): Proof of (III) and (IV).

Let $\phi = u_N^{n+\frac{1}{2}}$ in Eq (2.5). Taking the imaginary part, we have

which implies $\tilde{\mathcal M}^{n+1} = \tilde{\mathcal M}^{n}.$

Let $\phi = D_\tau{u_N^n}$ in Eq (2.5). Consider the real part, one has

where

Therefore, we have

This completes the proof.

Then we have the following unconditional convergence results.

Theorem 2. Suppose that the exact solutions of Schrödinger-Possion equations $(1.1)$ satisfy $(2.1)$. Then, there exists a positive constant $\tau_0$, such that when $\tau \leq \tau_0$, the fully-discrete system defined in Eqs $(2.5)–(2.7)$ has the numerical approximations $\{ u_N^m, \, v_N^m, \, w_N^m \}$, $m = 1, 2, 3, \ldots, N_t$, satisfying

where $C$ is a positive constant independent of $\tau$ and $N$.

3.   Numerical analysis

This section is concerned with the proof of convergence analysis of the fully conservative discrete scheme.

3.1.   Boundedness

Firstly, we present the boundedness of the numerical solutions $u_N^n, \, v_N^n$ and $w_N^n$ in $L^\infty$-norm and $H^1$-norm respectively.

Lemma 7. Suppose that the exact solutions satisfy $(2.1)$, then the following estimates hold:

where $K_i(i = 1, \, 2, \, \cdots, \, 5)$ are positive constants independent of $\tau$ and $N$.

Proof. By the discrete mass conservation $\tilde{M}^n = \Vert u_N^n \Vert^2 = \tilde{M}^0,$ we have $\Vert u_N^n \Vert\leq C_1$.

Let $\phi = v_N^n$ in Eq (2.6), we can obtain that

According to Lemma 3, we can get $\Vert v_N^n\Vert_{H^2}\le C_2\Vert|u_N^n|^2\Vert = C_2\Vert u_N^n\Vert_{L^4}^2.$

By using Gagliardo-Nirenberg inequalities, we can obtain

where $C_3$ is dependent on $\varepsilon$. Next, we will show the boundedness of $\Vert(u_N^n)_x\Vert$, $\Vert(v_N^n)_x\Vert$ and $w_N^n$ respectively. Firstly, if $\beta\le 0$, we can conclude from Theorem 1 that $\Vert(u_N^n)_x\Vert^2\le C_4$, $\Vert(v_N^n)_x\Vert^2\le C_5$ and $|w_N^n|\le C_6$.

Secondly, if $\beta > 0$, by Theorem 1, we have

Thus $\Vert(u_N^n)_x\Vert^2+(w_N^n)^2 = C_7+\frac{\beta}{2}\Vert(v_N^n)_x\Vert^2\le\frac{\beta}{2}(\varepsilon \Vert(u_N^n)_x\Vert^2+C_3\Vert u_N^n\Vert^6)+C_7,$ which can be rewritten as

Let $\varepsilon < 1/\beta$, we can get $\Vert(u_N^n)_x\Vert^2\le C_9$, $|w_N^n|\le K_1$. Then we can infer that $\Vert u_N^n\Vert_{H^1}\le K_2$, $\Vert v_N^n\Vert_{H^2}\le K_3$. By using the Soblev embedding theorem, we can obtain

The proof is completed.

We now present the following unconditionally optimal error estimates of numerical scheme (2.5)–(2.7).

3.2.   Proof of Theorem 2

Denote $\tilde{u}^n = \Pi_N u^n$, $e^n = \tilde{u}^n-u_N^n$, $e_0^{n+\frac{1}{2}} = \tilde{u}(t_{n+\frac{1}{2}})-{u}_N^{n+\frac{1}{2}}$, $\eta^n = v^n-v_N^n$ and $\zeta^n = w_n-w_N^n$. Let $u^n = u(x, t_n), v^n = (x, t_n), w^n = w(t_n)$. Subtracting Eqs (2.5)–(2.7) from Eqs (2.2)–(2.4), we have

where

Let $\phi = {e}^{n+\frac{1}{2}}$ in Eq (3.1), we can obtain

where we have used

Thus Eq (3.4) can be rewritten as

where $\Im(\cdot)$ means to take the imaginary part of a complex number. Taking the imaginary part of Eq (3.6), we have

We first estimate $\Im (G_1^{n+\frac{1}{2}}, \overline{e}^{n+\frac{1}{2}})$. Utilizing Cauchy-Schwarz inequality, we have

We know that

By Lemma 1, Lemma 4 and Theorem 7, we have $|G_1^{n+\frac{1}{2}}|\leq C(\tau^2+|{\eta}^{n+\frac{1}{2}}|+|{\zeta}^{n+\frac{1}{2}}|+ |{e}^{n+\frac{1}{2}}|)$. Using Lemma 4, we obtain $|R_1^{n+\frac{1}{2}}|\leq C\tau^2+C N^{-\sigma}$ and $\|(\tilde{u}(t_{n+\frac{1}{2}})-\tilde{u}^{n+\frac{1}{2}})_{xx}\|^2\leq C\tau^2$. Thus

Similarly, we can obtain

and $\Vert R_1^{n+\frac{1}{2}}\Vert_{H^1}\le C(N^{1-\sigma}+\tau^2)$.

Then, setting $\phi = e^{n+1}-e^n$ in Eq (3.1), we get

Taking the real part in Eq (3.10), we have

which can be written as

According to the Sobolev inequalities and the estimations of the truncation errors, we have

where $\epsilon > 0$ is a positive constant.

In order to estimate $\Vert D_\tau e^n\Vert_{H^{-1}}$ above, we consider Eq (3.1), from which we can get the following estimate for any test functions $\nu\in X_N\subseteq H_0^1(\Omega)$:

where $P_L\nu$ represents the $L^2$ projection operator of $\nu$ defined in Section 2.1.

By the duality between $H^{-1}(\Omega)$ and $H_0^1(\Omega)$, we can derive that

Then, multiplying (3.15) with $\epsilon\tau$ and summing up (3.13), (3.15), we get

Let $\phi = \eta^{n+1}$ in Eq (3.2), one has

where we have used Lemma 2. Then

which also implies

Further, taking the inner product of the Eq (3.3) with $2{\zeta}^{n+\frac{1}{2}}$, we can obtain

where

According to Eq (3.20), we can transform Eq (3.19) into

Similar to the analysis above, we can obtain

Based on inequality (3.18) and Theorem 7, we can conclude that

Combining the inequalities (3.9), (3.16), (3.18) and (3.22), we have

By using discrete Gronwall's inequality (Lemma 5), we can obtain

Then it follows that

Using the triangle inequality, we have

which completes the proof.

We further have the error estimates in $L^\infty$-norm.

According to Lemma 2 and (3.24), we can obtain that

Based on (3.17), we can get

By the triangle inequalities (3.28) and (3.29), we obtain the desired results (2.11).

4.   Numerical results

In this section, two numerical examples are presented to verify the theoretical results. In the numerical examples, the errors are defined as follows

where $u(x, t_n), \, v(x, t_n)$ represent the analytical solutions and $u_N^n, \, v_N^n$ represent the numerical solutions respectively. In the numerical simulations, we apply the iterative algorithms ^[36] to solve the nonlinear algebraic equations and the iterative tolerance is $10^{-15}$.

Example 1. Consider the following nonlinear Schrödinger equations

where $\Omega = (-1, 1)$ and $\psi_1(x)$ is a real potential function obtained by the following exact solutions

We set $N = 15$, $T = 1$ and test the convergence orders in the temporal directions. We show the maximum numerical errors at $T = 1$ in Table 1. The numerical results indicate that the numerical method is second-order accurate in temporal direction. Then, we set $\tau = 0.001$ and solve the problem with different $N$. We show the maximum numerical errors in Figure 1, respectively. One can see that when $N$ increases, the error decreases exponentially.

We test the structure-preserving properties of the numerical scheme(i.e., SAV-CN). We also discretized original equations (1.1) in space by Legendre-Galerkin spectral method and in temporal direction by the extrapolated Crank-Nicolson (ECN) method and the standard Crank-Nicolson (CN) method. We let $N = 20$ and $\tau = 0.1$ and show the discrepancies of the discrete mass and energy in Figure 2. One can see that the discrepancies of the discrete mass are sufficiently small for both the SAV-CN and CN methods. Furthermore, we also compare the discrepancies of the discrete energy between SAV-CN and the original energy (Origin) in Figure 2. We can find that the the discrepancies of the original energy is between $10^{-6}$ and $10^{-4}$. In contrast, the discrepancies of the discrete energy are of the order of the machine precision only for the SAV-CN methods.

Example 2. Consider the following nonlinear Schrödinger equations

where $\Omega = (0, 2)$, and $\psi_2(x)$ is a real potential function obtained by the following exact solutions

In this example, we firstly use the linear transformation $x = s+1(s\in(-1, 1))$ and solve the resulting problem by using the proposed method. To test the convergence orders in the temporal direction, we set $N = 15$ and refine the temporal stepsizes. The maximum numerical errors at $t = 1$ and convergence orders are shown in Table 2. The numerical results show that the convergence orders in the temporal directions are of $2$. Then, we set $\tau = 0.001$ and solve the problem with different $N$. We present the maximum numerical errors in Figure 3, respectively. Again, we find that when $N$ increases, the error decreases exponentially.

Next, we let $N = 20$, $\tau = 0.1$ and solve the problem by using previous mentioned methods. We show the discrepancies of the discrete mass and energy in Figure 4, respectively. Again, the discrepancies of the discrete mass are less than $10^{-15}$ for both the SAV-CN and CN methods. Besides, we also compare the discrepancies of the discrete energy between SAV-CN and the original energy (Origin) in Figure 4. We can conclude that the the discrepancies of the original energy is between $10^{-7}$ and $10^{-5}$. In contrast, the discrepancies of the discrete energy are sufficiently small only for the SAV-CN methods. These results further confirm the structure-preserving properties of the proposed method.

5.   Conclusions

In this paper, we consider numerical solutions of the nonlinear Schrödinger-Possion equations. The fully-discrete scheme is developed by combing the SAV approach and the Crank-Niclson Galerkin-Legendre spectral methods. The fully discrete scheme is proved structure-preserved. The unconditional stability and convergence results are obtained. It is shown the fully-discrete scheme is of order $2$ in the temporal direction and decreases exponentially in the spatial direction. Numerical results are given to confirm our theoretical findings.

Acknowledgments

The work was supported by NSFC (Grant Nos. 11771162, 62032023, 12002382, 11275269, 42104078, 61902411), Open Research Fund from State Key Laboratory of High Performance Computing of China (HPCL) (Grant No. 202101-01).

Conflict of interest

The authors declare no conflicts of interest.