Split-step quintic B-spline collocation methods for nonlinear Schrödinger equations

1.   Introduction

Efficient and reliable numerical methods are always urgently needed for nonlinear multi-dimensional problems due to the mass storage and the high computation cost. The split-step method, also known as the time-splitting method, is one kind of efficient skills for solving nonlinear parabolic or Schrödinger-type problems in multi-dimensions.

The split-step skill could be efficiently combined with the (compact) finite difference method ^[1,2], the finite element method ^[3,4], the spectral/pseudospectral method ^[5,6], et al. So the authors of ^[7] try to combine the skill with the cubic B-spline collocation (3BC) approach, and they formulate the split-step 3BC (SS3BC) method sucessfully. However, the SS3BC method is just second-order accuracy in space, and there is no numerical analysis in ^[7]. Genarally speaking, the 3BC method has second-order accuracy, and the quintic B-spline collocation (5BC) method is fourth-order ^[8]. Thus, split-step 5BC (SS5BC) methods are constructed in this paper to improve the accuracy in space, and some analyses are also given.

In this paper, we consider the nonlinear Schrödinger (NLS) equations as follows

with the initial condition

and the periodic boundary condition

where ${\bf L} = (L_x, L_y, L_z)$, and $L_x, L_y$ and $L_z$ are the periodic lengths respectively in the $x, y$ and $z$ direction. $u({\bf x}, t)$ is an unknown complex function, $V({\bf x}, t)$ and $u_0({\bf x})$ are given functions, $\alpha, \beta$ are real constants, $T > 0$ and $i^2 = -1$.

Computing the inner product of Eq (1.1) with $u$, and taking the imaginary part, one obtains the following conservation law

The rest of this paper is organized as follows. In Section 2, some preliminaries are introduced, and SS5BC methods are constructed in Section 3. In Section 4, solvability, conservation and linear stability of the schemes are discussed. Numerical experiments are carried out in Section 5, and the present methods are applied to study BECs in Section 6. Finally, some conclusions are given in Section 7.

2.   Preliminaries

We consider Eq (1.1) within ${\bf x}\in \Omega _d\subseteq {\bf R}^d$. Take $\Omega _3 = [x_L, x_R]\times [y_L, y_R]\times [z_L, z_R]$, and $x_R-x_L = L_x, y_R-y_L = L_y, z_R-z_L = L_z$. Let $\{x_k\}_{k = 0}^{N_x}\otimes \{y_l\}_{l = 0}^{N_y}\otimes\{z_m\}_{m = 0}^{N_z}$ be the partition of $\Omega _3$, such that

where $h_x = (x_R-x_L)/N_x, h_y = (y_R-y_L)/N_y$ and $h_z = (z_R-z_L)/N_z$ are the step sizes in space, and $N_x, N_y$ and $N_z$ are positive integers. Obviously, $\Omega_1$ and $\Omega_2$ are reduced forms of $\Omega_3$ which could also be partitioned. We divide the interval $[0, T]$ by the partition $\{t_n\}_{n = 0}^{N_t}$, where $t_n = n\tau$, $\tau = T/N_t$ is the step size in time, and $N_t$ is a positive integer.

For one dimension, let $\{B_j(x)\}_{j = -2}^{N_x+2}$ be the quintic B-spline basis functions ^[8] on the knots $\{x_j\}_{j = 0}^{N_x}$, such that

Similarly, we can define $\{B_k(y)\}_{k = -2}^{N_y+2}$ and $\{B_l(z)\}_{l = -2}^{N_z+2}$ as the basis functions on the knots $\{y_k\}_{k = 0}^{N_y}$ and $\{z_l\}_{l = 0}^{N_z}$, respectively.

A global approximation solution $U_{N}(x, t)$ of the exact solution $u(x, t)$ could be expressed in terms of the quintic B-spline as

where $\delta_j(t)$ are unknown time dependent parameters which should be determined. According to the property of the quintic B-spline (2.1), we can get the nodal values as follows:

where $j = 0, 1, 2, \cdots, N_x$.

For two dimensions, an approximation solution $U_{N}(x, y, t)$ could be expressed as

where $\delta_{jk}(t)$ are unknown time dependent parameters. According to (2.1), we can get the nodal value $U_N(x_j, y_k, t)$ as

