A numerical study of adaptive space and time discretisations for Gross-Pitaevskii equations.

As a basic principle, benefits of adaptive discretisations are an improved balance between required accuracy and efficiency as well as an enhancement of the reliability of numerical computations. In this work, the capacity of locally adaptive space and time discretisations for the numerical solution of low-dimensional nonlinear Schrödinger equations is investigated. The considered model equation is related to the time-dependent Gross-Pitaevskii equation arising in the description of Bose-Einstein condensates in dilute gases. The performance of the Fourier-pseudo spectral method constrained to uniform meshes versus the locally adaptive finite element method and of higher-order exponential operator splitting methods with variable time stepsizes is studied. Numerical experiments confirm that a local time stepsize control based on a posteriori local error estimators or embedded splitting pairs, respectively, is effective in different situations with an enhancement either in efficiency or reliability. As expected, adaptive time-splitting schemes combined with fast Fourier transform techniques are favourable regarding accuracy and efficiency when applied to Gross-Pitaevskii equations with a defocusing nonlinearity and a mildly varying regular solution. However, the numerical solution of nonlinear Schrödinger equations in the semi-classical regime becomes a demanding task. Due to the highly oscillatory and nonlinear nature of the problem, the spatial mesh size and the time increments need to be of the size of the decisive parameter [Formula: see text], especially when it is desired to capture correctly the quantitative behaviour of the wave function itself. The required high resolution in space constricts the feasibility of numerical computations for both, the Fourier pseudo-spectral and the finite element method. Nevertheless, for smaller parameter values locally adaptive time discretisations facilitate to determine the time stepsizes sufficiently small in order that the numerical approximation captures correctly the behaviour of the analytical solution. Further illustrations for Gross-Pitaevskii equations with a focusing nonlinearity or a sharp Gaussian as initial condition, respectively, complement the numerical study.

Introduction

General scope. In various fields, dynamical processes are modelled by time-dependent nonlinear partial differential equations, and, in most cases, the complexity of the considered applications requires the use of numerical methods. In order to enhance the reliability of numerical simulations and to facilitate the efficient computation of approximations which capture correctly the quantitative or qualitative solution behaviour, it is highly desirable to utilise algorithms involving an automatic selection of the space and time grid.

Main objective. In the present work, our main interest is to exploit the capacity of locally adaptive space and time discretisations. As scope of applications, we focus on low-dimensional nonlinear Schrödinger equations such as systems of coupled time-dependent Gross–Pitaevskii equations arising in the mathematical description of multi-component Bose–Einstein condensates in dilute gases at temperatures significantly below the critical temperature of the condensate, see [25,43,46]. Numerical illustrations include Gross–Pitaevskii equations with a defocusing or focusing nonlinearity, respectively, as well as a breathing-like solution with a sharp Gaussian as initial condition. Moreover, in regard to Schrödinger equations in the semi-classical regime, we study situations where an additional parameter reflecting the ratio of the kinetic and interaction energy reaches smaller values, see for instance [32] and references given therein.

Model problem. As model problem, we consider the following (dimensionless) time-dependent nonlinear Schrödinger equation for the complex-valued (macroscopic) wave function , subject to an initial condition in Wentzel-Kramers-Brillouin (WKB) form with real-valued (and sufficiently regular) functions and to asymptotic boundary conditions on the unbounded spatial domainThe partial differential equation involves a real-valued confining (sufficiently regular) potential , the coupling constant , and the parameter . Concerning existence and uniqueness results for time-dependent nonlinear Schrödinger equations, we refer to the monographs [17,50].

Space and time discretisation. For the space and time discretisation of the model problem (1) we study the Fourier-pseudo spectral method versus the finite element method as well as higher-order exponential operator splitting methods, see for instance [13,14,27]. Contrary to the Fourier-pseudo spectral method with basis functions supported on the entire domain, realised by fast Fourier transform techniques and thus constrained to uniform meshes, the finite element method is designed for a local adaptation of the spatial grid. In the context of exponential operator splitting methods recently exploited approaches based on a posteriori local error estimators [5] or embedded splitting pairs [34], respectively, enable the use of adaptive time stepsize selection strategies with relatively low additional computational effort.

