A numerical investigation of nonlinear Schrödinger equation using barycentric interpolation collocation method

1.   Introduction

In 1926, the Austrian physicist Schrödinger proposed the famous Schrödinger equation ^[1] in order to solve the uncertainties of the microscopic world and to reveal the fundamental laws of matter motion in the microphysical world, serving as a powerful tool to deal with all non-relativistic problems in atomic physics, which is equivalent to the role played by Newton's laws in classical mechanics. The Schrödinger equation plays a very wide and important role in the fields of atomic, molecular, solid state physics, nuclear physics, and chemistry, etc. However, it is usually difficult to find the exact solutions of the Schrödinger equation, which is replaced by the numerical solution. Therefore, the numerical method of nonlinear Schrödinger equation has attracted wide attention of scholars.

During the past two decades, enormous works have been devoted to numerical simulations and numerical analysis for the Schrödinger equation. In the literature ^[2], a local discontinuous galerkin finite element method was proposed to solve the nonlinear Schrödinger equation. Wang et al. ^[3] proposed two-grid finite element method to solve the nonlinear Schrödinger equation. Liao et al. ^[4] developed the finite difference explicit scheme and gave the error analysis in the sense of maximum norm. Sun ^[5] introduced an explicit difference scheme on an infinite domain using artificial boundary conditions. In the literature ^[6], a pseudo-spectral Fourier method was introduced to solve the Schrödinger equation in optical fibers. Mackenzie et al. ^[7] presented an hr-adaptive method for the cubic nonlinear Schrödinger equation. Zhai et al. ^[8] developed the strang splitting method for space fractional nonlinear Schrödinger equations. Feng et al. ^[9] focused on high-order conserving SAV-Gauss collocation finite element schemes for nonlinear Schrödinger equation. Hu et al. ^[10] introduced two-grid finite element algorithms for two-dimensional nonlinear Schrödinger equation with wave operator. Iqbal et al. ^[11] proposed a cubic B-spline Galerkin method for the coupled nonlinear Schrödinger equation. Cai et al. ^[12] developed the integral factor method, diagonal matrix technique and FFT algorithm to solve the high-dimensional Schrödinger equation with damping.

Based on previous research, different from the traditional element-based numerical method, we adopt a new meshless method called barycentric interpolation collocation method, which can effectively avoid the cumulative error caused by the difference scheme. {Compared with the spectral method, the barycentric interpolation collocation method can deal with irregular domains. As a novel meshless method, barycentric interpolation method is a good choice for dealing with polynomial interpolations since it is accurate highly, fast, stable and easy on program implementation.} Recently, barycentric interpolation collocation method has been widely used in elastic mechanics, microwave technology, fluid mechanics and many other applications. Many researchers extended this method to solve various types of differential equations, such as the stress separation in photoelasticity ^[13], Fredholm integral equation ^[14], heat conduction equation ^[15], two-dimensional Fractional Cable Equation ^[16], Allen-Cahn equation ^[17], Volterra integral equation ^[18], time-fractional order telegraph equation ^[19], Burgers equation ^[20], and so on.

Based on the above works, we focus on a fully discrete scheme for the Schrödinger equation, which is a second-order Crank-Nicolson scheme in time combined with the barycentric interpolation collocation method in space. Moreover, we will give consistency analysis of the semi-discretized scheme in space. This paper will further promote the development of high-precision numerical algorithms for solving nonlinear Schrödinger equations, and provide a new way for the study of nonlinear science.

The rest of the paper is structured as follows: In section 2, temporal discretization and the barycentric interpolation collocation method are introduced. In section 3, consistency analysis of semi-discretized system are carried out in detail for the Schrödinger equation. In section 4, fully discretized scheme based on Crank-Nicolson is presented. Numerical examples verify the accuracy and efficiency of the algorithm in section 5 and some conclusions are given in section 6.

2.   Description of the method

2.1.   Schrödinger equations

The initial boundary value problem of nonlinear Schrödinger equation is considered as follows.