where $j = 0, 1, 2, \cdots, N_x$ and $k = 0, 1, 2, \cdots, N_y$. Unfortunately, Eq (2.9) is complicated, and it is two-dimensional which goes against the dimensionality reduction by the splitting method.

For simplicity, new notations $\tilde\delta_{j, k}(t)$ and $\hat\delta_{j, k}(t)$ are introduced as follows

So Eq (2.9) could be rewritten as

or

Obviously, either Eq (2.12) or Eq (2.13) is more concise than Eq (2.9). The advantage of the new writing will also be shown in the following sections.

For three dimensions, an approximation solution $U_{N}(x, y, z, t)$ could be expressed as

where $\delta_{jkl}(t)$ are unknown time dependent parameters.

Denote

where the first subscript $j$ related to x direction is fixed. Similarly, $\hat\delta_{j, k, l}(t)$ related to y direction and $\check\delta_{j, k, l}(t)$ related to z direction could be defined, where the second subscript $k$ and the third subscript $l$ are respectively fixed. Therefore, the nodal value $U_N(x_j, y_k, z_l, t)$ could be written as

3.   Numerical methods

Firstly, the second-order Strang splitting ^[9,10] is applied, and Eq (1.1) could be separated as

So Eq (1.1) could be approximately solved within $t\in [t_n, t_{n+1}]$ by solving Eqs (3.1)–(3.3) in sequence.

Multipling Eq (3.1) by $\bar u$ and taking the imaginary part, it follows that

Thus $|u|^2$ in Eq (3.1) is independent of $t$. So one can take $|u({\bf x}, t)|^2 = |u({\bf x}, t_n)|^2$. Consequently, it follows from Eq (3.1) that

where $t\in [t_n, t_{n+1}]$. Equation (3.3) could be solved similarly.

3.1.   One-dimensional (1D) scheme

For $d = 1$, Eq (3.2) could be written as

Applying the Crank-Nicolson scheme within $t\in [t_n, t_{n+1}]$, one obtains

where $U_N^n(x)$ is the approximation of $u(x, t_n)$. Taking $x = x_j$, it follows from Eqs (2.3) and (2.5) that

where $r_x = \tau /h_x^2$.

Consequently, Eqs (3.1)–(3.3) for $d = 1$ could be solved in sequence as follows

where $j = 0, 1, 2, \cdots, N_x$, and $n = 0, 1, 2, \cdots, N_t-1$. $U_N^{n+1, 1}(x)$ and $U_N^{n+1, 2}(x) = \sum_{j = -2}^{N_x+2}\delta _j^{n+1, 2}B_j(x)$ are intermediate results.

3.2.   Two-dimensional (2D) scheme

For $d = 2$, Eq (3.2) could be written as

The first-order Lie splitting ^[10] is applied, and one has

As the operators $\frac {\partial ^2}{\partial x^2}$ and $\frac {\partial ^2u}{\partial y^2}$ are commutable, there is no splitting error here.

Similar to Eq (3.5), one can obtain from Eq (3.10) that

where $U_N^n(x, y)$ is the approximation of $u(x, y, t_n)$. Taking $(x, y) = (x_j, y_k)$, it follows from Eqs (2.13) and (2.5) that

The above equation is concide benefiting from the new notation $\tilde\delta_{j, k}^{n+1}$. Moreover, it forms a $(N_x+1)$ system for each value of $k$ which avoids solving a $(N_x+1)\times (N_y+1)$ system directly. That is the reason why the author applies the splitting method in this paper.

Similarly, it follows from Eq (3.11) that

where $r_y = \tau /h_y^2$.

Consequently, Eqs (3.1)–(3.3) for $d = 2$ could be solved approximately as follows

where $j = 0, 1, 2, \cdots, N_x, \ k = 0, 1, 2, \cdots, N_y$, and $n = 0, 1, 2, \cdots, N_t-1$. $U_N^{n+1, 1}(x, y)$ and

are intermediate results.

3.3.   Three-dimensional (3D) scheme

For $d = 3$, Eq (3.2) could be written as

It follows from the Lie splitting ^[10] that

Similar to the 2D case, Eqs (3.1)–(3.3) for $d = 3$ could be solved approximately as follows

where $j = 0, 1, 2, \cdots, N_x, \ k = 0, 1, 2, \cdots, N_y, \ l = 0, 1, 2, \cdots, N_z$, and $n = 0, 1, 2, \cdots, N_t-1$. $U_N^{n+1, 1}(x, y, z)$ and