Time-splitting approach. As standard for a time-stepping approach, numerical approximations to the values of the analytical solution to the initial-boundary value problem (1) at certain time grid points with associated time increments for are obtained from a recurrence relation of the forminvolving the numerical evolution operator . Exponential operator splitting methods utilise a natural decomposition of the function defining the right-hand side of the differential equation into (at least) two parts and the presumption that each of the resulting subproblems is solvable in an efficient (and accurate) manner. On the one hand, time-splitting methods applied to low-dimensional nonlinear Schrödinger equations such as (1) in most cases rely on the numerical solution of the linear subproblemIn the numerical computations we will utilise that the solution to the model problem (1) remains localised due to the confining potential so that the unbounded spatial domain can be restricted to a sufficiently large bounded domain . In particular, in connection with the spatial discretisation of (3a) by the Fourier pseudo-spectral or finite element method, respectively, we will choose the domain such that the perturbations caused by the imposed artificial periodic or homogeneous Dirichlet boundary conditions, respectively, are insignificant. In more involved situations, where this simplifying presumption is not justified, alternative though considerably more costly approaches rely on the Hermite pseudo-spectral method or the incorporation of suitable boundary conditions, see [3,4,16,56-58]. On the other hand, due to the invariance property of the associated analytical solution, the subproblemreduces to a linear initial value problem with solution given through pointwise multiplicationcollocated at the grid points corresponding to the space discretisation of (3a). Altogether, compositions of the (numerical) solutions to the subproblems (3) with suitably adjusted time increments in the substeps then yield approximations (2) of a certain order of consistency.

Error and efficiency behaviour of discretisations. For nonlinear Schrödinger equations of the form (1) with parameter , numerous contributions confirm the favourable behaviour of pseudo-spectral time-splitting methods regarding stability, accuracy, efficiency, and the preservation of conserved quantities; as a small excerpt, we mention the works [8,45]. Furthermore, numerical comparisons given in [16] show that higher-order exponential operator splitting methods and, in particular, a fourth-order scheme with small effective error constant proposed in [12] outperform standard time integrators such as explicit Runge–Kutta methods when lower tolerances are required or long-term computations are carried out. Numerical evidence is affirmed by a theoretical stability and convergence analysis of splitting methods for (non)linear evolution equations with sufficiently regular solutions and applications to Schrödinger equations [30,36,41,51]; global error estimates for full discretisations based on pseudo-spectral time-splitting methods applied to Gross–Pitaevskii systems are deduced in [23,52].

Error and efficiency behaviour of discretisations in the semi-classical regime. For parameter values , the numerical solution of Schrödinger equations becomes a demanding task. Numerical studies and theoretical results for pseudo-spectral time-splitting methods given in [9,10,19,20], e.g., confirm that due to the highly oscillatory and nonlinear nature of problems such as (1) the spatial mesh size and the time increments need to be of the size of the parameter ε, especially when it is desired to capture correctly the quantitative behaviour of the wave function itself. Consequently, for smaller parameter values, the required high resolution in space and time renders the numerical computation a challenge, already in two space dimensions. For a review of numerical methods that are taylored for linear Schrödinger equations in the semi-classical regime and their limitations when applied to nonlinear problems, we refer to [32] and the references given therein, see also [37,38].

Objectives. In this work, our concern is to provide a numerical study of space and time discretisations for low-dimensional nonlinear Schrödinger equations. Primarily, we intend to exploit the capacity of locally adaptive integration methods. As a basic principle, the main benefits of adaptive discretisations are an improved balance between required accuracy and efficiency as well as an enhancement of the reliability of the numerical computations. For this reason, we find it worthwhile to investigate the performance of the Fourier-pseudo spectral method constrained to uniform meshes versus the locally adaptive finite element method for the model problem (1) in different parameter regimes. Regarding the time integration by higher-order exponential operator splitting methods, we expect the local time stepsize control to be effective, in particular, in the following situations. On the one hand, for Schrödinger equations with mildly varying regular solutions, there is a potential gain in efficiency. On the other hand, for problems (1) with focusing nonlinearity, a breathing-like solution with a sharp initial Gaussian, or for smaller values of the parameter ε, especially when it is desired to quantitatively resolve the values of the analytical solution, an automatic time stepsize control facilitates to determine the time increments sufficiently small in order that the numerical approximation captures correctly the solution behaviour. We point out that the application of a standard time stepsize control is only justified in situations where the analytical solution to the considered problem is sufficiently regular such that the nonstiff order of convergence is retained. For this reason, we focus on situations where this regularity requirement is fulfiled and no order reduction is encountered due to the lack of regularity of the data of the problem or the effect of the imposed boundary conditions, see also [51,52] for numerical counter-examples.

Structure of the manuscript. The present manuscript is organised as follows. In Section 2, we state the considered space and time discretisations based on the Fourier pseudo-spectral versus the finite element method and higher-order exponential operator splitting methods. In particular, we detail the strategies for a local adaptation of the spatial mesh as well as the approaches of a posteriori local error estimators and embedded splitting pairs for an adaptive time stepsize control, used in our implementation. Numerical experiments for the Fourier-pseudo spectral versus the finite element method and a first-order versus a fourth-order time splitting method with variable time stepsizes are presented in Section 3. Conclusions and perspectives are finally given in Section 4.

Space and time discretisations

In this section, we recapitulate the necessary basics on the space and time discretisations under consideration, especially on the employed locally adaptive meshing strategies.

Space discretisation by Fourier pseudo-spectral and finite element method