with homogeneous Dirichlet boundary condition, {where $\Omega$ is a bounded domain in $R^{d}(d = 1, 2)$. It is easy to show that this equation conserves the total mass and the total energy ^[21].

2.2.   Temporal discretization

The time interval $[0, T]$ is discretized into $N$ subintervals of equal length $\tau = \frac{T}{N}$ such as $0 = t_0 < t_1 < \cdots < t_N = T$, where $t_n = n\tau$, $n = 0, 1, \cdots, N.$

Let $t_{n+\frac{1}{2}} = \dfrac{1}{2}\left(t_n+t_{n+1}\right)$, consider Eq (2.1) at $(x, t_{n+\tfrac{1}{2}})$, it can be modified as

where $0\le n\le N-1$.

From Eq (2.4), the first term can be modified as

The third term can be modified as

From Eqs (2.4–2.6), we get

By ignoring the high order term $o(\tau^2)$, we get

Thus we get Crank-Nicolson scheme, which is a nonlinear implicit scheme, and its numerical solution is second order in time.

First, we prove the lemma as follows:

Lemma 1. For $u^{n+\frac{1}{2}} = \frac{u^{n+1}+u^{n}}{2}\in C$, it holds that $\langle u^{n+\frac{1}{2}}, u^{n+\frac{1}{2}}\rangle \in {\rm R}$.

Proof.

Because $|u^{n+1}|^2+|u^{n}|^2\in {\rm R}$, so we only need to prove $u^{n+1}\overline{u}^{n}+u^{n}\overline{u}^{n+1}\in {\rm R}$. Suppose $u^{n+1} = a^{n+1}+ib^{n+1}$, $\overline{u}^{n+1} = a^{n+1}-ib^{n+1}$, then we can get

Hence, $\langle u^{n+\frac{1}{2}}, u^{n+\frac{1}{2}}\rangle \in {\rm R}.$

For the completeness of this paper, the mass and energy conservation of semi discrete scheme in time will be given here.

For the discrete inner product with $u^{n+\frac{1}{2}}$ on both sides of Eq (2.8)

Take the imaginary part of Eq (2.11), there are

Hence, from the second term of the Eq (2.12), using the integration and lemma 1, we can deduce

Using lemma 1, from the third term and the forth term of the Eq (2.12), we can get

Consequently, the first term$\dfrac{1}{2\tau}(||u^{n+1}||^2-||u^{n}||^2) = 0$. In other words, $||u^{n+1}||^2 = ||u^{n}||^2$. Hence, the total mass is conserved.

For the discrete inner product with $(u^{n+1}-u^n)$ on both sides of the Eq (2.8), and take the real parts on both sides of the equation

And we can deduce

Hence, we can obtain

In other words, the total energy is conserved.

2.3.   Barycentric Lagrange interpolation

Suppose there are $n+1$ distinct interpolation nodes $x_j$ in the direction of $x$, where corresponds to a set of real numbers $y_j$ $(j = 0, 1, \cdots, n)$ at the same time. Therefore, there exists a unique interpolation polynomial whose degree is not exceeding m, satisfying $p(x_j) = y_j$ $(j = 0, 1, \cdots, n)$. As it known to us, such polynomial $p(x)$ is unique and can be written in Lagrange form as

where $\pi_{j}(x)$ are Lagrange interpolation basis functions, satisfying the properties

Let

Defining the barycentric weights by

Then the interpolation basis function $\pi_{j}(x)$ can be expressed as

Inserting Eq (2.26) into Eq (2.21), we can get

The barycentric Lagrange interpolation formula can be achieved from Eqs (2.23–2.26), namely

The barycentric Lagrange interpolation is ill-formed for isometric nodes. However, choosing a node family with a density ratio of $(1-x^2)^{-\frac{1}{2}}$ can make it numerically stable. The simplest node distribution is the Chebyshev point family. Therefore, in order to ensure the numerical stability of the barycenter Lagrange interpolation, the second Chebyshev point family is adopted in this paper: $x_j = \cos(\frac{j}{m}\pi)$ ($j = 0, 1, \cdots, m$).

2.4.   Differential matrix

The derivative of $p(x)$ defined as Eq (2.28) with respect to $x$ as

The first and second order differentiation matrices can be obtained by the following formula

By mathematical induction, we have the following $v$-order differential matrix

3.   Semi-discretized system and consistency analysis

3.1.   Semi-discretized system based on barycentric interpolation collocation method

The spatial region $\Omega = [a, b]\times[c, d]$ of the Schrödinger equation is discretized into $(M+1)\times (N+1)$ second-type Chebyshev nodes, namely $a = x_{0} < x_{1} < \cdots < x_{M} = b$, $c = y_{0} < y_{1} < \cdots < y_{N} = d$. Then the $(M+1)\times (N+1)$ Chebyshev nodes on the region are $(x_{i}, t), i = 0, 1, \cdots, M$. Denoting $u(x_i, y, t) = u_i(y, t)$, $v(x_i, y) = v_i(y), i = 0, 1, \cdots, M$, barycentric interpolation approximation function can be modified as

where $Q_{l}(x)$ is the interpolation basis function of the barycentric interpolation. Substitute the above formula into the equation, and let the equation hold on nodes $x_{0}, x_{1}, \cdots, x_{M}$, the ordinary differential equation system can be obtained, namely

where $Q''_{l}(x_i) = \dfrac {{\mathrm{d}}^2 Q_{l}(x_i)}{\mathrm{d} x^2} = C_{il}^{(2)}$, $G(u) = \beta|u|^2u$. Eq (3.2) can be written in matrix form, namely

where C = $\left({\begin{array}{*{20}{c}} {C_{00}^{(2)}} & \cdots & {C_{0M}^{(2)}}\\ \vdots & {} & \vdots \\ {C_{M0}^{(2)}} & \cdots & {C_{MM}^{(2)}} \end{array}} \right)$.

Denote $u_i(y_j, t) = u_{ij}(t)$. Similarly, $u_i(y, t)$ can be written in barycentric interpolation form

where $\eta_{r}(y)$ is the basis function in barycentric interpolation on the direction of $y$.

Inserting Eq (3.4) into Eq (3.3) and making Eq (3.3) be identical at the nodes $y_{0}, y_{1}, \cdots, y_{N}$, we get

where ${\eta_r ''}(y_{j}) = \dfrac {{\mathrm{d}}^2 \eta_{r}(y_j)}{\mathrm{d} y^2} = D_{j r}^{(2)}$.

Let

Eq (3.5) can be rewritten in the following matrix form

In the above formula: the symbol $\otimes$ represents the Kronecker product of the matrix, $\boldsymbol{C}^{(2)}$ is the barycentric interpolation second-order differential matrix about the nodes $x_{0}, x_{1}, \cdots, x_{M}$; $\boldsymbol{D}^{(2)}$ is the barycentric interpolation second-order differential matrix about the nodes $y_{0}, y_{1}, \cdots, y_{N}$; $\boldsymbol{I}_{M}$ and $\boldsymbol{I}_{N}$ are the M-order identity matrix and N-order identity matrix.

3.2.   Consistency analysis

In this part, we present consistency estimates of the semi-discretized scheme with the collocation method. Let $p(x)$ is the Lagrange interpolation function approximating $u(x)$. According to theorem 3.1 in reference ^[19], we obtain the following approximation properties.

Lemma 2. For the error function defined above $u(x, t)$ and the function $u(x, t) \in {C^{m + 1}}[a, b]$

where function $e(x, t) = u(x, t)-p(x, t) = \frac{u^{(m+1)}(\xi)}{(m+1)!}\prod_{i = 0}^{m}(x-x_i)$, $C_{1}$, $C^{*}_{1}$, $C^{**}_{1}$ are constant independent of $x$, $e$ is natural logarithm and $h_{x}$ represents the length of the interval $\left[a, b\right]$.

Defining the error function as

In that way, we can get

In order to estimate the error of the barycentric interpolation Lagrange collocation method, an operator is defined

where$f(u)$ is equal to $\left| u\right| ^{2}$ in the equation.

Theorem 1. For $u(x, y, t) \in C^{m + 1}([a, b]\times [c, d]$, the error result is as follow

Proof.

where

For $R_1$, we have

Applying Lemma 2, we obtain

Similarly, for $R_{2}$ and $R_{3}$ we can get

Then

Assuming $u(x, y, t), u({x_m}, {y_n}, t)$ is bounded, we have

Substituting (3.12)–(3.17) into (3.11), this completes the proof.

It can be seen from the error estimate that the numerical scheme is exponentially convergent, and the order of the differential operator determines the order of convergence of the algebraic equation.

4.   Fully discretized scheme based on Crank-Nicolson method

In the section, we will solve the ODE system (3.6) by the Crank-Nicolson scheme. Partition the interval (0, T] into a uniform mesh with the time step $\tau = \dfrac{T}{l}$: $0 = t_{0} < t_{1} < \cdots < t_{l} = T$. Let ${\bf U}^k = {\bf U}(t_k), k = 0, 1, \cdots, l.$

Using Crank-Nicolson difference scheme in the time direction of the Eq (3.6), we can get

Where $\textbf{G}^{k} = \dfrac{\beta}{2}(|{\textbf{U}^{k+1}}|^{2}+|{\textbf{U}^{k}}|^{2})\textbf{U}^{k+\frac{1}{2}}$.

We can get the following scheme:

where $F^{k} = \left[iI- \dfrac{\tau}{2}({{\bf{C}}^{(2)}} \otimes {{\bf{I}}_N}) - \dfrac{\tau}{2}({{\bf{I}}_M} \otimes {{\bf{D}}^{(2)}}) -\dfrac{\tau}{2}diag(v({\bf{x}})) +\dfrac{\tau \beta}{4}(|{{\bf{U}}^{k}}|^{2}+|{{\bf{U}}^{k+1}}|^{2}) \right]{{\bf{U}}^{k}}$, $k = 0, 1, \cdots, l - 1$.

Using Newton iterative method for the nonlinear term $\dfrac{\tau}{2}|{{\bf{U}}^{k+1}}|^{2}U^{k+1}$, we can obtain

5.   Numerical experiments

In this section, several examples are provided to verify the high accuracy and numerical stability of the proposed scheme.

For convenience, we define the maximum absolute and relative error symbols as follows:

where $U_{i}$ denote the numerical solution, $u_{i}$ denote the exact solution respectively, $\parallel \cdot \parallel_{\infty}$ is the $L^{\infty}$ norm.

5.1.   Example 1

We consider the following problem on $\Omega = [0, 1] \times (0, 1]$ as

where

Fix $M = 16$, and time step $\tau$ is varied. Table 1 shows that the convergence order is second order in time.

Table 2 shows the maximum absolute error and the maximum relative error obtained by the second-order finite difference scheme and the barycentric interpolation collocation scheme. Second-order finite difference scheme ^[21] for the Schrödinger equation is as follows.

It can be seen that barycentric interpolation collocation method can achieve higher accuracy with less nodes. For example, when the temporal step $\tau = 10^{-4}$, barycentric Lagrange interpolation only needs to take 8 points in space to achieve the accuracy of $10^{-6}$, but the second-order finite difference scheme requires 256 points to achieve the accuracy of $10^{-5}$.

By applying the Fourier transform to both sides of (2.8) and the definition of the Laplacian given by (5.8), one gets

where $\hat{u}^{n}_{j}$ is the $j$ th Fourier coefficient.

Remark. Suppose the Laplacian $(-\Delta)$ has a complete set of orthonormal eigenfunctions $\{\varphi_{j}\}$ satisfying standard boundary conditions on a bounded region $\Omega \in R^{n}$, with corresponding eigenvalues $\lambda_{j}$, i.e., $(-\Delta)\varphi_{j} = \lambda_{j} \varphi_{j}$ on $\Omega$, and let

Then, for any $u \in \mathcal{U}$,

For detailed introduction of Fourier spectral method, please refer to reference ^[22].

Table 3 shows the maximum absolute error and the maximum relative error obtained by the spectral method. Compared with spectral method, we can find that the highest accuracy of the spectral method is only $10^{-9}$ in Table 3. However, the accuracy of barycentric interpolation collocation method is $10^{-10}$ when we take 12 points in space from Table 2. So our scheme is still superior.

5.2.   Example 2

The initial value condition is $u(x, 0) = sinx$, $\beta = 1,$ and the second term of Eq(2.1) is $\frac{1}{2} \Delta u$. The exact solution of this problem is $u(x, t) = exp(-\frac{3it}{2})sinx, 0 \leq x \leq 2\pi, 0 < t \leq 1.$

The errors in maximum norm and convergence rates are displayed in Table 4, we can find that the scheme is second order in time for 1D Schrödinger equation.

Table 5 gives the maximum error and relative error obtained by the barycentric interpolation collocation scheme and the second-order finite difference scheme. From Table 5, it can be found that the barycentric interpolation collocation scheme can achieve high accuracy by taking fewer points in space compared with the second-order finite difference scheme. Especially when the barycentric Lagrange interpolation collocation scheme takes 8 points, it achieves higher precision than the second-order finite difference scheme takes 128 points.

So we can draw images of numerical and exact solutions, as shown below. It can be seen from Figures 1 and 2 that the numerical solution obtained by the barycentric interpolation collocation scheme is consistent with the exact solution, indicating that the scheme is very stable.

The discrete mass and energy conservation properties of the Schrödinger equation are further verified. The space step is selected as $\tau = 10^{-4}$, the corresponding calculation interval is $[0, 2\pi]\times$ $[0, 1]$. It can be seen from Figures 3 and 4 that our collocation scheme preserves the total mass and energy in the discrete sense.

5.3.   Example 3

The initial value is $u_0(x) = 2sech(2x)exp(ix)$. Let the boundary conditions of the one-dimensional Schrödinger equation be $u(-10, t) = u(10, t) = 0$ and potential function is $v(x) = 0$, $\beta = 1.$

The solitary waveform under the initial condition $t = 0$ is shown in Figure 5. The time step is $\tau = 10^{-3}$ and the corresponding interval is $[-10, 10]\times[0, 5]$. The solitary waveforms of barycentric Lagrange interpolation at $t = 2$ and $t = 5$ are shown in Figure 6. Therefore, it can be seen that the solitary waveform of barycentric Lagrange interpolation is very stable, indicating that the numerical method is stable.

5.4.   Example 4

We consider the following problem on $\Omega = {[-1, 1]^2} \times (0, 1]$ as

where

The error comparison results of different discretization schemes in space can be obtained, as shown in Table 6. One can see that barycentric interpolation collocation method has high precision, as it only uses $9 \times 9$ mesh nodes can achieve the accuracy of second-order central difference approach with $80 \times 80$ mesh nodes.

Table 7 indicates our scheme is second order in time for 2D Schrödinger equation.

5.5.   Example 5

Consider 2D Schrödinger equation on $\Omega = {[0, 2\pi]} \times [0, 2\pi]$ with the following initial condition

where

The exact solution for the problem is $u = e^{-2ti} \sin x \sin y$.We take simulation parameters as

So we can draw images of numerical and exact solutions, as shown below. It can be seen from Figures 7 and 8 that the numerical solution obtained by the barycentric interpolation collocation scheme is consistent with the exact solution, indicating that the scheme is very stable. Then we need to verify again the conservation of energy and mass. It can be seen from the Figures 9 and 10 that our collocation scheme satisfies the discrete mass and energy conservation.

6.   Conclusions

In this paper, by a combination of the barycentric interpolation collocation method and Crank-Nicolson scheme, an efficient numerical scheme for Schrödinger equation is developed. In addition, the consistency analysis of the semi-discretized scheme is derived. Numerical examples show that our scheme is second order in time and convergent exponentially in space. The total mass and energy conservation in the discrete senses are also checked numerically. Compared with the finite difference method and spectral method, the barycenter interpolation collocation method can achieve high accuracy with fewer points. In future work, we plan to extend this method to coupled schrodinger equations, KdV equations and Klein-Gordon equations etc.

Acknowledgments

This work is in part supported by the NSF of China (No. 11701197), the Fundamental Research Funds for the Central Universities (No. ZQN-702), the Natural Science Foundation of Fujian Province (No. 2022J01308), the Key Laboratory of Intelligent Computing and Information Processing of Ministry of Education (Xiangtan University) (No. 2020ICIP03) and National Training Program of Innovation and Enterpreneurship for Undergraduates (No. 202110385029).

Conflict of interest

The authors declare that there are no conflicts of interest regarding the publication of this manuscript.