are intermediate results.

4.   Solvability, conservation and linear stability

In this section, solvability, conservation and linear stability are considered for the numerical methods. The 1D scheme is discussed in details, and the results could be extended to the 2D and 3D ones similarly.

From Eq (3.7), there are $N_x+1$ equations with $N_x+5$ unknows. So four additional equations are required. It follows from the periodic boundary condition (1.3) with $d = 1$ that

Taking $x = x_L$, one reaches

by using Eqs (2.3)–(2.7). It follows from Eqs (4.1)–(4.5) that

Therefore, Eqs (3.7) combined with (4.6) could be rewritten as

where

and

Lemma 4.1. The circulant matrices $A$ and $B$ respectively have eigenvalues as follows

where $j = 0, 1, 2, \cdots, N_x-1$.

Proof. The eigenvalues could be calculated directly (see ^[11] and the reference therein). □

Theorem 4.1. The solution of 1D SS5BC schemes (3.6)–(3.8) exists and is unique.

Proof. The coefficient matrix of scheme (3.7) is $M = iA+10\alpha r_xB$. Using Lemma 4.1, the eigenvalues of $M$ are

where $j = 0, 1, 2, \cdots, N_x-1$. All the eigenvalues are nonzero, so $M$ is invertible and the solution of scheme (3.7) exists and is unique. Moreover, Eqs (3.6) and (3.8) are respectively obtained from Eqs (3.1) and (3.3) exactly. Therefore, the theorem is proved. □

Similarly, one can prove that the solution of 2D SS5BC scheme (3.12)–(3.15) or 3D scheme (3.20)–(3.24) also exists and is unique.

Lemma 4.2. For any $N\times N$ real symmetric matrix $C$ and any complex vector $\delta = \left(\delta _1, \delta _2, \cdots, \delta _N\right)^\top$, $\delta ^{\mathrm{H}}C\delta$ is real, where $\delta^\mathrm{H}$ is the conjugate transpose of $\delta$ and $N$ is a positive integer.

Proof. By matrix operation, one has

where $c_{jk}$ is the element of the matrix lying on the intersection of the $j$th row and the $k$th column of $C$. Since $C$ is a real symmetric matrix, $\delta ^{\mathrm{H}}C\delta$ is real. □

Theorem 4.2. The 1D SS5BC schemes (3.6)–(3.8) conserves the discrete quantity

Proof. It follows from Eqs (3.6) and (3.8) that $\left|U_N^{n+1, 1}(x)\right| = \left|U_N^n(x)\right|$ and $\left|U_N^{n+1}(x)\right| = \left|U_N^{n+1, 2}(x)\right|$. So

where $n = 0, 1, 2, \cdots, N_t-1$.

Equation (4.7) is rewritten as

Multiplying Eq (4.12) by $\left[A\left(\delta ^{n+1, 2}+\delta ^{n+1, 1}\right)\right]^\mathrm{H}$ and taking the imaginary part, one has

where $A^\mathrm{H}$ is the conjugate transpose of $A$, and $A^\mathrm{H} = A$ since $A$ is real symmetric. It is obvious that

As $A^2$ is real symmetric and $\left(\delta ^{n+1, 2}\right)^\mathrm{H}A^2\delta ^{n+1, 1}$ and $\left(\delta ^{n+1, 1}\right)^\mathrm{H}A^2\delta ^{n+1, 2}$ are conjugate to each other, one gets

By applying Lemma 4.2, one obtains from Eq (4.14) that

Moreover, as $AB$ is a real symmetric matrix, one can obtain from Lemma 4.2 that

It follows from Eqs (4.13), (4.15) and (4.16) that

i.e.,

where

is used. Therefore, it follows from Eqs (4.11) and (4.17) that

Thus, Eq (4.10) is reached. □

Similarly, the 2D SS5BC schemes (3.12)–(3.15) and the 3D schemes (3.20)–(3.24) also respectively preserve the discrete quantity as follows

Next, the linear stability of the SS5BC scheme is considered.

Theorem 4.3. The 1D SS5BC schemes (3.6)–(3.8) is unconditionally stable.

Proof. Substituting $U_N^{n} = \xi ^{n}e^{i\beta _1x}$ and $U_N^{n+1, 1} = \xi ^{n+1, 1}e^{i\beta _1x}$ into Eq (3.6), one gets $\xi ^{n+1, 1} = G_1\xi ^{n}$, where