Fourier pseudo-spectral versus finite element method. For the spatial discretisation of the model problem (1) we consider two different approaches, namely, the Fourier pseudo-spectral versus the finite element method; for detailed information and instructive examples we also refer to the monographs [13,14,26,53]. Both approximation methods have in common that they rely on a representation of the numerical solution as a (finite) linear combination of basis functions; in the present situation, this affects the solution to the linear subproblem (3a). However, as one of their main differences, the Fourier basis functions are supported on the entire domain, contrary to the finite element basis functions that are supported only locally and thus designed for a local adaptation of the grid. Furthermore, the efficient realisation of the Fourier pseudo-spectral method through fast Fourier transform techniques is constrained to uniform meshes, whereas the implementation of local mesh adaptations is practicable in connection with the finite element method.

Fourier pseudo-spectral method

Basics (Fourier spectral method). A theoretical basis of the Fourier pseudo-spectral method is provided by the theory of selfadjoint operators and Sobolev spaces, see [54]. In particular, it is utilised that the associated basis functions , given byfor any multi-index , form a complete orthonormal system of the Lebesgue space and that they fulfil the following eigenvalue relation involving the Laplace operator here, the domain is assumed to be the cartesian product of bounded and, for simplicity, symmetric intervals with (sufficiently large) for . Furthermore, Parseval’s identity yieldsNote that we define the inner product in such that it is sesquilinear with respect to the second component; as standard, denotes the complex conjugate.

Numerical approximation. For the numerical realisation of the relations in (4) the infinite index set is restricted to a finite set and the spectral coefficients are approximated by means of the trapezoidal rule. More precisely, for some multi-index comprising positive and even integer numbers for , the finite setreplaces the index set . Moreover, imposing periodicity of the considered functions, approximations to the spectral coefficients are computed throughResults on the approximation rate of the Fourier pseudo-spectral method as well as instructive numerical examples are found in [13,26,53]. For instance, under suitable regularity requirements, the uniform boundholds for the truncation error, see [13, Ch. 2.12, Th. 9]. Furthermore, in [24, Lem. 4] it is shown that the approximation error of the Fourier pseudo-spectral method, caused by a truncation of the infinite series and an application of the trapezoidal rule, fulfils the estimatewhere is given by (5a) with for ; as usual, denotes the Sobolev space with exponent , see also [1]. Concerning the practical realisation of the Fourier pseudo-spectral method, fast Fourier transform techniques are utilised, see also [13,47] for detailed information; for instance, in a single space dimension, the number of required operations is .

Numerical solution of the linear subproblem. By means of the Fourier spectral decompositionand the corresponding eigenvalue relation for the Fourier basis functions, the representationfor the analytical solution to the linear subproblem (3a) follows, see (4); in particular, Parseval’s identity (4d) implies for if . As described before, truncating the infinite series and applying the trapezoidal rule for the numerical computation of the spectral coefficients further yields an approximation throughsee (5). We recall that for a function and its Fourier series equality generally holds in only. However, in the context of the numerical solution of the nonlinear Schrödinger equation (1), we tacitly assume that the analytical solution and in particular the considered initial condition satisfy suitable regularity requirements so that the arising pointwise evaluations (at the quadrature nodes) are well-defined; the investigation of such questions is part of a stability and convergence analysis of time-splitting methods, see for instance [36].

Finite element method

Preliminaries. As an alternative to the Fourier pseudo-spectral method, we apply the finite element method for the spatial discretisation of the linear subproblem (3a). In view of the considered model problem (1) involving a confining potential and asymptotic boundary conditions, we replace the unbounded domain by a (sufficiently large) bounded domain and impose homogeneous Dirichlet boundary conditions; in our implementation, for simplicity, we again let be the cartesian product of (sufficiently large and symmetric) bounded intervals, see also SubSection 2.1.1 Employing an equivalent formulation of the nonlinear Schrödinger equation (1) and especially of the linear subproblem (3a) as real-valued system, it meanwhile suffices to consider real-valued functions in only. In the following, in analogy to the corresponding relations for the Fourier-pseudo spectral method, we employ the same notation for analogical (though differing) quantities.

Basics (Galerkin approach). The finite element method relies on the Galerkin approach briefly described in the following, see also [13,14]. We consider a normed space and suppose that the family of finite dimensional subspaces satisfies the property of limited completeness . Further, we assume that any subspace is spanned by (suitably chosen) locally supported basis functions , that is, it holds ; clearly, the dimension of the linear space is equal to , the number of elements in . The above assumptions implywith denoting the projector on . For theoretical considerations it is often convenient to suppose that the basis functions form a Galerkin basis such that and . However, as this assumption does not hold in our realisation of the Galerkin approach, we indicate the dependence of the basis functions on K, unless it is a fixed number.

Galerkin approximation. For the following, we consider and, for notational simplicity, neglect the dependence of on K. For a given function its Galerkin approximation is defined by the requirementRepresenting as a linear combination of the (real-valued) basis functions spanning , this yields a system of linear equations for the unknown coefficients involving the (computable) quantities for .

