A unified framework of high order structure-preserving B-splines Galerkin methods for coupled nonlinear Schrödinger systems

Abstract

Using a general computational framework, we derive an optimal error estimate in the $L_2$ norm for a semi discrete method based on high order B-splines Galerkin spatial discretizations, applied to a coupled nonlinear Schrödinger system with cubic nonlinearity. A fully discrete method based on a conservative nonlinear splitting Crank-Nicolson time step is then proposed; and conservation of the mass and the energy is theoretically proven. To validate its accuracy in space and time, and its conservation properties, several numerical experiments are carried out with B-splines up to order 7.

Introduction

In this work we are interested in the numerical approximation of a general class of N-strongly coupled nonlinear Schrödinger systems using high order B-splines Galerkin methods as spatial discretization. The systems are of the form

$$i\frac{\partial {\psi }_m}{\partial t}=-{\alpha }_m\frac{{\partial }^2{\psi }_m}{\partial x^2}-f_m({\left|{\psi }_1\right|}^2,\dots ,{\left|{\psi }_N\right|}^2){\psi }_m,\text{ in }\mathbb{R}\times (0,T],$$

$${\psi }_m(x,0)={\psi }_m^o(x),\text{ in }\mathbb{R},$$

where the coefficients ${\alpha }_m\in \mathbb{R}$; and, there exists a differentiable function $F:{\mathbb{R}}^N\to \mathbb{R}$ such that for all $z\in {\mathbb{R}}^N$ we have $\mathrm{\nabla }F(z)={(f_1(z),f_2(z),\dots ,f_N(z))}^T$. Computationally, the system is solved on an open and bounded domain Ω; and to model the physical conditions ${\lim }_{|x|\to \infty }|{\psi }_m(x,\cdot )|=0$, periodic or homogeneous Dirichlet boundary conditions are used.

The model problem (1.1a), (1.1b) is strongly related to many important physical applications such as: interacting solitons in nonlinear optical fibers [27], electromagnetic wave propagation [21], high temperature plasma [3], etc. Furthermore, given the difficulty of obtaining closed form solutions for general nonlinear potentials, there is a high interest in the development of efficient and accurate numerical methods for these problems.

A B-splines Galerkin method of order p is a direct extension of the classical finite element method (FEM) that seeks an approximation in the space spanned by B-splines. For the maximum smoothness setting, this is the space of piecewise polynomials in ${\mathcal{C}}^{p-1}\left(\mathrm{\Omega }\right)$; consequently, they reproduce the smoothness of the exact solution, while retaining the same optimal order of convergence than the standard FEM but with a much smaller number of unknowns.

To put the present work into perspective let us briefly recall the previous work on B-splines Galerkin methods applied to nonlinear Schrödinger equations. Gardner and Gardner [14] proposed a quadratic B-splines Petrov-Galerkin method with piecewise constant weight functions and proved conservation of two invariants, for the semi-discrete formulation on a scalar nonlinear Schrödinger equation. Quadratic B-splines Galerkin methods were later presented by Dağ [10]; and, cubic and quintic B-splines were considered by Iqbal et al., [18], [17]; and more recently, in [16], [19] for 2-CNLS systems. In all these works, the classical Crank-Nicolson (CN) method was used as time discretization; however, it is well known that this method does not preserve the global energy of nonlinear Schrödinger equations with arbitrary nonlinear potentials.

Spline collocation methods have also been considered, e.g., [15], [23], [22]. These methods seek a spline approximation that satisfies the differential equation at some specific points, usually referred to as collocation points. Since the assembly of their discrete operators requires fewer number of operations, they have shown to be more efficient than standard Galerkin methods [25]; however, arbitrary distributions of collocation points can result in a reduction of their order of convergence: by one for p even; or by two for p odd. To circumvent this problem, one can apply the collocation method to a perturbation of the original differential equation [1], [9], technique known as modified collocation methods. Another approach consists of taking the collocation points as the zeros of some orthogonal polynomial [4], the so-called orthogonal collocation methods (OSC). However, the former methods are limited to uniform meshes; moreover, it is not clear if the perturbed equation preserves any invariants. To obtain optimal order of convergence, OSC methods are limited to polynomials in ${\mathcal{C}}^1$, as a result, they have a higher number of degrees of freedom (DOF). Specifically, the ratio between the DOF for spline approximations in ${\mathcal{C}}^1$ and ${\mathcal{C}}^{p-1}$ is $N(p-1)+2:N+p+1$ where N is the number of cells.

For periodic or homogeneous boundary conditions the model problem (1.1a)–(1.1b) exhibits the following invariants:

$$\begin{gathered} \mathcal{M}(t):=\sum _{m=1}^N\underset{\mathrm{\Omega }}{\int }|{\psi }_m{|}^2\text{dx}\text{ and } \\ \mathcal{E}(t):=\frac{1}{2}\underset{\mathrm{\Omega }}{\int }(\sum _{m=1}^N{\alpha }_m{\left|\frac{\partial {\psi }_m}{\partial x}\right|}^2-F(|{\psi }_1{|}^2,\dots ,|{\psi }_N{|}^2))\text{dx}, \end{gathered}$$

usually referred to as total mass $\mathcal{M}(t)$ and energy $\mathcal{E}(t)$. It is well known that conservation of these two particular invariants plays a crucial role in the design of numerical schemes to approximate the solution of nonlinear Schrödinger type of equations. For instance, as discussed in [24], the lack of conservation can lead to the undesirable phenomenon known as blow-up. Numerical stability in the $L_2$-norm is a direct consequence of mass conservation; and, in addition conservation of both invariants is a key ingredient for the analysis of error estimates.

The novelty and main contribution of the present work is the convergence analysis of a conservative semi discrete method based on B-splines of arbitrary order, applied to a well known Schrödinger system with cubic nonlinearity. Although, the main ideas to prove conservation of the physical invariants have been recently presented by the authors in [8], using the local discontinuous Galerkin method as spatial discretization, no error analysis was carried out. By considering the spaces spanned by B-splines not only we overcome some technical difficulties inherent in discontinuous polynomial spaces, but also, compared to the classical and discontinuous Galerkin methods, significantly fewer degrees of freedom are required, hence a less memory consuming and more efficient method is obtained. To the best of our knowledge, conservation of invariants has not been proven for fully discrete schemes based on B-splines Galerkin discretizations, and no error estimates for the B-splines Galerkin method applied to this kind of problems have been obtained before. Numerical experiments presented so far show an evident loss of conservation for at least one of the invariants, especially the energy which is a non quadratic quantity. Owing to [8], we derive a fully discrete conservative scheme for a general class of N-CNLS systems with the following properties:

• Using appropriate choices of knots depending on the boundary conditions of the problem, conservation of both, mass and energy, is obtained by a time advancing scheme based on a conservative splitting method consisting of a nonlinear Gauss-Seidel sweep that uses a modified Crank-Nicolson scheme on each component. As a result the global system is decomposed into a set of smaller nonlinear problems.

• High order accuracy in time is achieved by using the well known theory of composition methods.

• Spatial accuracy is inherited from the approximation properties of high order B-splines.

After a brief review of B-splines, Section 2 is devoted to the semi-discrete formulation of the B-splines Galerkin method applied to N-CNLS systems. Conservation of the mass and energy is discussed and an optimal error estimate is presented for a Schrödinger system with cubic nonlinearity. The description of a fully discrete method and its conservation properties are presented in Section 3. To validate the conservation properties and accuracy of the proposed scheme, a series of numerical experiments are carried out in Section 4.