Similarly, it follows from Eq (3.8) that $\xi ^{n+1} = G_3\xi ^{n+1, 2}$, where

Equation (3.7) chould be rewritten as

Plugging $U_N^{n+1, 2} = \xi ^{n+1, 2}e^{i\beta _1x}$ into the above equation, one has $\xi ^{n+1, 2} = G_2\xi ^{n+1, 1}$, where

So $\xi ^{n+1} = G\xi ^n$, where $G = G_3G_2G_1$. Thus, the 1D SS5BC scheme is unconditionally stable, since $|G| = 1$. □

Similarly, one can obtain that the 2D and 3D SS5BC schemes are also unconditionally stable.

5.   Numerical experiments

Denote

be the maximum error. For convenience, we take $N = N_x = N_y = N_z$ and $h = h_x = h_y = h_z$. The convergence order is approximated as

where $||e||_\infty(h_1)$ is the maximum error corresponding to the step size $h_1$.

5.1.   Test 1

Take $\alpha = \frac {1}{2}$ and $\beta = -1$. For $d = 1$, Eq (1.1) has an exact solution

with

For $d = 2$, Eq (1.1) has an exact solution

with

And Eq (1.1) with $d = 3$ has an exact solution

with

The initial condition (1.2) could be given according to the exact solution. Take ${\bf x}\in [0, 1]^d$ and $T = 1$. As $V({\bf x}, t) = V({\bf x})$ independent of $t$, we have

in the nonlinear schemes. And the linear schemes are solved by the Douple-Sweep method.

It follows from the results in ^[12] that

where $S_3$ and $S_5$ are respectively the cubic spline and quintic spline associated with the function $f$. So the accuracy order of the SS3BC scheme in ^[7] and the SS5BC scheme in this paper might be $O(\tau ^2+h^2)$ and $O(\tau ^2+h^4)$, respectively.

The above three examples are applied to verify the convergence order of the proposed SS5BC scheme and the SS3BC scheme in ^[7]. Thus we take $\tau = h^2$ for the SS5BC scheme, and $\tau = h$ for the SS3BC one. The numerical results are listed in Tables 1 and 2, respectively. Table 1 shows that all the 1D, 2D and 3D SS5BC schemes are convergent with second-order in time and fourth-order in space, and Table 2 shows that all the 1D, 2D and 3D SS3BC schemes are convergent with second-order both in time and in space. So the proposed SS5BC schemes do improve the accuracy order compared with the SS3BC schemes ^[7], which is the aim of this paper.

In the above tests, the SS5BC scheme possesses higher order of accuracy than the SS3BC one. However, the coefficient matrix is five diagonal for the former scheme while tridiagonal for the later one, and it seems that the SS5BC scheme might cost more. For further comparison, the computing time is compared between the two kinds of methods. For the above three examples, one attaines $||e||_\infty < 0.002$ when the step sizes are free ^[13]. Take ${\bf x}\in [0, 1]^d$ and $T = 0.5$. Table 3 shows that the SS5BC methods are more efficient than the SS3BC ones since the formers cost less time.

The conserved quantity $Q(t)$ is respectively simulated for $d = 1, 2, 3$, seeing Eqs (4.10), (4.18) and (4.19). Taking ${\bf x}\in [0, 1]^d, h = 0.05, \tau = 0.02$ and $T = 5$, the values of $Q^n$ at $t = 0, 1, 2, 3, 4, 5$ are listed in Table 4. It follows from the table that the 1D, 2D and 3D SS5BC schemes remain the conserved quantity $Q^n$ very well.

5.2.   Test 2

Take $\alpha = \frac {1}{2}$ and $\beta = -1$. For $d = 1$, Eq (1.1) has a soliton solution

with

For $d = 2$, Eq (1.1) has a soliton solution

with

For $d = 3$, Eq (1.1) has a soliton solution

with

The initial condition $u({\bf x}, 0)$ in Eq (1.2) could be given according to the exact solutions. Since $V({\bf x}, t)$ can not be integrated exactly, one may apply the following approximation

in the schemes. Take ${\bf x}\in[-5, 5]^d, t\in[0, 1]$, and $\tau = h^2$. The error and convergence order are listed in Table 5. It shows that the 1D, 2D and 3D SS5BC schemes are convergent with order $O(\tau ^2+h^4)$.