Error estimate for the Galerkin approximation. In the above situation, it is straightforward to deduce an estimate for the error of the Galerkin approximation in . An application of (8a) with and the inequality of Cauchy–Schwarz yieldsand thus ; recall that denotes the projector on . Consequently, by means of the triangle inequality, the boundwith c  = 2 follows, which relates the error of the Galerkin approximation to results from approximation theory.

Realisation of the Galerkin approach (Finite element method). In regard to our finite element implementation based on the library DEAL.II [6,7], developed by Wolfgang Bangerth and collaborators, we consider finite dimensional subspaces , spanned by (real-valued) piecewise polynomial functions on a regular grid. Especially, in two space dimensions, we use piecewise quadratic basis functions on a rectangular grid. More precisely, the basis functions are supported on a single rectangle (finite element) and defined as the quadratic interpolants through the values at the vertices and the midpoints of the edges; thus, the overall number of basis functions is the number of rectangles times 9. Due to the fact that the inner product of two basis functions supported on different rectangles vanishes, the coefficient matrix defining the system (8b) is sparse. In our implementation, using an interface to UMFPACK, the sparse linear systems are solved through LU-decompositions.

Discretisation of the linear subproblem. Regarding the implementation of the finite element method for the space discretisation of the model problem (1), it is convenient to reformulate the Schrödinger equation as system for the real-valued functions , defined by for and . In particular, the linear subproblem (3a) is rewritten as followsThe variational formulation of the above problem utilises partial integrationnote that the boundary terms vanish due to the imposed homogeneous Dirichlet boundary conditions and that the notation , which is convenient but slightly incorrect, is understod as given in (10). Moreover, replacing the functions by their Galerkin approximations , given bythe following system of linear differential equations for the time-dependent coefficients resultsas before, we set and further for . In compact form, the above initial value problem is rewritten asinvolving the unknowns and as well as the (computable) quantities and . As a first approach, to avoid time discretisation errors and stability issues, the solution of the spatially semi-discretised linear subproblem (11) is computed through the application of the matrix exponential to the starting values and .

Local mesh adaptation and the Dörfler marking strategy. As pointed out before, the finite element approach with locally supported basis functions is designed for a local adaptation of the spatial mesh. In the above formalism, iterative refinements of the mesh correspond to the generation of a sequence of subspaces with increasing dimension. In our implementation, we utilise a standard a posteriori estimator for the local error in space, see [33]; that is, for a problem in two space dimensions, the quantitywhich measures the jumps of the normal derivatives of the qudratic interpolant f along the faces of a rectangle R, is evaluated. In a mesh refinement step, rectangles with an error exceeding some prescribed threshold are subdivided into four rectangles. As the iterative refinement of a spatial mesh requires an adaptation of the threshold in order to obtain asymptotically efficient meshes, we further employ the Dörfler marking strategy [21] known to be asymptotically the best refinement strategy. For a set of rectangles defining the actual rectangular mesh the smallest possible subset is computed such that the conditionis fulfiled for some prescribed threshold c, which is independent of the actual refinement step; then, the rectangles contained in are refined subsequently. In a fully adaptive discretisation method a mesh refinement (or coarsening) is carried out at each time step by means of a residual-based error estimator; however, as in the numerical illustrations we study either the spatial approximation error or the behaviour of the time integration methods for a sufficiently accurate space discretisation, we do not further elaborate this point in the present work.

Error estimates. Various contributions provide error estimates for finite element space discretisations, see for instance [14] in the context of the Poisson equation or, more generally, time-independent equations involving elliptic differential operators. However, to our knowledge, it remains open to deduce a global error estimate for full discretisations of the nonlinear Schrödinger equation (1) based on a Galerkin approach combined with a higher-order time-splitting method as well as an appropriate time integration method for the numerical solution of the linear subproblem (3a).

Time discretisation by splitting methods

In this section, we introduce the general form of a higher-order exponential operator splitting method and further describe the two approaches of a posteriori local error estimators and embedded splitting pairs, which we use in our adaptive time stepsize control.

Exponential operator splitting methods

Abstract formulation and Lie-calculus. For the specification and theoretical analysis of higher-order exponential operator splitting methods applied to nonlinear Schrödinger equations, it is convenient to rewrite (1) as an abstract initial value problem of the formand to employ the formal calculus of Lie-derivatives, see for instance [22,42,35,48] in the context of numerical methods for Hamiltonian problems. The powerfulness of this abstract formulation and the calculus of Lie-derivatives arises in the similarities to ordinary differential equations and the formally straightforward extension of the less involved linear case. In connection with unbounded operators, we find it advantageous to employ the following defining relations for the evolution operator and the associated Lie-derivativewhere denotes the flow operator associated with F, that is, gives the value of the solution to the initial value problem (12) at time t; throughout, we tacitly suppose that the arising (un)bounded (non)linear operators such as and are defined on suitably chosen subspaces of the underlying function space so that compositions thereof remain well-defined.

General form of exponential operator splitting methods and examples. For the time integration of the nonlinear evolutionary problem (12) we apply an exponential operator splitting method based on the (approximate) solution of the subproblems and . We employ the following general form of a higher-order splitting schemewith real method coefficients ; here, we define the product downwards. The above scheme includes various example methods proposed in the literature, see for instance [12,27,39]. Presumably, one of the most widely-used splitting schemes is the second-order Strang [49] splitting method, given byunder suitable regularity requirements on the data of the problem, it can be proven that the classical order of convergence, reflected in the (nonstiff) order conditions, is also retained for linear and nonlinear Schrödinger equations, see also [23,30,36,51,52].

Practical realisation of time-splitting methods. As set forth in the introduction, the time integration of the model problem (1) and nonlinear Schrödinger equations of a similar form by exponential operator splitting methods relies on the numerical solution of suitably chosen subproblems. In connection with the Fourier pseudo-spectral method, it is evident to decompose the function defining the partial differential equation into a linear part comprising the Laplacian and the remaining nonlinearity, that is, the main computational effort is in the numerical solution of subproblem (3a) by fast Fourier transforms, whereas subproblem (3b) is solvable by rapid pointwise multiplications. Alternatively, numerical approximations to the solution of subproblem (3b) are obtained by a finite element spatial semi-discretisation combined with an appropriate time integration method.

Adaptive time stepsize control

Adaptive time stepsize control. According to [47, Ch.16.2], the purpose and benefit of an adaptive time stepsize control is to achieve some predetermined accuracy in the solution with minimum computational effort. In our implementation, we utilise a standard local error control, see for instance [28], adjusting the actual time stepsize τ by the relationwith p the order of the time integrator, the prescribed tolerance, and an estimate of the local error; further flexibility is provided by additional parameters adapted to the problem under consideration. In the following, we give a brief description of two recently exploited approaches for obtaining local error estimators for exponential operator splitting methods. We point out that the convergence analysis in [36,52] ensures that time-splitting methods retain their (non)stiff orders when applied to the model problem (1) with sufficiently regular data. However, for Schrödinger equations involving less regular data an order reduction is encountered constricting the applicability of a standard time stepsize control, see also [51,52] for numerical illustrations.

Local error estimation by embedded splitting pairs. The idea of embedded exponential operator splitting pairs is a rather intuitive one, see [34]. In order to keep the additional computational effort for the local error estimator low, for a given splitting method of order with coefficients , the basic integrator, another splitting scheme of order with coefficients , the local error estimator, is constructed in such a way that certain of the first subsequent compositions coincide; evidently, the computation of the remaining compositions is parallelisable. The value of the numerical solution computed by means of the higher-order splitting method then replaces the unknown analytical solution value, that is, an approximation to the local error is obtained through the difference of the two approximate solutionsIn practise, the roles of the two splitting schemes is often reversed, that is, the time integration is performed by the higher-order scheme and the lower-order scheme is utilised for the estimation of the local error.

Embedded splitting pair of orders 4(3). As an illustration, we include the coefficients of an embedded splitting pair of orders 4(3), see Table 1 . We choose a favourable fourth-order splitting scheme by Blanes and Moan [12] with real (and negative) coefficients , appropriate for the time integration of Hamiltonian systems, as basic integrator. The construction of an embedded third-order scheme involving compositions is then straightforward; in the present situation, an elegant and still practicable tool is the computation of a Gröbner basis of the corresponding order conditions stated for example in [51]. The given embedded splitting scheme is obtained under the assumptions that the first four compositions coincide with the basic integrator, i.e. , and that further .

Table 1

Coefficients of an embedded exponential operator splitting pair of orders 4(3) based on a fourth-order scheme by Blanes and Moan [12].

j	a_j	j	b_j
1	0	1,7	0.0829844064174052
2,7	0.245298957184271	2,6	0.3963098014983680
3,6	0.604872665711080	3,5	-0.0390563049223486
4,5	1/2-(a_2+a_3)	4	1-2(b_1+b_2+b_3)
j	aˆ_j	j	bˆ_j
1,2,3,4	a_j	1,2,3,4	b_j
5	0.3752162693236828	5	0.4463374354420499
6	1.4878666594737946	6	-0.0060995324486253
7	-1.3630829287974774	7	0

A posteriori local error estimators. A more systematic approach compared to embedded splitting pairs is the construction of asymptotically correct a posteriori local error estimators for splitting methods; for the derivation of a posteriori error estimators for the first-order Lie–Trotter splitting and the second-order Strang splitting method applied to a linear evolution equation and a theoretical error analysis in the context of linear Schrödinger equations, see [5]. The approach relies on the derivation of a suitable local error representation in integral form by utilising a differential equation for the numerical evolution operator and the variation-of-constants formula; then, quantities involving the unknown analytical solution are substituted by suitable approximations and the integral is replaced by an appropriate quadrature formula approximation.