Taking ${\bf x}\in [-5, 10]^d, h = 0.1, \tau = 0.02$ and $T = 5$, the values of $Q^n$ at $t = 0, 1, 2, 3, 4, 5$ are listed in Table 6. The table shows that the 1D, 2D and 3D SS5BC schemes also keep the conserved quantity $Q^n$ well in this test.

Numerical simulations of $|u|^2$ computed by the SS5BC schemes are plotted in Figures 1–3. In Figure 1, the 1D solitary wave transfers from the left to the right as time inceases which is in accordance with the exact one. In Figure 2, the 2D solitary waves at $t = 0, 1, 2, 3, 4, 5$ are plotted which also agree with the exact ones. For 3D, profiles of $|u|^2$ at $x = 0, y = 2$ and $z = 4$ are respectively plotted in Figure 3, which are still consistent with the exact ones.

6.   Numerical applications

Taking $\alpha = \frac {1}{2}$ and the trap potential

the NLS equation (1.1) is known as the Gross-Pitaevskii (GP) equation, which is usually used to model the properties of a Bose-Einstein condensate (BEC) at extremely low temperatures ^[1,14]. For normalization ^[14], the conserved quantity in Eq (1.4) is required as $Q(t) = 1$ for each $t\in [0, T]$.

For $d = 1$, the initial condition (1.2) is chosen as

The condensate width ^[14] of 1D BEC is numerically approximated as

where

Take $x\in[-8, 8], T = 10, h = 0.2, \tau = 0.02$, and $\beta = -3$ in Eq (1.1). The approximated condensate density $|u|^2$ and the condensate width $\sigma ^n$ are plotted in Figure 4. And the conserved quantity $Q^n$ listed in Table 7 is about 1, which simulates $Q(t) = 1$ very well.

For $d = 2$, the condensate widths ^[14] along the x- and y-axes are respectively approximated as

where

Taking $(x, y)\in[-8, 8]^2, T = 10, h = 0.2, \tau = 0.02$ and $\beta = -2$, two cases are considered as follows:

Case I. Let $\gamma _y = 1$. The initial condition (1.2) is chosen as

Case II. Let $\gamma _y = 2$. The initial condition is taken as

The approximated condensate widths of Cases I and II are plotted in Figure 5. $\gamma _y$ is a ratio of the trap frequencies in x- and y-direction ^[14]. As $\gamma _y = 1$ in Case I, the trap potential (6.1) and also the condensate are isotropic. So the condensate widths along the x- and y-axes should be the same for Case I, which is shown in Figure 5 numerically. The conserved quantities $Q^n\approx 1$ of the two cases are both given in Table 8.

For $d = 3$, the condensate widths ^[14] along the x-axis is numerically approximated as

where

Similarly, the condensate widths along the y- and z-axes could be approximated respectively. Take $(x, y, z)\in[-8, 8]^3, T = 5, h = 0.2$ and $\tau = 0.02$. The initial condition (1.2) is chosen as

and the following two cases are considered:

Case I. Let $\beta = -0.1, \gamma _y = 2, \gamma _z = 4$.

Case II. Let $\beta = -1, \gamma _y = 1, \gamma _z = 2$.

In Figure 6, the approximated condensate widths of Cases I and II are plotted. As $\gamma _y = 1$ in Case II, the condensate is symmetric in x- and y-directions, which is shown in Figure 6 that $\sigma _x^n$ equals $\sigma _y^n$. In Table 9, the conserved quantities $Q^n\approx 1$ are listed.

7.   Conclusions

In this paper, the SS5BC schemes are proposed for the one-dimensional and multi-dimensional NLS equations. The new notations are introduced for the 2D and 3D equations in order to make the schemes more brief and accomplishable. The solvability, conservation and linear stability are discussed for the methods. Variable numerical experiments and applications are carried out to prove that the present schemes are reliable and efficient. The convergence order, conserved quantity and solitary wave are verified numerically. Finally, the SS5BC methods are applied to study BECs. It is worth to say that the skills of analysis in this paper could also be applied to the SS3BC schemes, and advanced theoretical analyses are still open which are expected in the future work.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is supported by the National Natural Science Foundation of China under Grant No. 11701280.

Conflict of interest

The author declares no conflict of interest in this paper.