A posteriori local error estimator for Lie–Trotter splitting applied to linear problems. As illustration, we indicate the approach in the special case of the Lie–Trotter splitting method applied to the evolutionary problem (12) with linear operators A and B; in the present situation, the exact and numerical evolution operators are linear with respect to the initial value, given by for . On the one hand, inserting the Lie–Trotter splitting operator into the evolution equation, defines the defectOn the other hand, inserting the evolution operator into the Sylvester equation satisfied by the splitting operator a relation for the truncation error followsTaking the difference of the relations in (13), representing the solution in integral form, replacing the evolution operator by the splitting operator, which is a computable quantity, and employing a quadrature formula approximation finally yields the following a posteriori local error estimatorUnder the assumption that the analytical solution satisfies suitable regularity requirements, it is shown in [5] that the above a posteriori local error estimator, when applied to a linear Schrödinger equation involving a bounded potential, is asymptotically correct, that is, it holdswith ; recall also the a priori local error estimate for the Lie–Trotter splitting method.

A posteriori local error estimator for Lie–Trotter splitting applied to nonlinear problems. In the following, we indicate the extension of the a posteriori local error estimator (14) to nonlinear problems. A formal extension of the above relations by means of the Lie-calculus, namely, replacing the involved operators A and B by their Lie-derivatives and reversing the order of the compositions, results inThe specification to the model problem (1) yieldshere, the computation of corresponds to the numerical solution of the linear subproblem (3a), which is already part of the Lie–Trotter splitting solution. Consequently, in connection with the Fourier pseudo-spectral method, the realisation of the a posteriori local error estimator for the Lie–Trotter splitting methods requires additional fast Fourier transforms for the computation of and as well as further (inexpensive) pointwise multiplications.

Numerical experiments

Six months in the lab can save you a day in the library. 1 A picture is worth a thousand words.

In this section, we present numerical experiments for the previously introduced space and time discretisations applied to the model problem (1).

Fourier pseudo-spectral versus finite element approximation

In the following, we compare the error behaviour of the Fourier pseudo-spectral versus the finite element approximation for the WKB-type functionchosen in the lines of [10, Ex.4.6]; in one and two space dimensions and for different parameter values the profiles are displayed in Fig. 1 . We first include illustrations for the case of a single space dimension and then describe numerical examples for the more involved two-dimensional case. Throughout, we choose the spatial domain with . In the context of the Fourier pseudo-spectral method, we recall the error bounds (6) and (7) as well as the estimate (9) for the Galerkin approach.

Fig. 1

Row 1: (red line) and (dashed black line) for the WKB-type function (16) with and for (left to right). Row 2: for (16) with and for (left to right). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Significant Fourier pseudo-spectral coefficients. For the WKB-type function given by (16) with and parameter values for we determine the Fourier pseudo-spectral representationfor Fourier basis functions, see also (5a). In Fig. 2 , we display the absolute values of the significant Fourier pseudo-spectral coefficients satisfying versus the (shifted) index numbers (starting with one). In addition, the dependence of the numbers of significant coefficients on the inverse of the parameter is illustrated for further values of ε. The numerical results indicate that the number of significant Fourier pseudo-spectral coefficients is proportional to , that is, for a satisfactory approximation a spatial mesh size of the size of the parameter is needed.

Fig. 2

Columns 1–4: Significant Fourier pseudo-spectral coefficients of the WKB-type function (16) with and for (left to right). Absolute value versus (shifted) index k. Column 5: Dependence of the numbers of significant Fourier pseudo-spectral coefficients on .

Error of the Fourier pseudo-spectral approximation. In order to determine the error of the Fourier pseudo-spectral approximation for the WKB-type function (16) with and for , we evaluate at spatial meshes comprising , equidistant grid points; the function values at the finest mesh with grid points serve as reference solution values. We compute the corresponding Fourier pseudo-spectral coefficients and then evaluate the series for at the finest mesh, see (17); for the difference we determine the discrete maximum norm as well as the discrete -norm through Parseval’s identity. The obtained results are displayed in Fig. 3 ; as expected, for an accurate approximation is obtained for relatively few Fourier basis functions, whereas K has to be increased considerably for smaller parameter values. From the numerical observations in a single space dimension we draw the conclusion that the required high spatial resolution for smaller parameter values constricts the feasibility of numerical computations for the Fourier pseudo-spectral method in two and especially in three space dimensions, since spatial mesh sizes are required.

Fig. 3

Errors of the Fourier pseudo-spectral approximation for the WKBtype function (16) with and for (left to right).

Error of the finite element approximation. The corresponding numerical experiments on the error of the finite element approximation for the WKB-type function (16) with and for are given in Fig. 4 . We choose spatial meshes comprising , equidistant grid points; the values at a finer mesh with 3 K equidistant grid points are used as reference solution values. The basis functions are supported on each subinterval and defined as linear interpolants through the values at the two endpoints or as quadratic interpolants through the values at the two endpoints and the midpoints, respectively. In particular for the slopes reflect the expected orders two and three, respectively.

Fig. 4

Errors of the finite element approximation for the WKB-type function (16) with and for (left to right). Piecewise linear (row 1) and quadratic (row 2) basis functions.

Error of the Fourier pseudo-spectral versus finite element approximation. In Fig. 5 , we display the errors of the Fourier pseudo-spectral and finite element approximations for the WKB-type function (16) with and for ; for the ease of presentation, we omit the corresponding results with respect to the maximum norm, which are similar. Our locally adaptive finite element implementation based on the library DEAL.II [6,7] uses the meshing strategies described in Section 2.1.2. As a further illustration, we include the errors for the first Hermite basis functionas well as a comparison of the results obtained for by the finite element approximation applied with and without Dörfler marking. As expected, the Fourier spectral approximation is superiour to the finite element approximation; nevertheless, also the finite element approximation yields accurate results for refined spatial meshes.

Fig. 5

Row 1: errors of the Fourier pseudo-spectral and finite element approximations versus the total number of basis functions for the Hermite basis function (left). Errors of the finite element approximation for applied with and without Dörfler marking for (middle) and (right). Row 2: errors for the WKB-type function (16) with and for (left to right).

Conclusions. The numerical experiments confirm that for the Fourier pseudo-spectral approximation is favourable in accuracy, due to the spectral convergence rate, and in efficiency, due the applicability of fast Fourier transform techniques. For smaller values of the required high resolution constricts the feasibility of accurate numerical computations in higher space dimensions for both, the Fourier pseudo-spectral and the finite element method.

Time integration by variable stepsize splitting methods

In the following, we investigate variable stepsize splitting methods for the time integration of the model problem (1). If not stated otherwise, we consider (1) on the bounded domain , with coupling constant , subject to a scaled harmonic potential involving positive weights and discretised in space by the Fourier pseudo-spectral method; the initial condition is defined by (16). As usual, in a single space dimension we write and for short.

Time integration of the model equation (1D). As a first illustration, we perform the time integration of the model problem (1) with , and . The problem is discretised in space by the Fourier pseudo-spectral method involving basis functions. For the time discretisation of (1) we apply the first-order Lie–Trotter splitting method and a fourth-order splitting scheme by Blanes and Moan, see Table 1; our local time stepsize control with tolerances for is based on the a posteriori local error estimator (15a) and the embedded third-order scheme with coefficients given in Table 1, respectively. The solution profiles for , the sections at time , and the generated sequences of time stepsizes are displayed in Fig. 6 . In all cases, there is a good agreement of the obtained numerical solutions for the two splitting schemes and the two tolerances. Furthermore, we include a comparison of the solution profiles for at time computed by the Fourier pseudo-spectral with domain for with basis functions and by the Sine and Hermite pseudo-spectral method with functions (see remarks below) as well as by the finite element method with piecewise linear basis functions; in all cases, the time integration is performed by the embedded 4(3) pair for . For instance, for the finite element method with , although a slight difference is perceptible compared to the numerical solution obtained by the Fourier pseudo-spectral method, the full discretisation error is dominated by the temporal error, and compared to the Fourier pseudo-spectral method a similar progression of the time stepsize control is observed.

Fig. 6

Time integration of the model problem (1) with , and . Row 1: solution profiles for and corresponding sections at time for (columns 1–2) and (columns 3–4). Rows 2–3: generated sequences of time stepsizes for (row 2) and (row 3). Results for (left to right). Row 4: comparisons of the solution profiles obtained for by different pseudo-spectral (Fourier, Sine, Hermite, ) and piecewise linear finite element () space discretisations and the embedded 4(3) time-splitting pair for . Associated time stepsizes for the finite element method with grid points (right).

Time integration of the model equation (3D). As an illustration in three space dimensions, we perform the time integration of (1) with and for ; we discretise the problem in space using Fourier basis functions and apply the embedded splitting pair of orders 4(3) with . The profile at time and the associated sequence of time stepsizes are shown in Fig. 7 . For comparison, we also include the corresponding results for the cases , obtained with 128 equidistant grid points in each space direction. In order to validate the obtained approximations for , more precisely, to study the effect of undesired perturbations from the artificial periodic boundary conditions, we apply two simple strategies, see also Fig. 6: On the one hand, we restrict the spatial domain to with , which has no significant effect, whereas a further reduction to changes the solution significantly; we point out that in the latter case an order reduction in time occurs and thus the local error estimator cannot be expected to be asympotically correct. On the other hand, we compare the obtained solution with the numerical solution computed by the (more costly) Hermite pseudo-spectral method, which intrinsically captures the asymptotic boundary conditions on , and observe a good agreement of the solutions; for details on the Hermite pseudo-spectral method, we refer to [16] and references given therein. In the present situation, the investigations for the one-and two-dimensional cases are useful indications for the more laborious three-dimensional case.

Fig. 7

Row 1: solution profiles at for the model problem (1) with (left to right), , and , computed with the embedded pair of orders 4(3) for . Row 2: associated time stepsizes for .

Time integration of the model equation involving a smaller parameter (1D). Numerical illustrations for (1) with , and , computed by the Fourier pseudo-spectral method with basis functions and the adaptive fourth-order splitting scheme, are displayed in Fig. 8 . As a further comparison, we include the approximations obtained for and different tolerances for ; for the largest tolerance the numerical result is still rather poor, whereas there is a good qualitative and quantitative agreement of the obtained approximations for smaller tolerances.

Fig. 8

Time integration of the model problem (1) with , and by the embedded 4(3) pair. Row 1: solution profiles for and corresponding sections at time for (columns 1–2) and (columns 3–4). Row 2: generated sequences of time stepsizes for (columns 1–2) and columns 3–4). Results for and tolerances (columns 1/3) and (columns 2/4). Row 3: profile for at time for tolerances (left to right).

Time integration of the model equation with additional lattice potential, focusing nonlineariry, and sharp initial Gaussian (1D). As further numerical illustrations, we perform the time integration of the model problem (1) with , and for the following cases:For the space discretisation we apply the Fourier pseudo-spectral method with basis functions, where for (i) and for (ii)-(iii); for the time discretisation we apply the embedded 4(3) splitting pair with . In Fig. 9 , we display the solution profiles and the associated time stepsizes. In all cases, the local error control yields a reliable and highly accurate numerical solution.

Fig. 9

Time integration of (1) with , and by the embedded 4(3) pair for . Solution profiles for sections (first row) and associated time stepsizes (second row). Left: additional lattice potential with and defocusing nonlinearity with for . Middle: focusing nonlinearity with for . Right: defocusing nonlinearity with and sharp initial Gaussian with for .

Additional potential with , defocusing nonlinearity, and initial Gaussian with

Focusing nonlinearity with and initial Gaussian (18) with .

Defocusing nonlinearity and sharp initial Gaussian (18) with .

Conclusions. The numerical examples confirm the benefits of adaptive higher-order splitting methods in different situations, with an enhancement of the reliability of the numerical computations and a gain in efficiency, in particular, for problems with a mildly varying solution.

Conclusions

In the present work, we studied adaptive space and time discretisations for the numerical solution of Gross–Pitaevskii equations with regular solutions. We included numerical comparisons for the Fourier pseudo-spectral versus the locally adaptive finite element method, and we illustrated the capacity of variable stepsize higher-order exponential operator splitting methods for the time integration of the model problem (1) for the regime (additional lattice potential, focusing nonlinearity, sharp initial Gaussian) and for smaller parameter values . The given numerical examples confirm that the Fourier pseudo-spectral method combined with adaptive higher-order time-splitting schemes is favourable. Although constrained to uniform meshes, the Fourier pseudo-spectral approximation is superiour to the finite element approximation in accuracy and efficiency, due to the retained spectral convergence rate and the applicability of fast Fourier transform techniques, in particular, for Gross–Pitaevskii equations with a well localised and mildly varying solution. The use of an automatic time stepsize control is beneficial in all situations as it improves the balance between required accuracy and efficiency and enhances the reliability of the numerical computations; without an automatic time stepsize control, in particular, for problems involving smaller parameter values , it is most likely that the time stepsizes prescribed by the user are either underestimated, which is disadvantageous in view of efficiency, or overestimated, which will usually lead to an insufficient approximation and require to restart the computations.

For future work it is intended to proceed with the theoretical and numerical investigation of adaptive space and time discretisations for low-dimensional (non)linear Schrödinger equations. This also includes more involved situations, where in the absence of a confining potential the solution is not well localised and thus the incorporation of suitable boundary conditions [3,57,58] is needed. To our knowledge, it remains open to provide a stability and error analysis for full discretisations based on finite element and higher-order time-splitting methods. In view of various applications, it is also of interest to study blow-up phenomena arising for instance in Gross–Pitaevskii equations with focusing nonlinearities or solutions with sharp spatio-temporal gradients such as in the formation of shock waves, see [2,11,15,18,29,31,40,44,50,55]. However, for this type of problems the lack of regularity of the solutions (in time) implies order reductions, and, as a consequence, local error estimators are no longer asymptotically correct, and the use of a standard (time) stepsize control is not justified. For solutions of lower regularity a standard error local control will systematically overestimate the time stepsizes, which might lead to non-observance of the singularities or the rejection of many steps and potentially failure of the method. Strategies to detect order reductions and the investigation of a suitably modified local error control shall be part of future work.